用C语言编写FFT算法

用simulink建模的方法实现FFT的问题特别多，还不如手工编写，也不是很复杂。

以下是512点的单边谱FFT算法，在TMS320C6713的DSP开发板上调试通过。

顺便问一句，有人会编welch功率谱密度的C语言实现程序么？

#include \"math.h\"

#define PI 3.1415926

#define SAMPLENUMBER 512

void InitForFFT();

void MakeWave();

void FFT();

int INPUT[SAMPLENUMBER],DATA[SAMPLENUMBER];

float

fWaveR[SAMPLENUMBER],fWaveI[SAMPLENUMBER],w[SAMPLENUMBER];

float sin_tab[SAMPLENUMBER],cos_tab[SAMPLENUMBER];

main()

{

int i;

InitForFFT();

MakeWave();

for ( i=0;i

{

fWaveR[i]=INPUT[i];

fWaveI[i]=0.0f;

w[i]=0.0f;

}

FFT(fWaveR,fWaveI);

for ( i=0;i

— 精选文库

2

{

DATA[i]=w[i];

}

while ( 1 ); // break point

}

void FFT(float dataR[SAMPLENUMBER],float dataI[SAMPLENUMBER])

{

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;i

{

x0=x1=x2=x3=x4=x5=x6=x7=x8=0;

— 精选文库

3

精选文库

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;

dataI[xx]=dataR[i];

}

for ( i=0;i

{

dataR[i]=dataI[i]; dataI[i]=0;

}

for ( L=1;L<=9;L++ )

{

b=1; i=L-1;

while ( i>0 )

{

—

4

b=b*2; i--;

}

for ( j=0;j<=b-1;j++ )

{

p=1; i=9-L;

while ( i>0 )

{

p=p*2; i--;

}

p=p*j;

for ( k=j;k<512;k=k+2*b )

{

TR=dataR[k]; TI=dataI[k]; temp=dataR[k+b];

— 精选文库

5

dataR[k]=dataR[k]+dataR[k+b]*cos_tab[p]+dataI[k+b]*sin_tab[p];

dataI[k]=dataI[k]-dataR[k+b]*sin_tab[p]+dataI[k+b]*cos_tab[p];

dataR[k+b]=TR-dataR[k+b]*cos_tab[p]-dataI[k+b]*sin_tab[p];

dataI[k+b]=TI+temp*sin_tab[p]-dataI[k+b]*cos_tab[p];

}

}

}

for ( i=0;i

{

w[i]=sqrt(dataR[i]*dataR[i]+dataI[i]*dataI[i]);

}

}

void InitForFFT()

— 精选文库

6

{

int i;

for ( i=0;i

{

sin_tab[i]=sin(PI*2*i/SAMPLENUMBER);

cos_tab[i]=cos(PI*2*i/SAMPLENUMBER);

}

}

void MakeWave()

{

int i;

for ( i=0;i

{

— 精选文库

7

精选文库

INPUT[i]=sin(PI*2*i/SAMPLENUMBER*30)*1024;

}

}

— 8