ditr2fft_fixpt.c

用dsp解压mp3程序的算法

/*****************************************************************************
*    ditr2fft_fixpt.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(complex *X, unsigned int EXP, complex *W, unsigned int SCALE)
{
  lcomplex ltemp;         /* temporary storage of complex variable          */
  complex 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 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)(int)U.re)>>15) - (((long)(int)X[id].im*(long)(int)U.im)>>15))>>scale;
        ltemp.im = ((((long)(int)X[id].im*(long)(int)U.re)>>15) + (((long)(int)X[id].re*(long)(int)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)(int)U.re*(long)(int)W[L-1].re)>>15) - (((long)(int)U.im*(long)(int)W[L-1].im)>>15));
      U.im = (int)((((long)(int)U.re*(long)(int)W[L-1].im)>>15) + (((long)(int)U.im*(long)(int)W[L-1].re)>>15));
      U.re = (int)ltemp.re;
    }
  }
}