ditr2fft_fixpt_ccs.c

用dsp解压mp3程序的算法

/***********************************************************************************
* ditr2fft_fixpt_ccs.c - Fix-point complex radix-2 decimation-in-time FFT.
*                        Perform in place FFT, the output overwrite the input array
*
***********************************************************************************/
#include "def_complex_fixpt.h"        /* floating-point complex.h header file */
  
void ditr2fft_fixpt_ccs(complex *X, unsigned int EXP, complex *W, unsigned int SCALE)
{

  lcomplex ltemp;         /* temporary storage of complex variable */
  lcomplex U;             /* Twiddle factor W^k */
  unsigned int i,j;
  unsigned int id;        /* Index for lower point in butterfly */
  unsigned int N=1<<EXP;  /* Number of points for FFT */
  unsigned int L;         /* FFT stage */
  unsigned int LE;        /* Number of points in sub DFT at stage L
                             and offset to next DFT in stage */
  unsigned int LE1;       /* Number of butterflies in one DFT at
                             stage L.  Also is offset to lower point
                             in butterfly at stage L */
  int scale;
  scale = 1;     
  if (SCALE == 0)         
    scale = 0;
    
  for (L=1; L<=EXP; L++)  /* FFT butterfly */
  {
    LE=1<<L;        	  /* LE=2^L=points of sub DFT */
    LE1=LE>>1;      	  /* Number of butterflies in sub-DFT */
    U.re = 32767;
    U.im = 0;

    for (j=0; j<LE1;j++)
    {
      for(i=j; i<N; i+=LE) /* Do the butterflies */
      {
        id=i+LE1;
        ltemp.re = ((((long)(int)X[id].re*(long)U.re)>>15) - ((((long)(int)X[id].im*(long)U.im)>>15))>>scale);
        ltemp.im = ((((long)(int)X[id].im*(long)U.re)>>15) + ((((long)(int)X[id].re*(long)U.im)>>15))>>scale);

        X[id].re = (X[i].re>>scale) - (int)ltemp.re;
        X[id].im = (X[i].im>>scale) - (int)ltemp.im;

        X[i].re = (X[i].re>>scale) + (int)ltemp.re;
        X[i].im = (X[i].im>>scale) + (int)ltemp.im;
      }
            
      /* Recursive compute W^k as U*W^(k-1) */

      ltemp.re = ((((long)U.re*(long)(int)W[L-1].re)>>15) - (((long)U.im*(long)(int)W[L-1].im)>>15));
      U.im = ((((long)U.re*(long)(int)W[L-1].im)>>15) + (((long)U.im*(long)(int)W[L-1].re)>>15));
      U.re = ltemp.re;
    }
  }
}