fft.c

适合TMS320F2812芯片的嵌入式仪器

#include "DSP28_Device.h"       /*DSP28xx Headerfile Include File*/
#include "DSP281x_Examples.h"     /* DSP28xx Examples Include File*/
#include "fft.h"
#include "math.h"

#define NUM 1024
#define pi 3.1415926
#pragma DATA_SECTION(ipcb, "FFTipcb");

#pragma DATA_SECTION(mag, "FFTmag");/* Prototype statements for functions found within this file.
                                  Global variables used in this example*/

void dsp_fft(float *a, float *b, int N); 
float ipcb[NUM]; /*in place computation buffer*/
float mag[NUM];/*magnitude buffer*/
float a[NUM];
 

main()
{
  Uint16 i;
  for(i=0;i<NUM;i++)
{
   ipcb[i]= sin((2*pi*i)/64);
   mag[i]=0;
}

 dsp_fft(ipcb, mag, NUM);
 for(i=0; i<NUM; i++)
 {
	a[i] = sqrt(ipcb[i]*ipcb[i] + mag[i]*mag[i]);
 }
}

  void dsp_fft(float *a, float *b, int N)      //c为指向输入序列的指针,N为变换点数,inv用来确定是FFT还是IFFT
 {


	long m,l,k,ip;
	long i,j;	                     //循环变量
	float ur,ui,wr,wi;
	int le,lel;                               
	float ctempr,ctempi,tempr,tempi;                  //临时变量
	



/*****************************************************************************************/
//  得到碟形运算级数，m为级数
/*****************************************************************************************/
	long Ntemp=N;
	m=0;	
	while(Ntemp>1)
	{
		Ntemp=Ntemp/2;
		m++;
	} 


/*****************************************************************************************/
//  码位倒置
/*****************************************************************************************/
	for(i=1,j=1;i<N;++i)
	{
		if(i<j)
		{
			ctempr=a[j];
			ctempi=b[j];
			a[j]=a[i];
			b[j]=b[i];
			a[i]=ctempr;
			b[i]=ctempi;
			k=N/2;
			while(k<j)
			{
				j=j-k;
				k=k/2;
			}
			j=j+k;
		}
		else
		{
			k=N/2;
			while(k<j)
			{
				j=j-k;
				k=k/2;
			}
			j=j+k;
		}
	}
	

/*****************************************************************************************/
//  FFT变换
/*****************************************************************************************/
	for(l=1;l<=m;l++)          //m为分级运算的级数，在此N=1024，m=10;
	{
		le=pow(2,l)+0.5;           
		lel=le/2;			   	
		ur=1.0;                //旋转因子的实部，初始化为1
		ui=0.0;                //旋转因子的虚部，初始化为0
		wr=(float)cos(pi/lel);
		wi=(float)(-sin(pi/lel));
		for(j=1;j<=lel;++j)
		{
			i=j;
					
			while(i<N)         //蝶形运算
			{
				ip=i+(int)lel;
				ctempr=a[ip]*ur-b[ip]*ui;
				ctempi=b[ip]*ur+a[ip]*ui;
				a[ip]=a[i]-ctempr;
				b[ip]=b[i]-ctempi;
				a[i]=a[i]+ctempr;
				b[i]=b[i]+ctempi;
				i=i+le;
			}
			tempr=ur*wr-ui*wi;
			tempi=ur*wi+ui*wr;
			ur=tempr;     //同一级中，下两个点运算的旋转因子的实部
			ui=tempi;     //同一级中，下两个点运算的旋转因子的虚部
		}
	}
   
}