數(shù)字信號處理4 fft

第四章快速傅里葉變換（FFT)4.2、直接計算DFT的問題和改進(jìn)途徑4.3、按時間抽取的FFT算法(DIT—FFT)4.4、按頻率抽取的FFT算法(DIF—FFT)4.5、N為復(fù)合數(shù)的FFT算法—統(tǒng)一FFT算法4.6、分裂基FFT算法4.7、FFT應(yīng)用4.8、線性調(diào)頻z變換（CZT)算法頻譜分析和DFT運算很重要，但在很長一段時間里，由于DFT運算復(fù)雜，并沒有得到真正的運用。頻譜分析以前大多采用模擬信號濾波的方法解決，直到1965年首次提出DFT運算的一種快速算法以后，情況才發(fā)生了根本變化。4.1引言2有限長序列在數(shù)字技術(shù)中占有很重要的地位。有限長序列的一個重要特點是其頻域也可以離散化表示，即離散傅里葉變換（DFT）。

認(rèn)識到DFT運算的一些內(nèi)在規(guī)律，發(fā)展和完善了高速有效的運算方法——快速傅里葉變換（FFT）算法。FFT的出現(xiàn)，使DFT的運算大大簡化，運算時間縮短1～2個數(shù)量級，使DFT的運算在實際中得到廣泛應(yīng)用。4.1引言31

DFT運算的特點

分析有限長序列x(n)進(jìn)行一次DFT運算所需的運算量。

一般，和

都是復(fù)數(shù)，每計算一個值，要進(jìn)行N次復(fù)數(shù)相乘，N-1次復(fù)數(shù)相加，一共有N個點，完成全部DFT運算，需要N2次復(fù)數(shù)相乘和N(N-1)次復(fù)數(shù)相加。44.2

直接計算DFT的問題和改進(jìn)途徑乘法比加法復(fù)雜，需要的運算時間多，尤其是復(fù)數(shù)相乘；每個復(fù)數(shù)相乘包括4個實數(shù)相乘和2個實數(shù)相加：1

DFT運算的特點51

DFT運算的特點整個DFT運算（N點）需要4N2實數(shù)相乘和2N（2N-1）次實數(shù)相加。

每個復(fù)數(shù)相加包括2個實數(shù)相加；每計算一個要進(jìn)行4N次實數(shù)相乘和2N+2（N-1）=2（2N-1）次實數(shù)相加；61

DFT運算的特點

例，計算N=10點的DFT，需要100次復(fù)數(shù)相乘，而N=1024點時，需要1048576（一百多萬）次復(fù)數(shù)乘法。如果要求實時處理，則要求有很高的計算速度才能完成上述計算量。

在DFT計算中，不論是乘法和加法，運算量均與N2成正比。N較大時，運算量十分可觀。7

反變換IDFT與DFT的運算結(jié)構(gòu)相同，只是多乘一個常數(shù)1/N，所以二者的計算量相同。1

DFT運算的特點82FFT算法的基本思想考察DFT與IDFT的運算發(fā)現(xiàn)，利用以下兩個特性可減少運算量：1）系數(shù)是一個周期函數(shù)，它的周期性和對稱性可用來改進(jìn)運算，提高計算效率。92FFT算法的基本思想利用這些周期性和對稱性，使DFT運算中有些項可合并；2）把長度為N點的大點數(shù)的DFT運算依次分解為若干個小點數(shù)的DFT。因為DFT的計算量正比于N2，N小，計算量也就小。10 FFT算法是基于這樣的基本思想發(fā)展起來的。有多種形式，基本上可分為兩類：

時間抽取法（DIT）

頻率抽取法（DIF）

2FFT算法的基本思想114.3按時間抽取的FFT算法:DIT--FFT先從一個特殊情況開始，假定N是2的整數(shù)次方：

N=2M，M：正整數(shù)--基-2FFT首先將序列x(n)分解為兩組，一組為偶數(shù)項，一組為奇數(shù)項，

12將DFT運算也相應(yīng)分為兩組：131算法原理因為故其中1算法原理

注意：H(k)，G(k)只有N/2個點，即k=0，1，…，N/2-1，還必須應(yīng)用系數(shù)

的周期性和對稱性求X(k)的N/2～N-1點：1算法原理15得：一個N點的DFT被分解為兩個N/2點的DFT，這兩個N/2點的DFT再合成為一個N點DFT。1算法原理16得

依此類推，G(k)和H(k)可以繼續(xù)拆分下去，這種按時間抽取算法是在輸入序列分成越來越小的子序列上執(zhí)行DFT運算，最后再合成為Ｎ點的DFT。1算法原理17蝶形信號流圖：將G(k)和H(k)合成X(k)運算可歸結(jié)為：蝶形運算流圖符號181算法原理

流圖需一次乘法、兩次加減法。

運算結(jié)構(gòu)像一蝴蝶，稱作蝶形運算結(jié)構(gòu)，簡稱蝶形結(jié)；采用這種表示法，可以將FFT過程用流圖表示。19蝶形信號流圖：1算法原理N=23=8

20蝶形信號流圖：1算法原理N=23=8

21蝶形信號流圖：1算法原理N=23=8

22蝶形信號流圖：1算法原理N=23=8

23蝶形信號流圖：1算法原理N=23=8

24蝶形信號流圖：1算法原理N=23=8

25蝶形信號流圖：1算法原理N=23=8

26蝶形信號流圖：1算法原理N=23=8

27蝶形信號流圖：1算法原理

按照這個辦法，繼續(xù)把N/2用２除；

N=2M，N/2可以被２整除，可以對兩個N/2點的DFT再分別作進(jìn)一步的分解，即對{G(k)}和{H(k)}的計算；{G(k)}和{H(k)}又可以分別通過計算兩個長度為N/4=2點的DFT，進(jìn)一步節(jié)省計算量。

一個８點的DFT就可以分解為四個２點的DFT。28蝶形信號流圖：1算法原理29x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)第三級抽取過程：1算法原理301,3,5,70,2,4,6x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)第二級第三級抽取過程：1算法原理1,3,5,70,2,4,60,42,61,53,7x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)第一級第二級第三級31抽取過程：1算法原理

由四個2點DFT組成8點DFT32抽取過程：1算法原理

由四個2點DFT組成8點DFT34抽取過程：1算法原理

由四個2點DFT組成8點DFT35抽取過程：1算法原理最后剩下的是2點DFT，它可以用一個蝶形結(jié)表示：這樣，得到一個8點的完整的按時間抽取運算的流圖。36抽取過程：1算法原理由于這種方法每一步分解都是按輸入時間序列是屬于偶數(shù)還是奇數(shù)來抽取的，所以稱為“按時間抽取法”或“時間抽取法”。37抽取過程：1算法原理N點DIT―FFT運算流圖(N=8)38抽取過程：1算法原理2時間抽取法FFT的運算特點1）蝶形運算2）原位計算3）序數(shù)重排391）蝶形運算

對于N=2M，總是可以通過M次分解最后成為2點的DFT運算。這樣構(gòu)成從x(n)到X(k)的M級運算過程。從流圖可知，每一級運算都由N/2個蝶形運算構(gòu)成。每一級運算都需要N/2次復(fù)乘和N次復(fù)加。402時間抽取法FFT的運算特點1）蝶形運算經(jīng)過時間抽取后M級運算總共需要的運算：復(fù)乘復(fù)加直接運算時則與N2成正比。N=2048，N2=4194304，(N/2)log2N=11264，N2/[(N/2)log2N]=392.4。FFT顯然要比直接法快得多。412時間抽取法FFT的運算特點422時間抽取法FFT的運算特點432時間抽取法FFT的運算特點2）原位計算當(dāng)數(shù)據(jù)輸入到存儲器中以后，每一級運算的結(jié)果仍然儲存在同一組存儲器中，直到最后輸出，中間無需其它存儲器，這叫原位計算。

每一級運算均可在原位進(jìn)行，這種原位運算結(jié)構(gòu)可節(jié)省存儲單元，降低設(shè)備成本，還可節(jié)省尋址的時間。442時間抽取法FFT的運算特點452時間抽取法FFT的運算特點2）原位計算N點DIT―FFT運算流圖(N=8)382）原位計算2時間抽取法FFT的運算特點3）序數(shù)重排對按時間抽取FFT的原位運算結(jié)構(gòu)，當(dāng)運算完畢時，正好順序存放著X(0)，X(1)，X(2)，…，X(7)，因此可直接按順序輸出。但原位運算的輸入x(n)卻不是自然順序，而是按x(0)，x(4)，x(2)，x(6)，…，x(7)的順序存入存儲單元，這種順序看起來相當(dāng)雜亂，然而它也是有規(guī)律的。462時間抽取法FFT的運算特點3）序數(shù)重排當(dāng)用二進(jìn)制表示這個順序時，它正好是“碼位倒置”（位翻轉(zhuǎn)）的順序。

例如，原來的自然順序應(yīng)是x(1)的地方，現(xiàn)在放著x(4)，用二進(jìn)制碼表示這一規(guī)律時，則是在

x（001）處放著x（100），

x（011）處放著x（110）。472時間抽取法FFT的運算特點碼位倒置順序表（位翻轉(zhuǎn)）自然順序二進(jìn)碼表示碼位倒置碼位倒置順序00000000100110042010010230111106410000115101101561100103711111172時間抽取法FFT的運算特點3）序數(shù)重排在實際運算中，直接將輸入數(shù)據(jù)x(n)按碼位倒置的順序排好輸入很不方便，總是先按自然順序輸入存儲單元，然后再通過變址運算將自然順序的存儲轉(zhuǎn)換成碼位倒置順序的存儲，然后進(jìn)行FFT的原位計算。目前有許多通用DSP芯片支持這種碼位倒置的尋址功能。492時間抽取法FFT的運算特點倒序規(guī)律502時間抽取法FFT的運算特點4.4按頻率抽取[DIF]的FFT算法51DITFFT：52將

分成前后兩半：前半子序列：后半子序列：4.4按頻率抽取[DIF]的FFT算法53由定義：4.4按頻率抽取[DIF]的FFT算法按k的奇偶將X(k)分為兩部分：54則4.4按頻率抽取[DIF]的FFT算法令554.4按頻率抽取[DIF]的FFT算法4.4按頻率抽取[DIF]的FFT算法則56DIF―FFT一次分解運算流圖(N=8)4.4按頻率抽取[DIF]的FFT算法57DIF―FFT二次分解運算流圖(N=8)4.4按頻率抽取[DIF]的FFT算法58DIF―FFT運算流圖(N=8)4.4按頻率抽取[DIF]的FFT算法594.5N為復(fù)合數(shù)的FFT算法(統(tǒng)一的FFT算法)處理方法補零N為素數(shù)，直接DFT或CZTN為復(fù)合數(shù)，用統(tǒng)一的FFT算法，基-2是特例6061算法原理4.5N為復(fù)合數(shù)的FFT算法(統(tǒng)一的FFT算法)614.5N為復(fù)合數(shù)的FFT算法(統(tǒng)一的FFT算法)n1n001230012314567289101161算法原理4.5N為復(fù)合數(shù)的FFT算法(統(tǒng)一的FFT算法)61算法效率4.5N為復(fù)合數(shù)的FFT算法(統(tǒng)一的FFT算法)運算復(fù)數(shù)乘法復(fù)數(shù)加法M個L點DFTML2ML(L-1)乘因子NL個M點DFTLM2LM(M-1)N點復(fù)合數(shù)FFTN(M+L+1)N(M+L-2)直接計算N點DFTN2N(N-1)加速比例N=60=12*53.33.9原理：序列按模4同余數(shù)抽取為4個序列，一直分到僅剩4點。71基-4DIT-FFT算法4.5N為復(fù)合數(shù)的FFT算法(統(tǒng)一的FFT算法)展開：擴展：72基-4DIT-FFT算法4.5N為復(fù)合數(shù)的FFT算法(統(tǒng)一的FFT算法)蝶形圖734.5N為復(fù)合數(shù)的FFT算法(統(tǒng)一的FFT算法)基-4DIT-FFT算法規(guī)律小結(jié)輸入反序，輸出正序原位運算蝶形級數(shù)每級蝶形數(shù)復(fù)乘復(fù)加744.5N為復(fù)合數(shù)的FFT算法(統(tǒng)一的FFT算法)基-4DIT-FFT算法把x(n)按前后分作4段，求DFT。754.5N為復(fù)合數(shù)的FFT算法(統(tǒng)一的FFT算法)基-4DIF-FFT算法按模4同余數(shù)抽取X(k)764.5N為復(fù)合數(shù)的FFT算法(統(tǒng)一的FFT算法)基-4DIF-FFT算法蝶形圖

先不計算DFT，而繼續(xù)按先后分4段，直到剩4點為止。774.5N為復(fù)合數(shù)的FFT算法(統(tǒng)一的FFT算法)基-4DIF-FFT算法綜合基-2、基-4DIF而形成的一種高效、實用算法。784.6分裂基FFT算法4.6分裂基FFT算法794.6分裂基FFT算法814.6分裂基FFT算法一個N點的DFT計算可分解為一個N/2點DFT和兩個N/4點DFT。804.6分裂基FFT算法L-蝶形圖824.6分裂基FFT算法N/4個“L”形，每個L形包含2個乘法運算，其余為加法運算。蝶形圖中求和向上者按基-2，向下者按基-4(包含L形)進(jìn)行分解，直到4點為止。834.6分裂基FFT算法8點分裂基FFT流圖844.6分裂基FFT算法85規(guī)律結(jié)構(gòu)與基-2DIF相同，僅乘系數(shù)不同；運算量較基-2減少；每級包含L形個數(shù)：4.6分裂基FFT算法每個L形2次復(fù)乘，總共基-2復(fù)乘數(shù)基-4復(fù)乘數(shù)864.7FFT應(yīng)用1IDFT的運算方法

以上所討論的FFT算法可用于IDFT運算——簡稱為IFFT,比較IDFT的定義式：

87IDFT與DFT的差別：把DFT中的每一個系數(shù)

改為再乘以常數(shù)1/N

因此，以上所討論的時間抽取或頻率抽取的FFT運算均可直接用于IDFT運算，當(dāng)然，蝶形中的系數(shù)應(yīng)改為。881IDFT的運算方法

89第二種方法，完全不需要改動FFT程序，而是直接利用它作IFFT?？紤]到故

1IDFT的運算方法

90

IFFT計算分三步：①將X(k)取共軛，即虛部乘以-1

②對直接作FFT③對FFT的結(jié)果取共軛并乘以1/N，得x(n)1IDFT的運算方法

91FFT算法都是復(fù)數(shù)運算,包括序列x(n)也認(rèn)為是復(fù)數(shù);大多數(shù)場合,信號是實數(shù)序列,任何實數(shù)都可看成虛部為零的復(fù)數(shù);例如,求某實信號x(n)的復(fù)譜,可認(rèn)為是將實信號加上數(shù)值為零的虛部變成復(fù)信號(x(n)+j0),再用FFT求其離散傅里葉變換。2實數(shù)序列的FFT

這種作法很不經(jīng)濟(jì),把實序列變成復(fù)序列,存儲器要增加一倍,且計算機運行時,即使虛部為零,也要進(jìn)行涉及虛部的運算,浪費了運算量。

合理的解決方法是利用復(fù)數(shù)據(jù)FFT對實數(shù)據(jù)進(jìn)行有效計算,下面介紹兩種方法。2實數(shù)序列的FFT92用一個N點FFT同時計算兩個N點實序列的DFT

設(shè)x(n)、y(n)是彼此獨立的兩個N點實序列,且X

(k)=DFT[x

(n)]，Y

(k)=DFT[y(n)]

則X

(k)、Y(k)可通過一次FFT運算同時獲得。

2實數(shù)序列的FFT932實數(shù)序列的FFT

首先，將x

(n)、y(n)分別當(dāng)作一復(fù)序列的實部及虛部,令94對g(n)進(jìn)行N點FFT運算，得到則：2實數(shù)序列的FFT952實數(shù)序列的FFT用一個N點的FFT運算獲得一個2N點實序列的DFT

設(shè)x(n)是2N點的實序列,將x(n)分為偶數(shù)組x1(n)和奇數(shù)組x2(n)x1(n)=x(2n)n=0,1,…,N-1x2(n)=x(2n+1)n=0,1,…,N-1然后將x1(n)及x2(n)組成一個復(fù)序列：

y(n)=x1(n)+jx2(n)96通過N點FFT運算可得到：

Y(k)=X1(k)+jX2(k) 根據(jù)前面的討論，得到

972實數(shù)序列的FFT

為求2N點x(n)所對應(yīng)X(k)，需求出X(k)與X1(k)、X2(k)的關(guān)系：

2實數(shù)序列的FFT98所以2實數(shù)序列的FFT99計算步驟：1）將序列x(n)按奇偶序拆分成x1(n)及x2(n)，并組成復(fù)序列y(n)，經(jīng)FFT運算求得Y(k)，2）利用共軛對稱性求出X1(k)、X2(k)，3）最后利用上式求出X(k),實現(xiàn)用一個N點的FFT計算一個2N點的實序列的DFT。2實數(shù)序列的FFT100

線性卷積是求離散系統(tǒng)響應(yīng)的主要方法之一,許多重要應(yīng)用都建立在這一理論基礎(chǔ)上,如卷積濾波等。以前曾討論了用圓周卷積計算線性卷積的方法歸納如下:3線性卷積的FFT算法1013線性卷積的FFT算法

將長為N2的序列x(n)延長到L,補L-N2個零，將長為N1的序列h(n)延長到L,補L-N1個零，如果L≥N1+N2-1,則圓周(循環(huán))卷積與線性卷積相等,此時,可用FFT計算線性卷積，方法如下:

1023線性卷積的FFT算法a.計算X(k)=FFT[x(n)]b.求H(k)=FFT[h(n)]c.求Y(k)=H(k)X(k)k=0～L-1d.求y(n)=IFFT[Y(k)]n=0～L-1可見，只要進(jìn)行二次FFT,一次IFFT就可完成線性卷積計算。103計算表明,L>32時,上述計算線性卷積的方法比直接計算線卷積有明顯的優(yōu)越性,因此,也稱上述循環(huán)卷積方法為快速卷積法。上述結(jié)論適用于x(n)、h(n)兩序列長度比較接近或相等的情況。1043線性卷積的FFT算法

如果x(n)、h(n)長度相差較多,例如,h(n)為某濾波器的單位脈沖響應(yīng),長度有限,用來處理一個很長的輸入信號x(n),或者處理一個連續(xù)不斷的信號，按上述方法,h(n)要補許多零再進(jìn)行計算,計算量有很大的浪費，或者根本不能實現(xiàn)。為了保持快速卷積法的優(yōu)越性,可將x(n)分為許多段,每段的長度與h(n)接近,處理方法有兩種：3線性卷積的FFT算法1051)

重疊相加法——由分段卷積的各段相加構(gòu)成總的卷積輸出

h(n)x(n)3線性卷積的FFT算法106計算步驟：a.

事先計算濾波器參數(shù)H(k)=DFT[h(n)]，N點b.

用N點FFT計算Xi(k)=DFT[xi(n)]c.

Yi(k)=Xi(k)H(k)d.

用N點IFFT求yi(n)=IDFT[Yi(k)]e.

將重疊部分相加1073線性卷積的FFT算法3線性卷積的FFT算法1082)重疊保留

和第一種方法稍有不同，即將上面分段序列中補零的部分不是補零，而是保留原來的輸入序列值，這時，如利用DFT實現(xiàn)h(n)和xi(n)的循環(huán)卷積，則每段卷積結(jié)果中有N1-1個點不等于線性卷積值需舍去。

重疊保留法與重疊相加法的計算量差不多，但省去了重疊相加法最后的相加運算。

3線性卷積的FFT算法109110111相關(guān)的概念很重要，互相關(guān)運算廣泛應(yīng)用于信號分析與統(tǒng)計分析，如通過相關(guān)函數(shù)峰值的檢測測量兩個信號的時延差等。兩個長為N的實離散時間序列x（n）與y（n）的互相關(guān)函數(shù)定義為

4用FFT計算相關(guān)函數(shù)則可以證明，rxy(m)的離散傅里葉變換為

Rxy(k)=DFT[rxy(m)] =X*(k)Y

(k)0≤k≤N-1

其中X(k)=DFT[x(n)],Y(k)=DFT[y(n)],1124用FFT計算相關(guān)函數(shù)113證明：互相關(guān)函數(shù)定義為

x(n)及y(n)的卷積公式

相比較，可以得到相關(guān)和卷積的時域關(guān)系：

4用FFT計算相關(guān)函數(shù)114因得證畢。4用FFT計算相關(guān)函數(shù)115

當(dāng)x(n)=y(n)時，得到x(n)的自相關(guān)函數(shù)為：

維納--辛欽定理：

自相關(guān)函數(shù)與信號功率譜互為傅立葉變換對。

4用FFT計算相關(guān)函數(shù)116推導(dǎo)表明，互相關(guān)和自相關(guān)函數(shù)的計算可利用FFT實現(xiàn)。由于離散傅里葉變換隱含著周期性，所以用FFT計算離散相關(guān)函數(shù)也是對周期序列而言的。直接做N點FFT相當(dāng)于對兩個N點序列x(n)、y(n)作周期延拓，作相關(guān)后再取主值（類似圓周卷積）。

4用FFT計算相關(guān)函數(shù)117實際一般要求的是兩個有限長序列的線性相關(guān)，為避免混淆，需采用與循環(huán)卷積求線性卷積相類似的方法，先將序列延長補0后再用上述方法。4用FFT計算相關(guān)函數(shù)118設(shè)x(n)和y(n)的長均為N，求線性相關(guān)；為了使兩個有限長序列的線性相關(guān)可用其圓周相關(guān)代替而不產(chǎn)生混淆，選擇周期L≥2N-1，且L=2m，以使用FFT，將x(n)，y(n)補零至長為L。用FFT計算X(k)，Y(k)（k=0，1…，L-1）。利用FFT求兩個有限長序列的線性相關(guān)：4用FFT計算相關(guān)函數(shù)119R（k）=X*（k）Y（k）對R（k）作IFFT，取后N-1項，得取前N項，得4用FFT計算相關(guān)函數(shù)123

二維信號有圖像信號、時空信號、時頻信號等。二維離散傅里葉變換可用于處理二維離散信號。二維離散傅里葉變換的定義為：

5用FFT計算二維離散的傅里葉變換124

一個二維DFT可用二個一維DFT計算。

若兩次DFT都用快速算法求得,且M和N都是2的冪,則可使用基-2FFT算法,所需要乘法次數(shù)為

5用FFT計算二維離散的傅里葉變換125二維離散傅里葉變換所需的乘法次數(shù)為（r+q）MN，當(dāng)M和N比較大時用FFT運算