快速傅里叶变换的基2FFT算法的C++实现

[转自 http://hi.baidu.com/kxw102/blog/item/8672eaeddf54532e279791a4.html]

快速傅里叶变换的基2FFT算法的C++实现

2011-01-19 05:26

快速傅里叶变换的基本原理由于公式不好显示请读者参考其它文章或书籍，本文重点给出了时域抽取法FFT的C++实现。 下面是算法的流程图 倒序的流程图 C++实现代码： 1. FFT.h #pragma once #ifndef FFT_H #define FFT_H #include #include using namespace std; #define pi 3.141592 class FFT { private: void changeOrder(double* xr,double* xi,int n); public: FFT(void); void FFT_1D(double* ctxr,double* ctxi,double* cfxr,double* cfxi,int len); void rFFT_1D(double* ctxr,double* cfxr,double* cfxi,int len); void IFFT_1D(double* cfxr,double* cfxi,double* ctxr,double* ctxi,int len); void rIFFT_1D(double* cfxr,double* cfxi,double* ctxr,int len); ~FFT(void); }; #endif 2.FFT.cpp #include "./fft.h" FFT::FFT() { } /// //倒序实现 //xr实部，xi虚部，n为2的幂 /// void FFT::changeOrder(double *xr,double *xi,int n) { double temp; int k; int lh=n/2; int j=lh; int n1=n-2; for(int i=1;i<=n1;i++) { if(i { temp=xr[i]; xr[i]=xr[j]; xr[j]=temp; temp=xi[i]; xi[i]=xi[j]; xi[j]=temp; } k=lh; while(j>=k) { j=j-k; k=(int)(k/2+0.5); } j=j+k; } } // //复数FFT //ctxr和ctxi的长度为len //cfxr和cfxi的长度为2的幂 // void FFT::FFT_1D(double *ctxr,double *ctxi,double *cfxr,double *cfxi,int len) { int m=ceil(log((double)len)/log(2.0)); int l,b,j,p,k; double rkb,ikb; int n=1< double* rcos=new double[n/2]; double* isin=new double[n/2]; for(l=0;l { rcos[l]=cos(l*pi*2/n); isin[l]=sin(l*pi*2/n); } memcpy(cfxr,ctxr,sizeof(double)*len); memcpy(cfxi,ctxi,sizeof(double)*len); for(l=len;l { cfxr[l]=0; cfxi[l]=0; } changeOrder(cfxr,cfxi,n); //倒序 for(l=1;l<=m;l++) { b=(int)(pow(2,l-1)+0.5); for(j=0;j { p=j*(int)(pow(2,m-l)+0.5); for(k=j;k { rkb=cfxr[k+b]*rcos[p]+cfxi[k+b]*isin[p]; ikb=cfxi[k+b]*rcos[p]-cfxr[k+b]*isin[p]; cfxr[k+b]=cfxr[k]-rkb; cfxi[k+b]=cfxi[k]-ikb; cfxr[k]=cfxr[k]+rkb; cfxi[k]=cfxi[k]+ikb; } } } delete []rcos; delete []isin; } / //实数FFT //ctxr的长度为len //cfxr和cfxi的长度为2的幂 void FFT::rFFT_1D(double *ctxr,double *cfxr,double *cfxi,int len) { int m=ceil(log((double)len)/log(2.0)); int n=1< double* rcos=new double[n/2]; double* isin=new double[n/2]; for(int l=0;l { rcos[l]=cos(l*pi*2/n); isin[l]=sin(l*pi*2/n); } double* txr=new double[(len+1)/2]; double* txi=new double[(len+1)/2]; for(int i=0;i { txr[i]=ctxr[i*2]; txi[i]=ctxr[i*2+1]; } if(len%2==1) { txr[(len-1)/2]=ctxr[len-1]; txi[(len-1)/2]=0; } FFT_1D(txr,txi,cfxr,cfxi,(len+1)/2); double* x1r=new double[n/2]; double* x1i=new double[n/2]; double* x2r=new double[n/2]; double* x2i=new double[n/2]; x1r[0]=cfxr[0]; x1i[0]=0; x2r[0]=cfxi[0]; x2i[0]=0; for(int k=1;k { x1r[k]=(cfxr[k]+cfxr[n/2-k])/2.0; x1i[k]=(cfxi[k]-cfxi[n/2-k])/2.0; x2r[k]=(cfxi[k]+cfxi[n/2-k])/2.0; x2i[k]=(-cfxr[k]+cfxr[n/2-k])/2.0; } double rkb,ikb; for(i=0;i { rkb=x2r[i]*rcos[i]+x2i[i]*isin[i]; ikb=x2i[i]*rcos[i]-x2r[i]*isin[i]; cfxr[i+n/2]=x1r[i]-rkb; cfxi[i+n/2]=x1i[i]-ikb; cfxr[i]=x1r[i]+rkb; cfxi[i]=x1i[i]+ikb; } delete []txr; delete []txi; delete []x1r; delete []x1i; delete []x2r; delete []x2i; delete []rcos; delete []isin; } /// //复数IFFT //cfxr和cfxi的长度为n(2的幂) //ctxr和ctxi的长度为len /// void FFT::IFFT_1D(double *cfxr,double *cfxi,double *ctxr,double *ctxi,int len) { int m=ceil(log((double)len)/log(2.0)); int n=1< double* txr=new double[n]; double* txi=new double[n]; for(int i=0;i cfxi[i]=-cfxi[i]; FFT_1D(cfxr,cfxi,txr,txi,n); for(i=0;i { ctxr[i]=txr[i]/n; ctxi[i]=-txi[i]/n; } delete []txr; delete []txi; } // //实数IFFT //cfxr和cfxi的长度为n(2的幂) //ctxr的长度为len // void FFT::rIFFT_1D(double *cfxr,double *cfxi,double *ctxr,int len) { int m=ceil(log((double)len)/log(2.0)); int n=1< double* txr=new double[n]; double* txi=new double[n]; for(int i=0;i cfxi[i]=-cfxi[i]; FFT_1D(cfxr,cfxi,txr,txi,n); for(int i=0;i { ctxr[i]=txr[i]/n; } delete []txr; delete []txi; } FFT::~FFT(void) { } 3.test.cpp #include #include "FFT.h" using namespace std; void main() { double xr[10]={1,2,3,4,5,6,7,8,9,10}; //实部 double xi[10]={0,0,0,0,0,0,0,0,0,0}; //虚部 double *cxr,*cxi; cxr=new double[16]; cxi=new double[16]; FFT f; f.rFFT_1D(xr,cxr,cxi,10); for(int i=0;i<16;i++) { cout< <<"+j"< cout< } cout< double *fxr,*fxi; fxr=new double[10]; fxi=new double[10]; f.rIFFT_1D(cxr,cxi,fxr,10); for(i=0;i<10;i++) { cout< cout< } delete []fxr; delete []fxi; delete []cxr; delete []cxi; }