92数字信号处理讲义--第9章 离散傅里叶变换的计算

第９章离散傅里叶变换的计算教学目的1．了解直接计算DFT与序列长度N的关系，理解开发高效算法的的意义；2．理解高效算法的原理与实现方法，掌握按时间和按频率抽取的FFT算法原理和实现；3．掌握卷积实现DFT的方法；4．掌握复合数N的FFT算法；5．了解离散傅丽叶变换中有限寄存器长度的影响。教学重点与难点重点：本章是本课程的另一个重中这重。1．高效算法的原理与实现方法，掌握按时间和按频率抽取的FFT算法原理和实现；2．卷积实现DFT的方法．；3．离散傅丽叶变换中有限寄存器长度的影响。难点：1．掌握按时间和按频率抽取的FFT算法原理和实现；2．卷积实现DFT的方法。9.0引言离散傅里叶变换（DFT）有一种快速算法，我们平常称为快速傅里叶变换（FFT）。由于有限长序列在其频域也可离散化为有限长序列（DFT），因此离散傅里叶变换（DFT）在数字信号处理中是非常有用的。例如，在信号的频谱分析、系统的分析、设计和实现中都会用到DFT的计算。但是，在相当长的时间里，由于DFT的计算量太大，即使采用计算机也很难对问题进行实时处理，所以并没有得到真正的运用。直到1965年首次发现了DFT运算的一种快速算法以后，情况才发生了根本的变化。人们开始认识到DFT运算的一些内在规律，从而很快地发展和完善了一套高速有效的运算方法，这就是现在人们普遍称之为快速傅里叶变换（FFT）的算法。FFT出现后使DFT的运算大大简化，运算时间一般可缩短一二个数量级之多，从而使DFT的运算在实际中真正得到了广泛的应用。9.1离散傅里叶变换的高效计算9.1.1直接计算DFT的运算量问题设x(n)为N点有限长序列，其DFT为k=0,1,…,N-1（9-1）10)()(NnnkNWnxkX反变换（IDFT）为n=0,1,…,N-1（9-2）二者的差别只在于WN的指数符号不同，以及差一个常数乘因子1/N，所以IDFT与DFT具有相同的运算工作量。下面我们只讨论DFT的运算量。一般来说，x(n)和WNnk都是复数，X(k)也是复数，因此每计算一个X(k)值，需要N次复数乘法和N-1次复数加法。而X(k)一共有N个点（k从0取到N-1），所以完成整个DFT运算总共需要N2次复数乘法及N(N-1)次复数加法。在这些运算中乘法运算要比加法运算复杂，需要的运算时间也多一些。因为复数运算实际上是由实数运算来完成的，这时DFT运算式可写成由此可见，一次复数乘法需用四次实数乘法和二次实数加法；一次复数加法需二次实数加法。因而每运算一个X(k)需4N次实数乘法和2N+2(N-1)=2(2N-1)次实数加法。所以，整个DFT运算总共需要4N2次实数乘法和2N(2N-1)次实数加法。当然，上述统计与实际需要的运算次数稍有出入，因为某些WNnk可能是1或j，就不必相乘了，例如W0N=1，WNＮ／２=-1，WNN/4=-j等就不需乘法。但是为了便于和其他运算方法作比较，一般都不考虑这些特殊情况，而是把WNnk都看成复数，当N很大时，这种特例的影响很小。从上面的统计可以看到，直接计算DFT，乘法次数和加法次数都是和N2成正比的，当N很大时，运算量是很可观的，有时是无法忍受的。例9-1根据式（9-1），对一幅N×N点的二维图像进行DFT变换，如用每秒可做10万次复数乘法的计算机，当N=1024时，问需要多少时间（不考虑加法运算时间）？解直接计算DFT所需复乘次数为（N2)2≈1012次，因此用每秒可做10万次复数乘法的计算机，则需要近3000小时。10)(1)(NknkNWkXNnX])}Re[)](Im[]Im[)]((Re[][Im)](Im[]Re[)]({Re[]}Im[])]}{Re[(Im[)]({Re[)()(101010nkNnkNnkNNnnkNNnnkNnkNNnnkNWnxWnxjWnxWnxWjWnxjnxWnxkX这对实时性很强的信号处理来说，要么提高计算速度，而这样，对计算速度的要求太高了。另外，只能通过改进对DFT的计算方法，以大大减少运算次数。9.1.2改善途径能否减少运算量，从而缩短计算时间呢？仔细观察DFT的运算就可看出，利用系数WＮnk的以下固有特性，就可减少运算量：（1）WNnk的对称性2）WNnk的周期性（3）WNnk的可约性另外这样，利用这些特性，使DFT运算中有些项可以合并，并能使DFT分解为更少点数的DFT运算。而前面已经说到，DFT的运算量是与N2成正比的，所以N越小越有利，因而小点数序列的DFT比大点数序列的DFT的运算量要小。快速傅里叶变换算法正是基于这样的基本思想而发展起来的。它的算法形式有很多种，但基本上可以分成两大类，即按时间抽取(ＤecimationinTime，缩写为DIT）法和按频率抽取（Decimationinrequency，缩写为DIF）法。9.2按时间抽取（DIT）的FFT算法9.2.1算法原理设序列x(n)长度为N，且满足N=2M，M为正整数。按n的奇偶把x(n)分解为两个N/2点的子序列：（9-4）则可将DFT化为nkNnkNWW*)()()(NknNkNnNnkNmnkmNnkNnmkmNnkNkNNkNNNnkNknNNkNnN)2/(2/)()(,1,12,,1,0)()12()()2(21NrrxrxrxrxrkNNrkNrkNNrkrNNrrkNNrNnnnkNNnNnnnkNnkNWrxWWrxWrxWrxWnxWnxWnxnxDFTkX))(()()()12()2()()()()]([)(2120221201)12(1202120101010为奇数为偶数由于，故上式可表示成（9-5）式中，X1(k)与X2(k)分别是x1(r)及x2(r)的N/2点DFT：（9-6）（9-7）由此，我们可以看到，一个N点DFT已分解成两个N/2点的DFT这两个N/2点的DFT再按照式（9-5）组合成一个N点DFT。这里应该看到X1(k)，X2(k)只有N/2个点，即k=0,1,…,N/2-1。而X(k)却有N个点，即k=0,1,…,N-1，故用式（9-5）计算得到的只是X(k)的前一半的结果，要用X1(k)，X2(k)来表达全部的X(k)值，还必须应用系数的周期性，即这样可得到（9-8）同理可得（9-9）式（9-8）、式（9-9）说明了后半部分k值（N/2≤k≤N-1）所对应的X1(k)，Ｘ2(k)分别等于前半部分k值（0≤k≤N/2-1）所对应的X1(k)，X2(k)。再考虑到WkN的以下性质（9-10）2/22222NjNjNWeeW)()()()()(212/21202/1120kXWkXWrxWWrxkXkNrkNNrkNrkNNrrkNNrrkNNrrkNNrrkNNrWrxWrxkXWrxWrxkX2/1202/212022/1202/11201)12()()()2()()(rkNNkrNWW2/22/)()()(212/120122/12011kXWrxWrxkNXrkNNrkNrNNr)(222kXkNXkNkNNNkNN2/2这样，把式（9-8）、式（9-9）、式（9-10）代入式（9-5），就可将X(k)表达为前后两部分：（9-11）（9-12）图9-1时间抽取法蝶形运算流图符号图9-2按时间抽取将一个N点DFT分解为两个N/2点DFT（N=8）一个N点DFT分解为两个N/2点DFT，每一个N/2点DFT只需(N/2)2=N2/4次复数乘法，N/2(N/2-1)次复数加法。两个N/2点DFT共需2×(N/2)2=N2/2次复数乘法和N(N/2-1)次复数加法。此外，把两个N/2点DFT合成为N点DFT时，有N/2个蝶形运算，还需要N/2次复数乘法及2×N/2=N次复数加法。因而通过第一步分解后，总共需要(N2/2)+(N/2)=N(N+1)/2≈N2/2次复数乘法和N(N/2-1)+N=N2/2次复数加法。由此可见，通过这样分解后运算工作量差不多节省了一半。既然这样分解是有效的，由于N=2M，因而N/2仍是偶数，可以进一步把每个N/2点子序列再按其奇偶部分分解为两个N/4点的子序列。12,,1,0)()()(21NkkXWkXkXkN)()(22221221kXWkXNkXWNkXNkXkNNkN12,,1,0NkX2(k)X1(k)kNW－1X1(k)＋kNWX2(k)X1(k)－kNWX2(k)X1(0)X1(1)x1(0)＝x(0)X(0)X(1)X1(2)X1(3)－1－10NWDFT2点Nx1(1)＝x(2)x1(2)＝x(4)x1(3)＝x(6)X(2)X(3)X2(0)X2(1)x2(0)＝x(1)X(4)X(5)X2(2)X2(3)DFT2点Nx2(1)＝x(3)x2(2)＝x(5)x2(3)＝x(7)X(6)X(7)1NW2NW3NW－1－1（9-13）且式中：（9-14）（9-15）图9-3给出N=8时，将一个N/2点DFT分解成两个N/4点DFT，由这两个N/4点DFT组合成一个N/2点DFT的流图。图9-3N/2点DFT分解为两个N/4点DFTX2(k)也可进行同样的分解：式中：14,,1,0)()12()()2(4131Nllxlxlxlx14,,1,0)()()()()12()2()(42/31404/42/1404/3140)12(2/114022/11NkkXWkXWlxWWlxWlxWlxkXkNNllkNkNNllkNNlklNNllkN)()(442/31kXWkXkNXkN14,,1,0Nk1404/441404/33)()()()(NllkNNllkNWlxkXWlxkXDFT4点NX3(0)X3(1)x3(0)＝x1(0)＝x(0)x3(1)＝x1(2)＝x(4)X1(0)X1(1)DFT4点NX4(0)X4(1)x4(0)＝x1(1)＝x(2)x4(1)＝x1(3)＝x(6)X1(2)X1(3)－1－102/NW12/NW)()(4)()()(62/5262/52kXWkXkNXkXWkXkXkNkN14,,1,0Nk1404/21404/661404/21404/55)12()()()2()()(NllkNNllkNNllkNNllkNWlxWlxkXWlxWlxkX图9-4一个N=8点DFT分解为四个N/4点DFT根据上面同样的分析知道，利用四个N/4点的DFT及两级蝶形组合运算来计算N点DFT，比只用一次分解蝶形组合方式的计算量又减少了大约一半。最后，剩下的是2点DFT，对于此例N=8，就是四个N/4=2点DFT，其输出为X3(k)，X4(k)，X5(k)，X6(k)，k=0,1，这由式（9-14）～式（9-17）四个式子可以计算出来。例如，由式（9-14）可得：式中，,故上式不需要乘法。类似地,可求出X4(k)，X5(k)，X6(k)。这种方法的每一步分解，都是按输入序列在时间上的次序是属于偶数还是属于奇数来分解为两个更短的子序列，所以称为“按时间抽取法”。DFT4点NX3(0)X3(1)x3(0)＝x1(0)＝x(0)x3(1)＝x1(2)＝x(4)X1(0)X1(1)DFT4点NX4(0)X4(1)x4(0)＝x1(1)＝x(2)x4(1)＝x1(3)＝x(6)X1(2)X1(3)－1－10NW2NWX(0)X(