ditr2fft_fixpt_intr.c

用dsp解压mp3程序的算法

/************************************************************************************
*  ditr2fft_fixpt_intr.c - Fix-point complex radix-2 decimation-in-time FFT algorithm
*                          using C5000 intrinsics for CCS
*
************************************************************************************/

#include "def_complex_fixpt.h"        /* floating-point complex.h header file */

void ditr2fft_fixpt_intr(complex *X, unsigned int EXP, complex *W, unsigned int SCALE)
{
    lcomplex ltemp;         /* temporary storage of complex variable */
    complex temp;
    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 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 = _lsmpy(X[id].re, U.re);
                temp.re = (_smas(ltemp.re, X[id].im, U.im)>>16);  
                temp.re = _sadd(temp.re, 1)>>scale; /* Rounding & scale */
                ltemp.im = _lsmpy(X[id].im, U.re);
                temp.im = (_smac(ltemp.im, X[id].re, U.im)>>16);
                temp.im = _sadd(temp.im, 1)>>scale; /* Rounding & scale */
                X[id].re = _ssub(X[i].re>>scale, temp.re);
                X[id].im = _ssub(X[i].im>>scale, temp.im);
                X[i].re = _sadd(X[i].re>>scale, temp.re);
                X[i].im = _sadd(X[i].im>>scale, temp.im);
            
            }
            
            /* Recursive compute W^k as U*W^(k-1) */

			ltemp.re = _lsmpy(U.re, W[L-1].re);
            ltemp.re = _smas(ltemp.re, U.im, W[L-1].im);  
            ltemp.im = _lsmpy(U.re, W[L-1].im);
            ltemp.im = _smac(ltemp.im, U.im, W[L-1].re);
            U.re = ltemp.re>>16;            
            U.im = ltemp.im>>16;
        }
    }
}