用C語(yǔ)言實(shí)現(xiàn)FFT算法
用C語(yǔ)言實(shí)現(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 ;??? ??? /*變址運(yùn)算*/??? for(i=1;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??????? {??????????? 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??????????? {??????????????? /*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??????????????? {??????????????????? /*? 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??? {??????? s[i].real=sin(pp*i/32);??????? s[i].imag=0 ;??? }??? ??? FFT(s,Num);??? ??? for(i=1;i??? {??????? result[i]=sqrt(pow(s[i].real,2)+pow(s[i].imag,2));??? }}