基2FFT算法分析的介绍

1、快速傅里叶变换有限长序列可以通过离散傅里叶变换(DFT)将其频域也离散化成有限长序列.但其计算量太大，很难实时地处理问题，因此引出了快速傅里叶变换(FFT).1965年，Cooley和Tukey提出了计算离散傅里叶变换(DFT)的快速算法，将DFT的运算量减少了几个数量级。从此，对快速傅里叶变换(FFT)算法的研究便不断深入，数字信号处理这门新兴学科也随FFT的出现和发展而迅速发展。根据对序列分解与选取方法的不同而产生了FFT的多种算法，基本算法是基2DIT和基2DIF。FFT在离散傅里叶反变换、线性卷积和线性相关等方面也有重要应用。快速傅里叶变换(FFT)是计算离散傅里叶变换(DFT)的快速

2、算法。DFT的定义式为NknX(k)=xx(n)WNRN(k)n=0在所有复指数值W：n的值全部已算好的情况下，要计算一个X(k)需要N次复数乘法和N1次复数加法。算出全部N点X(k)共需N2次复数乘法和N(N-1)次复数加法。即计算量是与N2成正比的。FFT的基本思想：将大点数的DFT分解为若干个小点数DFT的组合，从而减少运算量。Wn因子具有以下两个特性，可使DFT运算量尽量分解为小点数的DFT运算：(1)周期性：W:咻1)n=W；n=W：nH)k对称性：wFn/2)_W：利用这两个性质，可以使DFT运算中有些项合并，以减少乘法次数。例子:求当N=4时，X(2)的值3X(2)=x(n)W：

3、n=x(0)W40xW42x(2)w44x(3)W：nz0=x(0)x(2)W40x(1)x(3)W42(周期性)=x(0)+x(2)x(1)+x(3)W：(对称性)通过合并，使乘法次数由4次减少到1次，运算量减少。FFT的算法形式有很多种，但基本上可以分为两大类：按时间抽取(DIT)和按频率抽取(DIF)。4.1按时间抽取(DIT)的FTT为了将大点数的DFT分解为小点数的DFT运算，要求序列的长度N为复合数，最常用的是N=2M的情况(M为正整数)。该情况下的变换称为基2FFT。下面讨论基2情况的算法。先将序列x(n)按奇偶项分解为两组,x(2r)=x1(r)x(2r+1)=x2(r)N川,

4、丁将DFT运算也相应分为两组N-1X(k)=DFTx(n)=、'x(n)W；nn=0N4N4=x(n)W；n+、x(n)W；nn=0n=0n为偶数n为奇数N/24N/24=''x(2r)W：k'x(2r1)WN2r伊r=0r=0N/2N/22rkk2rk=,、Xi(r)WnWn-X2()Wnr0r0N/2N/2=ZXi(r)wNk/2十w；工X2(r)wNk/2(因为W；rk=W：k/2)rOr4=Xi(k)W；X2(k)其中X1(k)、X2(k)分别是x1(n)、x2(n)的N/2点的dftN/2N/2rk一rkNX1(k)=ZXi(r)WN/2=£

5、x(2)Wn/2,0<k<-1r0r02N/2JN/2rk一rkNX2(k)=X2(r)WN/2=X(2r1)WN/2,0<k_-1r=0r=02至此，一个N点DFT被分解为两个N/2点的DFT。上面是否将全部N点的X(k)求解出来了？N分析：X(k)和X2(k)只有N/2个点(k=0,1，一一1),则由2kX(k)=X1(k)+WnX2«)只能求出X(k)的前N/2个点的DFT,要求出全部N点的X(k),需要找出X1(k)、X2(k)和X(k+N/2)的关系，其中k=0,1，£1。由式子X(k)=X1(k)+W；Xz(k)可得X(k+N/2)=X1(k+

6、N/2)+W；氏/2X2(k+N/2)化简得kNX(k+N/2)=X1(k)WNX2(k),k=0,1:，一12这样N点DFT可全部由下式确定出来：X(k)=X1(k)+W；Xz(k)N.kk=0,1,-1(*)X(k+N/2)=X1(k)WNX2(k)2上式可用一个专用的碟形符号来表示，这个符号对应一次复乘和两次复加运a-.j-bWNkaW；ba-W；b图蝶形运算符号2N2N2通过这样的分解以后，每一个N/2点的DFT只需要(一)2=次复数乘24NoN法，两个N/2点的DFT需要2()2=次复乘，再加上将两个N/2点22N2DFT合并成为N点DFT时有N/2次与W因子相乘，一共需要NN2一,

7、+也次复乘。可见，通过这样的分解，运算量节省了近一半。22因为N=2M,N/2仍然是偶数，因此可以对两个N/2点的DFT再分别作进一步的分解，将两个N/2点的DFT分解成两个N/4点的DFT。例如对x1(r),可以在按其偶数部分及奇数部分进行分解:xi(2l)=X3(l)_NI。1,一内(2l+1)=X4(l)4则的运算可相应分为两组：N/44Xi(k)=Zxi(2l)WN2l/k2l=0N/44xi(2l1)wN2l21)kl-0N/44N/44=x3(l)Wn/4WN/2x4(l)WN/4l=0l=0kN=X3(k)W：/2X4(k)k=0,1,14将系数统一为以N为周期，即W；/2=W；

8、k,可得,Xi(k)=X3(k)+WnX4(k)N.2kk=0,1,一-1Xi(kN/4)=X3(k)-Wr2kX4(k)4同样，对X2（k）也可进行类似的分解。一直分解下去，最后是2点的3一DFT,2点DFT的运算也可用碟形符号来表不。这样，对于一个N=2=8产*Mh.*1-出网化enp(3Hx(x(x(xlx(x(xX31)才司可可划刀xfxl(xxlx(xlx1x(的DFT运算，其按时间抽取的分解过程及完整流图如下图所示。x0bxHL8点DFTX(011Jux(?lQVroiQT口x30ftVfQlijy1xMl门OVfjfx510门f£lxMl»口x网沏°

9、X(G)°X(7)0)可©s'1)5)国71flflttfl-_l_tflXXXXXXXX*T*T，*TXXXXXXXX这种方法，由于每一步分解都是按输入序列在时域上的次序是属于偶数还是奇数来抽取的，故称为“时间抽取法”。分析上面的流图，N=2M,一共要进行M次分解，构成了从x(n)到X(k)的M级运算过程。每一级运算都是由N/2个蝶形运算构成，因此每一级运算都需要N/2次复乘和N次复加，则按时间抽取的M级运算后总共需要复数乘法次数：mF=M=；log2N复数加法次数：aF=NM=Nlog2N根据上面的流图，分析FFT算法的两个特点，它们对FFT的软硬件构成产生很大

10、的影响。(1)原位运算也称为同址运算，当数据输入到存储器中以后，每一级运算的结果仍然存储在原来的存储器中，直到最后输出，中间无需其它的存储器。根据运算流图分析原位运算是如何进行的。原位运算的结构可以节省存储单元，降低设备成本。(2)变址分析运算流图中的输入输出序列的顺序，输出按顺序，输入是“码位倒置”的顺序。见图。自然顺序一进制表小码位倒置码位倒置顺序0000000010011004201001023011110641000011510110156110011371111117X(0)123456X(4)X(2)X(6)X(1)X(5)X(3)码位倒置的变址处理X(7)码位倒置顺序在实际运算中

11、，直接将输入数据x(n)按码位倒置的顺序排好输入很不方便，一般总是先按自然顺序输入存储单元，然后通过变址运算将自然顺序的存储换成码位倒置顺序的存储，这样就可以进行FFT的原位运算。变质的功能如图所示。用软件实现是通用采用雷德(Rader)算法，算出I的倒序J以后立即将输入数据X(I)和X(J)对换。尽管变址运算所占运算量的比例很小，但对某些高要求的应用(尤其在实时信号处理中)，也可设法用适当的电路结构直接实现变址。例如单片数字信号处理器TMS320C25就有专用于FFT的二进制码变址模式。4.2按频率抽取（DIF）的FTT除时间抽取法外，另外一种普遍使用的FFT结构是频率抽取法。频率抽取法将输

12、入序列不是按奇、偶分组，而是将N点DFT写成前后两部分：N1knX(k)=DFTx(n)='、x(n)WNn=0(N/2)1N1=x(n)W；n+x(n)W；nnHn-N/2N/2N/2='、x(n)WN1k、x(nN/2)WNnN/2)kn=0n=SN/2="x(n)WNN/2)kx(nN/2)WN1kn=e因为W；/2=1,WNN/2)k=(1)k,k为偶数时(1)k=1,k为奇数时(-1)k=-1,由此可将X(k)分解为偶数组和奇数组：N/24X(k)=£x(n)+(-1)kx(n+N/2)W:n=0N/2X(2r)=vx(n)x(nN/2)W；nrn

13、=0N/24='x(n)x(nN/2)WNn；2n=0N/24X(2r1)="x(n)x(nN/2)WN2r1)nn=0N/24=、x(n)-x(nN/2)W；W；/2nz00,1,N/2-1x1(n)=x(n)+x(n+N/2)令<nn=x2(n)=x(n)-x(nN/2)Wn这两个序列都是N/2点的序列，对应的是两个N/2点的DFT运算:N/2JX(2r尸'、Xi(n)W%n卫N/2dX(2r1)='、X2(n)W；/2n=0这样，同样是将一个N点的DFT分解为两个N/2点的DFT。频率抽选法对应的碟形运算关系图如下:<X(0XlUJx1)xZ

14、)X(1)XR)x3oxMloMWxX团x(3)oxffl*x5)ox(6)oM小8点DFT网oX(4X(5)X(6>4点DFT4占DFTX0)X2)AX*oX6)+。Xl*。X3)zX5>ox(7)X0)XIX(2)oX7)X6)x(05xlh划x30x(5)ox6)0刈J。这种分组的办法由于每次都是按输出X(k)在频域的顺序上是属于偶数还是奇数来分组的，称为频率抽取法。与前面按时间抽取的方法相比，相同点问题：如何利用快速算法计算IDFT?分析IDFT的公式:1N1x(n)=IDFTX(k)二一'、X(k)WNTn=0,1,N-1N百比较DFT的公式：N-4X(k)=DF

15、Tx(n)='x(n)wNnk,k=0,1,N-1n=0得知可用两种方法来实现IDFT的快速算法：(1)只要把DFT运算中的每一个系数W:该为Wn,并且最后再乘以常数,就可以用时间抽取法N或频率抽取的FFT算法来直接计算IDFT。这种方法需要对FFT的程序和参数稍加改动才能实现。(2)因为1*nk*1._*、x(n)=£X(k)WN=DFTX(k),n=0,1,，N-1,也就Nk卫N是说，可先将X(k)取共轲变换，即将X(k)的虚部乘以1,就可直接调用FFT的程序，最后再对运算结果取一次共轲变换并乘以常数1/N即可得到x(n)的值。这种方法中，FFT运算和IFFT运算都可以共

16、用一个子程序块，在使用通用计算机或用硬件实现时比较方便。4.1.3混合基FFT算法以上讨论的是基2的FFT算法，即N=2M的情况，这种情况实际上使用得最多，这种FFT运算，程序简单，效率很高，用起来很方便。另外，在实际应用时，有限长序列的长度N到底是多少在很大程度时是由人为因素确定的，因此，大多数场合人们可以将N选定为N=2M,从而可以直接调用以2为基数的FFT运算程序。如果长度N不能认为确定，而N的数值又不是以2为基数的整数次方，一般可有以下两种处理方法：(1)将x(n)用补零的方法延长，使N增长到最邻近的一个N=2M数值。例如，N=30,可以在序列x(n)中补进x(30)=x(31)=0两

17、个零值点，使N=32。如果计算FFT的目的是为了了解整个频谱，而不是特定频率点，则此法可行。因为有限长序列补零以后并不影响其频谱X(ejw),只是频谱的采样点数增加而已。(2)如果要求特定频率点的频谱，则N不能改变。如果N为复合数，则可以用以任意数为基数的FFT算法来计算。快速傅里叶变换的基本思想就是要将DFT的运算尽量分小。例如，N=6时，可以按照N=3X2分解，将6点DFT分解为3组2点DFT。举例：N=9时的快速算法。4.2快速傅里叶变换的应用凡是可以利用傅里叶变换来进行分析、综合、变换的地方，都可以利用FFT算法及运用数字计算技术加以实现。FFT在数字通信、语音信号处理、图像处理、匹配

18、滤波以及功率谱估计、仿真、系统分析等各个领域都得到了广泛的应用。但不管FFT在哪里应用，一般都以卷积积分或相关积分的具体处理为依据，或者以用FFT作为连续傅里叶变换的近似为基础。4.2.1利用FFT求线性卷积一快速卷积在实际中常常遇到要求两个序列的线性卷积。如一个信号序列x(n)通过FIR滤波器时，其输出y(n)应是x(n)与h(n)的卷积：cdy(n)=x(n)h(n)=x(m)h(n-m)m二:二有限长序列x(n)与h(n)的卷积的结果y(n)也是一个有限长序列。假设x(n)与h(n)的长度分别为N1和N2,则y(n)的长度为N=N1+N2-1。若通过补零使x(n)与h(n)都加长到N点，

19、就可以用圆周卷积来计算线性卷积。这样得到用FFT运算来求y(n)值(快速卷积)的步骤如下：(1)对序列x(n)与h(n)补零至长为N,使N>N1+N2-1,并且N=2M(M为整数)，即x(n)=«x(n),n=0,1，N1-1。n=N1,N1+1,，N-1、|h(n),n=0,1，N2-1h(n)=«0,n=N2,N2+1，N1(2)用FFT计算x(n)与h(n)的离散傅里叶变换x(n)u(FFT)uX(k)(N点)h(n)u(FFT)uH(k)(N点)(3)计算丫(k)=X(k)H(k)(4)用IFFT计算Y(k)的离散傅里叶反变换得：y(n)=IFFTY(k)(N

20、点)4.2.2利用FFT求相关一快速相关互相关及自相关的运算已广泛的应用于信号分析与统计分析，应用于连续时间系统也用于离散事件系统。用FFT计算相关函数称为快速相关，它与快速卷积完全类似，不同的是一个应用离散相关定理，另一个应用离散卷积定理。同样都要注意到离散傅里叶变换固有的周期性，也同样用补零的方法来绕过这个障碍。设两个离散时间信号x(n)与y(n)为已知，离散互相关函数记作Rxy(n),定义为QORxy(n)='x(m)y(nm)m二：二如果x(n)与y(n)的序列长度分别为N1和N2,则用FFT求相关的计算步骤如下：(1)对序列x(n)与y(n)补零至长为N,使N>N1+N2-1,并且N=2m(M为整数)，即x(n),n=0,1，N1-1x(n)=«0,n=N1,N1+1,，N1,、y(n),n=0,1，N2-1y(n)=30,n=N2,N2+1，N-1(2)用FFT计算x(n)与y(n)的离散傅里叶变换x(n)u(FFT)uX(k)(N点)y(n)u(FFT)uY(k)(N点)(3)将X(k)的虚部ImX(k)改变符号，求得其共轲X*(k)(4)计算Rxy(k)