/*
 * SDSPRIDCT  IDCT S-function.
 * DSP Blockset S-Function to perform an IDCT on a real vector.
 *
 *  Copyright 1995-2000 The MathWorks, Inc.
 *  $Revision: 1.6 $  $Date: 2000/05/05 20:16:06 $
 */


#define S_FUNCTION_NAME sdspridct

#include <math.h>
#include <string.h>

#include "dsp_sim.h"

extern void fft(int_T n, real_T *xr, real_T *xi);

static void mdlInitializeSizes(SimStruct *S)
{
  ssSetNumInputs(        S, DYNAMICALLY_SIZED);
  ssSetNumOutputs(       S, DYNAMICALLY_SIZED);
  ssSetNumSampleTimes(   S, 1);
  ssSetDirectFeedThrough(S, 1);
  ssSetNumRWork(         S, DYNAMICALLY_SIZED);
  ssSetOptions(          S, SS_OPTION_EXCEPTION_FREE_CODE |
                            SS_OPTION_USING_ssGetUPtrs);  
}


static void mdlInitializeSampleTimes(SimStruct *S)
{

  ssSetSampleTimeEvent(S, 0, INHERITED_SAMPLE_TIME);
  ssSetOffsetTimeEvent(S, 0, 0.0);
}


static void mdlInitializeConditions(real_T *x0, SimStruct *S)
{
  
  int_T numinputs = ssGetNumInputs(S);   /* real inputs */
  
#ifdef MATLAB_MEX_FILE
  /*
   * Check to make sure input size is a power of 2:
   */
  
  if (numinputs != -1) {
    int_T dummy;
    if (frexp((real_T)numinputs, &dummy) != 0.5) {
      ssSetErrorStatus(S,"Width of input to IDCT must be a power of 2");
      return;
    }
    if (numinputs == 1) {
      mexWarnMsgTxt("Computing the IDCT of a scalar input.\n");
    }
  }
#endif /* MATLAB_MEX_FILE */
  
  /*  Compute complex weights  */
  {
    const real_T piN2 = 2.0*atan(1.0)/ (numinputs);  /* pi / (2*N) */
    const real_T fac  = sqrt(2*numinputs);
    real_T      *ww   = ssGetRWork(S);
    int_T        i;
    
    for (i = 0; i < numinputs; i++) {
      ww[i]             =  cos(i*piN2) * fac;
      ww[i + numinputs] =  sin(i*piN2) * fac;
    }
    
    ww[0] /= sqrt(2.0);

  }
}


static void mdlOutputs(real_T *y, const real_T *x, const real_T *u, 
                       SimStruct *S, int_T tid)
{
  if (ssIsMajorTimeStep(S)) {
    UPtrsType uptr      = ssGetUPtrs(S);
    const int_T N       = ssGetNumInputs(S); /* real inputs */
	const int_T N2      = N/2;
    real_T   *yr        = y;
    real_T   *yi        = y+N;
    real_T   *wr        = ssGetRWork(S);
    real_T   *wi        = wr+N;
    int_T     i;
        
    /* Multipy by weights */

    for (i = 0; i< N; i++) {
      *yr++ = **(uptr + i) * *(wr + i);
      *yi++ = **(uptr + i) * *(wi + i);
    }
    yr -= N; yi -= N;

    for (i=0; i<N; i++) {
      *yi++ = -y[i + N];
    }
    
    fft(N, y, y+N); 

    yi -= N;  
    for(i = 0; i<N; i++) {
      *yr++ =  y[i]   / N;
      *yi++ = -y[i+N] / N;
      
    }
    yr -= N; yi -= N;

	/* 
	 * Need to resort values
	 * Use yi to store values during shuffle
	 */

	/* Move second half of real (N/2 to N) to imag */
	for(i = 0; i<N2; i++) {
		*yi++ = y[i+N2];
	}
	yi -= N/2;
	
	/* Move elements in yr to correct order */
	for(i = N2 - 1; i >= 0; i--) {
		y[i*2] = y[i];
	}

	/* Move elements in yi to correct order in yr */
	for(i = 1; i <=N2; i++) {
		yr[2*i - 1]= yi[N2- 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)
{
}


#ifdef MATLAB_MEX_FILE

#define MDL_GET_INPUT_PORT_WIDTH
/* Function: mdlGetInputPortWidth =============================================
 *
 */
static int_T mdlGetInputPortWidth(SimStruct *S, int_T outputPortWidth)
{
  if (outputPortWidth % 2 != 0) {
    ssSetErrorStatus(S,"Invalid output port width");
    return(1);
  }
  return(outputPortWidth/2);
}

#define MDL_GET_OUTPUT_PORT_WIDTH
/* Function: mdlGetOutputPortWidth ============================================
 *
 */
static int_T mdlGetOutputPortWidth(SimStruct *S, int_T inputPortWidth)
{
  return(inputPortWidth*2);
}

#endif

/* Function: mdlSetWorkWidths =================================================
 * Abstract:
 *
 */
#ifdef MATLAB_MEX_FILE
#define MDL_SET_WORK_WIDTHS
static void mdlSetWorkWidths(SimStruct *S)
{
  ssSetNumRWork(S, (int_T)(ssGetNumInputs(S)*2));
}
#endif

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