FFT算法实现

FFT算法实现（转）
2007年10月17日星期三16:04
/**//***********************************
//DITradix-4FFTcomplex
//
//1.Maximiumpointsare2048becausethat
//functionmylog2(intN)hasthelimitofmaximalpoints
//
//2.Thelaststageis2DFT(for8,32,128,512,2048)
//or4DFT(for4,16,64,256,1024).
//
//
//24juillet2007
//purcharse*gmail.com
************************************/
#include
#include
#include
typedefdoublereal64;/**//*floatingpoint*/
typedeffloatreal32;/**//*32bitfixedpoint*/
structcomplex
{
real32real;
real32imag;
}complex;
staticstructcomplexmulticomplex(structcomplex,structcomplex);
staticintmylog2(int);
staticvoidDFT_2(structcomplex*,structcomplex*);
staticvoidDFT_4(structcomplex*,structcomplex*,structcomplex*,structcomplex*);
staticvoidRunFFT(structcomplex*,int);
staticvoidRunIFFT(structcomplex*,int);
staticvoidBitReverse(structcomplex*,int);
staticvoidFFT_R4(structcomplex*,int,int);
staticvoidFFT_L2(structcomplex*,int);
structcomplexs[257];
intNum=16;
constfloatPI=3.1415926535898;
main()
{
inti;
/**//*rectangle*/
/**//*
for(i=0;i{
s[i].real=0;
s[i].imag=0;
}
s[0].real=10;
s[0].imag=10;
*/
/**//*sinus*/
for(i=0;i{
s[i].real=sin(PI*i/Num);
s[i].imag=cos(PI*i/Num);
}
/**//**/
/**//*
for(i=0;i{
s[i].real=0;
s[i].imag=0;
}
for(i=Num*3/8;i{
s[i].real=(i-Num*3/8);
s[i].imag=0;
}
for(i=Num/2;i{
s[i].real=(Num*5/8-i);
s[i].imag=0;
}*/
printf("***********Donnees*****************\n");
for(i=0;i{
printf("%.4f",s[i].real);
printf("%.4f\n",s[i].imag);
}
RunFFT(s,Num);
printf("***********FFT*****************\n");
for(i=0;i{
printf("%.4f",s[i].real);
printf("%.4f\n",s[i].imag);
}
RunIFFT(s,Num);
printf("***********IFFT*****************\n");
for(i=0;i{
printf("%.4f",s[i].real);
printf("%.4f\n",s[i].imag);
}
}
/**//*****************functions********************/
staticstructcomplexmulticomplex(structcomplexb1,structcomplexb2)/**//*multiplicationofcomplex*/
{
structcomplexb3;
b3.real=b1.real*b2.real-b1.imag*b2.imag;
b3.imag=b1.real*b2.imag+b1.imag*b2.real;
return(b3);
}
staticintmylog2(intN)/**//*Max(N)=4098*/
{
intk=0;
if(N>>12){k+=12;N>>=12;}
if(N>>8){k+=8;N>>=8;}
if(N>>4){k+=4;N>>=4;}
if(N>>2){k+=2;N>>=2;}
if(N>>1){k+=1;N>>=1;}
returnk;
}

staticvoidBitReverse(structcomplex*xin,intN)
{
intLH,i,j,k;
structcomplextmp;
LH=N/2;
j=N/2;
for(i=1;i<=(N-2);i++)
{
if(i{
tmp=xin[j];
xin[j]=xin[i];
xin[i]=tmp;
}
k=LH;
while(j>=k)
{
j=j-k;
k=k/2;
}
j=j+k;
}
}
staticvoidDFT_2(structcomplex*b1,structcomplex*b2)
{
structcomplextmp;
tmp=*b1;
(*b1).real=(*b1).real+(*b2).real;
(*b1).imag=(*b1).imag+(*b2).imag;
(*b2).real=tmp.real-(*b2).real;
(*b2).imag=tmp.imag-(*b2).imag;
}
staticvoidDFT_4(structcomplex*b0,structcomplex*b1,structcomplex*b2,structcomplex*b3)
{
/**//*variableslocales*/
structcomplextemp[4];
/**//*calculx1*/
temp[0].real=(*b0).real+(*b1).real;
temp[0].imag=(*b0).imag+(*b1).imag;
/**//*calculx2*/
temp[1].real=(*b0).real-(*b1).real;
temp[1].imag=(*b0).imag-(*b1).imag;
/**//*calculx3*/
temp[2].real=(*b2).real+(*b3).real;
temp[2].imag=(*b2).imag+(*b3).imag;
/**//*calculx4+multiplicationwith-j*/
temp[3].imag=(*b3).real-(*b2).real;
temp[3].real=(*b2).imag-(*b3).imag;
/**//*thelaststage*/
(*b0).real=temp[0].real+temp[2].real;
(*b0).imag=temp[0].imag+temp[2].imag;
(*b1).real=temp[1].real+temp[3].real;
(*b1).imag=temp[1].imag+temp[3].imag;
(*b2).real=temp[0].real-temp[2].real;
(*b2).imag=temp[0].imag-temp[2].imag;
(*b3).real=temp[1].real-temp[3].real;
(*b3).imag=temp[1].imag-temp[3].imag;
}
staticvoidFFT_R4(structcomplex*xin,intN,intm)
{
inti,L,j;
doubleps1,ps2,ps3,p;
intle,B;
structcomplexw[4];
for(L=1;L<=m;L++)
{
le=pow(4,L);
B=le/4;/**//*thedistanceofbuttefly*/
p=N/le;
for(j=0;j<=B-1;j++)
{
//ps0=(2*pp/N)*0*j;
//w[0].real=cos(ps0);
//w[0].imag=-sin(ps0);
ps1=(2*PI/N)*p*2*j;
w[1].real=cos(ps1);
w[1].imag=-sin(ps1);
ps2=(2*PI/N)*p*1*j;
w[2].real=cos(ps2);
w[2].imag=-sin(ps2);
ps3=(2*PI/N)*p*3*j;
w[3].real=cos(ps3);
w[3].imag=-sin(ps3);
for(i=j;i<=N-1;i=i+le)/**//*controlethosesamebutteflies*/
{
/**//*multiplewithW*/
//xin[i]=multicomplex(xin[i],w[0]);
xin[i+B]=multicomplex(xin[i+B],w[1]);
xin[i+2*B]=multicomplex(xin[i+2*B],w[2]);
xin[i+3*B]=multicomplex(xin[i+3*B],w[3]);
/**//*DFT-4*/
DFT_4(xin+i,xin+i+B,xin+i+2*B,xin+i+3*B);
}
}
/**//*
printf("*****N°%d**********",L);
for(i=0;i{
printf("%.8f",xin[i].real);
printf("%.8f",xin[i].imag);
}
*/
}
}//finduFFT_R4
staticvoidFFT_L2(structcomplex*xin,intN)
{/**//*Forthelaststage2DFT*/
intj,B;
doublep,ps;
structcomplexw;
B=N/2;
for(j=0;j<=B-1;j++)
{
ps=(2*PI/N)*j;
w.real=cos(ps);
w.imag=-sin(ps);
/**//*multipleavecW*/
xin[j+B]=multicomplex(xin[j+B],w);
DFT_2(xin+j,xin+j+B);
}
}//finduFFT_L2
staticvoidRunFFT(structcomplex*xin,intN)
{
intm,i;
BitReverse(xin,N);
m=mylog2(N);
if((m%2)0)
{
/**//*Allthestagesareradix4*/
FFT_R4(xin,N,m/2);
}
else
{
/**//*thelaststageisradix2*/
FFT_R4(xin,N,m/2);
FFT_L2(xin,N);
}
}
staticvoidRunIFFT(structcomplex*xin,intN)
{/**//*inverseFFT*/
inti;
for(i=0;i{
xin[i].imag=-xin[i].imag;
}
RunFFT(xin,N);
for(i=0;i{
xin[i].real=xin[i].real/Num;
xin[i].imag=-xin[i].imag/Num;
}