FFT离散傅氏变换的快速算法

FFT(离散傅氏变换的快速算法)FFT(离散傅氏变换的快速算法)目录1算法简介2DFT算法3源码表步4MATLAB中FFT的使用方法1算法简介编辑FFT(FastFourierTransformation),即为快速傅氏变换，是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并没有新的FFT算法图(Bufferfly算法)发现,但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。设x(n)为N项的复数序列，由DFT变换，任一X(m)的计算都需要N次复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次''运算"(四次实数乘法和四次实数加法)，那么求出N项复数序列的X(m),即N点DFT变换大约就需要NT次运算。当N=1024点甚至更多的时候，需要N2=1048576次运算，在FFT中，利用WN的周期性和对称性,把一个N项序列(设N=2k,k为正整数)，分为两个N/2项的子序列，每个N/2点DFT变换需要(N/2)2次运算，再用N次运算把两个N/2点的DFT变换组合成一个N点的DFT变换。这样变换以后，总的运算次数就变成N+2*(N/2)a2=n+(Na2)/2°继续上面的例子^=1024时，总的运算次数就变成了525312次，节省了大约50%的运算量。而如果我们将这种"一分为二"的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N点的DFT变换就只需要Nlog2N次的运算,N在1024点时，运算量仅有10240次，是先前的直接算法的】％,点数越多，运算量的节约就越大,这就是FFT的优越性。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坏境下的源码源码⑴：注:亲测，这个版本无法运行，作者删改了重要内容[1]请参考源码(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 旨向时域数组的指针；complex*FD 才旨向频域数组的指针;r 2的幕数，即迭代次数LONGcount;//Fourier变换点数int“k;〃循环变量intbfsize,p;//中间变量doubleangle;//角度complex*W,*X1,*X2/*X;count=1<<r;//计算Fourier变换点数为1左移r位W=newcomplex[count/2];XI=newcomplex[count];X2=newcomplex[count];//分配运算所需存储器//计算加权系数(旋转因子w的i次^表)for(i=0;i<count/2;i++)(angle=-i*PI*2/count;W[i]=complex(cos(angle),sin(angle));}〃将时域点写入XImemcpy(Xl,TD,sizeof(complex)*count);//采用蝶形算法进行快速Fourier变换for(k=0;k<r;k++)(for(