用C語言編寫FFT算法

用C語言編寫FFT算法用simulink建模的方法實現(xiàn)FFT的問題特別多，還不如手工編寫，也不是很復(fù)雜。以下是512點的單邊譜FFT算法，在TMS320C6713的DSP開發(fā)板上調(diào)試通過。順便問一句，有人會編welch功率譜密度的C語言實現(xiàn)程序么？#include math.h#define PI 3.1415926#define SAMPLENUMBER 512void InitForFFT();void MakeWave();void FFT();int INPUTSAMPLENUMBER,DATASAMPLENUMBER;float fWaveRSAMPLENUMBER,fWaveISAMPLENUMBER,wSAMPLENUMBER;float sin_tabSAMPLENUMBER,cos_tabSAMPLENUMBER;main()int i;InitForFFT();MakeWave();for ( i=0;ifWaveRi=INPUTi;fWaveIi=0.0f;wi=0.0f;FFT(fWaveR,fWaveI);for ( i=0;iDATAi=wi;while ( 1 ); / break pointvoid FFT(float dataRSAMPLENUMBER,float dataISAMPLENUMBER)int x0,x1,x2,x3,x4,x5,x6,x7,x8,xx;int i,j,k,b,p,L;float TR,TI,temp;for ( i=0;ix0=x1=x2=x3=x4=x5=x6=x7=x8=0;x0=i&0x01; x1=(i/2)&0x01; x2=(i/4)&0x01; x3=(i/8)&0x01;x4=(i/16)&0x01; x5=(i/32)&0x01; x6=(i/64)&0x01;x7=(i/128)&0x01;x8=(i/256)&0x01;xx=x0*256+x1*128+x2*64+x3*32+x4*16+x5*8+x6*4+x7*2+x8;dataIxx=dataRi;for ( i=0;idataRi=dataIi; dataIi=0;for ( L=1;L0 )b=b*2; i-;for ( j=0;j0 )p=p*2; i-;p=p*j;for ( k=j;k512;k=k+2*b )TR=dataRk; TI=dataIk; temp=dataRk+b;dataRk=dataRk+dataRk+b*cos_tabp+dataIk+b*sin_tabp;dataIk=dataIk-dataRk+b*sin_tabp+dataIk+b*cos_tabp;dataRk+b=TR-dataRk+b*cos_tabp-dataIk+b*sin_tabp;dataIk+b=TI+temp*sin_tabp-dataIk+b*cos_tabp;for ( i=0;iwi=sqrt(dataRi*dataRi+dataIi*dataIi);void InitForFFT()int i;for ( i=0;isin_tabi=sin(PI*2*i/SAMPLENUMBER);cos_tabi=cos(PI*2*i/SAMPLENUMBER);