FFT算法的c++实现
- 格式：doc
- 大小：253.00 KB
- 文档页数：9
FFT算法源代码（C++）#include <iostream>#include <iomanip>#include "math.h"using namespace std;double pi = 3.1415926535897932;class complex{public://无参构造函数complex(){re=0;im=0;}//有参构造函数complex(double real,double imag){re=real;im=imag;}//加法complex operator + (complex& c){return complex( re + c.re , im + c.im );}//减法complex operator - (complex& c){return complex( re - c.re , im - c.im );}//乘法complex operator * (complex& c){return complex( (re * c.re)-(im * c.im) , (re * c.im)+(im * c.re) );}//除法complex operator / (complex& c){return complex( ( re*c.re + im*c.im )/( c.re*c.re + c.im*c.im ),((im * c.re)-(re * c.im))/((c.re*c.re)+(c.im*c.im)) );}//除2void half(){re=re/2;im=im/2;}//输出到ostreamvoid insert(ostream& out){if(re>=0) cout<<" ";out<<re<< ((im>=0)?'+':'-') <<'j'<< ((im>=0)?im:(0-im));}//显示void show(){cout<<re<< ((im>=0)?'+':'-') <<'j'<< ((im>=0)?im:(0-im)) ;}//设值void setValue(double real,double imag){if(real>0.00000001||real< -0.00000001) re=real;else re=0;if(imag>0.00000001||imag< -0.00000001) im=imag;else im=0;}friend complex c_pow(complex c,int y);private:double re,im;};//流输出ostream& operator << (ostream& out, complex c){c.insert(out);return out;}//乘方complex c_pow(complex c,int y){complex returnValue = c ;for(int i = 1 ; i < y ; i++ ){returnValue = returnValue * c ;}return returnValue;}//获得Wnvoid getW(complex w[],int len){for(int i=0 ; i<len ; i++ ){w[i].setValue( cos(0 - pi*2*(i)/len) , sin(0 - pi*2*(i)/len) );}}//获得mWnvoid getMW(complex w[],int len){for(int i=0 ; i<len ; i++ ){w[i].setValue( cos(pi*2*(i)/len) , sin(pi*2*(i)/len) );}}//倒序重排void resort(complex x[],int len){//初始化int half=len/2;int p1=0,p2;int k;complex temp;//重排for( p2=0 ; p2+1<len ; p2++ ){if(p2<p1 && p1+1<len){temp=x[p1];x[p1]=x[p2];x[p2]=temp;}k=half;while(k<p1+1 && p1+1<len){p1=p1-k;k=k/2;}p1=p1+k;}}//基2时间FFTvoid FFT2t(complex x[],complex y[],int len){//初始化int m=0,n,k,j,q;int tlen=len/2;while(tlen!=0){m++;tlen=tlen/2;}resort(x,len); //倒序重排for(j=0;j<len;j++){y[j]=x[j];}complex *xt = new complex[len];complex *w = new complex[len];getW(w,len);//开始迭代运算q=len;for(int i=1;i<=m;i++){q=q/2;n=len/q/2;for(j=0;j<len;j++){xt[j]=y[j];}if(i!=1){ for(j=0;j<q;j++){for(k=0;k<n;k++){xt[n+j*2*n+k]=xt[n+j*2*n+k]*w[k*q];}}for(j=0;j<n;j++){for(k=0;k<q;k++){y[j+k*2*n]=xt[j+k*2*n]+xt[j+n+k*2*n];}}for(;j<2*n;j++){for(k=0;k<q;k++){y[j+k*2*n]=xt[j-n+k*2*n]-xt[j+k*2*n];}}}delete[]xt,w;}//基2频率FFT和IFFTvoid FFT(complex x[],complex y[],int len,bool IFFT = false) {//初始化int m=0,n,k,j,q;int tlen=len/2;while(tlen!=0){m++;tlen=tlen/2;}for(j=0;j<len;j++){y[j]=x[j];}complex *xt = new complex[len];complex *w = new complex[len];if(IFFT==true) //若是IFFT则获取Wn序列{getMW(w,len);}else //若是FFT则获取mWn序列{getW(w,len);}//开始迭代运算for(int i=1;i<=m;i++){n=n/2;q=len/n/2;for(j=0;j<len;j++){xt[j]=y[j];}for(j=0;j<n;j++){for(k=0;k<q;k++){y[j+k*2*n]=xt[j+k*2*n]+xt[j+n+k*2*n];}}for(;j<2*n;j++){for(k=0;k<q;k++){y[j+k*2*n]=xt[j-n+k*2*n]-xt[j+k*2*n];}}if(i!=m){ for(j=0;j<q;j++){for(k=0;k<n;k++){y[n+j*2*n+k]=y[n+j*2*n+k]*w[k*q];}}}if(IFFT==true)//若是IFFT则除以2{for(j=0;j<len;j++){y[j].half();}}}resort(y,len); //倒序重排delete[]xt,w;}void main(){//初始化complex q(0.9,0.3);complex x[32],y[32],z[32],w[32],v[32],w1;w1.setValue(1,0);getW(w,32); //获取Wn序列x[0].setValue(1,0);for(int i=1;i<32;i++)//定原序列的值{x[i]=x[i-1]*q;}//输出原序列cout<<"原序列："<<endl;for(i=0;i<32;i+=2){cout<<setiosflags(ios::fixed)<<setprecision(6)<<x[i]<<"\t"<<x[i+1]<<endl;}for(i=0;i<32;i++) //进行公式计算{complex a=w1-(x[31]*q),b=w1-q*w[i];v[i]=a/b;}//输出公式计算值cout<<"公式计算值："<<endl;for(i=0;i<32;i+=2){cout<<v[i]<<"\t"<<v[i+1]<<endl;}FFT2t(x,y,32); //进行基2时间变换//输出基2时间变换后序列cout<<"基2时间变换后序列："<<endl;for(i=0;i<32;i+=2){cout<<y[i]<<"\t"<<y[i+1]<<endl;}FFT(y,z,32,1); //进行反变换//输出反变换后的序列cout<<"反变换后的序列："<<endl;for(i=0;i<32;i+=2)cout<<z[i]<<"\t"<<z[i+1]<<endl;}}运行结果。
