FFT離散傅氏變換的快速算法

FFT(離散傅氏變換的快速算法)FFT(離散傅氏變換的快速算法)目錄1算法簡(jiǎn)介2DFT算法3源碼表步4MATLAB中FFT的使用方法1算法簡(jiǎn)介編輯FFT(FastFourierTransformation),即為快速傅氏變換，是離散傅氏變換的快速算法，它是根據(jù)離散傅氏變換的奇、偶、虛、實(shí)等特性，對(duì)離散傅立葉變換的算法進(jìn)行改進(jìn)獲得的。它對(duì)傅氏變換的理論并沒有新的FFT算法圖(Bufferfly算法)發(fā)現(xiàn),但是對(duì)于在計(jì)算機(jī)系統(tǒng)或者說數(shù)字系統(tǒng)中應(yīng)用離散傅立葉變換，可以說是進(jìn)了一大步。設(shè)x(n)為N項(xiàng)的復(fù)數(shù)序列，由DFT變換，任一X(m)的計(jì)算都需要N次復(fù)數(shù)乘法和N-1次復(fù)數(shù)加法，而一次復(fù)數(shù)乘法等于四次實(shí)數(shù)乘法和兩次實(shí)數(shù)加法，一次復(fù)數(shù)加法等于兩次實(shí)數(shù)加法，即使把一次復(fù)數(shù)乘法和一次復(fù)數(shù)加法定義成一次''運(yùn)算"(四次實(shí)數(shù)乘法和四次實(shí)數(shù)加法)，那么求出N項(xiàng)復(fù)數(shù)序列的X(m),即N點(diǎn)DFT變換大約就需要NT次運(yùn)算。當(dāng)N=1024點(diǎn)甚至更多的時(shí)候，需要N2=1048576次運(yùn)算，在FFT中，利用WN的周期性和對(duì)稱性,把一個(gè)N項(xiàng)序列(設(shè)N=2k,k為正整數(shù))，分為兩個(gè)N/2項(xiàng)的子序列，每個(gè)N/2點(diǎn)DFT變換需要(N/2)2次運(yùn)算，再用N次運(yùn)算把兩個(gè)N/2點(diǎn)的DFT變換組合成一個(gè)N點(diǎn)的DFT變換。這樣變換以后，總的運(yùn)算次數(shù)就變成N+2*(N/2)a2=n+(Na2)/2°繼續(xù)上面的例子^=1024時(shí)，總的運(yùn)算次數(shù)就變成了525312次，節(jié)省了大約50%的運(yùn)算量。而如果我們將這種"一分為二"的思想不斷進(jìn)行下去，直到分成兩兩一組的DFT運(yùn)算單元，那么N點(diǎn)的DFT變換就只需要Nlog2N次的運(yùn)算,N在1024點(diǎn)時(shí)，運(yùn)算量?jī)H有10240次，是先前的直接算法的】％,點(diǎn)數(shù)越多，運(yùn)算量的節(jié)約就越大,這就是FFT的優(yōu)越性。2DFT算法編輯ForlengthNinputvectorx,theDFTisalengthNvectorX,withelementsNX(k)=sumx(n)*exp(-j*2*pi*(k-l)*(n-l)/N)/1<=k<=N.n=lTheinverseDFT(computedbyIFFT)isgivenbyNx(n)=(1/N)sumX(k)*exp(j*2*pi*(k-l)*(n-l)/N),1<=n<=N.k二13源碼表示編輯在C壞境下的源碼源碼⑴：注:親測(cè)，這個(gè)版本無法運(yùn)行，作者刪改了重要內(nèi)容[1]請(qǐng)參考源碼(2)(seepages507-508ofNumericalRecipesinC)Inputs:data[]:arrayofcomplex*datapointsofsize2*NFFT+1-data[0]isunused,*then'thcomplexnumberx(n),for0<二nv二length(x)-lzisstoredas:data[2*n+l]=real(x(n))data[2*n+2]=imag(x(n))iflength(Nx)<NFFT,theremainderofthearraymustbepaddedwithzerosnn:FFTorderNFFT.ThisMUSTbeapowerof2and>=length(x)-isign:ifsetto1,computestheforwardFFTifsetto-1,computesInverseFFT-inthiscasetheoutputvalueshavetobemanuallynormalizedbymultiplyingwith1/NFFT.Outputs:data[]:TheFFTorIFFTresultsarestoredindata,overwritingtheinput.Vvoidfourl(doubledata[]isign)intn,mmax,mJ,istep,i;doublewtemp,wrfwpr,wpi,wi,theta;doubletempr,tempi;n二nn<<1;j=1;for(i=1；ivn;i+=2){if(j>i){tempr=data[j];data(j]=data[i];data[i]=tempr;tempr=data[j+l];data[j+l]=data[i+l];data[i+l]=tempr;}m二n>>1;while(m>=2&&j>m){j-=m;m>>=1;}j+=m;}mmax=2;while(n>mmax){istep=2*mmax;theta=TWOPI/(isign*mmax);wtemp=sin*theta);wpr=*wtemp*wtemp;wpi=sin(theta);wr=;wi=;for(m=1;m<mmax;m+=2){for(i=m;iv二n;i+二istep){j=i+mmax;tempr=wr*data[j]-wi*data[j+l];tempi=wr*data[j+l]+wi*data[j];data(j]=data[i]-tempr;data(j+l]=data[i+l]-tempi;data[i]+二tempr;data[i+l]+=tempi;}wr=(wtemp=wr)*wpr-wi*wpi+wr;wi=wi*wpr+wtemp*wpi+wi;}mmax=istep;}}在c++壞境下的源碼boolFFT(complex*TDfcomplex*FD,intr)(〃一維快速Fourier變換。//complex*TD 旨向時(shí)域數(shù)組的指針；complex*FD 才旨向頻域數(shù)組的指針;r 2的幕數(shù)，即迭代次數(shù)LONGcount;//Fourier變換點(diǎn)數(shù)int“k;〃循環(huán)變量intbfsize,p;//中間變量doubleangle;//角度complex*W,*X1,*X2/*X;count=1<<r;//計(jì)算Fourier變換點(diǎn)數(shù)為1左移r位W=newcomplex[count/2];XI=newcomplex[count];X2=newcomplex[count];//分配運(yùn)算所需存儲(chǔ)器//計(jì)算加權(quán)系數(shù)(旋轉(zhuǎn)因子w的i次^表)for(i=0;i<count/2;i++)(angle=-i*PI*2/count;W[i]=complex(cos(angle),sin(angle));}〃將時(shí)域點(diǎn)寫入XImemcpy(Xl,TD,sizeof(complex)*count);//采用蝶形算法進(jìn)行快速Fourier變換for(k=0;k<r;k++)(for(