快速傅里葉變換的基2FFT算法的C實現(xiàn)

1、快速傅里葉變換的基2fft算法的c 實現(xiàn)快速傅里葉變換的基2fft算法的c+實現(xiàn)2011-01-19 05：26快速傅里葉變換的基本原理由于公式不好顯示請讀者參考其它文章或書籍，本文重點給出了時域抽取法fft的c+實現(xiàn)。下面是算法的流程圖倒序的流程圖c+實現(xiàn)代碼：1.fft.h#pragma once#ifndef fft_h#define fft_h#include cstring#include cmath using namespace std；#define pi 3.141592 class fftprivate：void changeorder(double*xr,double*x

2、i,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

3、.fft.cpp#include.fft.hfft：fft()/倒序?qū)崿F(xiàn)/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 j)temp=xri；xri=xrj；xrj=temp；temp=xii；xii=xij；xij=temp；k=lh；while(j=k)j=j-k；k=(int)(k/2+0.5)；j=j+k；/復(fù)數(shù)fft/ctxr和ctxi的長度為len/cfxr

4、和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 m；double*rcos=new doublen/2；double*isin=new doublen/2；for(l=0；l n/2；l+)rcosl=cos(l*pi*2/n)；isinl=sin(l*pi*2/n)；memcpy(cfxr,ctxr,sizeof(double)

5、*len)；memcpy(cfxi,ctxi,sizeof(double)*len)；for(l=len；l n；l+)cfxrl=0；cfxil=0；changeorder(cfxr,cfxi,n)；/倒序for(l=1；l=m；l+)b=(int)(pow(2,l-1)+0.5)；for(j=0；j b；j+)p=j*(int)(pow(2,m-l)+0.5)；for(k=j；k n；k+=(int)(pow(2,l)+0.5)rkb=cfxrk+b*rcos+cfxik+b*isin；ikb=cfxik+b*rcos-cfxrk+b*isin；cfxrk+b=cfxrk-rkb；cfxi

6、k+b=cfxik-ikb；cfxrk=cfxrk+rkb；cfxik=cfxik+ikb；delete rcos；delete isin；/實數(shù)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 m；double*rcos=new doublen/2；double*isin=new doublen/2；for(int l=0；l n/2；l+)rcosl=cos(l*pi

7、*2/n)；isinl=sin(l*pi*2/n)；double*txr=new double(len+1)/2；double*txi=new double(len+1)/2；for(int i=0；i len/2；i+)txri=ctxri*2；txii=ctxri*2+1；if(len%2=1)txr(len-1)/2=ctxrlen-1；txi(len-1)/2=0；fft_1d(txr,txi,cfxr,cfxi,(len+1)/2)；double*x1r=new doublen/2；double*x1i=new doublen/2；double*x2r=new doublen/2；d

8、ouble*x2i=new doublen/2；x1r0=cfxr0；x1i0=0；x2r0=cfxi0；x2i0=0；for(int k=1；k n/2；k+)x1rk=(cfxrk+cfxrn/2-k)/2.0；x1ik=(cfxik-cfxin/2-k)/2.0；x2rk=(cfxik+cfxin/2-k)/2.0；x2ik=(-cfxrk+cfxrn/2-k)/2.0；double rkb,ikb；for(i=0；i n/2；i+)rkb=x2ri*rcosi+x2ii*isini；ikb=x2ii*rcosi-x2ri*isini；cfxri+n/2=x1ri-rkb；cfxii+n

9、/2=x1ii-ikb；cfxri=x1ri+rkb；cfxii=x1ii+ikb；delete txr；delete txi；delete x1r；delete x1i；delete x2r；delete x2i；delete rcos；delete isin；/復(fù)數(shù)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 m；doub

10、le*txr=new doublen；double*txi=new doublen；for(int i=0；i n；i+)cfxii=-cfxii；fft_1d(cfxr,cfxi,txr,txi,n)；for(i=0；i len；i+)ctxri=txri/n；ctxii=-txii/n；delete txr；delete txi；/實數(shù)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

11、(2.0)；int n=1 m；double*txr=new doublen；double*txi=new doublen；for(int i=0；i n；i+)cfxii=-cfxii；fft_1d(cfxr,cfxi,txr,txi,n)；for(int i=0；i len；i+)ctxri=txri/n；delete txr；delete txi；fft：fft(void)3.test.cpp#include iostream#includefft.husing namespace std；void main()double xr10=1,2,3,4,5,6,7,8,9,10；/實部double xi10=0,0,0,0,0,0,0,0,0,0；/虛部double*cxr,*cxi；cxr=new double16；cxi=new double16；fft f；f.rfft_1d(xr,cxr,cxi,10)；for(int i=0；i 16；i+)cout cxri+jc