/*
 * SRIFFT  Real IFFT S-function.  Assumes input data is conjugate symmetric.
 * DSP Blockset S-Function to perform an inverse FFT on symmetric data.
 *	This SIMULINK S-function computes the IFFT of the input vector;
 *	it is called by the FFT blocks in the DSP Blockset.
 *	It computes the length N FFT of conjugate-symmetric data
 *	using a length N/2 complex FFT.
 *
 *  Note: This file must be linked with an FFT routine that does the core
 *        computation.  The FFT routine itself can be coded in either C or
 *        assembly language.
 *
 *  Copyright 1995-2000 The MathWorks, Inc.
 *  $Revision: 1.13 $  $Date: 2000/05/05 20:16:16 $
 */

#define S_FUNCTION_NAME srifft

#include <string.h>
#include <math.h>  /* frexp() */

#include "dsp_sim.h"
#include "fft_rt.h"

#define SAMPLE_TIME ssGetArg(S,0)
#define NUM_ARGS    1


#ifdef MATLAB_MEX_FILE
#define MDL_CHECK_PARAMETERS
static void mdlCheckParameters(SimStruct *S) {
	int_T M = mxGetM(SAMPLE_TIME);
	int_T N = mxGetN(SAMPLE_TIME);

	if ( (M != 1) || ((N != 1) && (N != 2)) ) {
		ssSetErrorStatus(S, "The sample time must be a "
							"scalar or a 2-element row vector");
	}
}
#endif


static void mdlInitializeSizes(SimStruct *S)
{
	ssSetNumSFcnParams(S, NUM_ARGS);

#if defined(MATLAB_MEX_FILE)
    if (ssGetNumSFcnParams(S) == ssGetSFcnParamsCount(S)) {
        mdlCheckParameters(S);
        if (ssGetErrorStatus(S) != NULL) {
            return;
        }
    } else {
        return; /* Simulink will report a parameter mismatch error */
    }
#endif

    ssSetNumInputs(        S, DYNAMICALLY_SIZED);
    ssSetNumOutputs(       S, DYNAMICALLY_SIZED);
    ssSetDirectFeedThrough(S, 1);
    ssSetNumSampleTimes(   S, 1);
    ssSetOptions(          S, SS_OPTION_EXCEPTION_FREE_CODE |
                              SS_OPTION_USING_ssGetUPtrs);
}


static void mdlInitializeSampleTimes(SimStruct *S)
{
    ssSetSampleTimeEvent(S, 0, mxGetPr(SAMPLE_TIME)[0]);
    ssSetOffsetTimeEvent(S, 0, (mxGetN(SAMPLE_TIME) == 2)
                              ? mxGetPr(SAMPLE_TIME)[1] : 0.0);
}


static void mdlInitializeConditions(real_T *x0, SimStruct *S)
{
#ifdef MATLAB_MEX_FILE
	/*
	 * Check to make sure input size is a power of 2
	 */
	int_T numinputs=ssGetNumInputs(S);
	if (numinputs != -1) {
		int_T dummy;
		if (!(frexp((real_T) numinputs,&dummy) == 0.5)) {
			ssSetErrorStatus(S,"Width of input to FFT must be a power of 2");
			return;
		}
	}
	  if (numinputs == 1) {
		mexWarnMsgTxt("Computing the IFFT of a scalar input.\n");
	  }
#endif /* MATLAB_MEX_FILE */
}


static void mdlOutputs(real_T *y, const real_T *x, const real_T *u, 
                       SimStruct *S, int_T tid)
{
    /* 
     * Use a length N complex IFFT to compute a length 2*N real IFFT.
     * Equivalent MATLAB code:
     *    function x = rifft(X)
     *    
     *	N = length(X)/2;
     *
     *	W=exp(j*2*pi*(0:N-1)/(2*N));
     *
     *	X0 = X(1:N) + X(1+N:2*N);
     *	X1 = W(:).*(X(1:N) - X(1+N:2*N));
     *
     *	Q = X0 + j * X1;
     *	V = imag(Q) + j*real(Q);
     *	v = fft(V);
     *	q = imag(v)/N + j*real(v)/N;
     *
     *	x(1:2:2*N-1) = 0.5 * real(q);
     *	x(2:2:2*N) = 0.5 * imag(q);
     */
    
	if (ssIsMajorTimeStep(S)) {
		int_T   numoutputs = ssGetNumOutputs(S);
		int_T   N2         = numoutputs/2;
		int_T   N          = N2/2;		      /* length of complex IFFT */
		real_T  theta      = 8 * atan(1.0) / N2;
		real_T  ct         = cos(theta);
		real_T  st         = sin(theta);
		real_T  Wr         = 1.0;
		real_T  Wi         = 0.0;
		real_T *yr         = y;
		real_T *yi         = y+N2;
		int_T   i;

		UPtrsType uptr = ssGetUPtrs(S);

		if (N2==1) {
			/* scalar input */
			*yr = **uptr++;
			*yi = **uptr;
			return;
		} else {

			/*
			 * First, generate X0 and X1
			 */
			for (i = 0; i < N; i++) {

			  real_T Xr  = *uptr[i];
			  real_T Xi  = *uptr[i + N2];
			  real_T XrN = *uptr[i + N];
			  real_T XiN = *uptr[i + N2 + N];
			  real_T X0_r = Xr + XrN;
			  real_T X0_i = Xi + XiN;
			  real_T X1_r = Wr * (Xr - XrN) - Wi*(Xi - XiN);
			  real_T X1_i = Wr * (Xi - XiN) + Wi*(Xr - XrN);
			  real_T Wr_t  = Wr;
			  yr[i] = X0_r - X1_i;
			  yi[i] = X0_i + X1_r;
			  Wr = Wr_t * ct - Wi * st;
			  Wi = Wr_t * st + Wi * ct;
			}

			/*
			 * now call fft; it does most of the work.  Note Re and Im
			 * parts are swapped to get ifft.
			 */
			
			fft(N, yi, yr);
			y += N2;

			{
				real_T recip2N = 1.0/(real_T)N2;
				for (i = N - 1; i >= 0; i--) {
				  *(--y) = recip2N * yi[i];
				  *(--y) = recip2N * yr[i];
				}
			}
		}
	}
}

static void mdlUpdate(real_T *x, const real_T *u, SimStruct *S, int_T tid)
{
}

static void mdlDerivatives(real_T *dx, const real_T *x, const real_T *u, 
                           SimStruct *S, int_T tid)
{
}


static void mdlTerminate(SimStruct *S)
{
}


#include "dsp_trailer.c"
 
/* [EOF]: srifft.c */