FFT/IFFT的完整VC实现

#include<iostream.h> #include<math.h> #include<stdio.h> #include<fstream.h> #include<memory.h> #define PI 3.1415926 #define LENGTH 16 /*--------------------------------*/ struct Complex { double real; double imag; }; /*--------------------------------*/ Complex Data_In1[LENGTH]; double Data_In2[LENGTH*2]; Complex FFT_Out1[LENGTH]; Complex FFT_Out2[LENGTH]; Complex FFT_Out3[LENGTH]; float DataI[LENGTH]; float DataQ[LENGTH]; float FFTMag_Out[LENGTH]; float FFTShift_Out[LENGTH]; bool FFT(Complex *data_in, Complex *fft_out, int n); bool FFT(double *data_in, Complex *fft_out, int n); bool FFT(float *datar_in, float *datai_in,Complex *fft_out, int datan,int fftn,bool flag); int InvBits(int B,int Bn); bool FFTMag(Complex *fft_in,float *fftmag_out,int fftn); bool FFTShift(float *fft_in,float *fftshift_out,int fftn); void AddWindow(float *datar_in,float *datai_in,int n,int flag);//加窗，flag选取加什么窗 void main() { for(int i=1;i<=8;i++) { DataI[i-1]=(float)i; DataQ[i-1]=(float)i; } //AddWindow(DataI,DataQ,10,1); FFT(DataI,DataQ,FFT_Out3,8,LENGTH,true); cout<<"FFT3的结果如下:\n"; for(i=0;i<LENGTH;i++) { cout<<FFT_Out3[i].real<<"+j"<<FFT_Out3[i].imag<<"\n"; } FFTMag(FFT_Out3,FFTMag_Out,LENGTH); FFTShift(FFTMag_Out,FFTMag_Out,LENGTH); } int InvBits(int B,int Bn) { int temp; temp=0; for(int i=0;i<Bn;i++) { temp=(temp<<1)|(B&1); B=B>>1; } return temp; } /*--------------------------------*/ //FFT算法 //实部与虚部交错存放 bool FFT(double *data_in, Complex *fft_out, int n) { int State,Group,Butterfly;//计数器 int State_Num;//总的级数 int Data_Space;//数据间隔 int Group_Num;//一级的蝶群数 int Butterfly_Num;//一个蝶群里的蝶数 int Rotate_Space;//旋转因子间隔 int Rotate_Addr;//取旋转因子的地址 int Length_Temp=n; for(State_Num=0;Length_Temp>=2;State_Num++,Length_Temp/=2);//根据长度求级数 //位反转 for(int i=0;i<n;i++) { int Addr_Temp; Addr_Temp=InvBits(i,State_Num); fft_out[i].real=data_in[Addr_Temp*2]; fft_out[i].imag=data_in[Addr_Temp*2+1]; } //FFT Data_Space=1; Group_Num=n/2; Butterfly_Num=1; Rotate_Space=n/2; Rotate_Addr=0; for(State=0;State<State_Num;State++) { for(Group=0;Group<Group_Num;Group++) { Rotate_Addr=0; for(Butterfly=0;Butterfly<Butterfly_Num;Butterfly++) { double Ror,Roi; Ror=cos(2*PI*Rotate_Addr/n); Roi=-sin(2*PI*Rotate_Addr/n); int x_addr=Butterfly+Group*Butterfly_Num*2; int y_addr=Butterfly+Group*Butterfly_Num*2+Data_Space; double tempr,tempi; tempr=Ror*fft_out[y_addr].real-Roi*fft_out[y_addr].imag; tempi=Ror*fft_out[y_addr].imag+Roi*fft_out[y_addr].real; fft_out[y_addr].real=fft_out[x_addr].real-tempr; fft_out[y_addr].imag=fft_out[x_addr].imag-tempi; fft_out[x_addr].real+=tempr; fft_out[x_addr].imag+=tempi; Rotate_Addr+=Rotate_Space; Rotate_Addr=Rotate_Addr%(n/2); } } Group_Num=Group_Num/2; Data_Space=Data_Space*2; Butterfly_Num=Butterfly_Num*2; Rotate_Space/=2; } return true; } /*--------------------------------*/ /*--------------------------------*/ //FFT算法 //实部与虚部分开存放 bool FFT(Complex *data_in,Complex *fft_out, int n) { int State,Group,Butterfly;//计数器 int State_Num;//总的级数 int Data_Space;//数据间隔 int Group_Num;//一级的蝶群数 int Butterfly_Num;//一个蝶群里的蝶数 int Rotate_Space;//旋转因子间隔 int Rotate_Addr;//取旋转因子的地址 int Length_Temp=n; for(State_Num=0;Length_Temp>=2;State_Num++,Length_Temp/=2);//根据长度求级数 //位反转 for(int i=0;i<n;i++) { int Addr_Temp; Addr_Temp=InvBits(i,State_Num); fft_out[i].real=data_in[Addr_Temp].real; fft_out[i].imag=data_in[Addr_Temp].imag; } //FFT Data_Space=1; Group_Num=LENGTH/2; Butterfly_Num=1; Rotate_Space=LENGTH/2; Rotate_Addr=0; for(State=0;State<State_Num;State++) { for(Group=0;Group<Group_Num;Group++) { Rotate_Addr=0; for(Butterfly=0;Butterfly<Butterfly_Num;Butterfly++) { double Ror,Roi; Ror=cos(2*PI*Rotate_Addr/n); Roi=-sin(2*PI*Rotate_Addr/n); int x_addr=Butterfly+Group*Butterfly_Num*2; int y_addr=Butterfly+Group*Butterfly_Num*2+Data_Space; double tempr,tempi; tempr=Ror*fft_out[y_addr].real-Roi*fft_out[y_addr].imag; tempi=Ror*fft_out[y_addr].imag+Roi*fft_out[y_addr].real; fft_out[y_addr].real=fft_out[x_addr].real-tempr; fft_out[y_addr].imag=fft_out[x_addr].imag-tempi; fft_out[x_addr].real+=tempr; fft_out[x_addr].imag+=tempi; Rotate_Addr+=Rotate_Space; Rotate_Addr=Rotate_Addr%(n/2); } } Group_Num=Group_Num/2; Data_Space=Data_Space*2; Butterfly_Num=Butterfly_Num*2; Rotate_Space/=2; } return true; } /*--------------------------------*/ bool FFT(float *datar_in, float *datai_in,Complex *fft_out, int datan,int fftn,bool flag) { int i; if(datan>fftn) return false; int State,Group,Butterfly;//计数器 int State_Num;//总的级数 int Data_Space;//数据间隔 int Group_Num;//一级的蝶群数 int Butterfly_Num;//一个蝶群里的蝶数 int Rotate_Space;//旋转因子间隔 int Rotate_Addr;//取旋转因子的地址 int Length_Temp=fftn; for(State_Num=0;Length_Temp>=2;State_Num++,Length_Temp/=2);//根据长度求级数 int Factor; if(flag) Factor=1; else Factor=fftn; //补零 float *datai_temp=new float[fftn]; float *datar_temp=new float[fftn]; memcpy(datai_temp,datar_in,datan*sizeof(float)); memcpy(datar_temp,datai_in,datan*sizeof(float)); for(i=datan;i<fftn;i++)//如果输入数据长度不够，则补0 { datai_temp[i]=0.0; datar_temp[i]=0.0; } //位反转 for(i=0;i<fftn;i++) { int Addr_Temp; Addr_Temp=InvBits(i,State_Num); fft_out[i].real=datar_temp[Addr_Temp]; fft_out[i].imag=datai_temp[Addr_Temp]; } //FFT Data_Space=1; Group_Num=fftn/2; Butterfly_Num=1; Rotate_Space=fftn/2; Rotate_Addr=0; for(State=0;State<State_Num;State++) { for(Group=0;Group<Group_Num;Group++) { Rotate_Addr=0; for(Butterfly=0;Butterfly<Butterfly_Num;Butterfly++) { float Ror,Roi; Ror=(float)cos(2*PI*Rotate_Addr/fftn); if(flag) Roi=(float)-sin(2*PI*Rotate_Addr/fftn); else Roi=(float)sin(2*PI*Rotate_Addr/fftn); int x_addr=Butterfly+Group*Butterfly_Num*2; int y_addr=Butterfly+Group*Butterfly_Num*2+Data_Space; float tempr,tempi; tempr=Ror*fft_out[y_addr].real-Roi*fft_out[y_addr].imag; tempi=Ror*fft_out[y_addr].imag+Roi*fft_out[y_addr].real; fft_out[y_addr].real=(fft_out[x_addr].real-tempr)/Factor; fft_out[y_addr].imag=(fft_out[x_addr].imag-tempi)/Factor; fft_out[x_addr].real=(fft_out[x_addr].real+tempr)/Factor; fft_out[x_addr].imag=(fft_out[x_addr].imag+tempi)/Factor; Rotate_Addr+=Rotate_Space; Rotate_Addr=Rotate_Addr%(fftn/2); } } Group_Num=Group_Num/2; Data_Space=Data_Space*2; Butterfly_Num=Butterfly_Num*2; Rotate_Space/=2; } delete[] datai_temp; delete[] datar_temp; return true; } bool FFTMag(Complex *fft_in,float *fftmag_out,int fftn) { if(fftn<=0) return false; for(int i=0;i<fftn;i++) { fftmag_out[i]=(float)sqrt(fft_in[i].real*fft_in[i].real+fft_in[i].imag*fft_in[i].imag); } return true; } bool FFTShift(float *fft_in,float *fftshift_out,int fftn) { if(fftn<=0) return false; if(fft_in!=fftshift_out) { for(int i=0;i<fftn/2-1;i++) { fftshift_out[i+fftn/2]=fft_in[i]; fftshift_out[i]=fft_in[i+fftn/2]; } } else { float *temp=new float[fftn/2]; memcpy(temp,fft_in,fftn*sizeof(float)/2); for(int i=0;i<fftn/2;i++) { fftshift_out[i]=fft_in[i+fftn/2]; fftshift_out[i+fftn/2]=temp[i]; } delete[] temp; } return true; } void AddWindow(float *datar_in,float *datai_in,int n,int flag) { int i; float *wn=new float[n]; int m=n-1; switch(flag) { case 1://blackman {