/*
 * SDSPIFFT  IFFT S-function.  
 * DSP Blockset S-Function to perform an IFFT on complex inputs. 
 * Output either real or complex data. In real output mode, the input
 * data is assumed to be conjugate symmetric. It computes the length N FFT
 * of conjugate-symmetric data using a length N/2 complex FFT.
 *
 *  Copyright 1995-2000 The MathWorks, Inc.
 *  $Revision: 1.20 $  $Date: 2000/06/11 23:24:13 $
 */
#define S_FUNCTION_NAME sdspifft
#define S_FUNCTION_LEVEL 2

#include "dsp_sim.h"
#include "dspfft_rt.h"
#include "dspifft_rt.h"

enum {INPORT=0, NUM_INPORTS}; 
enum {OUTPORT=0, NUM_OUTPORTS}; 
enum {BUFF_IDX=0, MAX_NUM_DWORKS};

enum {FCN_TYPE_ARGC=0, NCHANS_ARGC, NUM_ARGS};
#define FCN_TYPE_ARG (ssGetSFcnParam(S,FCN_TYPE_ARGC))
#define NCHANS_ARG   (ssGetSFcnParam(S,NCHANS_ARGC))

typedef enum {
    fcnOutputReal = 1,
    fcnOutputComplex,
    fcnOutputConjSym
} FcnType;

#define EDIT_OK(S, ARG) ((ssGetSimMode(S) != SS_SIMMODE_SIZES_CALL_ONLY) || !mxIsEmpty(ARG)) 


#define MDL_CHECK_PARAMETERS
static void mdlCheckParameters(SimStruct *S) 
{
#ifdef MATLAB_MEX_FILE
    real_T d;
    int_T  i;

    /* Function type: */
    if ( !mxIsDouble(FCN_TYPE_ARG)                  || 
         (mxGetNumberOfElements(FCN_TYPE_ARG) != 1) || 
         ((mxGetPr(FCN_TYPE_ARG)[0] < 1.0) && 
         (mxGetPr(FCN_TYPE_ARG)[0] > 3.0)) ) {
        THROW_ERROR(S,"Function type must be a scalar:"
            "1 (complex output), 2 (real output) or 3 (conjugate symmetric)");
    }

    /* Number of channels */
    if EDIT_OK(S, NCHANS_ARG) { 
        if ((mxGetNumberOfElements(NCHANS_ARG) != 1) || 
            !mxIsNumeric(NCHANS_ARG)  || 
            !mxIsDouble(NCHANS_ARG) || 
            mxIsEmpty(NCHANS_ARG)   ||
            mxIsComplex(NCHANS_ARG)) {
            THROW_ERROR(S, "Number of channels must be a scalar.");
        }

        d = mxGetPr(NCHANS_ARG)[0];
        i = (int_T)d;
        if ( (d != i) || (i < 1) ) {
            THROW_ERROR(S, "Number of channels must be an integer > 0");
        }
    }
#endif
}



static void mdlInitializeSizes(SimStruct *S)
{
    FcnType ftype;

    ssSetNumSFcnParams(S, NUM_ARGS);
#if defined(MATLAB_MEX_FILE)
    if (ssGetNumSFcnParams(S) != ssGetSFcnParamsCount(S)) return;
    mdlCheckParameters(S);
    if (ssGetErrorStatus(S) != NULL) return;
#endif

    ftype = (FcnType)((int_T)mxGetPr(FCN_TYPE_ARG)[0]);

    /* Inputs: */
    if (!ssSetNumInputPorts(S, NUM_INPORTS)) return;
    ssSetInputPortWidth(            S, INPORT, DYNAMICALLY_SIZED);
    ssSetInputPortDirectFeedThrough(S, INPORT, 1);
    ssSetInputPortComplexSignal(    S, INPORT, COMPLEX_YES);
    ssSetInputPortReusable(        S, INPORT, 1);

    /* Input is overwritable only if the user has requested a complex output
     * (i.e. real or complex mode, NOT conjugate symmetric mode): 
     */
    ssSetInputPortOverWritable(     S, INPORT, 1);

    /* Outputs: */
    if (!ssSetNumOutputPorts(S,NUM_OUTPORTS)) return;
    ssSetOutputPortWidth(S, OUTPORT, DYNAMICALLY_SIZED);
    {
        /* Note: Both real and complex modes generate complex results
         * Only conjugate symmetric generates real output 
         * In real mode, this block is followed by a "Real" block for efficiency
         */
        const FcnType ftype = (FcnType)((int_T)mxGetPr(FCN_TYPE_ARG)[0]);
        ssSetOutputPortComplexSignal(S, OUTPORT, (ftype == fcnOutputConjSym) ? COMPLEX_NO 
                                                                             : COMPLEX_YES);
        ssSetOutputPortReusable(    S, OUTPORT, 1);
    }

    if(!ssSetNumDWork(  S, DYNAMICALLY_SIZED)) return;

    ssSetNumSampleTimes(S, 1);
    ssSetOptions(       S, SS_OPTION_EXCEPTION_FREE_CODE |
                        SS_OPTION_USE_TLC_WITH_ACCELERATOR);
}


static void mdlInitializeSampleTimes(SimStruct *S)
{
    ssSetSampleTime(S, 0, INHERITED_SAMPLE_TIME);
    ssSetOffsetTime(S, 0, FIXED_IN_MINOR_STEP_OFFSET);
}



#define MDL_START
static void mdlStart(SimStruct *S)
{
#ifdef MATLAB_MEX_FILE
    const int_T   nchans = (int_T)mxGetPr(NCHANS_ARG)[0];
    const int_T   N      = ssGetInputPortWidth(S,INPORT) / nchans;
    int_T dummy;

    if (ssGetSampleTime(S,0) == CONTINUOUS_SAMPLE_TIME) {
        THROW_ERROR(S, "Continuous sample times not allowed for IFFT block.");
    }

    /* Check to make sure input size is a power of 2: */
    if (frexp((real_T)N, &dummy) != 0.5) {
        THROW_ERROR(S,"Width of input to IFFT must be a power of 2");
    }
    if (N == 1) {
      mexWarnMsgTxt("Computing the IFFT of a scalar input.\n");
    }
#endif
}


static void mdlOutputs(SimStruct *S, int_T tid)
{
    const FcnType ftype = (FcnType)((int_T)mxGetPr(FCN_TYPE_ARG)[0]);
    InputPtrsType uptr  = ssGetInputPortSignalPtrs(S,INPORT);
    const int_T   nchans = (int_T)mxGetPr(NCHANS_ARG)[0];
    const int_T   N      = ssGetInputPortWidth(S,INPORT) / nchans;

    switch(ftype) {
        case fcnOutputComplex:
        case fcnOutputReal:
        {
            creal_T *y  = (creal_T *)ssGetOutputPortSignal(S,OUTPORT);
            const boolean_T need_copy = (boolean_T)(ssGetInputPortBufferDstPort(S, INPORT) != OUTPORT);
            int_T f;

            for (f = 0; f++ < nchans; ) {

                if (N == 1) {
                    /* Scalar input - output is the input */
                    if (need_copy) {
                        *y++ = *((creal_T *)(*uptr++));
                    }

	        } else {
	          /* Copy inputs to outputs if not sharing buffer */
	          if (need_copy) {
                    int_T i;
                    for (i=0; i<N; i++) {
                        *y++ = *((creal_T *)(*uptr++));
                    }
                    y -= N;
	          }

	          dspifft(N, y);
                  y += N;
                }
            }
        }
        break;

        case fcnOutputConjSym:
        {
            real_T *y = ssGetOutputPortRealSignal(S,OUTPORT);
            int_T f;

            for (f = 0; f++ < nchans; ) {

                if (N == 1) {
                    /* Scalar input - output is just the real element */
                    const creal_T *u = (creal_T *)(*uptr++);
                    *y++ = u->re;

                } else {

                    creal_T *buff = (creal_T *)ssGetDWork(S, BUFF_IDX); 
                    creal_T  W    = {1.0, 0.0};
                    creal_T  twid;
                    real_T   theta = -8 * atan(1.0) / N;
                    int_T    N2 = N/2;
                    int_T    i;

                    /* Complex exponential: cos+jsin */
                    twid.re = cos(theta);
                    twid.im = sin(theta);

                    for (i = 0; i < N2; i++) {
                        creal_T a, b, c, d, e;

                        a = *((creal_T *)(uptr[i]));
                        b = *((creal_T *)(uptr[i+N2]));

                        c.re = a.re + b.re;
                        c.im = a.im + b.im;
    
                        d.re = a.re - b.re;
                        d.im = a.im - b.im;

                        e.re = CMULT_RE(d,W);
                        e.im = CMULT_IM(d,W);

                        /* Note: Re and Im parts are swapped to get ifft */
                        buff[i].im = c.re - e.im;
                        buff[i].re = c.im + e.re;

                        {
                            creal_T ctemp = W;
                            W.re = CMULT_RE(ctemp,twid);
                            W.im = CMULT_IM(ctemp,twid);
                        }
	            }

	            dspfft(N2, buff);
        
                    {
	                const real_T recip2N  = 1.0/(real_T)N;
                        const real_T last_val = buff->re;

                        for (i=0; i++ < N2-1;) {
                            *y++ = recip2N * (buff++)->im;
                            *y++ = recip2N * buff->re;
                        }
                        *y++ = recip2N * buff->im;
                        *y++ = recip2N * last_val;
                    }
                    /* Update input pointer for next frame */
                    uptr += N;
	        }
            }
        }
        break;
    }
}


static void mdlTerminate(SimStruct *S)
{
}


#if defined(MATLAB_MEX_FILE)
#define MDL_SET_WORK_WIDTHS
static void mdlSetWorkWidths(SimStruct *S)
{
    const FcnType ftype   = (FcnType)((int_T)mxGetPr(FCN_TYPE_ARG)[0]);
    const int_T   N       = ssGetInputPortWidth(S,INPORT);
    const int_T   N2      = N / 2;
    const int_T   nDWorks = ( (ftype == fcnOutputConjSym) && (N2 > 0) )
                             ? MAX_NUM_DWORKS : 0;

    if(!ssSetNumDWork(          S, nDWorks)) return;

    if (nDWorks > 0){
        ssSetDWorkWidth(        S, BUFF_IDX, N2);
        ssSetDWorkDataType(     S, BUFF_IDX, SS_DOUBLE);
        ssSetDWorkComplexSignal(S, BUFF_IDX, COMPLEX_YES);
    
    }
}
#endif

#ifdef	MATLAB_MEX_FILE
#include "simulink.c"
#else
#include "cg_sfun.h"
#endif

/* [EOF] sdspifft.c */