按时间抽取的基2FFT算法分析

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

2、傅里叶变换(DFT)的快速算法。DFT的定狡式为N-IX伙)=工心)呼心在所有复指数值Wf的值全部已算好的情况下，要计算一个X伙)需要N次 复数乘法和N-1次复数加法。算出全部N点X伙)共需N?次复数乘法和 N(N-I)次复数加法。即计算量是与N?成正比的。FFT的基本思想：将大点数的DFT分解为若干个小点数DFT的组合，从 而城少运算量。WN因子具有以下两个特性，可使DFT运算量尽量分解为小点数的DFT 运算：(1) 周期性：Wn = WN)k(2) 对称性：WN,2 =利用这两个性质，可以使DFT运算中有些项合并，以减少乘法次数。例子:求当N=4时，X(2)的值3X(2)=工 Xs)WJ

3、= x(O)VV4o + (1)W42 + (2)VV44 + x(3)46H-O=-r(0) + (2)VV40 +Wl) + x(3)VV42 (周期性)=(0) + x(2)-x(l) + X 网 (对称性)通过合并，使乘法次数由4次减少到1次，运算量减少。FFT的算法形式有很多种，但基本上可以分为两大类：按时间抽取(DlT)和按频率抽取(DlF)O按时间抽取(DlT)的FTT为了将大点数的DFT分解为小点数的DFT运算，要求序列的长度N为复 合数，最常用的是N = 2“的情况(M为正整数)。该情况下的变换称为基 2FFTo下面讨论基2情况的算法。先將序列X (n)按奇偶项分解为两组x(

4、2r) = XI (r)N 、x(2r + ) = x2(r)厂_，''勺"_将DFT运算也相应分为两组JV-IX(k) = DFTx(n)=x(n)W0NTNT=2>S)W(+;?=0H=OH为偶数n为奇数JV/2-1Af/2-1=Yjx(2r)Vrk ÷ £心 + 1)W严r-()r-()JV2-1W/2-1(r)W +; x2(r)r rOr-0/2-12-lX"叱2+叫2(M農(因为叱严=叱2)r-()r-0= X") + WjX2 伙)其中XI伙）、X2（Ar）分别是x1（n）s勺5）的"2点的DFTN

5、/2-172-lNX1W = X1 (r)VZ2 =工v(2r)VV2,0k-r-()r-0LN/2-1NJ 27NX2(k)=jx2(r)WZ2 = (2r + l)2,0Arl-l r-0r-0L至此，一个N点DFT被分解为两个N/2点的DFT。上面是否将全部N点的X伙）求解出来了N分析：X伙）和X?伙）只有N/2个点（&=0丄,一一1）,则由 2X伙）=Xl伙）+ WX2伙）只能求出X伙）的前N/2个点的DFT,要求出全部N点的X伙），需要找出Xl W > X2（k）和X伙+ N2）的关系，其中 =0,l, -,-Io 由式子 X伙)=Xl(k) + WX2(k)可得2X伙

6、 + N/2) = X1 (k + N2) + Wf2X2伙+ N2)化简得NX(k + N2) = = Xl(R)-W(X2(Q, R = (M,一_12这样N点DFT可全部由下式确定出来:X(k + N2) = Xk)-W.X2(k)2上式可用一个专用的碟形符号来表示，这个符号对应一次复乘和两次复加运算。a + Wba-W通过这样的分解以后，每一个N/2点、的DFT只需要()2 =次复数乘2NN法，两个N/2点的DFT需要2()2 =次复乘，再加上将两个N/2点 22NN NDFT合并成为N点DFT时有N/2次与W因子相乘，一共需要 + 222次复乘。可见，通过这样的分解，运算量节省了近一

7、半。因为N = 2'". N/2仍然是偶数，因此可以对两个N/2点的DFT再分别作进一步的分解，将两个N/2点的DFT分解成两个N/4点的DFTo例如对x1(r),可以在按其偶数部分及奇数部分进行分解：Fg心0,1,上1g + l) = x4()4則的运算可相应分为两组：yv4-lf4-1X")=工 m(2)Vc2 + > + 1)叱於“/-0/-()jV/4-1N/4-1=>3(W為+叱爲工兀叱為/-0Z-OX3(k) + W,2X4 伙)jt=0,l,-l一48 点 DFT×(0 (I) ×(2) 必(3) <×(

8、4) 啦5 叹6ZXXXXXXX(XXXXXXXX (0)M)(21问 £ x(×(×(×(x(×(×(X2dftXXXXXXXXOOOOoOoOO-OOOQO OoX1W = X3伙) + w,i伙),Al N I<k = O _ 1X、(k + N4) = X3(k)-WtkX4 伙) 4同样，对X2(Zr)也可进行类似的分解。一直分解下去，最后是2点的DFT, 2点DFT的运算也可用碟形符号来表示。这样，对于一个N = 23 =S 的DFT运算，其按时间抽取的分解过程及完整流图如下图所示。x(0- XrIb ×(

9、2kX (3Io- × M)O- x(5J0.X (6卩 ×Rk(O)同(2)问(1)(5)Pl ×(x(x(×(×(x(x(X这种方法，由于每一步分解都是按输入序列在时域上的次序是属于偶数 还是奇数来抽取的，故称为"时间抽取法S分析上面的流图，N = 2“，一共要进行M次分解，构成了从x(n)到 X(k)的M级运算过程。毎一级运算都是由N/2个蝶形运算构成，因此每一级 运算都需要N/2次复乘和N次复加，则按时间抽取的M级运算后总共需要NN复数乘法次数：HIF =M=-IosNf 22 复数加法次数：UIt =N-M= /Vlog2

10、N根据上面的流图，分析FFT算法的两个特点，它们对FFT的软硬件 构成产生很大的影响。(1 )原位运算也称为同址运算，当数据输入到存储器中以后，毎一级运算的结果仍然存储 在原来的存储器中，直到我后输出，中间无需其它的存储器。根据运算流图 分析原位运算是如何进行的。原位运算的结构可以节省存储单元，降低设备 成本。(2) 变址分析运算流图中的输入输出序列的顺序，输出按顺序，输入是“码位倒置” 的顺序。见图。自然顺序二进制表示码位倒置码位倒置顺序OOOO000010011004201001023011110641000011510110156110011371111117在实际运算中，直接舟输入数据

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

12、输入序列不是按奇、偶分组，而是将N点DFT写成前后两部分:N-IX(k) = DFnx(n)= x(n)WnH-O(T2)-1N-I= DS)W+ 2>(")wj0/r-4V2JV/2-1N/2-1=+ 2>(n + N2)Wyz2M-0r0N/2-1=yjx(n) +JV/2)V(H + / 2)W-O因为Wf" =_1,Wf 2)* =(_i)x , k为偶数时(_1)女=1, k为奇数时(一1)*=一1,由此可将X(k)分解为偶数纽和奇数纽：72-lXg= 2>() + (-1)饭+ N/2)卩曾n-0jV/2-1X(2r)=工x()+心+ N2)W

13、"H-Of/2-1=x(n) + x(n + N 2)V2rU)JV/2-1X (2r +1)=工x() - x(n + N / 2) Wi"-0N/2-1 =2>U)7G + N2)W.W篇 2片-0Xl (；0 = X(II) + x(zz + yv2)令 <H = O 丄，N/2 12() = x() 一 x(n + N2) ;这两个序列都是N/2点的序列，对应的是两个N/2点的DFT运算：2-lX(2J= x,(nM'2H-O2-lX（2 厂+ 1）=J>2陀;2-4）这样，同样是将一个N点的DFT分解为两个N/2点的DFT 了。频率抽选法

14、对应的碟形运算关系图如下：对于N=8吋频率抽取法的FFT流图如下:* IJ -IJ IJ J -IJ -IJ IJ 卩 1 234 567 fl¾ fl¾ fl¾ ( ( H ( ××××××××OOOAVO : : 厂 ； 二O歹，矽OAVOo-Or*cb77-.2占C-八“DFT送2占DFT-1厂F>o2占C-八"DFTo-12占厶八、DFTr U 一 一/ “ 卩g尼他 P R XXXX XXXX 06OOO600 kkk-IJ -IJ U -IJ -IJ IJ

15、-IJ 0426 1537«11 ll «11 «11 H »11 <l Xxxx Xxxx这种分组的办法由于每次都是按输出X(k)在频域的顺序上是属于偶数还是 奇数来分组的，称为频率抽取法。与前面按时间抽取的方法相比，相同点 问题：如何利用快速算法计算IDFT分析IDFT的公式：1 AMX(H) = IDFT X 伙)=若 X 伙)叱卩","=0,1,N 1比较DFT的公式：,-lXa) = DFTx(n) = YXS)W化 & =0,屮一1n0得知可用两种方法来实现IDFT的快速算法：(1 )只要把DFT运算中的每一

16、 个系数Wy该为Wy 并且最后再乘以常数丄，就可以用时间抽取法或N频率抽取的FFT算法来直接计算IDFTo这种方法需要对FFT的程序和参数稍 加 改 动 才 能 实 现。 (2) 因 为1 N-1W0 = 77工(QW严=77MT(灯M = O,1,N-1,也就 N A：-oN是说，可先将X(k)取共純变换，即将X(k)的虚部乘以一 1 ,就可直接调用 FFT的程序，罠后再对运算结果取一次共純变换并乘以常数1/N即可得到X (n) 的值。这种方法中，FFT运算和IFFT运算都可以共用一个子程序块，在使 用通用计算机或用硬件实现吋比较方便。4.1.3混合基FFT算法以上讨论的是基2的FFT算法，

17、即N = 2M的情况，这种情况实际上 使用得最多，这种FFT运算，程序简单，效率很鬲，用起来很方便。另外， 在实际应用时，有限长序列的长度N到底是多少在很大程度时是由人为因素 确定的，因此，大多数场合人们可以将N选定为N = 2M 9从而可以直接调 用以2为基数的FFT运算程序。如果长度N不能认为确定，而N的数值又不是以2为基数的整数次方， 一般可有以下两种处理方法：(1 ) 将x(n)用补零的方法延长，使N增长到最邻近的一个 N = 2M数值。例如,N=30,可以在序列x(n)中补进x(30)=x(31) =O两个零值点，使N=32。如果计算FFT的目的是为了 了解整 个频谱，而不是特定频率

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

19、。但不管FFT在哪里应用，一般都以卷积积分或相关积分的具体 处理为依据，或者以用FFT作为连续傅里叶变换的近似为基础。利用FFT求线性卷积一快速卷积在实际中常常遇到要求两个序列的线性卷积。如一个信号序列x(n)通 过FlR滤波器时，其输出y(n)应是x(n)与h(n)的卷积：Xy(n) = x(n) * h(n)= 工 x(ni)h(n 一 In)“10C有限长序列x(n)与h(n)的卷积的结果y(n)也是一个有限长序列。假设x (n) 与h(n)0长度分别为Nl和N2,则y(n)的长度为N=N1+N2T°若通过补零使 x(n)与h(n)都加长到N点，就可以用圆周卷积来计算线性卷枳。

20、这样得到 用FFT运算来求y(n)值(快速卷积)的步骤如下：(1 )对序列x(n)与h(n)补零至长为N,使NN1+N2-1 ,并且N = 2M (M为整数),即Lr(n), H = 0,1,9 /Vl -1X02)=O, H = NhNz ，N-If/?(/?), n = O,1,7V2-1 Il(H)= <O, n = V2,N2 + l,N-1(2) 用FFT计算x(n)与h(n)的离散傅里叶变换X(U) <=> (FFT) o X伙)(N 点)h(n) <=> (FFT) o H伙)(N 点)(3) 计算 Y(k)=X(k)H(k)(4) 用IFFT计算Y

21、(k)的离散傅里叶反变换得：y(n) = IFFTY(k) (N 点)4. 2. 2利用FFT求相关一快速相关互相关及自相关的运算已广泛的应用于信号分析与统计分析，应用于连 续时间系统也用于离散事件系统。用FFT计算相关函数称为快速相关，它与快速卷积完全类似，不同的是 一个应用离散相关定理，另一个应用离散卷积定理。同样都要注意到离散傅 里叶变换固有的周期性，也同样用补零的方法来绕过这个障碍。设两个离散时间信号x(n)与y(n)为已知，离散互相关函数记作&、.(”)，定义为OCRan= YJX(If j)y(n + m)J?!-X如果x(n)与y(n)的序列长慶分别为NI和N2,则用FFT求相关的计算步驟 如下：(1 )对序列x(n)与y(n)补零至长为N,使NMNI+N2T,并且N = 2"(M为整数),即x(n),2 = 0丄，Nl-IXW =O, n = Nl,M + l,，N-1b(n), =0,，N2-1y(n)=, n = N2N2 + X 、N-(2) 用FFT计算x(n)与y(n)的离散傅里叶变换X(H) O (FFT) o X伙)(N 点)y(n) 0 (FFT) 0 Yg(N 点)(3) 将X(k)的虚部lmX