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

第四章離散傅里葉變換快速算法

FastFourierTransform

(FFT)1

快速傅里葉變換（FFT）是計算DFT的一種快速有效方法，

并不是新的傅立葉變換式。

從前面的討論中看到，有限長序列在數(shù)字技術(shù)中占有很重要的地位。有限長序列的一個重要特點是其頻域也可以離散化，即離散傅里葉變換（DFT）。2DFT是信號分析與處理中的一種重要變換。但直接計算DFT的計算量與變換區(qū)間長度N的平方成正比，當(dāng)N較大時，計算量太大，直接用DFT算法進行譜分析和信號的實時處理是不切實際的。1965年庫里和圖基發(fā)現(xiàn)了DFT的一種快速算法，使DFT的運算效率提高1－2個數(shù)量級，為數(shù)字信號處理技術(shù)應(yīng)用于各種信號的實時處理創(chuàng)造了條件，推動了數(shù)字處理技術(shù)的發(fā)展。1984年，提出了分裂基快速算法，使運算效率進一步提高；3一、直接計算DFT的問題及改進的方法DFT的定義兩者形式類似，差別只在于的指數(shù)符號不同，及常數(shù)因子運算量是相同的41、運算量(a+jb)(c+jd)=(ac-bd)+j(ad+bc)∑

-=10)(NnnkNWnx562、減少運算量的途徑之一具有如下特性：使DFT運算中的有些項可以合并。對稱性周期性可約性knNNkNnNWW=--)()(=kNNkNWW-=+)2/(NNW-=2/,17減少運算量的途徑之二可將長序列的DFT分解為短序列的DFT。一個長序列DFT乘加運算比多個短序列DFT乘加運算多長序列DFT可由多個短序列DFT組合而成8二、按時間抽取(DIT)的基-2FFT算法

設(shè)M為自然數(shù)將長度為N的序列x(n)按n的奇偶分成兩組)12/(,1,0),12()()12/(,1,0),2()(21-=+=-==NrrxrxNrrxrx……1、算法原理9則x(n)的DFT為)()(21kXWkXkN+=∑∑

--+=12/2/212/2/1)()(NkrNkNNkrNWrxWWrxr=0r=01、算法原理10其中是x(2r)與x(2r+1)的N/2點DFT。1、算法原理11由于得12,,1,0-=Nk…1、算法原理12信號流圖

將X1

(k)和X2

(k)合成X(k)運算可歸結(jié)為：蝶形運算的簡化Wa+bWa-bW-W-1a+bWa-bWWabab(a)(b)X1(k)X2(k)X1(k)+WNkX2(k)X1(k)WNkX2(k)WNk圖：蝶形運算符號134點基2時間抽取FFT算法流圖x[0]x[2]x[1]x[3]X1[0]X1[1]X2[0]X2[1]2點DFT2點DFTX

[0]X

[1]X

[2]X

[3]14最后剩下的是2點DFT，它可以用一個蝶形結(jié)表示：x(0)x(1)X(0)X(1)158點基2FFT算法N/2點DFTN/2點DFTx(0)x(2)x(4)x(6)x(1)x(3)x(5)x(7)X1(0)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)WN0WN1WN2WN3X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)｛N=8點的DIT-2FFT(時域抽取圖)一次分解圖16

(3)第二次分解：

將x1(r)按r取奇、偶可分解成2個長度為N/4的子序列

x3(l)=x1(2l)、

x4(l)=x1(2l+1)，根據(jù)上面推導(dǎo)可得：X1(k)=X3(k)+WN/2kX4(k)，k=0,1,…,N/2-1將x2(r)按r取奇、偶可分解成2個長N/4的子序列

x5(l)=x2(2l)，l=0,1,…，N/4x6(l)=x2(2l+1)；同理得l=0,1,…，N/4-1；X1(k+N/2)=X3(k)WN/2kX4(k)，k=0,1,…,N/4-1；X1(k)=X3(k)+WN/2kX4(k)，k=0,1,…,N/4-1；X2(k)=X5(k)+WN/2kX6(k)，k=0,1,…,N/4-1；X2(k+N/2)=X5(k)WN/2kX6(k),k=0,1,…,N/4-1；17N/4點DFTx(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)X1(0)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)WN0WN1WN2WN3X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)N/4點DFTN/4點DFTN/4點DFTX3(0)X3(1)X4(0)X4(1)X5(0)X5(1)X6(0)X6(1)WN/20WN/21WN/21WN/20N=8點DFT的二次時域抽取分解圖X1(k+N/2)=X3(k)WN/2kX4(k)X1(k)=X3(k)+WN/2kX4(k)X2(k)=X5(k)+WN/2kX6(k)X2(k+N/2)=X5(k)WN/2kX6(k)k=0,1,…,N/4-1；18再次分解，對N=8點，可分解三次。X(5)N=8點DIT-FFT運算流圖x(0)x(4)x(2)x(6)x(1)x(5)x(3)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)WN0WN1WN2WN3X(0)X(1)X(2)X(3)X(4)X(6)X3(1)X4(0)X4(1)X5(0)X5(1)X6(0)X6(1)WN/20WN/21WN/20WN/40WN/40WN/40x(7)WN/21WN/40L=1級L=2L=3X(7)X3(0)X1(0)192點DFT的蝶形結(jié)表示：x(0)x(1)X(0)X(1)208點基2DIT-FFT算法N=8點DIT-FFT運算流圖x(0)x(4)x(2)x(6)x(1)x(5)x(3)WN0WN1WN2WN3X(0)X(1)X(2)X(3)X(4)X(5)X(6)WN0WN2WN0WN0WN0WN0x(7)WN2WN0L=1級L=2L=3X(7)21基2

DIT-FFT算法

由于這種方法每一步分解都是按輸入時間序列是屬于偶數(shù)還是奇數(shù)來抽取的，所以稱為“按時間抽取法”或“時間抽取法”。22

DIT―FFT算法與直接計算DFT運算量的比較1、直接DFT運算N點運算：復(fù)數(shù)乘次數(shù)：N×N

復(fù)數(shù)加次數(shù)：N×(N-1)2、用DIT-FFT作N點運算：復(fù)數(shù)乘次數(shù)：M×N/2=N/2×log2N；復(fù)加次數(shù)：2×N/2×M=N×log2N。整個運算流圖中有M級蝶形，每一級運算流圖中有N/2個蝶形，每個蝶形需一次復(fù)乘和兩次復(fù)數(shù)加運算。2324時間抽取法FFT的運算特點（1）蝶形運算（2）原位計算（3）序數(shù)重排（4）蝶形類型隨迭代次數(shù)成倍增加25（1）蝶形運算對于N=2M，總是可以通過M次分解最后成為2點的DFT運算。這樣構(gòu)成從x(n)到X(k)的M級運算過程。從上面的流圖可看到，每一級運算都由N/2個蝶形運算構(gòu)成。因此每一級運算都需要N/2次復(fù)乘和N次復(fù)加，這樣經(jīng)過時間抽取后M級運算總共需要的運算：復(fù)乘復(fù)加

而直接運算時則與N2成正比。26（2）原位計算當(dāng)數(shù)據(jù)輸入到存儲器中以后，每一級運算的結(jié)果仍然儲存在同一組存儲器中，直到最后輸出，中間無需其它存儲器，這叫原位計算。每一級運算均可在原位進行，這種原位運算結(jié)構(gòu)可節(jié)省存儲單元，降低設(shè)備成本，還可節(jié)省尋址的時間。27（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ī)律的。當(dāng)用二進制表示這個順序時，它正好是“碼位倒置”的順序。例如，原來的自然順序應(yīng)是x(1)的地方，現(xiàn)在放著x(4)，用二進制碼表示這一規(guī)律時，則是在

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

x（011）處放著x（110）。28自然順序二進碼表示碼位倒置碼位倒序0000000010011004201001023011110641000011510110156110010371111117

在實際運算中，一般直接將輸入數(shù)據(jù)x（n）按碼位倒置的順序排好輸入很不方便，總是先按自然順序輸入存儲單元，然后再通過變址運算將自然順序的存儲轉(zhuǎn)換成碼位倒置順序的存儲，然后進行FFT的原位計算。目前有許多通用DSP芯片支持這種碼位倒置的尋址功能。29造成倒位序的原因是輸入序列x(n)，按標號(序號)n的奇偶不斷地分組造成的。倒位序的形成30（4）蝶形類型隨迭代次數(shù)成倍增加觀察8點FFT的三次迭代運算:第一級迭代，有一種類型的蝶形運算系數(shù)W80，兩個數(shù)據(jù)點間隔為1第二級迭代，有二種類型的蝶形運算系數(shù)W80

、W82

，參加運算的兩個數(shù)據(jù)點間隔為2。第三級迭代，有四類蝶形運算系數(shù)W80

、W81

、W82

、W83

，參加運算的兩個數(shù)據(jù)點間隔為4。結(jié)論：每迭代一次，蝶形類型增加一倍，數(shù)據(jù)點間隔也增大一倍。每一級的取數(shù)間隔和蝶形類型種類均為2i-1，i=1,2,…M。31N點DIT—FFT運算流圖中，每級都有N/2個蝶形。每個蝶形都要乘以因子，稱其為旋轉(zhuǎn)因子，p稱為旋轉(zhuǎn)因子的指數(shù)。但各級的旋轉(zhuǎn)因子和循環(huán)方式都有所不同，為了編寫計算程序，應(yīng)先找出旋轉(zhuǎn)因子與運算級數(shù)的關(guān)系。用L表示從左到右的運算級數(shù)（L=1，2，、、，M）第L級共有個不同的旋轉(zhuǎn)因子。N=8時各級旋轉(zhuǎn)因子表示如下：32對的一般情況，第L級的旋轉(zhuǎn)因子為從而可以確定第L級運算的旋轉(zhuǎn)因子。33FFT算法流圖3435倒位序的實現(xiàn)（a）只要將順序二進制數(shù)（n2n1n0）的二進制位倒置，得（n0n1n2）。根據(jù)這種規(guī)律，容易用硬件電路和匯編語言產(chǎn)生倒位序數(shù)。（b）用高級語言程序?qū)崿F(xiàn)時，必須找出產(chǎn)生倒序數(shù)的十進制運算規(guī)律。36000000001

001

4

100201020103

011

6

1104

100

1

001510151016

110

3

01171117111左邊為按自然順序排列的二進制數(shù)，下面的一個數(shù)是上面一個數(shù)的最低位上加上1，且向高位進位。右邊為倒位序數(shù)，下面的一個數(shù)是上面一個數(shù)的最高位上加上1，且由高位向低位進位。稱為反向進位加法順序與倒序二進制對照表可由當(dāng)前任一倒序值求得下一個倒序值37

N=2M，用M位二進制數(shù)表示，則從左至右的十進制權(quán)值為：

N/2、N/4、N/8,…、2,1

對倒序數(shù)J，其下一個序數(shù)是在該序數(shù)J的二進制首位碼加1，相當(dāng)于十進制運算J+N/2。計算機上倒序數(shù)的實現(xiàn)過程為：J>N/2？J>N/4輸入當(dāng)前倒序數(shù)十進制數(shù)值JNYNYJ>N/2MNY結(jié)束J=J-N/2倒序數(shù)的十進制運算規(guī)律38N=8倒位序的變址處理分析上圖N=8點倒序規(guī)律，順序數(shù)I與倒序數(shù)J

存在關(guān)系：

I=J時，不交換；

I<J時，交換存儲器內(nèi)容。

I>J時，不交換，直接更新序數(shù)值；39LH=N/2J=LHN1=N-2I=1,N1I≥JT=A(I)A(I)=A(J)A(J)=TK=LHJ<KJ=J+KJ=J-KK=K/2NN

YY倒序程序框圖計算倒序值的程序流圖4041三、按頻率抽?。―IF）的基-2FFT算法1、算法原理仍設(shè)序列點數(shù)為N=2M，M為正整數(shù)。將X(k)按k的奇偶分組之前，先將輸入序列按前一半、后一半分開。k=0,1,…,N-142

k=0,1,…,N-1按k的奇偶可將X(k)分為兩部分:12,,1,0-=Nr…4312,,1,0-=Nr…令則x(n)x1(n)=x(n)+x(n+N/2)x2(n)=[x(n)x(n+N/2)]WNnWNnx(n+N/2)DIF-FFT蝶形運算流圖符號44X(3)N/2點DFTN/2點DFTx(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)x1(0)x1(1)x1(2)x1(3)x2(0)x2(1)x2(2)x2(3)WN0WN1WN2WN3X(0)X(2)X(4)X(6)X(1)X(5)X(7)N=8時第1次分解的運算流圖45N=2M，N/2仍是偶數(shù)，繼續(xù)將N/2點進行分解。將輸入序列x1(n)、x2(n)分別按前、后對半分解成4個長N/4的子序列,其n=0,1,…，N/4-1｛x3(n)=x1(n)+x1(n+N/4)

x4(n)=[x1(n)-x1(n+N/4)]WN/2n

x5(n)=x2(n)+x2(n+N/4)

x6(n)=[x2(n)-x2(n+N/4)]WN/2n｛N/4點DFTWN0WN1WN2WN3x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)X(0)X(4)X(2)X(6)X(1)X(5)X(3)X(7)WN0WN2WN0WN2N/4點DFTN/4點DFTN/4點DFTx1(0)x1(1)x1(2)x1(3)x2(0)x2(1)x2(2)x2(3)x3(0)x3(1)x4(0)x4(1)x5(0)x5(1)x6(0)x6(1)DIF―FFT二次分解運算流圖(N=8)

46經(jīng)過三次分解后，DIF―FFT運算流圖(N=8)為：x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)x1(0)x1(1)x1(2)x1(3)x2(0)x2(1)x2(2)x2(3)WN0WN1WN2WN3X(0)X(4)X(2)X(6)X(1)X(5)X(3)X(7)x3(0)x3(1)x4(0)x4(1)x5(0)x5(1)x6(0)x6(1)WN0WN2WN0WN2WN0WN0WN0WN0472、DIF-FFT運算規(guī)律

DIF-FFT算法也可采用原位計算；N=2M時，共有M級運算，每級共有N/2個蝶形，DIT與DIF算法的運算次數(shù)相同。DIT與DIF不同的是：

DIF-FFT算法輸入序列為自然序列，輸出為倒序序列。因此，在M級運算完成后，需對輸出數(shù)據(jù)進行倒序才能得到自然順序的X(k)。

蝶形運算符號不同：DIT-FFT蝶形是先相乘，后加/減；而DIF-FFT蝶形是先加/減，后相乘。DIT-FFT蝶形運算符號X1(k)X2(k)X1(k)+WNkX2(k)X1(k)WNkX2(k)WNkx2(n)=[x(n)x(n+N/2)]WNnx(n)x1(n)=x(n)+x(n+N/2)WNnx(n+N/2)DIF-FFT蝶形運算流圖符號48FFT的變形運算流圖(一)49FFT的變形運算流圖(二)50FFT減少運算量的方法旋轉(zhuǎn)因子的幾個特殊點取值存儲旋轉(zhuǎn)因子表提高運算速度51五、FFT應(yīng)用中的幾個問題1)IDFT的運算方法

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

DFT：X(k)=DFT[x(n)]=52IDFT與DFT的差別：1）把DFT中的每一個系數(shù)改為，2）再乘以常數(shù)1/N，則以上所討論的時間抽取或頻率抽取的FFT運算均可直接用于IDFT運算，當(dāng)然，蝶形中的系數(shù)應(yīng)改為。53DIT―IFFT運算流圖實際中，為了避免發(fā)生溢出，將1/N分配到每一級蝶形運算中，因為，所以每一級的每個蝶形輸出到乘以1/2.

54

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

IFFT計算分三步：①將X（k）取共軛（虛部乘以-1）

②對直接作FFT③對FFT的結(jié)果取共軛并乘以1/N，得x（n）。雖然2次用了取共軛運算，但可以和FFT共用一個子程序，實現(xiàn)方便。5556求某實信號x(n)的復(fù)譜,可認為是將實信號加上數(shù)值為零的虛部變成復(fù)信號(x(n)+j0),再用FFT求其離散傅里葉變換。把實序列變成復(fù)序列,存儲器要增加一倍,且計算機運行時,即使虛部為零,也要進行涉及虛部的運算。合理的解決方法是利用復(fù)數(shù)據(jù)FFT對實數(shù)據(jù)進行有效計算。2）實數(shù)序列的FFT

57

（1）用

一個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運算同時獲得。首先將x

(n)、y(n)分別當(dāng)作一復(fù)序列的實部及虛部,令g(n)=x

(n)+jy(n)對g(n)進行N點FFT運算，得到則：58

(2)用一個N點的FFT運算一個2N點實序列的DFT設(shè)x(n)是2N點的實序列,現(xiàn)人為地將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)=