數(shù)字信號處理第4章

第4章快速傅里葉變換(FFT)

4.1引言4.2基2FFT算法4.3進一步減少運算量的措施4.4其他快速算法簡介習題與上機題4.1引言

DFT是數(shù)字信號分析與處理中的一種重要變換

直接計算N點DFT：與變換區(qū)間長度N的平方成正比

N較大時，計算量太大——難以進行實時處理快速傅里葉變換FFT(FastFourierTransform)

高效計算DFT的算法

使DFT運算效率提高1~2個數(shù)量級4.2基2FFT算法4.2.1直接計算DFT的特點及減少運算量基本途徑一、直接計算DFT

有限長序列x(n)的N點DFT為

x(n)為復數(shù)序列的一般情況：計算X(k)的1個值：需要

次復數(shù)乘法和次復數(shù)加法；計算X(k)的所有N個值：共需N2次復數(shù)乘法和N(N－1)次復數(shù)加法運算

（N>>1時，N(N－1)=N2

）。(4.2.1)N(N－1)N較大時，運算量很大：例如N=1024時，N2=1048576實時處理信號難以實現(xiàn)，對運算速度要求高二、減少運算量的基本途徑：1、N較大，運算量大減小N，把N點DFT分解為幾個較短的DFT

2、旋轉(zhuǎn)因子具有周期性和對稱性周期性：

對稱性：（4.2.2）或者(4.2.3a)(4.2.3b)2：利用的周期性和對稱性來減少DFT的運算次數(shù)減少運算量的基本途徑：1：不斷地把長序列的DFT分解成幾個短序列的DFT4.2.2時域抽取法基2FFT基本原理基2FFT算法：

頻域抽取法FFT(DIF－FFT)時域抽取法：設(shè)序列x(n)的長度為N，且滿足N=2M，M為自然數(shù)。

（若不滿足可補零）按n的奇偶分解x(n)為兩個N/2點的子序列如下：時域抽取法FFT(DIT－FFT)則x(n)的N點DFT為而(4.2.4)其中Ｘ1(k)和X2(k)分別為x1(r)和x2(r)的N/2點DFT：(4.2.5)(4.2.6)所以由于X1(k)和X2(k)均以N/2為周期，且(4.2.7)(4.2.8)1、X1(k)和X2(k)兩個N/2點DFT2、上式的運算圖4.2.1蝶形運算符號因此X(k)又可表示為:計算N點DFT分解為：

——可用蝶形運算符號表示圖4.2.28點DFT一次時域抽取分解運算流圖要完成一個蝶形運算，需要和運算經(jīng)過一次分解后，計算1個N點DFT共需要計算：兩個N/2點DFT；

計算N點DFT時：總共需要的復數(shù)乘法次數(shù)為復數(shù)加法次數(shù)為結(jié)論：經(jīng)過一次分解，運算量減少近一半一次復數(shù)乘法兩次復數(shù)加法N/2個蝶形運算因為N=2M，N/2仍然是偶數(shù)，可對N/2點DFT再作進一步分解。同上，將x1(r)按奇偶分解成兩個N/4點的子序列x3(l)和x4(l)：X1(k)又可表示為(4.2.9)同理，X3(k)和X4(k)以N/4為周期，，可得：(4.2.10)用同樣的方法可計算出(4.2.11)式中：其中：經(jīng)過第二次分解，N/2點DFT分解為：2個N/4點DFT和N/4個蝶形運算，如圖4.2.3依此類推，經(jīng)過M次分解，N點DFT分解成：N個1點DFT和M級蝶形運算，而1點DFT就是時域序列本身完整的8點DIT－FFT運算流圖如圖4.2.4所示圖4.2.38點DFT二次時域抽取分解運算流圖圖中：1、用到關(guān)系式：2、輸入序列不是順序排列，但有排列規(guī)律3、數(shù)組A用于存放輸入序列和每級運算結(jié)果圖4.2.48點DIT-FFT運算流圖4.2.3DIT-FFT算法與直接計算DFT運算量的比較由DIT-FFT算法的分解過程及圖4.2.4可見，N=2M

時，其運算流圖應有M級蝶形，每一級都由N/2個蝶形運算構(gòu)成。因此，每一級運算都需要N/2次復數(shù)乘和N次復數(shù)加(每個蝶形需要兩次復數(shù)加法)。所以，M級運算總共需要的復數(shù)乘次數(shù)為復數(shù)加次數(shù)為而直接計算DFT的復數(shù)乘為N2次，復數(shù)加為N(N－1)次。當N>>1時，N2>>(N/2)lbN，所以，DIT-FFT算法比直接計算DFT的運算次數(shù)大大減少。例如，N=210=1024時，這樣，就使運算效率提高200多倍。圖4.2.5為FFT算法和直接計算DFT所需復數(shù)乘法次數(shù)CM與變換點數(shù)N的關(guān)系曲線。由此圖更加直觀地看出FFT算法的優(yōu)越性，顯然，N越大時，優(yōu)越性就越明顯。圖4.2.5DIT-FFT算法與直接計算DFT所需復數(shù)乘法次數(shù)的比較曲線4.2.4DIT-FFT的運算規(guī)律及編程思想

1．原位計算DIT-FFT的運算規(guī)律：(1)N=2M，共進行M級運算，每級由N/2個蝶形運算組成；(2)同一級中，每個蝶形的兩個輸入數(shù)據(jù)只對該蝶形有用，

每個蝶形的輸入、輸出數(shù)據(jù)在一條水平線上——輸入、輸出可用同一存儲單元即：第一級：N個存儲單元存儲N個序列值x(n)

M級蝶形運算后，該N個存儲單元中依次存放X(k)的N個值原位計算：利用同一存儲單元存儲蝶形計算輸入、輸出數(shù)據(jù)——可節(jié)省大量內(nèi)存，從而使設(shè)備成本降低。

2．旋轉(zhuǎn)因子的變化規(guī)律旋轉(zhuǎn)因子：N點DIT-FFT運算流圖中，每個蝶形都要乘以因子，稱其為旋轉(zhuǎn)因子，p為旋轉(zhuǎn)因子的指數(shù)。各級旋轉(zhuǎn)因子和循環(huán)方式都有所不同，需找出旋轉(zhuǎn)因子與運算級數(shù)的關(guān)系。用L表示從左到右的運算級數(shù)(L=1，2，…，M)。第L級共有2L－1個不同的旋轉(zhuǎn)因子。N=23=8時的各級旋轉(zhuǎn)因子：因為所以

(4.2.12)

(4.2.13)——按(4.2.12)和(4.2.13)式可確定第L級運算的旋轉(zhuǎn)因子(實際編程序時，L為最外層循環(huán)變量)。一般的，N=2M時，第L級的旋轉(zhuǎn)因子為

3．蝶形運算規(guī)律設(shè)序列x(n)經(jīng)時域抽選(倒序)后，按圖4.2.4所示的次序(倒序)存入數(shù)組A中。如果蝶形運算的兩個輸入數(shù)據(jù)相距B個點，應用原位計算，則蝶形運算可表示成如下形式：式中下標L：第L級運算AL(J)：第L級運算后A(J)的值（第L級蝶形輸出）AL－1(J)：第L級運算前A(J)的值（第L級蝶形輸入）

4.編程思想及程序框圖

歸納運算規(guī)律：第L級中，每個蝶形的兩個輸入數(shù)據(jù)相距B=2L－1個點；每級有B個不同的旋轉(zhuǎn)因子；同一旋轉(zhuǎn)因子對應著間隔為2L點的2M－L個蝶形。總結(jié)上述運算規(guī)律，可采用下述運算方法:(1)從輸入端(第1級)開始，逐級進行，共進行M級運算;(2)進行第L級運算時，依次求出B個不同的旋轉(zhuǎn)因子，同時計算其對應的2M－L個蝶形。圖4.2.6DIT-FFT運算和程序框圖L為最外層循環(huán)變量求B個不同的旋轉(zhuǎn)因子計算每個旋轉(zhuǎn)因子所對應的2M－L個蝶形DIT-FFT算法運算流圖：輸出X(k)為自然順序，但輸入序列不是按x(n)的自然順序排列。

這種經(jīng)過M次偶奇抽選后的排序稱為序列x(n)的倒序(倒位)。因此，在運算M級蝶形之前應先對序列x(n)進行倒序。

5．序列的倒序DIT-FFT算法輸入序列的排序：倒序的規(guī)律：N=2M，順序數(shù)可用M位二進制數(shù)(nM－1nM－2…n1n0)表示。例如N=8，用n2n1n0表示x(0)-x(7)：000-111M次偶奇時域抽選過程如圖4.2.7所示。圖4.2.7形成例序的樹狀圖（N=23）(1)第一次按最低位n0的0和1將x(n)分解為偶奇兩組；(2)第二次按次低位n1的0、1值分別對偶奇組分組；(3)依次類推，第M次按nM－1位分解，最后所得二進制倒

序數(shù)如圖4.2.7所示。表4.2.1順序和倒序二進制數(shù)對照表表4.2.1列出了N=8時以二進制數(shù)表示的順序數(shù)和倒序數(shù)，由表可見，只要將順序數(shù)(n2n1n0)的二進制位倒置，則得對應的二進制倒序值(n0n1n2)?！布娐泛蛥R編語言程序容易實現(xiàn)(1)自然順序數(shù)I增加1，二進制數(shù)最低位加1；(2)相應的倒序數(shù)J則在M位二進制數(shù)最高位加1。例如，自然順序數(shù)0：000最低位加1：001

對應的倒序數(shù)0：000最高位加1：100自然順序數(shù)1：001最低位加1：010

對應的倒序數(shù)4：100最高位為1：010——最高位加1要向次高位進位，其實質(zhì)是將最高位變?yōu)?，再在次高位加1。用這種算法，可以從當前任一倒序值求得下一個倒序值但用有些高級語言程序?qū)崿F(xiàn)倒序時，直接倒置二進制數(shù)位是不行的，因此必須找出產(chǎn)生倒序數(shù)的十進制運算規(guī)律。由表4.2.1可見：若用J表示當前倒序數(shù)的十進制數(shù)值，對于N=2M，M位二進制數(shù)的權(quán)值從左向右依次為N/2，N/4，N/8，…，2，1。二進制最高位加1：十進制：J+N/2——最高位是0(J<N/2)，則直接加1(J

J+N/2)；最高位是1(J≥N/2),先將最高位變成0(J

J－N/2)

然后次高位加1(J+N/4)，(次高位加1時，同樣要判斷0、1值)

J+N/4——次高位為0(J<N/4)，則直接加1(JJ+N/4);

次高位為1(J≥N/4),先將次高位變成0(JJ－N/4)

然后下一位加1(J+N/8)，

再判斷下一位的0或1值，依此類推，直到完成最高位加1，逢2向右進位的運算。(圖4.2.9所示的虛線框為倒序的程序框圖)

形成倒序J后，將原數(shù)組A中存放的輸入序列重新按倒序排列。(1)原輸入序列x(I)先按自然順序存入數(shù)組A中。

對N=8，A(0)，A(1)，…，A(7)中依次存放著x(0)，x(1)，x(2)，…，x(7)(2)對x(n)重新排序(倒序)：規(guī)律如圖4.2.8所示：①第一個序列值x(0)和最后一個序列值x(N－1)不需要重排；②每計算出一個倒序值J，便與循環(huán)語句自動生成的順序I

比較：當I=J時，不需要交換數(shù)據(jù)；

當I≠J時,A(I)與A(J)交換數(shù)據(jù)。③為了避免再次調(diào)換前面已調(diào)換過的一對數(shù)據(jù)，只對I<J的

情況調(diào)換A(I)和A(J)的內(nèi)容。圖4.2.8倒序規(guī)律圖4.2.9倒序程序框圖4.2.5頻域抽取法FFT(DIF-FFT)設(shè)序列x(n)長度為N=2M，首先將x(n)前后對半分開，得到兩個子序列，其DFT如下：式中將X(k)分解成偶數(shù)組與奇數(shù)組：當k取偶數(shù)(k=2m,m=0,1,…,N/2－1)時(4.2.14)當k取奇數(shù)(k=2m+1,m=0,1,…,N/2－1)時，令將x1(n)和x2(n)分別代入(4.2.14)和(4.2.15)式，可得

(4.2.16)——X(k)按k值的奇偶分為兩組：偶數(shù)組是x1(n)的N/2點DFT，

奇數(shù)組是x2(n)的N/2點DFT。x1(n)、x2(n)和x(n)之間的關(guān)系也可用蝶形運算流圖符號表示圖4.2.11DIF-FFT第一次分解運算流圖（N=8）由于N=2M，N/2仍然是偶數(shù)，繼續(xù)將N/2點DFT分成偶數(shù)組合奇數(shù)組，這樣每個N/2點DFT又可由兩個N/4點DFT形成，其輸入序列分別是x1(n)和x2(n)按上下對半分開形成的四個子序列。圖4.2.12示出了N=8時第二次分解運算流圖。

這樣繼續(xù)分解下去，經(jīng)過M－1次分解，最后分解為2M－1個兩點DFT，兩點DFT就是一個基本蝶形運算流圖。當N=8，經(jīng)兩次分解，便分解為四個兩點DFT。N=8的完整DIF-FFT運算流圖如圖4.2.13所示。圖4.2.12DIF-FFT第二次分解運算流圖（N=8）圖4.2.13DIF-FFT運算流圖（N=8）這種算法是對X(k)進行奇偶抽取分解的結(jié)果，所以稱之為頻域抽取法FFT。觀察圖4.2.13可知：DIF-FFT算法與DIT-TTF算法類似，可以原位計算；DIF-FFT算法共有M級運算，每級共有N/2個蝶形運算，所以兩種算法的運算次數(shù)亦相同；不同的是DIF-FFT算法輸入序列為自然順序，而輸出為倒序排列；M級運算完后，要對輸出數(shù)據(jù)進行倒序才能得到自然順序的X(k)；蝶形運算略有不同，DIT-FFT蝶形先乘后加(減)，而DIF-FFT蝶形先加(減)后相乘。4.2.6IDFT的高效算法上述FFT算法流圖也可以用于計算IDFT。比較DFT和IDFT的運算公式：只要將DFT運算式中的系數(shù)

改變?yōu)?，最后乘?/N，就是IDFT運算公式；現(xiàn)在流圖的輸入是X(k)，輸出就是x(n)；原來的DIT-FFT→DIF-IFFT；DIF-FFT→DIT-IFFT。 如果希望直接調(diào)用FFT子程序計算IFFT，則可用下面的方法：由于1、將X(k)取復共軛；2、直接調(diào)用FFT子程序，或者送入FFT專用硬件設(shè)備進行DFT運算；3、最后取復共軛并乘以1/N得到序列x(n)。這種方法雖然用了兩次取共軛運算，但可以與FFT共用同一子程序，因而用起來很方便。4.3進一步減少運算量的措施4.3.1多類蝶形單元運算DIT-FFT運算：N=2M點FFT共需要MN/2次復數(shù)乘法由(4.2.12)式，當L=1時，只有一種旋轉(zhuǎn)因子，第一級不需要乘法運算；當L=2時，共有兩個旋轉(zhuǎn)因子：和，第二級也不需要乘法運算；依此類推在DFT中，又稱其值為±1和±j的旋轉(zhuǎn)因子為無關(guān)緊要的旋轉(zhuǎn)因子，如，，等。綜上所述：(1)先除去第一、二兩級后，所需復數(shù)乘法次數(shù)應是

(4.3.1)(2)當L=3時，有兩個無關(guān)緊要的旋轉(zhuǎn)因子和，同一旋轉(zhuǎn)因子對應的碟形運算個數(shù)：2M－L=N/2L第三級不需要進行復數(shù)乘法運算的蝶形：2N/23=N/4(3)依此類推，當L≥3時，第L級的2個無關(guān)緊要的旋轉(zhuǎn)因子

減少復數(shù)乘法的次數(shù)為：2N/2L=N/2L－1因此，DIT-FFT的復乘次數(shù)降至

(4.3.3)

下面再討論FFT中特殊的復數(shù)運算，以便進一步減少復數(shù)乘法次數(shù)。(4.3.2)這樣，從L=3至L=M共減少復數(shù)乘法次數(shù)為：一般的，實現(xiàn)一次復數(shù)乘法運算：四次實數(shù)乘；

兩次實數(shù)加。但對這一特殊復數(shù)，任一復數(shù)(x+jy)與其相乘時，只需要兩次實數(shù)加和兩次實數(shù)乘就可實現(xiàn)：——這樣，對應的每個蝶形節(jié)省兩次實數(shù)乘。在DIT-FFT運算流圖中，從L=3至L=M級，每級都包含旋轉(zhuǎn)因子，第L級中，對應N/2L個蝶形運算。因此從第三級至最后一級，旋轉(zhuǎn)因子節(jié)省的實數(shù)乘次數(shù)與(4.3.2)式相同。所以從實數(shù)乘運算考慮，計算N=2M點DIT-FFT所需實數(shù)乘法次數(shù)為(4.3.4)在基2FFT程序中：1、若包含了所有旋轉(zhuǎn)因子，則稱該算法為一類碟形單元運算；2、若去掉的旋轉(zhuǎn)因子，則稱之為二類碟形單元運算；3、若再去掉的旋轉(zhuǎn)因子，則稱為三類碟形單元運算；4、若再判斷處理，則稱之為四類碟形運算。后三種運算稱為多類碟形單元運算。顯然，碟形單元類型越多，編程就越復雜，但當N較大時，乘法運算的減少量是相當可觀的。例如，N=4096時，三類碟形單元運算的乘法次數(shù)為一類碟形單元運算的75%。4.3.2旋轉(zhuǎn)因子的生成在FFT運算中，旋轉(zhuǎn)因子，求正弦和余弦函數(shù)值的計算量是很大的。所以編程時，產(chǎn)生旋轉(zhuǎn)因子的方法直接影響運算速