详解FFT快速傅里叶变换FFT

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

2、)是计算离散傅里叶变换(DFT)的快速算法。DFT 的定义式为N- 1knX (k )=汇 x(n)W rn (k)n=0kn在所有复指数值 WN的值全部已算好的情况下，要计算一个X (k )需要N2次复数乘法和N 1 次复数加法。算出全部N 点 X (k ) 共需 N 次复数乘法2和 N ( N - 1) 次复数加法。即计算量是与N 2 成正比的。FFT 的基本思想：将大点数的DFT 分解为若干个小点数DFT 的组合，从而减少运算量。WN 因子具有以下两个特性，可使 DFT 运算量尽量分解为小点数的DFT运算：=Wkn = W (n+N )k( 1) 周期性：W (k+N )n N NN(

3、2) 对称性：W (k+N / 2) = - W kNN利用这两个性质，可以使 DFT 运算中有些项合并，以减少乘法次数。例子：求当n = 4时，X(2)的值3汇 44444X (2) =x(n)W2n = x(0)W 0 + x(1W 2 + x(2)W 4 + x(3)W 6n=0= x(0) +x(2)W 0 +x(1) + x(3)W2(周期性)44=x(0) + x(2)-x(1) + x(3)W(4(对称性)通过合并，使乘法次数由4次减少到1次，运算量减少。FFT的算法形式有很多种，但基本上可以分为两大类：按时间抽取(DIT)和按频率抽取(DIF)。4.1按时间抽取(DIT )的F

4、TT为了将大点数的 DFT分解为小点数的复合数，最常用的是 N = 2 M的情况(M基2FFT。下面讨论基2情况的算法。先将序列x(n)按奇偶项分解为两组?x(2r ) = x(r)DFT运算，要求序列的长度N为为正整数)。该情况下的变换称为?N?x(2r +1) = x2(r )r = 0，1L，2- 1将DFT运算也相应分为两组N - 1Xz knX (k ) = DFT x(n) =x(n)WNn=0N- 1N- 1汇 x(n)Wkn 汇 x(n)Wkn+n=0n为偶数n=0n为奇数N / 2- 1E (2 )2rk二 xr=0r WnN / 2- 1E(2 x r=0N / 2- 12

5、rkN / 2- 1+ 1)(2 r+1 )kWn2rkx2(r)WN=Xi (r)WNWnr=0N / 2- 1k N / 2- 1Nrk/ 2 +WNrk (因为 W 2 rk W rk )=2_j X12_J X2N / 2WN - WN / 2)(r)W(r)Wr =0r=0= X 1(k ) + WNkX 2(k )其中 Xi(k)、X 2(k)分别是 Xi(n)、X2 (n)的 N/2 点的 dftN / 2- 1N / 2- 1 rkrkXi(k )= E Xi(r)WN / 2 = E x(2r)WN / 2 ,0 <k <- 1r=0r=0N / 2- 1N /

6、2- 1rkrkX 2(k)= E X2 (r)WN / 2 = E x(2r +1)Wn / 2 ,0 <k<1r=0r=0至此，一个N点DFT被分解为两个N/2 点的 DFT。上面是否将全部 N点的X (k )求解出来了 ?分析：X1(k )和JN- 1)，则由X 2 (k )只有 N/2 个点(k = 0,1,L ,2X (k) = X 1 (k ) +W ：X 2(k )只能求出X (k )的前N/2个点的DFT,要求出全部N点的X (k )，需要找出X(k )、X 2 (k )和X (k + N / 2)的关系，其kN(k ) +WNX2(k )可得中 k = 0,1L

7、, - 1。由式子 X (k ) = X 1X (k + N / 2) = X (k + N / 2) +W k+N / 2 X 2(k + N / 2)化简得1NkX (k + N / 2) = = X 1(k ) - Wn X 2 (k ) ， k = 0,1,L这样N点DFT可全部由下式确定出来:kk = 0,1L N2*)?X (k) = X 1 (k)+WN X 2(k)?k?X (k + N / 2) = X 1(k) - W X2 (k)上式可用一个专用的碟形符号来表示，这个符号对应一次复乘和两次复加运N图蝶形运算符号通过这样的分解以后,每一个N/2点的DFT只需要(N)22N2

8、次复数乘法，两个N/2点的DFT需要2(N ) 2 = N一次复乘，再加上将两个 N/2点22DFT合并成 为N点DFT时有 N /2次与 W 因子相乘，一共需要N 2一次复乘。可见，通过这样的分解，运算量节省了近一半。2因为N = 2M , N/2仍然是偶数，因此可以对两个N/2点的DFT再分别作进一步的分解，将两个N/2点的DFT分解成两个N/4点的DFT。例如对xi(r),可以在按其偶数部分及奇数部分进行分解:?Xi(2l ) = X3 (l )?Xi (2l +1) = X4(l )则的运算可相应分为两组：1=0,1L，4N / 4- 121kX1(k )=三 X1 (21)WN /

9、2l=0N / 4- 1+ E xQ1W1=0(2l+1)kN / 2N / 4- 1E 3()lk/ 4x 1 Wn X 1=0N / 4 1E 4()x l Wl=0lkN / 4一k将系数统一为以N为周期，即 W N / 2k=X 3(k ) +Wn /2 X 4(k )k = 0,1L ,4 - 1=W 2 k,可得N2kk = 0,1L-N - 1?X1 (k) = X 3 (k) +Wn X 4 (k)22kXi (k + N / 4) = X3 (k) - Wn X 4 (k)同样，对X2 (k)也可进行类似的分解。一直分解下去，最后是2点的DFT , 2点DFT的运算也可用碟形

10、符号来表示。这样，对于一个N = 23 = 8DFT 运算，其按时间抽取的分解过程及完整流图如下图所示。X0)KX1*oX2声。X3&ft4点T/小DFT-*ft.F0F 号 X(7)4点OFT0)制|)印 引力 xlHMSTHxfRxf01234567-Il fl -I ( ( ( -Il -Il xxxxxxxx1J 1J. u 1J 0 4 2 6X X X XJ u 口出p R X X X Xx0)oxnb1日点DFTX(0 ill 111 JuX 2n°xnjn V (9-1I1-J x310°xkJo 乂 fOlx 4。r 41x|51o0A 1*1rt

11、W 不 1| "JI UxfGlo°XPIft VJF £!xQo“x网(71x0<> 乂网。*团 乂呵 xl) 乂冏 x3) x(7)”XD Xl) XRl X网X同 ?X5 小X同 M凶7 -1这种方法，由于每一步分解都是按输入序列在时域上的次序是属于偶数 还是奇数来抽取的，故称为“时间抽取法”分析上面的流图， N = 2 M , 一共要进行 M次分解，构成了从 x(n)到X(k)的M级运算过程。每一级运算都是由 N/2个蝶形运算构成，因此每一 级运算都需要N/2次复乘和N次复加，则按时间抽取的 M级运算后总共需 要复数乘法次数：mF = N ?

12、= N log n M 222复数加法次数：aF = N ? M = N 10g2 N根据上面的流图，分析F FT算法的两个特点，它们对F FT的软硬件构成产生很大的影响。(1) 原位运算 也称为同址运算，当数据输入到存储器中以后，每一级运算的结果仍然存储 在原来的存储器中， 直到最后输出，中间无需其它的存 储器。根据运算流图 分析原位运算是如何进行的。原位运算的结构可以节省存储单元，降低设备成本。(2) 变址 分析运算流图中的输入输出序列的顺序，输出按顺序，输入是 “码位倒置”的顺序。见图。自然顺序一进制表小码位倒置码位倒置顺序00000000100110042010010230111106

13、41000011510110156110011371111117码位倒置顺序X(0) X(4) X(2)X(6) X(1) X(5)X(3)X(7)码位倒置的变址处理在实际运算中，直接将输入数据x(n)按码位倒置的顺序排好输入很不方便，一般总是先按自然顺序输入存储单元，然后通过变址运算将自然顺序的存储换成码位倒置顺序的存储，这样就可以进行FFT 的原位运算。变质的功能如图所示。用软件实现是通用采用雷德(Rader)算法，算出I的倒序J以后立即将输入数据X(I)和X(J)对换。尽管变址运算所占运算量的比例很小，但对某些高要求的应用(尤其在实时信号处理中)，也可设法用适当的电路结构直接实现变址。例

14、如单片数字信号处理器TMS320C25 就有专用于 FFT 的二进制码变址模式。4. 2 按频率抽取(DIF )的FTT除时间抽取法外，另外一种普遍使用的FFT 结构是频率抽取法。频率抽取法将输入序列不是按奇、偶分组，而是将N点DFT写成前后两部分：N -1knx(n)WNn =0N -1( N / 2)-1汇 x(n)Wkn+ 汇 x(n)Wknn=0Nn=N / 2N / 2- 1( )N / 2- 1nk +(n+N / 2)kE ( +/ 2)x n WN x n NWNn=0n=0nkN / 2- 1(N / 2)k x(n + N / 2)WN=E x(n) +Wnn =0W N

15、/ 2 = - 1,W( N / 2) k = (- 1)k ， k 为 偶 数时 (- 1)k =1 ， k 为奇数时 NN(- 1)k = - 1 ，由此可将X(k) 分解为偶数组和奇数组：N / 2- 1X (k)= n=0 x(n) + (- 1)k x(n + N / 2)WnkN / 2- 1n=0X (2r)=N / 2- 1Ex(n) + x(n + N / 2)W 2nrx(n) + x(n + N / 2)W n=0N / 2 nrN / 2-1X (2r +1>n=)x(n) - x(n + N / 2)W(2r+1)nN / 2- 1x(n) - x(n + N

16、/ 2)W nWnr n=0n = 0,1L , N / 2 - 1N?x1(n) = x(n) + x(n + N / 2) 令?X2 (n) = x(n) - x(n + N / 2)W n这两个序列都是 N/2点的序列，对应的是两个 N/2点的DFT运算:N / 2- 1二1N / 2n=0nrX (2r) x (n)WN / 2- 1k 2 N / 2rnX (2r+1)= n=0 x (n)W这样，同样是将一个N点的DFT分解为两个N/2点的DFT 了。频率抽选法 对应的碟形运算关系图如下：-1Wn对于N=8时频率抽取法的 FFT流图如下:x(Dj0xnio8 点 DFTr V fQ

17、lI。”x3L"X ZrtVniX门q 八 1JJnV Mlx51oX°x网乂C(7)X X(Q)*。X|2) xM) x|6| XU)a。X|3)"X(5)x(0). x(1)o 刈4 x(30x(5)0 x(6)0 孙2点DFT2点DFT2占DFTX(0) or口>o2占 匚八% DFTX X(7)这种分组的办法由于每次都是按输出X(k)在频域的顺序上是属于偶数还是奇数来分组的，称为频率抽取法。与前面按时间抽取的方法相比，相同点问题：如何利用快速算法计算IDFT ?分析IDFT的公式：N - 1X=IDFTX (k) 二三 X (k)WN-nk,n =

18、0,1L ,N -N k=o比较DFT的公式:N - 1X (k) = DFTx(n) = £； x(n)WNnk ,k = 0,1L , N - 1n=0得知可用两种方法来实现 IDFT的快速算法：(1)只要把 DFT运算中的每个系数Wnk该为W-nk,并且最后再乘以常数NN1、，一，、一.一，就可以用时间抽取法N或频率抽取的FFT算法来直接计算IDFT。这种方法需要对 FFT的程序和参数稍加 改动才能实现。 (2) 因 为/、_ 1 N-11*x(n) # X * (k)W/nk * >DFTX (k), n = 0,1L ,N - 1,也就N k=0N是说，可先将 X(k

19、)取共轲变换，即将 X(k)的虚部乘以1 ,就可直接调用FFT的程序，最后再对运算结果取一次共轲变换并乘以常数1/N即可得到x(n)的值。这种方法中，FFT运算和IFFT运算都可以共用一个子程序块， 在使用通用计算机或用硬件实现时比较方便。4.1.3混合基FFT算法以上讨论的是基 2的FFT算法，即N = 2 M的情况，这种情况实际上使用得最多，这种 FFT运算，程序简单，效率很高，用起来很方便。另 外，在实际应用时，有限长序列的长度N到底是多少在很大程度时是由人为因素确定的，因此，大多数场合人们可以将N选定为N = 2 M ,从而可以直接调用以2为基数的 FFT运算程序。如果长度N不能认为确

20、定，而 N的数值又不是以2为基数的整数次 方，一般可有以下两种处理方法：(1) 将x(n)用补零的方法延长，使 N增长到最邻近的一个N = 2 M数值。仞0口， N=30,可以在序列 x(n)中补进x(30) = x(31) = 0两个零值点， 使N=32。如果af算 FFT 的目的是为了了解整个频谱， 而不是特定频率点， 则此 法可行。因为有限长序列补零以后并不影响其频谱X (ejw ),只是频谱的采样点数增加而已。(2) 如果要求特定频率点的频谱，则 N不能改变。如果 N为复合数，则可以用以任意数为基数的FFT算法来计DFT的运算N=3X 2分解，算。快速傅里叶变换的基本思想就是要将 尽量

21、分小。例如，N=6时，可以按照DFT。将6点DFT分解为3组2点举例： N=9 时的快速算法。4.2 快速傅里叶变换的应用凡是可以利用傅里叶变换来进行分析、综合、 变换的地方，都可以利用 FFT 算法及运用数字计算技术加以实现。FFT 在数字通信、语音信号处理、 图像处理、匹配滤波以及功率谱估计、仿真、系统分析等各个领域都得到了广泛的应用。但不管FFT 在哪里应用，一般都以卷积积分或相关积分的具体处理为依据，或者以用FFT 作为连续傅里叶变换的近似为基础。4.2.1 利用 FFT 求线性卷积快速卷积在实际中常常遇到要求两个序列的线性卷积。如一个信号序列x(n)通过FIR滤波器时，其输出y(n)

22、应是x(n)与h(n)的卷积：y(n) = x(n) ? h(n) =x(m)h(n - m)有限长序列x(n)与h(n)的卷积的结果y(n)也是一个有限长序列。假设 x(n)与 h(n)的长度分别为 N1和N2,则y(n)的长度为 N=N1+N2-1。若通过补零使 x(n)与h(n)都加长到N点，就可以用圆周卷积来计算线性卷积。这样得到用FFT 运算来求y(n) 值(快速卷积)的步骤如下：(1) 对序列 x(n)与h(n)补零至长为 N，使 N > N1+N2-1 ,并且 N =皆(M为整数)，即?x(n), n = 01,L , N1 - 1 x(n) = ?0, n = N1, N

23、1 +1,L , N - 1?h(n), n = 01,L , N 2 - 1 h(n) = ?0, n = N2, N 2 +1,L , N - 1(2)用FFT计算x(n)与h(n)的离散傅里叶变换x(n) ? (FFT ) ? X (k )( N 点)h(n) ? (FFT ) ? H (k )N 点4.2 快速傅里叶变换的应用(3 ) 计算 Y(k)=X(k)H(k)(4) 用IFFT计算Y(k)的离散傅里叶反变换得：y(n)=IFFTY(k) ( N点)4.2.2 利用 FFT 求相关快速相关互相关及自相关的运算已广泛的应用于信号分析与统计分析，应用于连续时间系统也用于离散事件系统。

24、用 FFT 计算相关函数称为快速相关，它与快速卷积完全类似，不同的是一个应用离散相关定理，另一个应用离散卷积定理。同样都要注意到离散傅里叶变换固有的周期性，也同样用补零的方法来绕过这个障碍。设两个离散时间信号 x(n)与y(n)为已知，离散互相关函数记作Rxy (n),定义为Rxy (n)=汇 x(m)y(n + m)m=- °°如果x(n)与y(n)的序列长度分别为 N1和N2,则用FFT求相关的计算步骤如下：(1)对序列x(n)与y(n)补零至长为 N，使 N > N1+N2-1 ,并且N =2" (M为整数)，即?x(n), n = 01,L , N1 - 1 x(n) = ?0, n = N1, N1 +1,L , N - 1?y(n), n =01,L , N 2 - 1 y(n) = ?0, n = N2, N 2 +1,L , N - 1(2 )用FFT计算x(n)与y(n)的离散傅里叶变换x(n) ? (FFT ) ?X (k )( N 点)y(n) ? ( FF