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

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

FastFourierTransform

(FFT)1

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

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

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

-=10)(NnnkNWnx562、減少運(yùn)算量的途徑之一具有如下特性：使DFT運(yùn)算中的有些項(xiàng)可以合并。對(duì)稱性周期性可約性knNNkNnNWW=--)()(=kNNkNWW-=+)2/(NNW-=2/,17減少運(yùn)算量的途徑之二可將長序列的DFT分解為短序列的DFT。一個(gè)長序列DFT乘加運(yùn)算比多個(gè)短序列DFT乘加運(yùn)算多長序列DFT可由多個(gè)短序列DFT組合而成8二、按時(shí)間抽取(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點(diǎn)DFT。1、算法原理11由于得12,,1,0-=Nk…1、算法原理12信號(hào)流圖

將X1

(k)和X2

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

[0]X

[1]X

[2]X

[3]14最后剩下的是2點(diǎn)DFT，它可以用一個(gè)蝶形結(jié)表示：x(0)x(1)X(0)X(1)158點(diǎn)基2FFT算法N/2點(diǎn)DFTN/2點(diǎn)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點(diǎn)的DIT-2FFT(時(shí)域抽取圖)一次分解圖16

(3)第二次分解：

將x1(r)按r取奇、偶可分解成2個(gè)長度為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個(gè)長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點(diǎn)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點(diǎn)DFTN/4點(diǎn)DFTN/4點(diǎn)DFTX3(0)X3(1)X4(0)X4(1)X5(0)X5(1)X6(0)X6(1)WN/20WN/21WN/21WN/20N=8點(diǎn)DFT的二次時(shí)域抽取分解圖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再次分解，對(duì)N=8點(diǎn)，可分解三次。X(5)N=8點(diǎn)DIT-FFT運(yùn)算流圖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級(jí)L=2L=3X(7)X3(0)X1(0)192點(diǎn)DFT的蝶形結(jié)表示：x(0)x(1)X(0)X(1)208點(diǎn)基2DIT-FFT算法N=8點(diǎn)DIT-FFT運(yùn)算流圖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級(jí)L=2L=3X(7)21基2

DIT-FFT算法

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

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

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

而直接運(yùn)算時(shí)則與N2成正比。26（2）原位計(jì)算當(dāng)數(shù)據(jù)輸入到存儲(chǔ)器中以后，每一級(jí)運(yùn)算的結(jié)果仍然儲(chǔ)存在同一組存儲(chǔ)器中，直到最后輸出，中間無需其它存儲(chǔ)器，這叫原位計(jì)算。每一級(jí)運(yùn)算均可在原位進(jìn)行，這種原位運(yùn)算結(jié)構(gòu)可節(jié)省存儲(chǔ)單元，降低設(shè)備成本，還可節(jié)省尋址的時(shí)間。27（3）序數(shù)重排對(duì)按時(shí)間抽取FFT的原位運(yùn)算結(jié)構(gòu)，當(dāng)運(yùn)算完畢時(shí)，正好順序存放著x(0)，x(1)，x(2)，…x(7)，因此可直接按順序輸出，但這種原位運(yùn)算的輸入x（n）卻不能按這種自然順序存入存儲(chǔ)單元中，而是按x(0)，x(4)，x(2)，x(6)，…，x(7)的順序存入存儲(chǔ)單元，這種順序看起來相當(dāng)雜亂，然而它也是有規(guī)律的。當(dāng)用二進(jìn)制表示這個(gè)順序時(shí)，它正好是“碼位倒置”的順序。例如，原來的自然順序應(yīng)是x(1)的地方，現(xiàn)在放著x(4)，用二進(jìn)制碼表示這一規(guī)律時(shí)，則是在

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

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

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

、W82

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

、W81

、W82

、W83

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

001

4

100201020103

011

6

1104

100

1

001510151016

110

3

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

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

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

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

存在關(guān)系：

I=J時(shí)，不交換；

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

I>J時(shí)，不交換，直接更新序數(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倒序程序框圖計(jì)算倒序值的程序流圖4041三、按頻率抽?。―IF）的基-2FFT算法1、算法原理仍設(shè)序列點(diǎn)數(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蝶形運(yùn)算流圖符號(hào)44X(3)N/2點(diǎn)DFTN/2點(diǎn)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時(shí)第1次分解的運(yùn)算流圖45N=2M，N/2仍是偶數(shù)，繼續(xù)將N/2點(diǎn)進(jìn)行分解。將輸入序列x1(n)、x2(n)分別按前、后對(duì)半分解成4個(gè)長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點(diǎn)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點(diǎn)DFTN/4點(diǎn)DFTN/4點(diǎn)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二次分解運(yùn)算流圖(N=8)

46經(jīng)過三次分解后，DIF―FFT運(yùn)算流圖(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運(yùn)算規(guī)律

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

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

蝶形運(yùn)算符號(hào)不同：DIT-FFT蝶形是先相乘，后加/減；而DIF-FFT蝶形是先加/減，后相乘。DIT-FFT蝶形運(yùn)算符號(hào)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蝶形運(yùn)算流圖符號(hào)48FFT的變形運(yùn)算流圖(一)49FFT的變形運(yùn)算流圖(二)50FFT減少運(yùn)算量的方法旋轉(zhuǎn)因子的幾個(gè)特殊點(diǎn)取值存儲(chǔ)旋轉(zhuǎn)因子表提高運(yùn)算速度51五、FFT應(yīng)用中的幾個(gè)問題1)IDFT的運(yùn)算方法

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

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

54

b）第二種方法，完全不需要改動(dòng)FFT程序，而是直接利用它作IFFT。考慮到故

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

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

57

（1）用

一個(gè)N點(diǎn)FFT同時(shí)計(jì)算兩個(gè)N點(diǎn)實(shí)序列的DFT

設(shè)x

(n)、y

(n)是彼此獨(dú)立的兩個(gè)N點(diǎn)實(shí)序列,且X

(k)=DFT[x

(n)]，Y

(k)=DFT[y(n)]則X

(k)、Y(k)可通過一次FFT運(yùn)算同時(shí)獲得。首先將x

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

(n)+jy(n)對(duì)g(n)進(jìn)行N點(diǎn)FFT運(yùn)算，得到則：58

(2)用一個(gè)N點(diǎn)的FFT運(yùn)算一個(gè)2N點(diǎn)實(shí)序列的DFT設(shè)x(n)是2N點(diǎn)的實(shí)序列,現(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)組成一個(gè)復(fù)序列：y(n)=