1、第4章 快速傅里叶变换(FFT)第第4章章 快速傅里叶变换快速傅里叶变换(FFT)4.1 引言引言 4.2 基2FFT算法 4.3 进一步减少运算量的措施 4.4 其他快速算法简介习题与上机题第4章 快速傅里叶变换(FFT)4.1 引引 言言DFT是数字信号分析与处理中的一种重要变换。但直接计算DFT的计算量与变换区间长度N的平方成正比，当N较大时，计算量太大，所以在快速傅里叶变换FFT(Fast Fourier Transform)出现以前，直接用DFT算法进行谱分析和信号的实时处理是不切实际的。直到1965年提出DFT的一种快速算法以后，情况才发生了根本的变化。第4章 快速傅里叶变换(FF
2、T)自从1965年库利(T.W.Cooley)和图基(J.W.Tuky)在计算数学(Math.Computation,Vol.19,1965)杂志上发表了著名的机器计算傅里叶级数的一种算法论文后，桑德(G.Sand)图基等快速算法相继出现，又经人们进行改进，很快形成一套高效计算方法，这就是现在的快速傅里叶变换（FFT）。这种算法使DFT的运算效率提高了1 2个数量级，为数字信号处理技术应用于各种信号的实时处理创造了条件，大大推动了数字信号处理技术的发展。第4章 快速傅里叶变换(FFT)人类的求知欲和科学的发展是永无止境的。多年来，人们继续寻求更快、更灵活的好算法。1984年，法国的杜哈梅尔(P
3、.Dohamel)和霍尔曼(H.Hollmann)提出的分裂基快速算法，使运算效率进一步提高。本章主要讨论基2FFT算法及其编程思想。第4章 快速傅里叶变换(FFT)4.2 基基2FFT算法算法4.2.1 直接计算直接计算DFT的特点及减少运算量的基本途径的特点及减少运算量的基本途径有限长序列x(n)的N点DFT为 (4.2.1)考虑x(n)为复数序列的一般情况，对某一个k值，直接按(4.2.1)式计算X(k)的1个值需要N次复数乘法和(N1)次复数加法。因此，计算X(k)的所有N个值，共需N2次复数乘法和N(N1)次复数加法运算。10110 )()(NnknNNkWnxkX，第4章 快速傅里
4、叶变换(FFT)1N当时，N(N1)N2。由上述可见，N点DFT的乘法和加法运算次数均为N2。当N较大时，运算量相当可观。例如N=1024时，N2=1 048 576。这对于实时信号处理来说，必将对处理设备的计算速度提出难以实现的要求。所以，必须减少其运算量，才能使DFT在各种科学和工程计算中得到应用。如前所述，N点DFT的复乘次数等于N2。显然，把N点DFT分解为几个较短的DFT，可使乘法次数大大减少。另外，旋转因子具有明显的周期性和对称性。其周期性表现为 （4.2.2）mNmNlNmNlNmNWW2j)(2jee第4章 快速傅里叶变换(FFT)(4.2.3a)(4.2.3b)FFT算法就是
5、不断地把长序列的DFT分解成几个短序列的DFT，并利用的周期性和对称性来减少DFT的运算次数。算法最简单最常用的是基2FFT。其对称性表现为mNNmNWW或者 mNmNNWW*mNNmNWW2knNW第4章 快速傅里叶变换(FFT)4.2.2 时域抽取法基时域抽取法基2FFT基本原理基本原理基2FFT算法分为两类：时域抽取法FFT(DecimationIn Time FFT，简称DITFFT)；频域抽取法FFT(Decimation In Frequency FFT，简称DIFFFT)。本节介绍DITFFT算法。设序列x(n)的长度为N，且满足N=2M，M为自然数。按n的奇偶把x(n)分解为两
6、个N/2点的子序列12()(2)0112()(21)0112Nx rxrrNx rxrr，第4章 快速傅里叶变换(FFT)则x(n)的DFT为/2 1/2 12(21)00/2 1/2 1221200()()()(2)(21)()()knknNNnnNNkrkrNNrrNNkrkkrNNNrrX kx n Wx n Wxr WxrWx r WWx r W偶数奇数第4章 快速傅里叶变换(FFT)因为所以(4.2.4)22jj222/2eekrkrkrkrNNNNWW/2 1/2 11/22/20012()()()()()0,1,2,-1 NNkrrkrNNNrrkNX kx r WWx r WX
7、 kW XkkN第4章 快速傅里叶变换(FFT)其中1(k)和X2(k)分别为x1(r)和x2(r)的N/2点DFT,即(4.2.5)(4.2.6)由于X1(k)和X2(k)均以N/2为周期，且 ，因此X(k)又可表示为/2 111/2120()()DFT()NkrNNrX kx r Wx r/2 122/2220()()DFT()NkrNNrXkx r Wx rkNNkNWW2第4章 快速傅里叶变换(FFT)(4.2.7)(4.2.8)这样，就将N点DFT分解为两个N/2点DFT和(4.2.7)式以及(4.2.8)式的运算。(4.2.7)和(4.2.8)式的运算可用图4.2.1所示的流图符号
8、表示，称为蝶形运算符号。采用这种图示法，经过一次奇偶抽取分解后，N点DFT运算图可以用图4.2.2表示。图中，N=23=8，X(0)X(3)由(4.2.7)式给出，而X(4)X(7)则由(4.2.8)式给出。1210)()()(21NkkXWkXkXkN，1210)()()2(21NkkXWkXNkXkN，第4章 快速傅里叶变换(FFT)图4.2.1 蝶形运算符号第4章 快速傅里叶变换(FFT)图4.2.2 8点DFT一次时域抽取分解运算流图第4章 快速傅里叶变换(FFT)由图4.2.1可见，要完成一个蝶形运算，需要一次复数乘法和两次复数加法运算。由图4.2.2容易看出，经过一次分解后，计算1
9、个N点DFT共需要计算两个N/2点DFT和N/2个蝶形运算。而计算一个N/2点DFT需要(N/2)2次复数乘法和N/2(N/21)次复数加法。所以，按图4.2.2计算N点DFT时，总共需要的复数乘法次数为 221(1)22222NNNN NN第4章 快速傅里叶变换(FFT)复数加法次数为由此可见，仅仅经过一次分解，就使运算量减少近一半。既然这样分解对减少DFT的运算量是有效的，且N=2M,N/2仍然是偶数，故可以对N/2点DFT再作进一步分解。与第一次分解相同，将x1(r)按奇偶分解成两个N/4点的子序列x3(l)和x4(l)，即222122NNNN1410)12()()2()(1423Nll
10、xlxlxlx，第4章 快速傅里叶变换(FFT)X1(k)又可表示为(4.2.9)1210)()()()()12()2()(42/314/04/42/14/04/314/0)12(2/114/022/11NkkXWkXWlxWWlxWlxWlxkXkNNlklNkNNlklNNllkNNlklN，第4章 快速傅里叶变换(FFT)式中同理，由X3(k)和X4(k)的周期性和的对称性最后得到：(4.2.10)/4 133/4340()()DFT()NklNNlXkx l Wx l/4 144/4440()()DFT()NklNNlXkx l Wx lmNW2/kNNkNWW2/4/2/14/10)
11、()()4/()()()(42/3142/31NkkXWkXNkXkXWkXkXkNkN，第4章 快速傅里叶变换(FFT)用同样的方法可计算出 (4.2.11)其中1410 )()(4)()()(62/5262/52NkkXWkXNkXkXWkXkXkNkN，/4 155/4540/4 166/46405262()()DFT()()()DFT()()(2)01/4 1()(21)NklNNlNklNNlXkx l Wx lXkx l Wx lx lxllNx lxl，第4章 快速傅里叶变换(FFT)这样，经过第二次分解，又将N/2点DFT分解为2个N/4点DFT和(4.2.10)式或(4.2.
12、11)式所示的N/4个蝶形运算，如图4.2.3所示。依次类推，经过M次分解，最后将N点DFT分解成N个1点DFT和M级蝶形运算，而1点DFT就是时域序列本身。一个完整的8点DITFFT运算流图如图4.2.4所示。图中用到关系式。图中输入序列不是顺序排列，但后面会看到，其排列是有规律的。图中的数组A用于存放输入序列和每级运算结果，在后面讨论编程方法和倒序时要用到。mkNkmNWW/第4章 快速傅里叶变换(FFT)图4.2.3 8点DFT二次时域抽取分解运算流图第4章 快速傅里叶变换(FFT)图4.2.4 8点DIT-FFT运算流图第4章 快速傅里叶变换(FFT)4.2.3 DIT-FFT算法与直
13、接计算算法与直接计算DFT运算量的比较运算量的比较由DIT-FFT算法的分解过程及图4.2.4可见，N=2M 时，其运算流图应有M级蝶形，每一级都由N/2个蝶形运算构成。因此，每一级运算都需要N/2次复数乘和N次复数加(每个蝶形需要两次复数加法)。所以，M级运算总共需要的复数乘次数为复数加次数为lb22MNNCMNlbACN MN N第4章 快速傅里叶变换(FFT)而直接计算DFT的复数乘为N2次，复数加为N(N1)次。当N1时，N2(N/2)lbN，所以，DIT-FFT算法比直接计算DFT的运算次数大大减少。例如，N=210=1024时，这样，就使运算效率提高200多倍。图4.2.5为FFT
14、算法和直接计算DFT所需复数乘法次数CM与变换点数N的关系曲线。由此图更加直观地看出FFT算法的优越性，显然，N越大时，优越性就越明显。8.2045120576 048 1lb22NNN第4章 快速傅里叶变换(FFT)图4.2.5 DIT-FFT算法与直接计算DFT所需复数乘法次数的比较曲线第4章 快速傅里叶变换(FFT)4.2.4 DIT-FFT的运算规律及编程思想的运算规律及编程思想1 原位计算原位计算由图4.2.4可以看出，DIT-FFT的运算过程很有规律。N=2M点的FFT共进行M级运算，每级由N/2个蝶形运算组成。同一级中，每个蝶形的两个输入数据只对计算本蝶形有用，而且每个蝶形的输入
15、、输出数据结点又同在一条水平线上，这就意味着计算完一个蝶形后，所得输出数据可立即存入原输入数据所占用的存储单元(数组元素)。这样，经过M级运算后，原来存放输入序列数据的N个存储单元(数组A)中便依次存放X(k)的N个值。这种利用同一存储单元存储蝶形计算输入、输出数据的方法称为原位(址)计算。原位计算可节省大量内存，从而使设备成本降低。第4章 快速傅里叶变换(FFT)2 旋转因子的变化规律旋转因子的变化规律如上所述，N点DIT-FFT运算流图中，每级都有N/2个蝶形。每个蝶形都要乘以因子，称其为旋转因子，p为旋转因子的指数。但各级的旋转因子和循环方式都有所不同。为了编写计算程序，应先找出旋转因子
16、与运算级数的关系。用L表示从左到右的运算级数(L=1，2，M)。观察图4.2.4不难发现，第L级共有2L1个不同的旋转因子。N=23=8时的各级旋转因子表示如下：pNWpNW第4章 快速傅里叶变换(FFT)对N=2M的一般情况，第L级的旋转因子为3,2,1,0 31,0 20 1222/24/JWWWLJWWWLJWWWLJJNpNJJNpNJJNpNLLL时时时第4章 快速傅里叶变换(FFT)因为所以(4.2.12)(4.2.13)这样，就可按(4.2.12)和(4.2.13)式确定第L级运算的旋转因子(实际编程序时，L为最外层循环变量)。L-12,0,1,2,21LpJNWWJMLMLML
17、N222212,2,1,0122LJNJNpNJWWWLMMLLMJp2第4章 快速傅里叶变换(FFT)3 蝶形运算规律蝶形运算规律设序列x(n)经时域抽选(倒序)后，按图4.2.4所示的次序(倒序)存入数组A中。如果蝶形运算的两个输入数据相距B个点，应用原位计算，则蝶形运算可表示成如下形式：式中下标L表示第L级运算，AL(J)则表示第L级运算后的数组元素A(J)的值（即第L级蝶形的输出数据）。而AL1(J)表示第L级运算前A(J)的值（即第L级蝶形的输入数据）。11()()()pLLLNAJAJAJB W11()()()pLLLNA JBAJAJB W12 0,1,21;1,2,MLLpJJ
18、LM第4章 快速傅里叶变换(FFT)4.编程思想及程序框图编程思想及程序框图 仔细观察图4.2.4，还可以归纳出一些对编程有用的运算规律：第L级中，每个蝶形的两个输入数据相距B=2L1个点；每级有B个不同的旋转因子；同一旋转因子对应着间隔为2L点的2ML个蝶形。总结上述运算规律，便可采用下述运算方法。先从输入端(第1级)开始，逐级进行，共进行M级运算。在进行第L级运算时，依次求出B个不同的旋转因子，每求出一个旋转因子，就计算完它对应的所有2ML个蝶形。这样，我们可用三重循环程序实现DIT-FFT运算，程序框图如图4.2.6所示。第4章 快速傅里叶变换(FFT)图4.2.6 DIT-FFT运算和
19、程序框图第4章 快速傅里叶变换(FFT)另外，DIT-FFT算法运算流图的输出X(k)为自然顺序，但为了适应原位计算，其输入序列不是按x(n)的自然顺序排列，这种经过M次偶奇抽选后的排序称为序列x(n)的倒序(倒位)。因此，在运算M级蝶形之前应先对序列x(n)进行倒序。下面介绍倒序算法。5 序列的倒序序列的倒序DIT-FFT算法的输入序列的排序看起来似乎很乱，但仔细分析就会发现这种倒序是很有规律的。由于N=2M，因此顺序数可用M位二进制数(nM1nM2n1n0)表示。M次偶奇时域抽选过程如图4.2.7所示。第4章 快速傅里叶变换(FFT)第一次按最低位n0的0和1将x(n)分解为偶奇两组，第二
20、次又按次低位n1的0、1值分别对偶奇组分组；依次类推，第M次按nM1位分解，最后所得二进制倒序数如图4.2.7所示。表4.2.1列出了N=8时以二进制数表示的顺序数和倒序数，由表显而易见，只要将顺序数(n2n1n0)的二进制位倒置，则得对应的二进制倒序值(n0n1n2)。按这一规律，用硬件电路和汇编语言程序产生倒序数很容易。但用有些高级语言程序实现时，直接倒置二进制数位是不行的，因此必须找出产生倒序数的十进制运算规律。第4章 快速傅里叶变换(FFT)由表4.2.1可见，自然顺序数I增加1，是在顺序数的二进制数最低位加1，逢2向高位进位。而倒序数则是在M位二进制数最高位加1，逢2向低位进位。例如
21、，在(000)最高位加1，则得（100），而(100)最高位为1，所以最高位加1要向次高位进位，其实质是将最高位变为0，再在次高位加1，得到(010)。用这种算法，可以从当前任一倒序值求得下一个倒序值。第4章 快速傅里叶变换(FFT)图4.2.7 形成例序的树状图（N=23）第4章 快速傅里叶变换(FFT)表表4.2.1 顺序和倒序二进制数对照表顺序和倒序二进制数对照表第4章 快速傅里叶变换(FFT)为了叙述方便，用J表示当前倒序数的十进制数值。对于N=2M，M位二进制数最高位的十进制权值为N/2，且从左向右二进制位的权值依次为N/4，N/8，2，1。因此，最高为加1相当于十进制运算J+N/2
22、。如果最高位是0(JN/2)，则直接由J+N/2得下一个倒序值；如最高位是1(JN/2),则先将最高位变成0(JJN/2)，然后次高位加1(J+N/4)。但次高位加1时，同样要判断0、1值，如果为0(JN/4)，则直接加1(JJ+N/4)，否则将次高位变成0(JJN/4)，再判断下一位；依此类推，直到完成最高位加1，逢2向右进位的运算。图4.2.9所示的倒序的程序框图中的虚线框内就是完成计算倒序值的运算流程图。第4章 快速傅里叶变换(FFT)形成倒序J后，将原数组A中存放的输入序列重新按倒序排列。设原输入序列x(I)先按自然顺序存入数组A中。例如，对N=8，A(0)，A(1)，A(7)中依次存
23、放着x(0)，x(1)，x(2)，x(7)。对x(n)的重新排序(倒序)规律如图4.2.8所示。倒序的程序框图如图4.2.9所示。由图4.2.8可见，第一个序列值x(0)和最后一个序列值x(N1)不需要重排，每计算出一个倒序值J，便与循环语句自动生成的顺序I比较，当I=J时，不需要交换，当IJ时,A(I)与A(J)交换数据。另外，为了避免再次调换前面已调换过的一对数据，框图中只对IJ的情况调换A(I)和A(J)的内容。第4章 快速傅里叶变换(FFT)图4.2.8 倒序规律第4章 快速傅里叶变换(FFT)图4.2.9 倒序程序框图第4章 快速傅里叶变换(FFT)第3章介绍的MATLAB函数fft
24、是一个计算DFT的智能程序，如果计算点数N=2M，则自动按DIT-FFT快速算法计算，否则，直接计算DFT。所以，调用该函数计算DFT时，最好选取N=2M，使处理速度大大提高。第4章 快速傅里叶变换(FFT)4.2.5 频域抽取法频域抽取法FFT(DIF-FFT)在基2FFT算法中，频域抽取法FFT也是一种常用的快速算法，简称DTF-FFT。设序列x(n)长度为N=2M，首先将x(n)前后对半分开，得到两个子序列，其DFT可表示为如下形式：10/2 110/2/2 1/2 1(/2)00/2 1/20()DFT()()()()()2()2NknNnNNknknNNnn NNNknk n NNN
25、nnNkNknNNnX kx nx n Wx n Wx n WNx n Wx nWNx nWx nW第4章 快速傅里叶变换(FFT)式中将X(k)分解成偶数组与奇数组，当k取偶数(k=2m,m=0,1,N/21)时/21(1)1kNkNkWk 偶数奇数，/220/2 1/20(2)()2()(4.2.14)2NmnNnNmnNnNXmx nx nWNx nx nW第4章 快速傅里叶变换(FFT)当k取奇数(k=2m+1,m=0,1,N/21)时，/2 1(21)0/2 1/20(21)()2()2 (4.2.15)NnmNnNnnmNNnNXmx nx nWNx nx nWW令122102)(
26、)(2)()(21NnWNnxnxnxNnxnxnxnN，第4章 快速傅里叶变换(FFT)将x1(n)和x2(n)分别代入(4.2.14)和(4.2.15)式，可得 (4.2.16)(4.2.16)式表明，X(k)按奇偶k值分为两组，其偶数组是x1(n)的N/2点DFT，奇数组则是x2(n)的N/2点DFT。x1(n)、x2(n)和x(n)之间的关系也可用图4.2.10所示的蝶形运算流图符号表示。图4.2.11表示N=8时第一次分解的运算流图。/2 11/20/2 12/20(2)()(21)()NmnNnNmnNnXmx n WXmx n W第4章 快速傅里叶变换(FFT)图4.2.10 D
27、TFFFT蝶形运算流图符号第4章 快速傅里叶变换(FFT)图4.2.11 DIF-FFT第一次分解运算流图（N=8）第4章 快速傅里叶变换(FFT)由于N=2M，N/2仍然是偶数，继续将N/2点DFT分成偶数组合奇数组，这样每个N/2点DFT又可由两个N/4点DFT形成，其输入序列分别是x1(n)和x2(n)按上下对半分开形成的四个子序列。图4.2.12示出了N=8时第二次分解运算流图。这样继续分解下去，经过M1次分解，最后分解为2M1个两点DFT，两点DFT就是一个基本蝶形运算流图。当N=8，经两次分解，便分解为四个两点DFT。N=8的完整DIF-FFT运算流图如图4.2.13所示。第4章
28、快速傅里叶变换(FFT)图4.2.12 DIF-FFT第二次分解运算流图（N=8）第4章 快速傅里叶变换(FFT)图4.2.13 DIF-FFT运算流图（N=8）第4章 快速傅里叶变换(FFT)这种算法是对X(k)进行奇偶抽取分解的结果，所以称之为频域抽取法FFT。观察图4.2.13可知，DIF-FFT算法与DIT-TTF算法类似，可以原位计算，共有M级运算，每级共有N/2个蝶形运算，所以两种算法的运算次数亦相同。不同的是DIF-FFT算法输入序列为自然顺序，而输出为倒序排列。因此，M级运算完后，要对输出数据进行倒序才能得到自然顺序的X(k)。另外，蝶形运算略有不同，DIT-FFT蝶形先乘后加
29、(减)，而DIF-FFT蝶形先加(减)后相乘。最后要说明的是，上述两种FFT的算法流图形式不是唯一的。只要保证各支路传输比不变，改变输入与输出点以及中间结点的排列顺序，就可以得到其他变形的FFT运算流图，各种运算流图各有特点1。第4章 快速傅里叶变换(FFT)4.2.6 IDFT的高效算法的高效算法上述FFT算法流图也可以用于计算IDFT。比较DFT和IDFT的运算公式：1010()DFT()()1()IDFT()()NknNnNknNkX kx nx n Wx nx nX k WN第4章 快速傅里叶变换(FFT)只要将DFT运算式中的系数改变为，最后乘以1/N，就是IDFT运算公式。所以，只
30、要将上述的DIT-FFT与DIF-FFT算法中的旋转因子改为，最后的输出再乘以1/N就可以用来计算IDFT。只是现在流图的输入是X(k)，输出就是x(n)。因此，原来的DIT-FFT改为IFFT后，称为DIF-IFFT更合适；DIF-FFT改为IFFT后,应称为DIT-IFFT。knNWpNWpNWknNW第4章 快速傅里叶变换(FFT)如果希望直接调用FFT子程序计算IFFT，则可用下面的方法：由于所以，可以先将X(k)取复共轭，然后直接调用FFT子程序，或者送入FFT专用硬件设备进行DFT运算，最后取复共轭并乘以1/N得到序列x(n)。这种方法虽然用了两次取共轭运算，但可以与FFT共用同一
31、子程序，因而用起来很方便。*10*)(1)(1)(kXDFTNWkXNnxNkknN第4章 快速傅里叶变换(FFT)4.3 进一步减少运算量的措施进一步减少运算量的措施4.3.1 多类蝶形单元运算多类蝶形单元运算由DIT-FFT运算流图已得出结论，N=2M点FFT共需要MN/2次复数乘法。由(4.2.12)式，当L=1时，只有一种旋转因子 ，所以，第一级不需要乘法运算。当L=2 时，共有两个旋转因子：和 ，因此，第二级也不需要乘法运算。在DFT中，又称其值为1和j的旋转因子为无关紧要的旋转因子，如 ，等。10NW10NWjWNN4/0NW2/NNW4/NNW第4章 快速傅里叶变换(FFT)综上
32、所述，先除去第一、二两级后，所需复数乘法次数应是(4.3.1)进一步考虑各级中的无关紧要旋转因子。当L=3时，有两个无关紧要的旋转因子和，因为同一旋转因子对应着2ML=N/2L个碟形运算，所以第三级共有2N/23=N/4 个碟形不需要复数乘法运算。依此类推，当L3时，第L级的2个无关紧要的旋转因子减少复数乘法的次数为2N/2L=N/2L1。这样，从L=3至L=M共减少复数乘法次数为(2)2MNCM0NW4/NNW第4章 快速傅里叶变换(FFT)因此，DIT-FFT的复乘次数降至 (4.3.3)下面再讨论FFT中特殊的复数运算，以便进一步减少复数乘法次数。一般实现一次复数乘法运算需要四次实数乘，
33、两次实数加。但对这一特殊复数，任一复数(x+jy)与其相乘时，(4.3.2)222122331NNNMLLMLL(2)2(3)2222MNNNCMM2/2)1(8/jWNN第4章 快速傅里叶变换(FFT)IRyxyxyxyxyxj)(j)(22)jj(22)j)(j1(22def)(22)(22)(22xyyxIyxR第4章 快速傅里叶变换(FFT)只需要两次实数加和两次实数乘就可实现。这样，对应的每个蝶形节省两次实数乘。在DIT-FFT运算流图中，从L=3至L=M级，每级都包含旋转因子 ，第L级中，对应N/2L个蝶形运算。因此从第三级至最后一级，旋转因子 节省的实数乘次数与(4.3.2)式相
34、同。所以从实数乘运算考虑，计算N=2M点DIT-FFT所需实数乘法次数为(4.3.4)8/NNW8/NNW8/NNW8/NNW134(3)22210222MNNRMNM第4章 快速傅里叶变换(FFT)在基2FFT程序中，若包含了所有旋转因子，则称该算法为一类碟形单元运算；若去掉的旋转因子，则称之为二类碟形单元运算；若再去掉的旋转因子，则称为三类碟形单元运算；若再判断处理，则称之为四类碟形运算。我们将后三种运算称为多类碟形单元运算。显然，碟形单元类型越多，编程就越复杂，但当N较大时，乘法运算的减少量是相当可观的。例如，N=4096时，三类碟形单元运算的乘法次数为一类碟形单元运算的75%。1mNW
35、 jWrN(1j)2/2mNW第4章 快速傅里叶变换(FFT)4.3.2 旋转因子的生成旋转因子的生成在FFT运算中，旋转因子，求正弦和余弦函数值的计算量是很大的。所以编程时，产生旋转因子的方法直接影响运算速度。一种方法是在每级运算中直接产生；另一种方法是在FFT程序开始前预先计算出，m=0，1，N/21，存放在数组中，作为旋转因子表，在程序执行过程中直接查表得到所需旋转因子值，不再计算。这样使运算速度大大提高，其不足之处是占用内存较多。)/2cos(NmWmN)/2sin(jNmmNW第4章 快速傅里叶变换(FFT)4.3.3 实序列的实序列的FFT算法算法在实际工作中，数据x(n)常常是实
36、数序列。如果直接按FFT运算流图计算，就是把x(n)看成一个虚部为零的复序列进行计算，这就增加了存储量和运算时间。处理该问题的方法有三种。早期提出的方法是用一个N点FFT计算两个N点实序列的FFT，一个实序列作为x(n)的实部，另一个作为虚部，计算完FFT后，根据DFT的共轭对称性，用例3.2.2 所述的方法由输出X(k)分别得到两个实序列的N点DFT。第二种方法是用N/2点FFT计算一个N点实序列的DFT。第三种方法是用离散哈特莱变换（DHT）1。下面简要介绍第二种方法。第4章 快速傅里叶变换(FFT)设x(n)为N点实序列，取x(n)的偶数点和奇数点分别作为新构造序列y(n)的实部和虚部，
37、即对y(n)进行N/2点FFT，输出Y(k)，则1210)(j)()(1210)12()()2()(2121NnnxnxnyNnnxnxnxnx，1122()DFT()()011()DFT()()2epopX kx nYkNkXkx njYk，第4章 快速傅里叶变换(FFT)根据DIT-FFT的思想及式(4.2.7)和(4.2.8)，可得到X(k)的前 个值：(4.3.5)式中，。由于x(n)为实序列，因此X(k)具有共轭对称性，X(k)的另外N/2点的值为12N210)()()(21NkkXWkXkXkN，)0(2),0(22211XNXXNX1210)()(*NkkXkNX，第4章 快速傅
38、里叶变换(FFT)计算 点FFT的复乘次数为 ，计算式(4.3.5)的复乘次数为，所以用这种算法,计算X(k)所需复数乘法次数为。相对一般的N点FFT算法，上述算法的运算效率为，当N=2M=210时，=20/11，运算速度提高近1倍。2N)1(4MN2N)1(42)1(4MNNMN)1/(2)1(4/2MMMNMN第4章 快速傅里叶变换(FFT)4.4 4.4 其他快速算法简介其他快速算法简介快速傅里叶变换算法是信号处理领域重要的研究课题。由于教材篇幅和教学大纲所限，本章仅介绍算法最简单、编程最容易的基2FFT算法原理及其编程思想，使读者建立快速傅里叶变换的基本概念，了解研究FFT算法的主要途
39、径和编程思路。其他高效快速算法请读者参考文献1、3、12。例如，分裂基FFT算法、离散哈特莱变(换DHT)、基4FFT、基8FFT、基rFFT、混合基FFT，以及进一步减少运算量的途径等内容，对研究新的快速算法都是很有用的。本节简要介绍其他几种快速算法的运算量及其主要特点，以便读者选择快速算法时参考。第4章 快速傅里叶变换(FFT)从理论上讲，不同基数的FFT算法的运算效率不同，实际中最常用的是基2FFT、基4 FFT、分裂基FFT和DHT 1。为此，下面简要介绍后三种FFT算法的特点和运算效率，以扩展读者的视野。其具体算法请参考文献1、12。在基rFFT算法中，基4FFT算法运算效率与基8F
40、FT很接近，但基4FFT算法实现程序简单，且判断开销少。可以证明，当FFT的基大于4时，不会明显降低计算量。基4FFT要求N=4M，M为自然数。其复数乘法次数为12（4.4.1）3lb()8NNCM(基4)第4章 快速傅里叶变换(FFT)其中未计入乘以j和1的计算。比较基2FFT的复数乘法次数 ，基4FFT的复数乘法次数减少25%。1984年，法国的杜梅尔（P.Dohamel）和霍尔曼（H.Hollmann）将基2分解和基4分解糅合，提出了分裂基FFT算法，其复数乘法次数接近FFT理论最小值，但其运算流图却与基2FFT很相似，编程简单，运算程序也很短，是一种很实用的高效算法。分裂基FFT算法复
41、数乘法次数为2 CM(分裂基)=（4.4.2）NNCMlb21)2(22lb()1399MNNN 第4章 快速傅里叶变换(FFT)只考虑（4.4.2）式的第一项，分裂基FFT算法的复数乘法次数就比基2FFT减少33%，比基4FFT减少11%。应当说明，在比较时，未考虑（4.4.2）式后2项减少的运算量，所以分裂基FFT算法的效率更高。由以上比较可见，分裂基FFT算法的效率最高，所以得到广泛应用。但是，对实序列x(n)，上述各种FFT算法仍将其看成虚部为零的复序列存储和计算。而一次复数乘法需要四次实数乘法和二次实数加法。所以，必然浪费存储资源和增加多余的运算量。我们知道，实序列的N点DFT具有共
42、轭对称性，即*()(),0,1,2,1X NkXkkN第4章 快速傅里叶变换(FFT)所以，只要计算出X（k）的前面N/2个值，则其后面的N/2个值可以由对称性求得。因此，FFT算法得到的N个X（k）值有一半是多余的。由以上分析可见，对实序列一定存在更高效的快速算法。离散哈特莱变换（DHT）就是针对实序列的一种高效变换算法，相对一般的FFT算法，DHT的快速算法FHT可以减少近一半的计算量1。N点基2时域抽取快速DHT（基2DIT-FHT）算法的实数乘法次数为1 （4.4.3）FHT34MNMN第4章 快速傅里叶变换(FFT)N点基2DIT-FHT算法的实数加法次数为1 （4.4.4）由式（4
43、.4.3）可见，基2DIT-FHT算法的实数乘法次数约为基2DIT-FFT算法的一半。与前面三种FFT算法比较，对实序列，基2DIT-FHT算法的实数乘法次数最少。应当说明，DHT是与DFT不同的变换，所以要想得到实序列DFT，还要根据二者的关系式进行转换。下面会看到该关系非常简单。FHT1322MAN第4章 快速傅里叶变换(FFT)DHT具有以下主要优点：(1)DHT是实数变换，在对实信号进行处理时避免了复数运算，运算效率高，且实现硬件简单经济。(2)DHT的正、逆变换（除了因子1/N外）具有相同的形式，所以实现硬件或程序亦相同。N点DHT定义如下：（4.4.5）1H02()DHT()()c
44、os()0,1,2,1NNnXkx nx nknkNN第4章 快速傅里叶变换(FFT)N点逆DHT变换定义为 （4.4.6）(3)DHT满足循环卷积定理，所以，可以直接用FHT实现实序列的快速卷积，大大提高处理速度1，并使处理硬件简化。1HH012()IDHT()()cos()0,1,2,1NNkx nXkXkknnNNN第4章 快速傅里叶变换(FFT)(4)DHT与DFT之间的关系非常简单，容易实现二者之间的转换，关系式如下：（4.4.7）所以对实信号x(n)进行谱分析时，可以先对x(n)进行FHT,得到XH(k)=DHTx(n)N，然后再将H(k)转换成X(k)=DFTx(n)N，这样可以
45、提高分析速度，减少存储空间。HHHH11()()()j()()22X kXkXNkXkXNk第4章 快速傅里叶变换(FFT)习题与上机题习题与上机题1 如果某通用单片计算机的速度为平均每次复数乘需要4 s，每次复数加需要1 s，用来计算N=1024点DFT，问直接计算需要多少时间。用FFT计算呢？照这样计算，用FFT进行快速卷积来对信号进行处理时，估计可实现实时处理的信号最高频率。2 如果将通用单片机换成数字信号处理专用单片机TMS320系列，计算复数乘和复数加各需要10 ns。请重复做上题。第4章 快速傅里叶变换(FFT)3 已知X(k)和Y(k)是两个N点实序列x(n)和y(n)的DFT，
46、希望从X(k)和Y(k)求x(n)和y(n),为提高运算效率，试设计用一次N点IFFT来完成的算法。4 设x(n)是长度为2N的有限长实序列，X(k)为x(n)的2N点DFT。(1)试设计用一次N点FFT完成计算X(k)的高效算法；(2)若已知X(k)，试设计用一次N点IFFT实现求X(k)的2N点IDFT运算。第4章 快速傅里叶变换(FFT)5 分别画出16点基2DIT-FFT和DIF-FFT运算流图，并计算其复数乘次数，如果考虑三类碟形的乘法计算，试计算复乘次数。6*按照下面的IDFT算法编写MATLAB语言IFFT程序，其中的FFT部分不用写出清单，可调用fft函数。并分别对单位脉冲序列、矩形序列、三角序列和正弦序列进行FFT和IFFT，验证所编程序。*1()IDFT()DFT()x nX kXkN
