/*
 * SDSPIDCT2  IDCT S-function.
 * DSP Blockset S-Function to perform an IDCT.
 *
 *  Copyright 1995-2000 The MathWorks, Inc.
 *  $Revision: 1.13 $  $Date: 2000/06/14 14:05:01 $
 */
#define S_FUNCTION_NAME sdspidct2
#define S_FUNCTION_LEVEL 2

#include "dsp_sim.h"

enum {INPORT=0, NUM_INPORTS}; 
enum {OUTPORT=0, NUM_OUTPORTS}; 
enum {NUM_ARGS=0};
enum {WEIGHTS_IDX=0, BUFF_IDX, NUM_DWORK};

extern void dspfft(const int_T n, creal_T *y);
extern void dspifft(const int_T n, creal_T *y);


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

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

    /* Inputs: */
    if (!ssSetNumInputPorts(S, NUM_INPORTS)) return;
    if (!ssSetInputPortDimensionInfo(S, INPORT, DYNAMIC_DIMENSION)) return;
    ssSetInputPortFrameData(         S, INPORT, FRAME_INHERITED);
    ssSetInputPortDirectFeedThrough( S, INPORT, 1);
    ssSetInputPortComplexSignal(     S, INPORT, COMPLEX_INHERITED);
    ssSetInputPortReusable(          S, INPORT, 1);
    ssSetInputPortOverWritable(      S, INPORT, 0);

    /* Outputs: */
    if (!ssSetNumOutputPorts(S,NUM_OUTPORTS)) return;
    if (!ssSetOutputPortDimensionInfo(S, OUTPORT, DYNAMIC_DIMENSION)) return;
    ssSetOutputPortFrameData(         S, OUTPORT, FRAME_YES);
    ssSetOutputPortComplexSignal(     S, OUTPORT, COMPLEX_INHERITED);
    ssSetOutputPortReusable(          S, OUTPORT, 1);

    /* DWorks: */
    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)
{
    const boolean_T cmplx   = (boolean_T)(ssGetInputPortComplexSignal(S,INPORT) == COMPLEX_YES);
    const int_T  framebased = ssGetInputPortFrameData(S,INPORT);

    const int_T  numDims   = ssGetInputPortNumDimensions(S, INPORT);
    const int_T *inDims    = ssGetInputPortDimensions(S, INPORT);
    const int_T  inRows    = inDims[0];
    const int_T  inCols    = (numDims == 2) ? inDims[1] : 1;
    const int_T  frameSize = (inRows==1 && !framebased) ? inCols : inRows;
    
    const real_T piN2       = 2.0*atan(1.0)/frameSize;

    real_T       fac  = (cmplx) ? sqrt(2*frameSize) : sqrt(2*frameSize)/frameSize;
    creal_T     *ww   = (creal_T *)ssGetDWork(S, WEIGHTS_IDX);

#ifdef MATLAB_MEX_FILE
    if (ssGetSampleTime(S,0) == CONTINUOUS_SAMPLE_TIME) {
        THROW_ERROR(S, "Continuous sample times not allowed for IDCT block.");
    }
    
    /* Check to make sure input size is a power of 2: */
    {
        int_T  dummy;
        if (frexp((real_T)frameSize, &dummy) != 0.5) {
            THROW_ERROR(S,"Width of input to IDCT must be a power of 2");
        }
        if (frameSize == 1) {
            mexWarnMsgTxt("Computing the IDCT of a scalar input.\n");
        }
    }
#endif

  /*
   * Need to compute weights to multiply DFT coefficients.
   * By computing the real and imag parts
   *
   * ww = (exp(-i*(0:n-1)*pi/(2*n))/sqrt(2*n)).';
   * ww(1) = ww(1) / sqrt(2);
   *
   * exp(-i*x) = cos(x) - i*sin(x) 
   */
    
   /* Only need FrameSize number of coefficients */
    {
        int_T i;

        for (i = 0; i < frameSize; i++) {
            ww[i].re = cos(i*piN2) * fac;
            ww[i].im = sin(i*piN2) * fac;
        }
        
        if (!cmplx){
            ww[0].re /= sqrt(2.0);
        } else {
            ww[0].re *= sqrt(2.0);
        }
    }
}


static void mdlOutputs(SimStruct *S, int_T tid)
{
    const boolean_T cmplx = (boolean_T)(ssGetInputPortComplexSignal(S,INPORT) == COMPLEX_YES);

    const int_T  framebased = ssGetInputPortFrameData(S,INPORT);
    const int_T  numDims    = ssGetInputPortNumDimensions(S, INPORT);
    const int_T *inDims     = ssGetInputPortDimensions(S, INPORT);
    const int_T  inRows     = inDims[0];
    const int_T  inCols     = (numDims == 2) ? inDims[1] : 1;
    const int_T  frameSize  = (inRows==1 && !framebased) ? inCols : inRows;
    const int_T  nChans     = (inRows==1 && !framebased) ? 1      : inCols;
    int_T     ch;
    
    if (!cmplx) { /* REAL Input */
        
        InputRealPtrsType uptr = ssGetInputPortRealSignalPtrs(S, INPORT);
        real_T           *y    = ssGetOutputPortRealSignal(S,OUTPORT);
        
        
        if (frameSize == 1) {
            for(ch=0; ch<nChans; ch++) {
                *y++ = **uptr++;
            }
            
        } else {
            const int_T  N2   = frameSize/2; /* length of complex FFT       */
            int_T        i;
            
            for(ch=0; ch<nChans; ch++) {
                creal_T *ww   = ssGetDWork(S, WEIGHTS_IDX);
                creal_T *buff = ssGetDWork(S, BUFF_IDX);
                
                /* Multipy by weights */
                for (i = 0; i< frameSize; i++) {
                    buff[i].re =   **uptr * ww->re;
                    buff[i].im = -(**uptr++ * (ww++)->im);
                }
                
                dspfft(frameSize, buff); 
                
                for (i=0; i<N2; i++) {
                    *y++ = buff[i].re;
                    *y++ = buff[frameSize-1-i].re;
                }
            }
        }
        
    } else {  /* COMPLEX Input */
        
        InputPtrsType uptr = ssGetInputPortSignalPtrs(S, INPORT);
        creal_T      *y    = (creal_T *)ssGetOutputPortRealSignal(S,OUTPORT);
        
        for(ch=0; ch<nChans; ch++) {
            
            creal_T      *ww   = (creal_T *)ssGetDWork(S, WEIGHTS_IDX);
            creal_T      *buff = (creal_T *)ssGetDWork(S, BUFF_IDX);
            int_T         i;
            
            /* Multiply by weights */
            for (i = 0; i< frameSize; i++) {
                const creal_T *u = (creal_T *)uptr[i];
                buff[i].re = CMULT_RE(*u,ww[i]);
                buff[i].im = CMULT_IM(*u,ww[i]);
            }
            
            buff[frameSize].re = 0.0;
            buff[frameSize].im = 0.0;
            
            for(i = 1; i < frameSize; i++) {
                const creal_T *u = (creal_T *)uptr[frameSize-i];
                buff[frameSize+i].re = CMULT_IM(*u, ww[i]);
                buff[frameSize+i].im = -CMULT_XYCONJ_RE(*u, ww[i]);
            }
            
            dspifft(2*frameSize, buff);
            
            for(i=0; i<frameSize; i++) {
                *y++ = *buff++;
            }
            
            uptr += frameSize;  /* Move to next channel input */
            
        } /* for(nChans) */
    }
}


static void mdlTerminate(SimStruct *S)
{
}


#if defined(MATLAB_MEX_FILE)
#define MDL_SET_INPUT_PORT_DIMENSION_INFO
static void mdlSetInputPortDimensionInfo(SimStruct *S, 
                                      int_T port,
                                      const DimsInfo_T *dimsInfo)
{
    if(!ssSetInputPortDimensionInfo(S, port, dimsInfo)) return;
    
    {
        const int_T outputPortWidth = ssGetOutputPortWidth(S,OUTPORT);
        const int_T numDims = dimsInfo->numDims;
        const int_T inRows  = dimsInfo->dims[0];
        const int_T inCols  = (numDims == 2) ? dimsInfo->dims[1] : 1;

        /* If the input is a sample-based row, then the output is a 
         * single channel frame-based column. (convenience mode)
         */
        const int_T framebased = ssGetInputPortFrameData(S,INPORT);
        const int_T outRows    = (framebased || inRows != 1) ? inRows : inCols;
        const int_T outCols    = (framebased || inRows != 1) ? inCols : inRows;
        
        if (outputPortWidth ==  DYNAMICALLY_SIZED) {
            if(!ssSetOutputPortMatrixDimensions(S, OUTPORT, outRows, outCols)) return;
        }
    }
}


#define MDL_SET_OUTPUT_PORT_DIMENSION_INFO
static void mdlSetOutputPortDimensionInfo(SimStruct        *S, 
                                          int_T            port,
                                          const DimsInfo_T *dimsInfo)
{
    if(!ssSetOutputPortDimensionInfo(S, port, dimsInfo)) return;
}
#endif


#if defined(MATLAB_MEX_FILE)
#define MDL_SET_WORK_WIDTHS
static void mdlSetWorkWidths(SimStruct *S)
{
    const boolean_T cmplx = (boolean_T)(ssGetInputPortComplexSignal(S,INPORT) == COMPLEX_YES);

    const int_T  framebased = ssGetInputPortFrameData(S,INPORT);
    const int_T  numDims    = ssGetInputPortNumDimensions(S, INPORT);
    const int_T *inDims     = ssGetInputPortDimensions(S, INPORT);
    const int_T  inRows     = inDims[0];
    const int_T  inCols     = (numDims == 2) ? inDims[1] : 1;
    const int_T  frameSize  = (inRows==1 && !framebased) ? inCols : inRows;
    
    if(!ssSetNumDWork(      S, NUM_DWORK)) return;

    ssSetDWorkWidth(        S, WEIGHTS_IDX, frameSize);
    ssSetDWorkDataType(     S, WEIGHTS_IDX, SS_DOUBLE);
    ssSetDWorkComplexSignal(S, WEIGHTS_IDX, COMPLEX_YES);

    ssSetDWorkWidth(        S, BUFF_IDX, (cmplx) ? 2*frameSize : frameSize); /* Length of IFFT */
    ssSetDWorkDataType(     S, BUFF_IDX, SS_DOUBLE);
    ssSetDWorkComplexSignal(S, BUFF_IDX, COMPLEX_YES);
}
#endif

#include "dsp_cplxhs11.c"

#include "dsp_trailer.c"

/* [EOF] sdspidct2.c */