用C語言實現(xiàn)FFT算法

/*****************fft programe*********************/
#include "typedef.h"
#include "math.h"

struct compx EE(struct compx b1,struct compx b2)
{
??? struct compx b3 ;
??? b3.real=b1.real*b2.real-b1.imag*b2.imag ;
??? b3.imag=b1.real*b2.imag+b1.imag*b2.real ;
??? return(b3);
}

void FFT(struct compx*xin,int N)
{
??? int f,m,nv2,nm1,i,k,j=1,l ;
??? /*int f,m,nv2,nm1,i,k,j=N/2,l;*/
??? struct compx v,w,t ;
??? nv2=N/2 ;
??? f=N ;
??? for(m=1;(f=f/2)!=1;m++)
??? {
??????? ;
??? }
??? nm1=N-1 ;
???
??? /*變址運算*/
??? for(i=1;i<=nm1;i++)
??? {
??????? if(i??????? {
??????????? t=xin[j];
??????????? xin[j]=xin[i];
??????????? xin[i]=t ;
??????? }
??????? k=nv2 ;
??????? while(k??????? {
??????????? j=j-k ;
??????????? k=k/2 ;
??????? }
??????? j=j+k ;
??? }
???
??? {
??????? int le,lei,ip ;
??????? float pi ;
??????? for(l=1;l<=m;l++)
??????? {
??????????? le=pow(2,l);
??????????? // 這里用的是L而不是1? !!!!
??????????? lei=le/2 ;
??????????? pi=3.14159 ;
??????????? v.real=1.0 ;
??????????? v.imag=0.0 ;
??????????? w.real=cos(pi/lei);
??????????? w.imag=-sin(pi/lei);
??????????? for(j=1;j<=lei;j++)
??????????? {
??????????????? /*double p=pow(2,m-l)*j;
????????????????? double ps=2*pi/N*p;
????????????????? w.real=cos(ps);
????????????????? w.imag=-sin(ps);*/
??????????????? for(i=j;i<=N;i=i+le)
??????????????? {
??????????????????? /*? w.real=cos(ps);
????????????????????? w.imag=-sin(ps);*/
??????????????????? ip=i+lei ;
??????????????????? t=EE(xin[ip],v);
??????????????????? xin[ip].real=xin[i].real-t.real ;
??????????????????? xin[ip].imag=xin[i].imag-t.imag ;
??????????????????? xin[i].real=xin[i].real+t.real ;
??????????????????? xin[i].imag=xin[i].imag+t.imag ;
??????????????? }
??????????????? v=EE(v,w);
??????????? }
??????? }
??? }
??? return ;
}

/*****************main programe********************/

#include
#include
#include
#include "typedef.h"

float result[257];
struct compx s[257];
int Num=256 ;
const float pp=3.14159 ;

main()
{
??? int i=1 ;
??? for(;i<0x101;i++)
??? {
??????? s[i].real=sin(pp*i/32);
??????? s[i].imag=0 ;
??? }
???
??? FFT(s,Num);
???
??? for(i=1;i<0x101;i++)
??? {
??????? result[i]=sqrt(pow(s[i].real,2)+pow(s[i].imag,2));
??? }
}