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

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

4.1引言4.2基2FFT算法4.3進(jìn)一步減少運(yùn)算量的措施4.4其他快速算法簡介習(xí)題與上機(jī)題4.1引言

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

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

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

高效計(jì)算DFT的算法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(4.2.12)

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

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

4.編程思想及程序框圖

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

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

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

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

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

對應(yīng)的倒序數(shù)4：100最高位為1：010——最高位加1要向次高位進(jìn)位，其實(shí)質(zhì)是將最高位變?yōu)?，再在次高位加1。用這種算法，可以從當(dāng)前任一倒序值求得下一個(gè)倒序值但用有些高級語言程序?qū)崿F(xiàn)倒序時(shí)，直接倒置二進(jìn)制數(shù)位是不行的，因此必須找出產(chǎn)生倒序數(shù)的十進(jìn)制運(yùn)算規(guī)律。由表4.2.1可見：若用J表示當(dāng)前倒序數(shù)的十進(jìn)制數(shù)值，對于N=2M，M位二進(jìn)制數(shù)的權(quán)值從左向右依次為N/2，N/4，N/8，…，2，1。二進(jìn)制最高位加1：十進(jìn)制：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時(shí)，同樣要判斷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向右進(jìn)位的運(yùn)算。(圖4.2.9所示的虛線框?yàn)榈剐虻某绦蚩驁D)

形成倒序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所示：①第一個(gè)序列值x(0)和最后一個(gè)序列值x(N－1)不需要重排；②每計(jì)算出一個(gè)倒序值J，便與循環(huán)語句自動(dòng)生成的順序I

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

當(dāng)I≠J時(shí),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)前后對半分開，得到兩個(gè)子序列，其DFT如下：式中將X(k)分解成偶數(shù)組與奇數(shù)組：當(dāng)k取偶數(shù)(k=2m,m=0,1,…,N/2－1)時(shí)(4.2.14)當(dāng)k取奇數(shù)(k=2m+1,m=0,1,…,N/2－1)時(shí)，令將x1(n)和x2(n)分別代入(4.2.14)和(4.2.15)式，可得

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

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

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

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

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

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

(4.3.3)

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

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