/*
 * SDSPDCT  DCT S-function.
 * DSP Blockset S-Function to perform a DCT on a complex vector.
 *
 *  Copyright 1995-2000 The MathWorks, Inc.
 *  $Revision: 1.9 $  $Date: 2000/05/05 20:16:00 $
 */

#define S_FUNCTION_NAME sdspdct

#include <math.h>
#include "dsp_sim.h"
#include "fft_rt.h"

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


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

static void mdlInitializeConditions(real_T *x0, SimStruct *S)
{
  /*
   * Check to make sure input size is a power of 2:
   */
  int_T numinputs = ssGetNumInputs(S);

#ifdef MATLAB_MEX_FILE
  if (numinputs != -1) {
    int_T dummy;
    if (frexp((real_T)numinputs, &dummy) != 0.5) {
      THROW_ERROR(S,"Width of input to DCT must be a power of 2");
    }
    if (numinputs == 1) {
      mexWarnMsgTxt("Computing the DCT 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) 
   */
  
  {
    real_T  pi = 4.0*atan(1.0);
    int_T   n;
    int_T   twice_numcplxinputs = numinputs;
    int_T   numcplxinputs       = numinputs/2;
    real_T  den                 = sqrt(twice_numcplxinputs);
    real_T *ww                  = ssGetRWork(S);
    
    for (n = 0; n < numcplxinputs; n++){
      /* first part is real */
      ww[n] = cos( n*pi / twice_numcplxinputs ) / den;
      
      /* second part is imag */
      ww[n + numcplxinputs] = -sin( n*pi / twice_numcplxinputs ) / den;
    }
    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);
    int_T     numcplxinputs = ssGetNumInputs(S)/2;
    int_T     Nfft          = 2*numcplxinputs;
    int_T     j;
  
    
    /* Suppose 10 complex inputs (10 real followed by 10 imag = 20 inputs)
     * 20 complex outputs (20 real followed 20 imag = 40 outputs     
     * 
     * Copy the reals to the first 10 elements in output
     *
     * Reverse copy the real input to the second 10 elements of the output 
     * so that 1 2 3 becomes 3 2 1
     *
     * Copy the imag elements to first 10 elements of output
     *
     * Reverse copy the imag input to the second 10 elements of the output 
     */
      
    for (j = 0; j<2; j++) {
		int_T i;

		/* Copy the reals to the first 10 elements in output */
		for (i = 0; i<numcplxinputs; i++){
		  *y++ = **(uptr++);
		}

		/*Reverse copy the real input to the second 10 elements of the output*/ 
		for(i = 0; i<numcplxinputs; i++){
		  *y++ = **(--uptr);
		}
		/*
		 * Input pointer was left at beginning of real input. 
		 * Move it to begining of imag before doing copies for imag 
		 */
		uptr += numcplxinputs;
	}

    /* Move y pointer to top of output vector before doing FFT */
    y -= 4*numcplxinputs;

    /* Do in-place double-length FFT: */
    fft(Nfft, y, y+Nfft);


    /* Multiply by weights: */
    {
      real_T *yr = y;
      real_T *yi = y + Nfft;
      real_T *wr = ssGetRWork(S);
      real_T *wi = wr + numcplxinputs;
      int_T   i;

      for (i=0; i<numcplxinputs; i++) {
		real_T new_re = *yr * *wr   - *yi * *wi;
		real_T new_im = *yr * *wi++ + *yi * *wr++;
		
		*yr++ = new_re;
		*yi++ = new_im;
      }
    }
  }
}


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


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