fftfunction.c

TMS320VC5502的FFT算法程序验证

/********************************************************************************
*FFT函数                                                                    *
*********************************************************************************
*********************************************************************************
*																				*
* 函数定义：void fft1024(COMPLEX *Y, int N)                                     *
*           void fft512(COMPLEX *Y, int N)                                      *
*           void fft256(COMPLEX *Y, int N)                                      *
*  					                                                            *
* 功    能：基2的频率抽取FFT算法                								*
*																				*
* 入口参数：Y －－ 所需进行FFT变换数组的首地址 								    *
*          	N －－ FFT变换点数  	                                            *
* 出口参数：无										                            *
*                                                                               *
*********************************************************************************/
#include <csl_std.h>
#include "i_cmplx.h"     	/* definition of the complex type */  
#include "twiddle1024.h" 	/* quantised and scaled Twiddle factors */
#define LL 1024          	/* Maximum length of FFT */ 
extern unsigned int SampleLong;
/* fft radix-2 funtion using Decimation In Frequency */

//#pragma CODE_SECTION(fft, "mycode"); // Makes the program run from internal memory

void fft1024(COMPLEX *Y, int N) /* FFT(input sample array, # of points)     */
	{
   	Int32 temp1R, temp1I, temp2R,temp2I;  /* 32 bits temporary storage for */
   									   /* intermediate results			*/					
   	short tempR, tempI, c, s;      /* 16 bits temporary storages    */
   	/* variables */
   	Int32 TwFStep,  /* Step between twiddle factors */
   		  TwFIndex, /* Index of twiddle factors */
   		  BLStep,    /* Step for incrementing butterfly index */
   		  BLdiff,   /* Difference between upper and lower butterfly legs */
   		  upperIdx,
   		  lowerIdx, /* upper and lower indexes of buterfly leg */ 
   		  i, j, k;  /* loop control variables */
   		
   	BLdiff=N;
   	TwFStep=1;
   	for(k=N;k>1;k=(k>>1)) /* Do Log(base 2)(N) Stages */
   		{
      	BLStep=BLdiff;
      	BLdiff=BLdiff>>1; 
      	TwFIndex=0;
      	for(j=0;j<BLdiff;j++)/* Nbr of twiddle factors to use=BLDiff  */
      		{
      		c=w[TwFIndex].real;
      		s=w[TwFIndex].imag;
      		TwFIndex=TwFIndex+TwFStep;
      		for(upperIdx=j;upperIdx<N;upperIdx+=BLStep)
         		{                              
/* Calculations inside this loop avoid overflow by shifting left once
   the result of every adittion/substration and by shifting left 15 
   places the result of every multiplication. Double precision temporary
   results (32-bit) are used in order to avoid losing information because
   of overflow. Final DFT result is scaled by N (number of points), i.e.,
   2^(Nbr of stages) =2^(log(base 2) N) = N								*/
             		
            	lowerIdx=upperIdx+BLdiff;
            	temp1R     = (Y[upperIdx].real - Y[lowerIdx].real)>>1;
            	temp2R     = (Y[upperIdx].real + Y[lowerIdx].real)>>1;
            	Y[upperIdx].real  =  (short) temp2R;
            	temp1I     = (Y[upperIdx].imag - Y[lowerIdx].imag)>>1;
            	temp2I     = (Y[upperIdx].imag + Y[lowerIdx].imag)>>1;
            	Y[upperIdx].imag  =  (short) temp2I;
            	temp2R     = (c*temp1R - s*temp1I)>>15;
            	Y[lowerIdx].real  = (short) temp2R;
            	temp2I     = (c*temp1I + s*temp1R)>>15;
            	Y[lowerIdx].imag  =  (short) temp2I;
            	}
         	}
         	TwFStep = TwFStep<<1; /* update separation of twiddle factors)*/
       	}

/* bit reversal for resequencing data */

	j=0;
    for (i=1;i<(N-1);i++)
   	{
      k=N/2;
      while (k<=j)
      	{
         j = j-k;
         k=k/2;
         }
      j=j+k;
      if (i<j)
      	{
         tempR=Y[j].real;
         tempI=Y[j].imag;
         Y[j].real=Y[i].real;
         Y[j].imag=Y[i].imag;
         Y[i].real=tempR;
         Y[i].imag=tempI;
         }
      }
  	return;
	}

void fft512(COMPLEX *Y, int N) /* FFT(input sample array, # of points)     */
	{
   	Int32 temp1R, temp1I, temp2R,temp2I;  /* 32 bits temporary storage for */
   									   /* intermediate results			*/					
   	short tempR, tempI, c, s;      /* 16 bits temporary storages    */
   	/* variables */
   	Int32 TwFStep,  /* Step between twiddle factors */
   		  TwFIndex, /* Index of twiddle factors */
   		  BLStep,    /* Step for incrementing butterfly index */
   		  BLdiff,   /* Difference between upper and lower butterfly legs */
   		  upperIdx,
   		  lowerIdx, /* upper and lower indexes of buterfly leg */ 
   		  i, j, k;  /* loop control variables */
   		
   	BLdiff=N;
   	TwFStep=1;
   	for(k=N;k>1;k=(k>>1)) /* Do Log(base 2)(N) Stages */
   		{
      	BLStep=BLdiff;
      	BLdiff=BLdiff>>1; 
      	TwFIndex=0;
      	for(j=0;j<BLdiff;j++)/* Nbr of twiddle factors to use=BLDiff  */
      		{
      		c=w512[TwFIndex].real;
      		s=w512[TwFIndex].imag;
      		TwFIndex=TwFIndex+TwFStep;                 
      		/* Now do N/BLStep butterflies */  
      		for(upperIdx=j;upperIdx<N;upperIdx+=BLStep)
         		{                              
/* Calculations inside this loop avoid overflow by shifting left once
   the result of every adittion/substration and by shifting left 15 
   places the result of every multiplication. Double precision temporary
   results (32-bit) are used in order to avoid losing information because
   of overflow. Final DFT result is scaled by N (number of points), i.e.,
   2^(Nbr of stages) =2^(log(base 2) N) = N								*/
             		
            	lowerIdx=upperIdx+BLdiff;
            	temp1R     = (Y[upperIdx].real - Y[lowerIdx].real)>>1;
            	temp2R     = (Y[upperIdx].real + Y[lowerIdx].real)>>1;
            	Y[upperIdx].real  =  (short) temp2R;
            	temp1I     = (Y[upperIdx].imag - Y[lowerIdx].imag)>>1;
            	temp2I     = (Y[upperIdx].imag + Y[lowerIdx].imag)>>1;
            	Y[upperIdx].imag  =  (short) temp2I;
            	temp2R     = (c*temp1R - s*temp1I)>>15;
            	Y[lowerIdx].real  = (short) temp2R;
            	temp2I     = (c*temp1I + s*temp1R)>>15;
            	Y[lowerIdx].imag  =  (short) temp2I;
            	}
         	}
         	TwFStep = TwFStep<<1; /* update separation of twiddle factors)*/
       	}

/* bit reversal for resequencing data */

	j=0;
   for (i=1;i<(N-1);i++)
   	{
      k=N/2;
      while (k<=j)
      	{
         j = j-k;
         k=k/2;
         }
      j=j+k;
      if (i<j)
      	{
         tempR=Y[j].real;
         tempI=Y[j].imag;
         Y[j].real=Y[i].real;
         Y[j].imag=Y[i].imag;
         Y[i].real=tempR;
         Y[i].imag=tempI;
         }
      }
  	return;
	}

void fft256(COMPLEX *Y, int N) /* FFT(input sample array, # of points)     */
	{
   	Int32 temp1R, temp1I, temp2R,temp2I;  /* 32 bits temporary storage for */
   									   /* intermediate results			*/					
   	short tempR, tempI, c, s;      /* 16 bits temporary storages    */
   	/* variables */
   	Int32 TwFStep,  /* Step between twiddle factors */
   		  TwFIndex, /* Index of twiddle factors */
   		  BLStep,    /* Step for incrementing butterfly index */
   		  BLdiff,   /* Difference between upper and lower butterfly legs */
   		  upperIdx,
   		  lowerIdx, /* upper and lower indexes of buterfly leg */ 
   		  i, j, k;  /* loop control variables */
   		
   	BLdiff=N;
   	TwFStep=1;
   	for(k=N;k>1;k=(k>>1)) /* Do Log(base 2)(N) Stages */
   		{
      	BLStep=BLdiff;
      	BLdiff=BLdiff>>1; 
      	TwFIndex=0;
      	for(j=0;j<BLdiff;j++)/* Nbr of twiddle factors to use=BLDiff  */
      		{
      		c=w256[TwFIndex].real;
      		s=w256[TwFIndex].imag;
      		TwFIndex=TwFIndex+TwFStep;                 
      		/* Now do N/BLStep butterflies */  
      		for(upperIdx=j;upperIdx<N;upperIdx+=BLStep)
         		{                              
/* Calculations inside this loop avoid overflow by shifting left once
   the result of every adittion/substration and by shifting left 15 
   places the result of every multiplication. Double precision temporary
   results (32-bit) are used in order to avoid losing information because
   of overflow. Final DFT result is scaled by N (number of points), i.e.,
   2^(Nbr of stages) =2^(log(base 2) N) = N								*/
             		
            	lowerIdx=upperIdx+BLdiff;
            	temp1R     = (Y[upperIdx].real - Y[lowerIdx].real)>>1;
            	temp2R     = (Y[upperIdx].real + Y[lowerIdx].real)>>1;
            	Y[upperIdx].real  =  (short) temp2R;
            	temp1I     = (Y[upperIdx].imag - Y[lowerIdx].imag)>>1;
            	temp2I     = (Y[upperIdx].imag + Y[lowerIdx].imag)>>1;
            	Y[upperIdx].imag  =  (short) temp2I;
            	temp2R     = (c*temp1R - s*temp1I)>>15;
            	Y[lowerIdx].real  = (short) temp2R;
            	temp2I     = (c*temp1I + s*temp1R)>>15;
            	Y[lowerIdx].imag  =  (short) temp2I;
            	}
         	}
         	TwFStep = TwFStep<<1; /* update separation of twiddle factors)*/
       	}

/* bit reversal for resequencing data */

	j=0;
   for (i=1;i<(N-1);i++)
   	{
      k=N/2;
      while (k<=j)
      	{
         j = j-k;
         k=k/2;
         }
      j=j+k;
      if (i<j)
      	{
         tempR=Y[j].real;
         tempI=Y[j].imag;
         Y[j].real=Y[i].real;
         Y[j].imag=Y[i].imag;
         Y[i].real=tempR;
         Y[i].imag=tempI;
         }
      }
  	return;
	}