数字信号第4章

第4章快速傅里叶变换(FFT)第4章快速傅里叶变换(FFT)4.1引言4.2基2FFT算法4.3进一步减少运算量的措施4.4其他快速算法简介习题与上机题第4章快速傅里叶变换(FFT)4.1引言DFT是数字信号分析与处理中的一种重要变换。但直接计算DFT的计算量与变换区间长度N的平方成正比，当N较大时，计算量太大，所以在快速傅里叶变换FFT(FastFourierTransform)出现以前，直接用DFT算法进行谱分析和信号的实时处理是不切实际的。直到1965年提出DFT的一种快速算法以后，情况才发生了根本的第4章快速傅里叶变换(FFT)自从1965年库利(T.W.Cooley)和图基(J.W.Tuky)在《计算数学》(Math.Computation,Vol.19,1965)杂志上发表了著名的《机器计算傅里叶级数的一种算法》论文后，桑德(G.Sand)—图基等快速算法相继出现，又经人们进行改进，很快形成一套高效计算方法，这就是现在的快速傅里叶变换（FFT）。这种算法使DFT的运算效率提高了1~2个数量级，为数字信号处理技术应用于各种信号的实时处理创造了条件，大大推动了数字信号处理技术的发展。第4章快速傅里叶变换(FFT)人类的求知欲和科学的发展是永无止境的。多年来，人们继续寻求更快、更灵活的好算法。1984年，法国的杜哈梅尔(P.Dohamel)和霍尔曼(H.Hollmann)提出的分裂基快速算法，使运算效率进一步提高。本章主要讨论基2FFT第4章快速傅里叶变换(FFT)4.2基2FFT4.2.1直接计算DFT的特点及减少运算量的有限长序列x(n)的N点DFT为(4.2.1)考虑x(n)为复数序列的一般情况，对某一个k值，直接按(4.2.1)式计算X(k)的1个值需要N次复数乘法和(N－1)次复数加法。因此，计算X(k)的所有N个值，共需N2次复数乘法和N(N－1)次复数加法运算。10110)()(NnknNNkWnxkX，，，第4章快速傅里叶变换(FFT)1N当时，N(N－1)≈N2。由上述可见，N点DFT的乘法和加法运算次数均为N2。当N较大时，运算量相当可观。例如N=1024时，N2=1048576。这对于实时信号处理来说，必将对处理设备的计算速度提出难以实现的要求。所以，必须减少其运算量，才能使DFT在各如前所述，N点DFT的复乘次数等于N2。显然，把N点DFT分解为几个较短的DFT，可使乘法次数大大减少。另外，旋转因子具有明显的周期性和对称性。其周期性表现为（4.2.2）mNmNlNmNlNmNWWπ2j)(π2jee第4章快速傅里叶变换(FFT)(4.2.3a)(4.2.3b)FFT算法就是不断地把长序列的DFT分解成几个短序列的DFT，并利用的周期性和对称性来减少DFT的运算次数。算法最简单最常用的是基2FFT。其对称性表现为mNNmNWW或者mNmNNWW*][mNNmNWW2knNW第4章快速傅里叶变换(FFT)4.2.2时域抽取法基2FFT基本原理基2FFT算法分为两类：时域抽取法FFT(DecimationInTimeFFT，简称DIT－FFT)；频域抽取法FFT(DecimationInFrequencyFFT，简称DIF－FFT)。本节介绍DIT－FFT设序列x(n)的长度为N，且满足N=2M，M为自然数。按n的奇偶把x(n)分解为两个N/2点的子序列12()(2)0112()(21)0112NxrxrrNxrxrr，，，，，，，，第4章快速傅里叶变换(FFT)则x(n)的DFT为/21/212(21)00/21/21221200()()()(2)(21)()()knknNNnnNNkrkrNNrrNNkrkkrNNNrrXkxnWxnWxrWxrWxrWWxrW偶数奇数第4章快速傅里叶变换(FFT)因为所以(4.2.4)2π2πjj222/2eekrkrkrkrNNNNWW/21/211/22/20012()()()()()0,1,2,,-1NNkrrkrNNNrrkNXkxrWWxrWXkWXkkN第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)又可表示为/2111/2120()()DFT[()]NkrNNrXkxrWxr/2122/2220()()DFT[()]NkrNNrXkxrWxrkNNkNWW2第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所示的流图符号表示，称为蝶形运算符号。采用这种图示法，经过一次奇偶抽取分解后，N点DFT运算图可以用图4.2.2表示。图中，N=23=8，X(0)~X(3)由(4.2.7)式给出，而X(4)~X(7)则由(4.2.8)式给出。1210)()()(21NkkXWkXkXkN，，，，1210)()()2(21NkkXWkXNkXkN，，，，第4章快速傅里叶变换(FFT)图4.2.1蝶形运算符号第4章快速傅里叶变换(FFT)图4.2.28点DFT一次时域抽取分解运算流图第4章快速傅里叶变换(FFT)由图4.2.1可见，要完成一个蝶形运算，需要一次复数乘法和两次复数加法运算。由图4.2.2容易看出，经过一次分解后，计算1个N点DFT共需要计算两个N/2点DFT和N/2个蝶形运算。而计算一个N/2点DFT需要(N/2)2次复数乘法和N/2(N/2－1)次复数加法。所以，按图4.2.2计算N点DFT时，总共需要的复数乘法次数为221(1)22222NNNNNN第4章快速傅里叶变换(FFT)复数加法次数为由此可见，仅仅经过一次分解，就使运算量减少近一半。既然这样分解对减少DFT的运算量是有效的，且N=2M,N/2仍然是偶数，故可以对N/2点DFT再作进一步分解。与第一次分解相同，将x1(r)按奇偶分解成两个N/4点的子序列x3(l)和x4(l)，即222122NNNN1410)12()()2()(1423Nllxlxlxlx，，，第4章快速傅里叶变换(FFT)X1(k)又可表示为(4.2.9)1210)()()()()12()2()(42/314/04/42/14/04/314/0)12(2/114/022/11NkkXWkXWlxWWlxWlxWlxkXkNNlklNkNNlklNNllkNNlklN，，，第4章快速傅里叶变换(FFT)式中同理，由X3(k)和X4(k)的周期性和的对称性最后得到：(4.2.10)/4133/4340()()DFT[()]NklNNlXkxlWxl/4144/4440()()DFT[()]NklNNlXkxlWxlmNW2/kNNkNWW2/4/2/14/10)()()4/()()()(42/3142/31NkkXWkXNkXkXWkXkXkNkN，，，，第4章快速傅里叶变换(FFT)用同样的方法可计算出(4.2.11)其中1410)()(4)()()(62/5262/52NkkXWkXNkXkXWkXkXkNkN，，，/4155/4540/4166/46405262()()DFT[()]()()DFT[()]()(2)01/41()(21)NklNNlNklNNlXkxlWxlXkxlWxlxlxllNxlxl，，，，第4章快速傅里叶变换(FFT)这样，经过第二次分解，又将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点DIT－FFT运算流图如图4.2.4所示。图中用到关系式。图中输入序列不是顺序排列，但后面会看到，其排列是有规律的。图中的数组A用于存放输入序列和每级运算结果，在后面讨论编程方法和倒序时要用到。mkNkmNWW/第4章快速傅里叶变换(FFT)图4.2.38点DFT二次时域抽取分解运算流图第4章快速傅里叶变换(FFT)图4.2.48点DIT-FFT运算流图第4章快速傅里叶变换(FFT)4.2.3DIT-FFT算法与直接计算DFT运算量的比较由DIT-FFT算法的分解过程及图4.2.4可见，N=2M时，其运算流图应有M级蝶形，每一级都由N/2个蝶形运算构成。因此，每一级运算都需要N/2次复数乘和N次复数加(每个蝶形需要两次复数加法)。所以，M级运算总共需要的复数乘次数为复数加次数为lb22MNNCMNlbACNMNN第4章快速傅里叶变换(FFT)而直接计算DFT的复数乘为N2次，复数加为N(N－1)次。当N1时，N2(N/2)lbN，所以，DIT-FFT算法比直接计算DFT的运算次数大大减少。例如，N=210=1024时，这样，就使运算效率提高200多倍。图4.2.5为FFT算法和直接计算DFT所需复数乘法次数CM与变换点数N的关系曲线。由此图更加直观地看出FFT算法的优越性，显然，N8.20451205760481lb22NNN第4章快速傅里叶变换(FFT)图4.2.5DIT-FFT算法与直接计算DFT所需复数乘法次数的比较曲线第4章快速傅里叶变换(FFT)4.2.4DIT-FFT的运算规律及编程思想1．由图4.2.4可以看出，DIT-FFT的运算过程很有规律。N=2M点的FFT共进行M级运算，每级由N/2个蝶形运算组成。同一级中，每个蝶形的两个输入数据只对计算本蝶形有用，而且每个蝶形的输入、输出数据结点又同在一条水平线上，这就意味着计算完一个蝶形后，所得输出数据可立即存入原输入数据所占用的存储单元(数组元素)。这样，经过M级运算后，原来存放输入序列数据的N个存储单元(数组A)中便依次存放X(k)的N个值。这种利用同一存储单元存储蝶形计算输入、输出数据的方法称为原位(址)计算。第4章快速傅里叶变换(FFT)2．旋转因子的变化规律如上所述，N点DIT-FFT运算流图中，每级都有N/2个蝶形。每个蝶形都要乘以因子，称其为旋转因子，p为旋转因子的指数。但各级的旋转因子和循环方式都有所不同。为了编写计算程序，应先找出旋转因子与运算级数的关系。用L表示从左到右的运算级数(L=1，2，…，M)。观察图4.2.4不难发现，第L级共有2L－1个不同的旋转因子。N=23=8时的各级旋转因子表示pNWpNW第4章快速傅里叶变换(FFT)对N=2M的一般情况，第L级的旋转因子为3,2,1,031,0201222/24/J章快速傅里叶变换(FFT)因为所以(4.2.12)(4.2.13)这样，就可按(4.2.12)和(4.2.13)式确定第L级运算的旋转因子(实际编程序时，L为最外层循环变量)L-12,0,1,2,,21LpJNWWJMLMLMLN222212,2,1,0122LJNJNpNJ2第4章快速傅里叶变换(FFT)3．设序列x(n)经时域抽选(倒序)后，按图4.2.4所示的次序(倒序)存入数组A中。如果蝶形运算的两个输入数据相距B个点，应用原位计算，则蝶形运算可表示成如下形式：式中下标L表示第L级运算，AL(J)则表示第L级运算后的数