1、第4章 快速傅里叶变换(FFT),4.1 引言 4.2 基2FFT算法 4.3 进一步减少运算量的措施 4.4 其他快速算法简介,4.1 引 言DFT是数字信号分析与处理中的一种重要变换。但直接计算DFT，当N较大时，计算量太大，所以在快速傅里叶变换FFT(Fast Fourier Transform)出现以前，直接用DFT算法进行谱分析和信号的实时处理是不切实际的。直到1965年提出DFT的一种快速算法以后，情况才发生了根本的变化。,自从1965年库利和图基在计算数学杂志上发表了著名的机器计算傅里叶级数的一种算法论文后，桑德图基等快速算法相继出现，又经人们进行改进，很快形成一套高效计算方法，
2、这就是现在的快速傅里叶变换（FFT）。这种算法使DFT的运算效率提高了1 2个数量级，为数字信号处理技术应用于各种信号的实时处理创造了条件，大大推动了数字信号处理技术的发展。,人类的求知欲和科学的发展是永无止境的。多年来，人们继续寻求更快、更灵活的好算法。1984年，法国的杜哈梅尔(P. Dohamel)和霍尔曼(H. Hollmann)提出的分裂基快速算法，使运算效率进一步提高。本章主要讨论基2FFT算法。,4.2 基2FFT算法 4.2.1 直接计算DFT的特点及减少运算量的基本途径有限长序列x(n)的N点DFT为 考虑x(n)为复数序列的一般情况，对某一个k值，直接按(4.2.1)式计算
3、X(k)的1个值需要N次复数乘法和 (N1)次复数加法。因此，计算X(k)的所有N个值，共需N2次复数乘法和N(N1)次复数加法运算。,（4.2.1）,当 时，N(N1)N2。由上述可见，N点DFT的乘法和加法运算次数均为N2。当N较大时，运算量相当可观。例如N=1024时，N2=1 048 576。这对于实时信号处理来说，必将对处理设备的计算速度提出难以实现的要求。所以，必须减少其运算量，才能使DFT在各种科学和工程计算中得到应用。 如前所述，N点DFT的复乘法次数等于N2。显然，把N点DFT分解为几个较短的DFT，可使乘法次数大大减少。,FFT算法就是不断地把长序列的DFT分解成几个短序列
4、的DFT，并利用 的周期性和对称性来减少DFT的运算次数。算法最简单最常用的是基2FFT。,其对称性表现为,(4.2.3b),另外，旋转因子具有明显的周期性和对称性。其周期性表现为,4.2.2 时域抽取法基2FFT基本原理基2FFT算法分为两类：时域抽取法FFT(Decimation In Time FFT，简称DITFFT )； 频域抽取法FFT (Decimation In Frequency FFT，简称DIFFFT)。本节介绍DITFFT算法。设序列x(n)的长度为N，且满足N=2M，M为自然数。按n的奇偶把x(n)分解为两个N/2点的子序列 ,则x(n)的DFT为,因为,所以,其中1
5、(k)和X2(k)分别为x1(r)和x2(r)的N/2点DFT, 即,由于X1(k)和X2(k)均以N/2为周期，且 ，因此X(k)又可表示为,(4.2.5),(4.2.6),这样，就将N点DFT分解为两个N/2点DFT和(4.2.7)式以及(4.2.8)式的运算。(4.2.7)和(4.2.8)式的运算可用图4.2.1所示的流图符号表示，称为蝶形运算符号。采用这种图示法，经过一次奇偶抽取分解后，N点DFT运算图可以用图4.2.2表示。图中，N=23=8， X(0)X(3)由(4.2.7)式给出，而X(4) X(7)则由(4.2.8)式给出。,(4.2.7),(4.2.8),图4.2.1 蝶形运
6、算符号,偶数点的N/2 DFT,奇数点的N/2 DFT,序列DFT的N/2个点,序列DFT的后N/2个点,图4.2.2 8点DFT一次时域抽取分解运算流图,由图4.2.1可见，要完成一个蝶形运算，需要一次复数乘法和两次复数加法运算。由图4.2.2容易看出，经过一次分解后，计算1个N点DFT共需要计算两个N/2点DFT和N/2个蝶形运算。而计算一个N/2点DFT需要(N/2)2次复数乘法和N/2(N/21)次复数加法。所以，按图4.2.2计算N点DFT时，总共需要的复数乘法次数为,复数加法次数为,由此可见，仅仅经过一次分解，就使运算量减少近一半。既然这样分解对减少DFT的运算量是有效的，且N=2
7、M, N/2仍然是偶数，故可以对N/2点DFT再作进一步分解。与第一次分解相同，将x1(r)按奇偶分解成两个N/4点的子序列x3(l)和x4(l)，即,X1(k)又可表示为,(4.2.9),式中 同理，由X3(k)和X4(k)的周期性和 的对称性 最后得到：,(4.2.10),用同样的方法可计算出 ,(4.2.11),其中：,这样，经过第二次分解，又将N/2点DFT分解为2个N/4点DFT和(4.2.10)式或(4.2.11)式所示的N/4个蝶形运算，如图4.2.3所示。依次类推，经过M次分解，最后将N点DFT分解成N个1点DFT和M级蝶形运算，而1点DFT就是时域序列本身。一个完整的8点DI
8、TFFT运算流图如图4.2.4所示。图中用到关系式 。图中输入序列不是顺序排列，但后面会看到，其排列是有规律的。图中的数组A用于存放输入序列和每级运算结果。,图4.2.3 8点DFT二次时域抽取分解运算流图,图4.2.4 8点DIT-FFT运算流图,4.2.3 DIT-FFT算法与直接计算DFT运算量的比较由DIT-FFT算法的分解过程及图4.2.4可见，当N=2M 时，其运算流图应有M级蝶形，每一级都由N/2个蝶形运算构成。因此，每一级运算都需要N/2次复数乘和N次复数加(每个蝶形需要两次复数加法)。所以，M级运算总共需要的复数乘次数为,复数加次数为,而直接计算DFT的复数乘为N2次，复数加
9、为N(N1)次。当N1时，N2 (N/2) log2N，所以，DIT-FFT算法比直接计算DFT的运算次数大大减少。 例如，N=210=1024时， 这样，就使运算效率提高200多倍。图4.2.5为FFT算法和直接计算DFT所需复数乘法次数CM与变换点数N的关系曲线。由此图更加直观地看出FFT算法的优越性，显然，N越大时，优越性就越明显。,图4.2.5 DIT-FFT算法与直接计算DFT所需复数乘法次数的比较曲线,4.2.4 DIT-FFT的运算规律及蝶形画法1 原位计算 由图4.2.4可以看出，DIT-FFT的运算过程很有规律。 N=2M点的FFT共进行M级运算，每级由N/2个蝶形运算组成。
10、 同一级中，每个蝶形的两个输入数据只对计算本蝶形有用，而且每个蝶形的输入、输出数据结点又同在一条水平线上，这就意味着计算完一个蝶形后，所得输出数据可立即存入原输入数据所占用的存储单元(数组元素)。 这样，经过M级运算后，原来存放输入序列数据的N个 存储单元(数组A)中便依次存放X(k)的N个值。,8点DIT-FFT运算流图的画法,2 旋转因子的变化规律如上所述，N点DIT-FFT运算流图中，每级都有N/2个蝶形。每个蝶形都要乘以因子 ，称其为旋转因子，p为旋转因子的指数。但各级的旋转因子都有所不同。为了画出蝶形图，应先找出旋转因子 与运算级数的关系。用L表示从左到右的运算级数(L=1，2，M)
11、。观察图4.2.4不难发现，第L级共有2L1个不同的旋转因子。N=23=8时的各级旋转因子表示如下：,对N=2M的一般情况，第L级的旋转因子为：,因为 所以 ,（4.2.12）,（4.2.13）,这样，就可按(4.2.12)和(4.2.13)式确定第L级运算的旋转因子。,3 序列的倒序 DIT-FFT算法的输入序列的排序看起来似乎很乱，但仔细分析就会发现这种倒序是很有规律的。由于N=2M，因此顺序数可用M位二进制数(nM1 nM2n1n0)表示。 表4.2.1列出了N=8时以二进制数表示的顺序数和倒序数，由表显而易见，只要将顺序数(n2n1n0)的二进制位倒置，则得对应的二进制倒序值(n0n1
12、n2),表4.2.1 顺序和倒序二进制数对照表,4.2.5 频域抽取法FFT(DIF-FFT)在基2FFT算法中，频域抽取法FFT也是一种常用的快速算法，简称DTF- FFT。设序列x(n)长度为N=2M，首先将x(n)前后对半分开，得到两个子序列，其DFT可表示为如下形式：,式中 将X(k)分解成偶数组与奇数组，当k取偶数(k=2m, m=0, 1, , N/21)时 ,（4.2.14）,当k取奇数(k=2m+1, m=0, 1, , N/21)时，,令,(4.2.15), 将x1(n)和x2(n)分别代入(4.2.14)和(4.2.15)式，可得 (4.2.16)式表明，X(k)按奇偶k值
13、分为两组，其偶数组是x1(n)的N/2点DFT，奇数组则是x2(n)的N/2点DFT。x1(n)、x2(n)和x(n)之间的关系也可用图4.2.10所示的蝶形运算流图符号表示。图4.2.11表示N=8时第一次分解的运算流图。,（4.2.16）,图4.2.10 DTFFFT蝶形运算流图符号,序列的前半部分,序列的后半部分,图4.2.11 DIF-FFT第一次分解运算流图（N=8）,由于N=2M，N/2仍然是偶数，继续将N/2点DFT分成偶数组和奇数组，这样每个N/2点DFT又可由两个N/4点DFT形成，其输入序列分别是x1(n)和x2(n)按上下对半分开形成的四个子序列。 图4.2.12示出了N
14、=8时第二次分解运算流图。以这种方式分解下去，经过M1次分解，最后分解为2M1个两点DFT，两点DFT就是一个基本蝶形运算流图。 当N=8，经两次分解，便分解为四个两点DFT。N = 8的完整DIF-FFT运算流图如图4.2.13所示。,图4.2.12 DIF-FFT第二次分解运算流图（N = 8）,图4.2.13 DIF-FFT运算流图（N =8）,这种算法是对X(k)进行奇偶抽取分解的结果，所以称之为频域抽取法FFT。 观察图4.2.13可知，DIF-FFT算法与DIT-TTF算法类似，共有M级运算，每级共有N/2个蝶形运算，所以两种算法的运算次数亦相同。 不同的是DIF-FFT算法输入序
15、列为自然顺序，而输出为倒序排列。因此，M级运算完后，要对输出数据进行排序才能得到自然顺序的X(k)。 另外，蝶形运算略有不同，DIT-FFT蝶形先乘后加(减)，而DIF-FFT蝶形先加(减)后相乘。,4.2.6 IDFT的高效算法上述FFT算法流图也可以用于计算IDFT。比较DFT和IDFT的运算公式： ,只要将DFT运算式中的系数 改为 ，最后乘以1/N，就是IDFT运算公式。 所以，只要将上述的DIT-FFT与DIF-FFT算法中的旋转因子 改为 ，最后的输出再乘以1/N就可以用来计算IDFT。只是现在流图的输入是X(k)，输出就是x(n)。 因此，原来的DIT-FFT改为IFFT后，称为
16、DIF-IFFT更合适；DIF-FFT改为IFFT后, 应称为DIT-IFFT。 ,如果希望直接调用FFT子程序计算IFFT，则可用下面的方法： 由于 所以，可以先将X(k)取复共轭，然后直接调用FFT子程序，最后取复共轭并乘以1/N得到序列x(n)。这种方法虽然用了两次取共轭运算，但可以与FFT共用同一子程序，因而用起来很方便。,图4.2.4 8点DIT-FFT运算流图,L=1,L=2,L=3,4.3 进一步减少运算量的措施 4.3.1 多类蝶形单元运算由DIT-FFT运算流图已得出结论，N=2M点FFT共需要MN/2次复数乘法。由(4.2.12)式，当L=1时，只有一种旋转因子 所以，第一
17、级不需要乘法运算。当L=2 时，共有两个旋转因子： 和 ，因此，第二级也不需要乘法运算。在DFT中，又称其值为1和j的旋转因子为无关紧要的旋转因子如 ， ， 等。,综上所述，先除去第一、二两级后，所需复数乘法次数应是 进一步考虑各级中的无关紧要旋转因子。当L=3时，有两个无关紧要的旋转因子 和 ，因为同一旋转因子对应着2ML=N/2L个碟形运算，所以第三级共有2N/23=N/4 个碟形不需要复数乘法运算。依此类推，当L3时，第L级的2个无关紧要的旋转因子减少复数乘法的次数为2N/2L=N/2L1。这样，从L=3至L=M共减少复数乘法次数为,(4.3.1),因此，DIT-FFT的复乘次数降至 下
18、面再讨论FFT中特殊的复数运算，以便进一步减少复数乘法次数。一般实现一次复数乘法运算需要四次实数乘，两次实数加。但对 这一特 殊复数，任一复数(x+jy)与其相乘时，,(4.3.2),（4.3.3）,只需要两次实数加和两次实数乘就可实现。这样， 对应的每个蝶形节省两次实数乘。,在DIT-FFT运算流图中，从L=3至L=M级，每级都包含旋转因子 ，第L级中， 对应N/2L个蝶形运算。因此从第三级至最后一级，旋转因子 节省的实数乘次数与(4.3.2)式相同。所以从实数乘运算考虑，计算N=2M点DIT-FFT所需实数乘法次数为,(4.3.4),在基2FFT程序中，若包含了所有旋转因子，则称该算法为一
19、类碟形单元运算；若去掉 的旋转因子，则称之为二类碟形单元运算；若再去掉 的旋转因子，则称为三类碟形单元运算；若再判断处理 ，则称之为四类碟形运算。我们将后三种运算称为多类碟形单元运算。显然，碟形单元类型越多，编程就越复杂，但当N较大时，乘法运算的减少量是相当可观的。例如，N=4096时，三类碟形单元运算的乘法次数为一类碟形单元运算的75%。,4.3.2 旋转因子的生成在FFT运算中，旋转因子 ，求正弦和余弦函数值的计算量是很大的。所以编程时，产生旋转因子的方法直接影响运算速度。一种方法是在每级运算中直接产生；另一种方法是在FFT程序开始前预先计算出 ，m = 0，1，N/21， 存放在数组中，
20、作为旋转因子表，在程序执行过程中直接查表得到所需旋转因子值，不再计算。这样使运算速度大大提高，其不足之处是占用内存较多。,4.3.3 实序列的FFT算法在实际工作中，数据x(n)常常是实数序列。如果直接按FFT运算流图计算，就是把x(n)看成一个虚部为零的复序列进行计算，这就增加了存储量和运算时间。处理该问题的方法有两种。早期提出的方法是用一个N点FFT计算两个N点实序列的FFT，一个实序列作为x(n)的实部，另一个作为虚部，计算完FFT后，根据DFT的共轭对称性，由输出X(k)分别得到两个实序列的N点DFT（例题3.2.2）。第二种方法是用N/2点FFT计算一个N点实序列的DFT。下面简要介
21、绍第二种方法。,设x(n)为N点实序列，取x(n)的偶数点和奇数点分别作为新构造序列y(n)的实部和虚部，即 对y(n)进行N/2点FFT，输出Y(k)，则,根据DIT-FFT的思想及式(4.2.7)和(4.2.8)，可得到 X(k)的前 个值： 式中， 。由于x(n)为实序列， 因此X(k)具有共轭对称性，X(k)的后N/2点的值为,（4.3.5）,计算 N/2点FFT的复乘次数为N(M-1)/4，计算式(4.3.5)的复乘次数为N/2，所以用这种算法, 计算X(k)所需复数乘法次数为 。相对一般的N点FFT算法，上述算法的运算效率为当N=2M=210时，=20/11，运算速度提高近1倍。,
22、4.4 其他快速算法简介快速傅里叶变换算法是信号处理领域重要的研究课题。本章仅介绍算法最简单、编程最容易的基2FFT算法原理及其编程思想，使读者建立快速傅里叶变换的基本概念，了解研究FFT算法的主要途径和编程思路。例如，分裂基FFT算法、离散哈特莱变换(DHT)、基4FFT、基8FFT、基rFFT、混合基FFT，以及进一步减少运算量的途径等内容，对研究新的快速算法都是很有用的。本节简要介绍其他几种快速算法的运算量及其主要特点，以便读者选择快速算法时参考。,从理论上讲，不同基数的FFT算法的运算效率不同，实际中最常用的是基2FFT、基4 FFT、分裂基FFT和DHT。为此，下面简要介绍后三种FF
23、T算法的特点,以扩展读者的视野。 在基rFFT算法中，基4FFT算法运算效率与基8FFT很 接近，但基4FFT算法实现程序简单，且判断开销少。可以证明，当FFT的基大于4时，不会明显降低计算量。,分裂基FFT算法是将基2分解和基4分解糅合提出的，其复数乘法次数接近FFT理论最小值，但其运算流图 却与基2FFT很相似，编程简单，运算程序也很短，是 一种很实用的高效算法。 但是，对实序列x(n)，上述各种FFT算法仍将其看成虚部为零的复序列存储和计算。而一次复数乘法需要四次实数乘法和二次实数加法。所以，必然浪费存储资源和增加多余的运算量。 ,我们知道，实序列的N点DFT具有共轭对称性，即,所以，只
24、要计算出X（k）的前面N/2个值，则其后面的N/2个值可以由对称性求得。因此，FFT算法得到的N个X（k）值有一半是多余的。,所以，对实序列一定存在更高效的快速算法。离散哈特莱变换（DHT）就是针对实序列的一种高效变换算法，相对一般的FFT算法，DHT的快速算法FHT可以减少近一半的计算量。,应当说明，DHT是与DFT不同的变换，所以要想得到实序列DFT，还要根据二者的关系式进行转换。下面会看到该关系非常简单。,DHT具有以下主要特点： (1) DHT是实数变换，在对实信号进行处理时避免了复数运算，运算效率高。 (2) DHT的正、逆变换（除了因子1/N外）具有相同的形式。N点DHT定义如下：,N点逆DHT变换定义为 (3) DHT满足循环卷积定理，所以，可以直接用FHT实现实序列的快速卷积，大大提高处理速度。,（4.4.6）,(4) DHT与DFT之间的关系非常简单，容易实现二者之间 的转换，关系式如下：,（4.4.7）,对实信号x(n)进行谱分析时，可以先对x(n)进行FHT， 然后再将H(k)转换成X(k)，这样可以提高分析速度,
