《數(shù)字信號(hào)處理-基于數(shù)值計(jì)算》課件-第3章

第3章離散傅立葉變換及其應(yīng)用3.1引言3.2離散傅立葉變換3.3快速傅立葉變換FFT3.4離散傅立葉變換的實(shí)際應(yīng)用問(wèn)題3.5快速傅立葉變換FFT典型用法

3.1引言

離散傅立葉變換(DFT)是數(shù)字信號(hào)處理領(lǐng)域中功能最強(qiáng)大的常用方法之一，它與數(shù)字濾波器理論構(gòu)成了數(shù)字信號(hào)處理的最基本內(nèi)容。DFT可以看成是一般傅立葉變換的數(shù)字化形式。二百多年前，法國(guó)科學(xué)家傅立葉(Fourier)在他的一篇著名論文中證明了“任一周期函數(shù)都可以表示成正弦函數(shù)和的形式，其中正弦函數(shù)的頻率是周期函數(shù)頻率的整數(shù)倍”這一理論，對(duì)科學(xué)發(fā)展和工程實(shí)踐產(chǎn)生了巨大而深遠(yuǎn)的影響。后來(lái)人們?cè)谛盘?hào)分析領(lǐng)域提出了基于輸入與輸出都是離散形式的傅立葉變換算法(DFT)，它巧妙地契合了數(shù)字計(jì)算機(jī)的工作方式，伴隨著計(jì)算機(jī)運(yùn)算速度的極大提高，也得益于1965年庫(kù)利(Cooley)和圖基(Tukey)提出的高效率計(jì)算DFT的快速方法，使得DFT的各種實(shí)際應(yīng)用獲得蓬勃發(fā)展，尤其是在卷積計(jì)算、調(diào)制、解調(diào)、離散余弦變換等數(shù)字信號(hào)處理算法中。我們?cè)诘?章講到周期序列的傅立葉級(jí)數(shù)展開(kāi)，即DFS算法，它的輸入與輸出都呈現(xiàn)離散形式，不過(guò)由于它們是周期序列，無(wú)始無(wú)終，盡管全部的信息都只蘊(yùn)藏在有限點(diǎn)的一個(gè)周期內(nèi)，但還是不能直接使用計(jì)算機(jī)進(jìn)行計(jì)算。幸運(yùn)的是，現(xiàn)實(shí)中要處理的信號(hào)都可以認(rèn)為是有限長(zhǎng)度的，哪怕是持續(xù)很長(zhǎng)時(shí)間的信號(hào)，我們也可以合理地將其分成許多小段，再頭尾相接來(lái)逼近。因此，研究一個(gè)有限長(zhǎng)序列的計(jì)算顯得至關(guān)重要。另一方面，我們把有限長(zhǎng)序列(假設(shè)N點(diǎn))想像成是某個(gè)周期為N點(diǎn)的周期序列的一個(gè)周期，這種思路稱(chēng)為有限長(zhǎng)序列的周期延拓，從而可以借用周期序列的DFS進(jìn)行計(jì)算，最后在計(jì)算結(jié)果中截取其一個(gè)周期。這個(gè)過(guò)程就是以下要介紹的DFT算法。

3.2離散傅立葉變換

3.2.1DFT定義

為了引用周期序列的概念，我們先來(lái)討論周期序列和有限長(zhǎng)序列之間的相互表達(dá)。假定一個(gè)周期序列x(n)，它是由長(zhǎng)為N點(diǎn)的有限長(zhǎng)序列x(n)經(jīng)周期延拓而成，即~(3.2.1)(3.2.2)周期序列x(n)的第一個(gè)周期n=0~(N－1)定義為x(n)的主值區(qū)間,也就是說(shuō)，x(n)是x(n)的主值序列，即

x(n)=x(n)RN(n)

(3.2.3)式(3.2.3)中RN(n)為N點(diǎn)因果矩形序列。如圖3.2.1所示是一個(gè)N=6點(diǎn)的有限長(zhǎng)序列的周期延拓序列。~~~~圖3.2.1有限長(zhǎng)序列N=6的周期延拓同樣，對(duì)于DFS的周期性的系數(shù)X(k)，與主值頻譜序列X(k)也有如下關(guān)系：~(3.2.4)(3.2.5)為了便于對(duì)照，將第2章的DFS變換對(duì)重寫(xiě)如下：式中WN=e－j(2π/N)稱(chēng)為旋轉(zhuǎn)因子，幅度為1，相角為－2π/N，任何復(fù)數(shù)與其相乘將被附加一個(gè)負(fù)相角，在極坐標(biāo)圖中相當(dāng)于被順時(shí)針旋轉(zhuǎn)2π/N弧度。

WN是離散傅立葉變換中重要的單位復(fù)指數(shù)，它是序列點(diǎn)數(shù)N的函數(shù)。注意到DFS的求和運(yùn)算區(qū)間就是在周期序列的主值區(qū)間進(jìn)行的，因此，我們定義有限長(zhǎng)N點(diǎn)序列的DFT為式(3.2.6)和(3.2.7)分別稱(chēng)為離散傅立葉正變換DFT和離散傅立葉逆變換IDFT。再次強(qiáng)調(diào)的是，雖然以上兩個(gè)式子都是針對(duì)有限長(zhǎng)序列的，卻都應(yīng)該將其理解為周期序列中的一個(gè)周期。

x(n)與X(k)都是長(zhǎng)度為N

的序列(復(fù)序列)，都有N個(gè)獨(dú)立值，因而具有等量的信息。已知x(n)就能唯一地確定X(k)，反之亦然。改寫(xiě)式(3.2.6)如下：可以看出，從k=0，1，2，…，N－1，相當(dāng)于頻率從直流DC到(N－1)fs的N個(gè)頻率等分點(diǎn)。這些頻率點(diǎn)上的余弦序列和正弦序列稱(chēng)之為頻率單元或分析頻點(diǎn)。也就是說(shuō)，輸入時(shí)域序列x(n)與頻率單元做序列點(diǎn)積運(yùn)算而得到頻譜的實(shí)部和虛部，即該頻率點(diǎn)所分解到的復(fù)系數(shù)X(k)。

為了方便DFT的使用，特別是用MATLAB軟件編程時(shí)，經(jīng)常將其表示成矩陣形式。令

x={x(0),x(1)，x(2),…,x(N－1)}T構(gòu)成時(shí)域序列的列矩陣

X={X(0),X(1)，X(2),…,X(N－1)}T構(gòu)成頻域序列的列矩陣

那么式(3.2.6)可以寫(xiě)成

X=WNx

(3.2.8)式中WN是N×N的方陣，而且關(guān)于主對(duì)角線對(duì)稱(chēng)，即(3.2.9)同理，式(3.2.7)的IDFT式可寫(xiě)成(3.2.10)DFT隱含的周期性可以從定義中看出來(lái)，對(duì)任意整數(shù)m，總有

【例3.2.1】

某復(fù)合正弦信號(hào)x(t)由幅度為1個(gè)單位、頻率為2kHz和幅度為0.5個(gè)單位、頻率為6kHz且滯后135°相角的兩個(gè)正弦分量復(fù)合構(gòu)成：

x(t)=sin(2π×2000×t)+0.5sin(2π×6000×t－3π/4)

若對(duì)其用fs=16kHz采樣率進(jìn)行離散化(即16000個(gè)樣點(diǎn)數(shù)據(jù)/秒)，獲得一串x(nT)。x(nT)=x(t)|t=nT，T=1/16000。現(xiàn)在我們僅取其16個(gè)數(shù)據(jù)(相當(dāng)于觀察了1ms的復(fù)合正弦信號(hào))進(jìn)行分析，可以獲得這些數(shù)據(jù)所提供的頻譜信息。試計(jì)算這16個(gè)數(shù)據(jù)的DFT變換X(k)值，繪制X(k)的幅度序列和相位序列圖并加以仔細(xì)考察。解先用MATLAB計(jì)算出16個(gè)采樣數(shù)據(jù)(復(fù)合正弦序列)。

t=0:1/16000:15/16000；%t=0開(kāi)始，時(shí)間增量1/16ms，觀察16點(diǎn)，即1ms

xt=sin(2*pi*2000*t)+0.5*sin(2*pi*6000*t-3*pi/4)；%復(fù)合正弦序列

figure(1);stem(t,xt);xlabel(′n′);ylabel(′x(n)′);

16點(diǎn)采樣序列值x(n)={－0.35，0.71，1.35，0.21，0.35，－0.71，－1.35，－0.21，－0.35，0.71，1.35，0.21，0.35，－0.71，－1.35，－0.21}；繪出的序列桿圖如圖3.2.2所示。圖3.2.2復(fù)合正弦的采樣序列構(gòu)造16×16的變換矩陣WN，并計(jì)算出頻譜X(k)。

n=0:15;k=0:15；%兩個(gè)序號(hào)行向量

WN=exp((－j*2*pi/16)).^(n′*k)；%構(gòu)造變換矩陣，注意是群運(yùn)算

X=xt*WN;Xa=abs(X);%進(jìn)行DFT運(yùn)算，獲得頻譜序列X(k)，求幅度

Xb=(angle(X))*180/pi；%求相角，單位從弧度換成角度

figure(2);

subplot(2,1,1);stem(k,Xa);xlabel(′k′);ylabel(′X(k)′)；%幅度譜

subplot(2,1,2);stem(k,Xb);xlabel(′k′);ylabel(′φ(k)deg′);%相位譜

可以得到頻譜序列值：

X(k)=｛0，0，－8i，0，0， 0，－2.8+2.8i，0，0，0，－2.8－2.8i，0，0，0，+8i，0｝；

程序運(yùn)行后分別繪出幅頻和相頻特性如圖3.2.3所示。(注意MATLAB計(jì)算結(jié)果中近似于0的數(shù)值的表示方式，同時(shí)，對(duì)于非常小的幅度信號(hào)分量所對(duì)應(yīng)的相位值可以不必考慮。)圖3.2.3復(fù)合正弦序列的頻譜圖3.2.2DFT性質(zhì)

下面討論DFT的性質(zhì)。假設(shè)有限長(zhǎng)序列x1(n)和x2(n)的長(zhǎng)度分別為N1和N2，取N=max[N1，N2]，即短的序列在其后補(bǔ)相應(yīng)個(gè)0。設(shè)x1(n)和x2(n)的N點(diǎn)DFT分別為：

X1(k)=DFT[x1(n)]

X2(k)=DFT[x2(n)]

(1)線性性質(zhì)。假設(shè)y(n)=ax1(n)+bx2(n)，則有

Y(k)=DFT[y(n)]=aX1(k)+bX2(k)，0≤k≤N－1

(3.2.11)

式中a、b為任意常數(shù)。

(2)循環(huán)移位。N點(diǎn)序列x(n)的循環(huán)移位m點(diǎn)后的序列y(n)定義為

y(n)=x((n+m))NRN(n)

(3.2.12)式(3.2.12)的x((n+m))N符號(hào)表示這樣的操作：先以x(n)當(dāng)主值周期做周期延拓成為x(n)，然后進(jìn)行移位，最后由RN(n)截取0~N－1區(qū)間成為N點(diǎn)的y(n)。這說(shuō)明，N點(diǎn)序列循環(huán)移位后仍然是N點(diǎn)序列。如圖3.2.4所示，圖(a)為N=6點(diǎn)的序列；圖(b)為進(jìn)行周期延拓的序列；圖(c)為周期序列右移3點(diǎn)的序列；圖(d)為截取n=0~5主值周期。

圖3.2.5說(shuō)明了同一個(gè)序列進(jìn)行線性移位和循環(huán)移位操作的結(jié)果差異，左邊是線性位移，0從左邊右移入主值區(qū)而序列值被移出；右邊是循環(huán)移位，6個(gè)值都在主值區(qū)內(nèi)。請(qǐng)仔細(xì)體會(huì)。~圖3.2.4有限長(zhǎng)N點(diǎn)序列循環(huán)移位圖3.2.5有限長(zhǎng)6點(diǎn)序列線性位移和循環(huán)移位的區(qū)別考慮到，我們有

上式說(shuō)明，經(jīng)過(guò)循環(huán)移位后的序列的DFT與移位前的DFT差別僅是復(fù)乘以旋轉(zhuǎn)因子而發(fā)生相位上的改變。這個(gè)改變量

與移位置m以及離散頻率k成正比，對(duì)幅頻特性沒(méi)有影響。同樣，根據(jù)DFT變換的對(duì)偶性，有(3.2.13)(3.2.14)

(3)共軛對(duì)稱(chēng)性。設(shè)x*(n)是x(n)的復(fù)共軛序列，長(zhǎng)度為N。若X(k)=DFT[x(n)],則

DFT[x*(n)]=X*(N－k)，0≤k≤N－1(3.2.15)注意，當(dāng)k=0時(shí)，X*(N－k)=X*(N)，已經(jīng)超出主值區(qū)間，我們?cè)俅螐闹芷谘油氐母拍钊ダ斫猓?/p>

就有X*(N)=X*(0)。

證明如下：有了式(3.2.15)，就可以推導(dǎo)出許多有用的對(duì)稱(chēng)性質(zhì)。復(fù)序列的實(shí)部Re｛x(n)｝的DFT：

(3.2.16)Xe(k)稱(chēng)為X(k)的共軛偶對(duì)稱(chēng)分量。復(fù)序列的虛部Im{x(n)}的DFT：(3.2.17)

Xo(k)稱(chēng)為X(k)的共軛奇對(duì)稱(chēng)分量(也稱(chēng)共軛反對(duì)稱(chēng))。復(fù)序列x(n)的共軛偶對(duì)稱(chēng)分量與共軛奇對(duì)稱(chēng)分量如圖3.2.6所示。圖3.2.6x(n)的共軛偶對(duì)稱(chēng)分量與共軛奇對(duì)稱(chēng)分量現(xiàn)在討論Xe(k)與Xo(k)兩個(gè)分量本身的對(duì)稱(chēng)性。對(duì)于將式中R以N－k代入，得(3.2.18)得到說(shuō)明Xe(k)具有共軛偶對(duì)稱(chēng)特點(diǎn)。類(lèi)似地，可得到(3.2.19)即Xo(k)具有共軛反對(duì)稱(chēng)特點(diǎn)，稱(chēng)其為X(k)的共軛奇對(duì)稱(chēng)分量。現(xiàn)實(shí)中采集的時(shí)域序列都是純實(shí)數(shù)序列x(n)，即x(n)=x*(n)；x(n)=Re{x(n)}，那么，頻率X(k)只有共軛偶對(duì)稱(chēng)部分，即X(k)=Xe(k)，而Xo(k)=0。這表明實(shí)數(shù)序列的DFT具有共軛對(duì)稱(chēng)性。利用這一特性，意味著只要知道一半數(shù)目的X(k)，k=0,1，2，…，N/2－1，就可得到另一半的X(k)，k=N/2，…，N－1。這一特點(diǎn)在求DFT時(shí)可以加以利用，省去一半運(yùn)算量,以提高運(yùn)算效率。根據(jù)DFT的對(duì)偶特性，我們也可以找到頻譜X(k)的實(shí)部Re{X(k)}、虛部Im{X(k)}與序列x(n)的共軛偶部xe(n)與共軛奇部xo(n)的關(guān)系。

分別以xe(n)及xo(n)表示序列x(n)的圓周共軛偶部與圓周共軛奇部，即應(yīng)把N點(diǎn)序列看成是周期延拓后的序列，顯然，x(N)=x(0)=x(N－0)，也就是從所謂的圓周意義上來(lái)理解?？勺C明：(3.2.20)(3.2.21)

(4)選頻特性。

對(duì)復(fù)指數(shù)函數(shù)x(t)=ejrω0t(r為整數(shù))，進(jìn)行采樣得復(fù)指數(shù)序列(余弦和正弦)

x(n)=ejrω0n

當(dāng)任意頻率取ω0=2π/N時(shí)，x(n)的離散傅立葉變換為(3.2.22)現(xiàn)在來(lái)深入理解這個(gè)現(xiàn)象。我們把【例3.2.1】中的2kHz頻率分量改成2.3kHz，且為了看得更清楚，去掉6kHz分量。即x(t)=1×sin(2π×2300×t)信號(hào)，經(jīng)過(guò)T=1/16ms采樣，用N=16個(gè)數(shù)據(jù)計(jì)算出的頻譜幅度結(jié)果，如圖3.2.7所示。虛線是2.3kHz信號(hào)被截取后的連續(xù)頻譜圖。

本來(lái)期望于2.3kHz的地方出現(xiàn)幅度為1×N/2=8的頻譜，可是沒(méi)有?，F(xiàn)在由于DFT的選頻特點(diǎn)，因?yàn)閒s/N=1kHz，它只會(huì)表達(dá)0~15kHz整數(shù)頻率點(diǎn)信號(hào)，對(duì)2.3kHz是沒(méi)有頻率單元與之對(duì)應(yīng)的，最接近的是2kHz單元，它獲得最大的輸出約為6.5；并且所有各頻率分析點(diǎn)都不位于連續(xù)頻譜的過(guò)零點(diǎn)而觀察到了泄漏的情況，特別是第一、第二旁瓣的幅度。要注意的是，這條虛線不是sinc(x)或Sa(x)函數(shù)，它已經(jīng)包含了周期化后的旁瓣的互相串?dāng)_，也可以說(shuō)是“混疊的Sa(x)”。圖3.2.7頻譜泄漏現(xiàn)象的揭示再看該X(k)所對(duì)應(yīng)的采樣數(shù)據(jù)x(n)。如圖3.2.8所示，注意觀察那些小幅值點(diǎn)，可以發(fā)現(xiàn)數(shù)據(jù)并沒(méi)有出現(xiàn)完整的周期性規(guī)律，換句話說(shuō)，周期信號(hào)采樣應(yīng)該采集到一個(gè)完整周期里的N個(gè)數(shù)據(jù)，然后用這個(gè)序列片段通過(guò)周期延拓后，恰好與真實(shí)的周期信號(hào)序列相同，這才合乎邏輯。正因?yàn)镈FT本身隱含著周期化處理，圖中這個(gè)16點(diǎn)的x(n)周期化后已不是x(nT)=sin(2π×2300×nT)，

－∞＜n＜∞。也可以理解成另外某個(gè)信號(hào)，當(dāng)然頻譜就肯定不同。請(qǐng)讀者比較圖3.2.9(a)和圖3.2.9(b)兩種情況。同一個(gè)信號(hào)，截取的點(diǎn)數(shù)(長(zhǎng)度)不同，(a)圖截取2個(gè)完整周期，經(jīng)過(guò)周期化后跟原周期序列是一樣的，而(b)圖只截取1.75個(gè)周期，周期化時(shí)在“接頭”的地方將出現(xiàn)跳變，已非原來(lái)的正弦序列。圖3.2.82.3kHz正弦信號(hào)以16kHz采樣率取樣的數(shù)據(jù)圖3.2.9正弦序列的兩種截?cái)嗲闆r如何才能更真實(shí)地獲得結(jié)果呢？假如我們繼續(xù)增長(zhǎng)對(duì)2.3kHz信號(hào)的觀察時(shí)間，以1/16ms的間隔獲得更多的采樣數(shù)據(jù)，直到序列有一個(gè)完整周期。比如N=160，那么做DFT時(shí)，其頻率分析點(diǎn)間隔是fs/N=16kHz/160=0.1kHz，即每點(diǎn)遞增100Hz，第23點(diǎn)就恰好準(zhǔn)確地觀察到2.3kHz信號(hào)最高幅度，其他都為0，避開(kāi)了泄漏現(xiàn)象引起的幅度誤判。(盡管泄漏減小了些，但始終存在著，只是沒(méi)有觀察到，考慮一下有什么辦法能看到第一旁瓣的大致幅度？)我們還可以在圖3.2.8的16點(diǎn)序列后面添加0直到160點(diǎn)，再計(jì)算DFT，這實(shí)際上是增加并調(diào)整DFT頻率分析點(diǎn)位置和密度，使得2.3kHz上有一個(gè)分析點(diǎn)，這樣也能看到頻譜在第23點(diǎn)處有一個(gè)最大幅度，但它將比前面那種方法要小些。還有一點(diǎn)不同，這個(gè)方法在其他的分析頻率點(diǎn)也都有輸出，并不為0，反映的正是泄漏的那些旁瓣大小。值得指出的是，我們?cè)缫阎婪侵芷谛蛄行盘?hào)的DTFT頻譜是連續(xù)頻率的，在其頻帶范圍里分布著無(wú)限多的頻率成分。顯然用DFT計(jì)算這種序列時(shí)，無(wú)論N取多大，總有些頻率不能剛好對(duì)應(yīng)到有限個(gè)輸出頻率單元上，結(jié)果一定會(huì)受到上述泄漏現(xiàn)象的影響，它的擴(kuò)散到其他頻率點(diǎn)的特征很容易引起我們的誤判。實(shí)際的信號(hào)序列頻譜都是如此，雖然有許多辦法能減少泄漏以提高DFT計(jì)算精確度，但頻譜泄漏卻是不可避免的。泄漏的問(wèn)題歸根結(jié)底是無(wú)始無(wú)終的信號(hào)被截取成一段，即乘以RN(n)=U(n)－U(n－N)而造成的。

(5)DFT與Z變換的關(guān)系。

有限長(zhǎng)的序列總存在Z變換：與DFT定義對(duì)比，就可以看出有如下關(guān)系成立：式中z=ej(2π/N)k，說(shuō)明是z平面單位圓上相角為(2π/N)k(k=0，1，…，N－1)的N個(gè)點(diǎn)，X(k)就是x(n)在這些點(diǎn)上的Z變換。如果結(jié)合序列的DTFT頻譜來(lái)理解，就是連續(xù)譜在頻率域被離散化。頻率間隔是2π/N，如圖3.2.10所示是頻率點(diǎn)分布。幅頻則如圖3.2.11所示，虛線表示的包絡(luò)線就是有限長(zhǎng)序列的周期頻譜的主周期圖形，即(3.2.23)(3.2.24)圖3.2.10單位圓上等分點(diǎn)處的Z變換即DFT3.2.3頻率域采樣

N點(diǎn)時(shí)域序列x(n)，其DTFT是ω的連續(xù)函數(shù)，即頻譜X(ejω)；而我們用DFT計(jì)算x(n)時(shí)是N點(diǎn)的X(k)，它是連續(xù)頻譜X(ejω)的頻率域N點(diǎn)等間隔采樣，如圖3.2.11所表達(dá)的。

現(xiàn)在要討論的問(wèn)題是：設(shè)頻率域0~2π間的均勻采樣點(diǎn)數(shù)為M，它可不可以比N小？或者比N大？我們?cè)僖淮螐腄FT對(duì)偶特性來(lái)分析。圖3.2.11序列的DTFT連續(xù)頻譜與DFT離散頻譜關(guān)系時(shí)域里連續(xù)信號(hào)被采樣成離散信號(hào)時(shí)，會(huì)使得頻譜發(fā)生周期化。時(shí)域采樣間隔T決定了頻譜周期化的周期大小，即fs=1/T，為防止頻譜混疊發(fā)生，fs應(yīng)足夠大，大到超過(guò)信號(hào)帶寬。同樣，頻域里連續(xù)頻譜被采樣成等間隔離散頻率點(diǎn)，即彼此呈諧波關(guān)系，而使得時(shí)域?qū)?yīng)表現(xiàn)為周期化。頻率采樣間隔F決定了時(shí)域周期化的周期大小，即ts=1/F，為防止時(shí)域混疊發(fā)生，ts應(yīng)足夠大，大于信號(hào)長(zhǎng)度。顯然，頻域采樣點(diǎn)數(shù)M=2π/F越大，其頻率間隔F越小，時(shí)域周期化的周期長(zhǎng)度ts越長(zhǎng)，它換算成時(shí)域序列點(diǎn)數(shù)ts/T=Ns就越大，只要原連續(xù)頻譜對(duì)應(yīng)的時(shí)域序列點(diǎn)數(shù)N少于Ns，就不會(huì)發(fā)生時(shí)域序列的混疊。這就是頻率域采樣定理。圖3.2.12X(ejω)的32個(gè)頻譜樣點(diǎn)的幅度和相位圖

【例3.2.2】

頻率域取樣的例子。

一個(gè)連續(xù)的頻譜X(ejω)在一周期里等間隔取樣了32個(gè)頻率數(shù)據(jù)X(k)，X(k)={40，－32.6－j30.1,－14.3+j38.1,20.6+j7.0，0.3－j0.7,…,－32.6+j30.1}，如圖3.2.12所示。經(jīng)過(guò)IDFT逆變換后得到對(duì)應(yīng)的時(shí)域序列：

x(n)=｛2，－1，－3，－5，－2，2.1204e－016，1，2，4，5，7，9，8，6，3，1，1，2，1.2204e－016，

0，1.304e－016，0，0，0，0，0，0，0，0，0，0，0｝

如圖3.2.13所示，注意x(0)=2，x(1)=－1，…，x(16)=1，x(17)=2，從x(18)起都為0，說(shuō)明頻譜所對(duì)應(yīng)的x(n)是有限長(zhǎng)N=18點(diǎn)的時(shí)間序列。

若對(duì)X(ejω)在一周期里等間隔取樣了16個(gè)頻率數(shù)據(jù)X1(k)，就是等效從32點(diǎn)的X(k)中抽取偶數(shù)序號(hào)的頻譜點(diǎn)構(gòu)成，如圖3.2.14所示。

X1(k)={40，－14.3+j38.1，0.3－j0.7，…，－14.3－j38.1}圖3.2.13X(k)的IDFT后的序列x(n)幅度圖圖3.2.14X(ejω)的16個(gè)頻譜樣點(diǎn)的幅度和相位圖對(duì)應(yīng)的IDFT后得到的時(shí)間序列x1(n)為16點(diǎn)：

x1(n)={3，1，－3，－5，－2，2.1204e-016，1，2，4，

5，7，9，8，6，3，1}3.2.4循環(huán)卷積定理

所謂循環(huán)卷積的概念是為了配合有限長(zhǎng)序列的DFT性質(zhì)和應(yīng)用而引入的，即頻率域相乘操作對(duì)應(yīng)于時(shí)間域的卷積運(yùn)算。我們知道，兩個(gè)周期相同的序列可以在一個(gè)周期內(nèi)進(jìn)行卷積運(yùn)算，稱(chēng)為周期卷積，其卷積結(jié)果依然保有周期性?，F(xiàn)在對(duì)于兩個(gè)一樣長(zhǎng)度的有限長(zhǎng)序列先周期延拓后再進(jìn)行卷積運(yùn)算，就稱(chēng)為循環(huán)卷積。它類(lèi)同于周期卷積，只不過(guò)它的周期性是延拓想像出來(lái)的，循環(huán)卷積的結(jié)果仍然是個(gè)同長(zhǎng)度的有限長(zhǎng)序列。

假設(shè)x(n)、y(n)是兩個(gè)長(zhǎng)度為N的有限長(zhǎng)序列，它們的N點(diǎn)DFT分別為X(k)、Y(k)，若F(k)=X(k)Y(k)，那么

(3.2.25)運(yùn)算成立。

首先把x(n)、y(n)周期性延拓成周期序列x(n)、y(n)，再把式(3,2.25)看做是周期序列x(n)和y(n)的周期卷積后再取其主值序列。將F(k)周期延拓，簡(jiǎn)便記為F((k))N=X((k))NY((k))N,對(duì)應(yīng)的周期卷積式為~~~~因?yàn)樵?≤m≤N－1內(nèi)，x((m))N=x(m)，可以看到式(3.2.25)的運(yùn)算過(guò)程與周期卷積是一樣的，只是最后僅取結(jié)果的主值序列。由于卷積過(guò)程只在主值區(qū)間0≤m≤N－1內(nèi)進(jìn)行，因此對(duì)于y((n－m))N實(shí)際上就是y(m)的循環(huán)移位，稱(chēng)為“循環(huán)卷積”，以區(qū)別于線性卷積及周期卷積。循環(huán)卷積習(xí)慣上用一個(gè)數(shù)字外加一個(gè)圓圈來(lái)表示，數(shù)字表示參與卷積的序列長(zhǎng)度。兩個(gè)序列的循環(huán)卷積計(jì)算方法：

(1)由兩個(gè)有限長(zhǎng)序列x(n)、y(n)延拓構(gòu)造出周期序列x((n))N和y((n))N。

(2)計(jì)算x((n))N和y((n))N的周期卷積f((n))N。

(3)取卷積結(jié)果f((n))N主值，即f(n)=x(n)○y(n)=f((n))NRN(n)。

我們用一個(gè)例子來(lái)驗(yàn)證DFT的循環(huán)卷積的性質(zhì)。

【例3.2.3】

有兩個(gè)長(zhǎng)度都為6點(diǎn)的序列x(n)和y(n)，其頻譜分別記X(k)和Y(k)。要求驗(yàn)證DFT循環(huán)卷積性質(zhì)，x(n)=｛－2，5，－1，3，4，7｝和y(n)=｛1，2，7，3，4，6｝。N解分析：把y(n)周期化，有

y((n))6=｛…,4,6,1,2,7,3,4,6,1,2,7,3,4,…｝

↑

對(duì)y(0)=1處左右翻轉(zhuǎn)后成為

y((－m))6={…，6，4，3，7，2，1，6，4，3，7，2，1，…}

↑

按照定義式子，計(jì)算得

f((n))6={…,75,68,50,82,52,41,75,68,50,82,52,41,75,…}

↑

取主值序列，得：

f(n)=f((n))6R6(n)={75，68，50，82，52，41}由例3.2.3說(shuō)明，DFT的卷積性質(zhì)是屬于循環(huán)卷積的?，F(xiàn)在用表3.2.1總結(jié)一下我們所講的三種卷積的異同點(diǎn)。

【例3.2.4】

線性卷積的數(shù)據(jù)來(lái)源于例3.2.3。不妨將y(n)看成是某離散系統(tǒng)的單位脈沖響應(yīng)，x(n)是其輸入，那么系統(tǒng)的零狀態(tài)響應(yīng)就是二者的線性卷積，即x(n)*y(n)=f(n)。調(diào)用MATLAB信號(hào)處理的內(nèi)部函數(shù)conv(x,y)，它計(jì)算兩個(gè)序列線性卷積，即可得到6+6－1=11個(gè)點(diǎn)的輸出響應(yīng)數(shù)據(jù)：

f(n)=｛－2，1，－5，30，10，41，77，67，55，52，42｝

現(xiàn)我們?cè)趚(n)后面添加5個(gè)0，使得序列成為11個(gè)點(diǎn)，即

x(n)=｛－2，5，－1，3，4，7，0，0，0，0，0｝；n=0~10

然后DFT求出X(k)，k=0~10。用同樣辦法構(gòu)造y(n)，也在其后面添5個(gè)0，再經(jīng)過(guò)DFT得到Y(jié)(k)。

最后求IDFT｛X(k)Y(k)｝而得到輸出響應(yīng)f(n)=x(n)*y(n)。結(jié)果與直接卷積conv(x,y)一樣，從而實(shí)現(xiàn)了用DFT求取系統(tǒng)響應(yīng)的目的。圖3.2.15用DFT計(jì)算線性卷積為什么序列經(jīng)過(guò)添0處理后再進(jìn)行循環(huán)卷積就會(huì)有線性卷積的結(jié)果呢？對(duì)照?qǐng)D3.2.15，如果把添0后的y(n)周期化成y((n))11再翻折，然后開(kāi)始移位乘加運(yùn)算?？梢园l(fā)現(xiàn)，所添的0在這個(gè)過(guò)程中恰好把原有的y(n)的前后周期都隔離開(kāi)，加大了循環(huán)移位的周期，配合上x(chóng)(n)所添的0，效果就好比做線性卷積時(shí)有值序列以外都用0替代一樣，從而避免了y(n)延拓后的主值周期以外的數(shù)據(jù)被循環(huán)進(jìn)入乘加運(yùn)算，達(dá)到利用循環(huán)卷積計(jì)算線性卷積的目的。

一般地，一個(gè)N點(diǎn)的序列x(n)和一個(gè)M點(diǎn)的序列h(n)的線性卷積結(jié)果是長(zhǎng)度為L(zhǎng)=N+M－1的序列，那么，進(jìn)行x(n)后添M－1個(gè)0而y(n)后添N－1個(gè)0的預(yù)處理，再進(jìn)行L點(diǎn)循環(huán)卷積就可以得到線性卷積的結(jié)果。

3.3快速傅立葉變換FFT

3.3.1減少運(yùn)算量的思路

我們從DFT定義來(lái)分析有限長(zhǎng)N點(diǎn)序列x(n)進(jìn)行一次DFT運(yùn)算所需的運(yùn)算量。如k=5時(shí)的X(5)，(3.3.1)一般地，x(n)和WnkN都是復(fù)數(shù)，因此，計(jì)算X(5)的值時(shí)，要進(jìn)行N次復(fù)數(shù)相乘和N－1次復(fù)數(shù)相加。X(k)有k=0，1，2，…，N－1總共N個(gè)點(diǎn)，故完成全部的DFT運(yùn)算，就需要N×N次復(fù)數(shù)相乘和N×(N－1)次復(fù)數(shù)相加。進(jìn)一步，從(a+jb)(c+jd)=(ac－bd)+j(ad+cb)可知,一次復(fù)數(shù)相乘包含了4次實(shí)數(shù)乘法和2次實(shí)數(shù)加法,而一次復(fù)數(shù)加法包含了2次實(shí)數(shù)加法。所以，完成一次N點(diǎn)DFT需要4×N2次實(shí)數(shù)乘法和2×N2+2×N×(N－1)=4×N2－2×N次實(shí)數(shù)加法。我們知道，在計(jì)算機(jī)中，加法運(yùn)算要比乘法運(yùn)算快得多，相比之下加法可以忽略不計(jì)，那么，DFT的運(yùn)算量就基本上約等于4N2次實(shí)數(shù)乘法。當(dāng)N=1024點(diǎn)時(shí)，該乘法量將達(dá)到四百多萬(wàn)的驚人次數(shù)。

慶幸的是，算式中乘數(shù)之一——的序列值具有周期重復(fù)和對(duì)稱(chēng)特點(diǎn)，令q=nk，那么整數(shù)q的范圍是0~(N－1)2，只有q的值在0~N－1之間時(shí)，WnkN實(shí)部或虛部才是獨(dú)立的數(shù)，q為其他整數(shù)時(shí)，都將被重復(fù)。例如N=8的情況，整數(shù)q的取值為0，1，2，3，…，49。WnkN的實(shí)部和虛部如圖3.3.1所示。圖3.3.1N=8時(shí)WNnk的實(shí)部和虛部序列可以看出，q=0~7是正弦主值周期，值為｛0，0.707，1，0.707，0，-0.707，-1，-0.707｝才是獨(dú)立的，其余的q值所對(duì)應(yīng)的正弦值都已在主值周期里出現(xiàn)過(guò)。由于對(duì)稱(chēng)性，正弦主值區(qū)間里q=1和q=3的序列值是相等的，q=4和q=6的序列值也是相等的，并且還與q=1的序列值僅差一負(fù)號(hào)，剩下特殊序列值是±1和0。余弦也類(lèi)似，主值為{1，0.707，0，-0.707，-1，-0.707，0，0.707}。那么x(n)與這些相同的值做乘法運(yùn)算時(shí)就可以化簡(jiǎn)合并，例如式(3.3.1)的X(5)中的某2項(xiàng)：

…+x(1)W58+…+x(5)W258+…

=…+x(1)(－0.707+j0.707)+…+x(5)(0.707-j0.707)+…

=[x(1)－x(5)](－0.707+j0.707)+…付出1次x(1)與x(5)的加法而節(jié)省了1次復(fù)數(shù)乘法。實(shí)際上，正弦函數(shù)只要知道1/4周期的值就可全部確定，因?yàn)槠浜瘮?shù)值在0~π/2和π/2~π范圍內(nèi)是關(guān)于π/2點(diǎn)左右對(duì)稱(chēng)的，而在0~π和π~2π范圍內(nèi)又關(guān)于橫軸正負(fù)半周對(duì)稱(chēng)。恰恰是這些特性的巧妙利用，成功地減少了DFT運(yùn)算量，才造就了FFT。3.3.2基2-FFT算法

1.基2時(shí)間抽取(DIT)的FFT

假定序列x(n)的長(zhǎng)度為N點(diǎn)，N=2M，M為正整數(shù)。

我們將x(n)依照序號(hào)分解為兩組子序列，一個(gè)為偶數(shù)號(hào)項(xiàng)組成，另一個(gè)為奇數(shù)號(hào)項(xiàng)組成:(3.3.2)因此，求x(n)的DFT如下：(3.3.3)而X(k)的后半段k=N/2，N/2+1，N/2+2，…，N－1，可由周期性和對(duì)稱(chēng)性求得。

因?yàn)閄1(k+N/2)=X1(k)，X2(k+N/2)=X2(k)，所以(3.3.4)我們由此得出結(jié)論：一個(gè)N點(diǎn)的DFT運(yùn)算任務(wù)X(k)被分解為兩個(gè)N/2點(diǎn)的DFT運(yùn)算，這兩個(gè)較少點(diǎn)的X1(k)和X2(k)可以再進(jìn)一步組合，獲得一個(gè)N點(diǎn)的X(k)。當(dāng)然隨后的這個(gè)組合運(yùn)算也是要花時(shí)間的，不過(guò)這個(gè)開(kāi)銷(xiāo)很小。因此，基本上可以認(rèn)為，這個(gè)按奇偶分組的方法，其計(jì)算效率提高了將近一倍。如式(3.3.3)和式(3.3.4)的組合運(yùn)算過(guò)程，可以用如圖3.3.2或圖3.3.3所示的蝶形信號(hào)流圖表達(dá)。圖3.3.2蝶形信號(hào)流圖圖3.3.3另一形式的蝶形圖對(duì)于X1(k)或X2(k)的求解，自然會(huì)想到繼續(xù)如法炮制，比如將其原序列x1(r)再次按奇偶分組，計(jì)算2個(gè)N/4點(diǎn)的更短DFT，然后組合得到X1(k)。這個(gè)過(guò)程可以一直下去，最后到2個(gè)點(diǎn)的DFT運(yùn)算：這也可以歸結(jié)為x(0)和x(1)兩點(diǎn)輸入經(jīng)過(guò)如圖3.3.2所示蝶形運(yùn)算后輸出的兩個(gè)點(diǎn)值，這恰是基-2名稱(chēng)的由來(lái)。我們以圖3.3.4所示N=8為例子，說(shuō)明基-2的DIT算法過(guò)程。如圖3.3.4所示，將x(n)按奇偶序號(hào)分成2個(gè)4點(diǎn)的子序列，分別用4點(diǎn)DFT算出中間的值A(chǔ)(k)和B(k)；經(jīng)過(guò)蝶形組合運(yùn)算得到8個(gè)點(diǎn)的DFT值X(k)。注意，A(0)與B(0)是一雙對(duì)偶節(jié)點(diǎn)，下節(jié)點(diǎn)B(0)要乘以WNk，k的值與上節(jié)點(diǎn)A(0)序號(hào)一致，其他類(lèi)推。我們現(xiàn)在對(duì)虛線框里的4點(diǎn)DFT進(jìn)行再次分組，它的輸入序列｛x(0)，x(2)，x(4)，x(6)｝按其位置奇偶分開(kāi)成偶序號(hào)組｛｛x(0)，x(4)｝和奇序號(hào)組｛x(2)，x(6)｝，對(duì)｛x(1)，x(3)，x(5)，x(7)｝也同樣進(jìn)行，如圖3.3.5所示。圖3.3.4兩個(gè)4點(diǎn)DFT組合成一個(gè)8點(diǎn)DFT圖3.3.54個(gè)2點(diǎn)DFT組合成2個(gè)4點(diǎn)DFT全部用流圖繪出如圖3.3.6所示。圖3.3.6基2-DIT的8點(diǎn)DFT運(yùn)算流圖相對(duì)于序號(hào)的自然規(guī)律稱(chēng)為正序排列，我們把這種看似混亂的輸入序列秩序稱(chēng)為倒序，英文為bitreversal，其本意是“按位倒轉(zhuǎn)”。它指的是這樣一種過(guò)程：一個(gè)自然序列x(n)，它的值依照序號(hào)n=0，1，2，3，…，N－1順序地排列?，F(xiàn)在將序號(hào)用二進(jìn)制表達(dá)，根據(jù)長(zhǎng)度N的大小，顯然需要不同位數(shù)(bit)的二進(jìn)制，N=2M，則需要M位才能表達(dá)，比如N=32，就需要5位二進(jìn)制。以圖3.3.6為例，序號(hào)0~7，可用3個(gè)bit的二進(jìn)制表達(dá)。如表3.3.1所示，把原來(lái)的自然序x(n)的每一個(gè)值，根據(jù)其序號(hào)的對(duì)應(yīng)倒序位置重新排列，所得到的新序列就是倒序規(guī)律的。值得注意的是，同一個(gè)整數(shù)在不同長(zhǎng)度N里，其倒序位置是不同的。排列過(guò)程并不復(fù)雜，比如，x(0)的序號(hào)0，其倒序號(hào)也是0，則x(0)還是放在n=0的位置；而x(1)的序號(hào)1對(duì)應(yīng)的倒序號(hào)是4，那么x(1)應(yīng)該和x(4)號(hào)互相調(diào)換位置，如此類(lèi)推。應(yīng)該注意，這樣調(diào)換的工作按自然序逐個(gè)進(jìn)行，只需完成到N/2序號(hào)就結(jié)束，否則，會(huì)發(fā)生重復(fù)調(diào)換又回到自然序的情況。MATLAB的信號(hào)處理工具箱里的bitrevorder(x)函數(shù)就是用來(lái)完成這個(gè)工作的，它把自然序的數(shù)據(jù)x排列調(diào)整成倒序排列。讀者可以試一下。

因此，我們按照流圖3.3.6進(jìn)行運(yùn)算之前，需要對(duì)數(shù)據(jù)進(jìn)行一次倒序排列。不過(guò)，也可以對(duì)流圖進(jìn)行重新整理，通過(guò)升降流圖中的某些水平橫向的整條支路，把輸入x(n)調(diào)整成自然序，輸出頻譜序列X(k)則將成為倒序排列。通過(guò)整理后如圖3.3.7所示，它是輸入正序輸出倒序的基-2時(shí)間抽取FFT的另一流圖。圖3.3.7輸入正序輸出倒序的基-2時(shí)間抽取8點(diǎn)FFT流圖圖中的中間變量C、D、A、B實(shí)際上是無(wú)所謂的，只是個(gè)標(biāo)記。最后得到的輸出頻譜是倒序排列的，它可以通過(guò)倒序整理環(huán)節(jié)，成為自然序。如果現(xiàn)在僅調(diào)整第3級(jí)的輸出點(diǎn)位置也成為自然序，可以想像，蝶形圖就會(huì)失去同址存儲(chǔ)的規(guī)律，不利于編程。我們?cè)谶@里又一次體會(huì)了相同功能可以用不同的流程實(shí)現(xiàn)的概念。以下再介紹另一種常用的FFT算法。

2.基2頻率抽取的FFT(DIF)

與DIT將輸入數(shù)據(jù)按序號(hào)奇偶分組不同，頻率抽取方法DIF則是將輸入序列x(n)按前后對(duì)半斷開(kāi)，成為兩個(gè)短序列部分。這樣便將N點(diǎn)DFT做成前后兩截:奇頻組，當(dāng)k=2r+1時(shí)，得(3.3.7)顯然，式(3.3.6)和式(3.3.7)都是N/2點(diǎn)的DFT，被變換的短序列a(n)是原x(n)的前后兩段直接對(duì)加,而短序列b(n)稍微復(fù)雜些,是x(n)的前后兩段對(duì)減再乘以WnN,這可以用流圖3.3.8表示。如果對(duì)a(n)或b(n)也分別進(jìn)行前后截兩段及相應(yīng)的處理，那么就是重復(fù)上述過(guò)程，而成為4個(gè)N/4點(diǎn)的DFT。如此類(lèi)推，最后是2點(diǎn)的蝶形運(yùn)算。仍然以N=8為例子，如圖3.3.9所示，輸入序列是正序，輸出頻譜是倒序排列的，使用時(shí)可對(duì)結(jié)果再行排序處理。圖3.3.8頻率抽取的蝶形圖圖3.3.9N=8時(shí)的DIF頻率抽取流圖這里的WqN類(lèi)型也是N/2個(gè)，全部計(jì)算的乘法次數(shù)和DIT相同。比如，直流分量X(0)的計(jì)算，X(0)=c(0)+c(1)=[a(0)+a(2)]+[a(1)+a(3)]=x(0)+x(4)+x(1)+x(5)+x(2)+x(6)+x(3)+x(7)。事實(shí)上，對(duì)于直流，DFT就是把所有N點(diǎn)數(shù)據(jù)直接相加，

。

同樣地，對(duì)圖3.3.9進(jìn)行調(diào)整，為了輸出X(k)按自然序排列，流圖中整條水平橫線支路上下調(diào)換位置，就可以得到輸入倒序而輸出正序的DIF另一種形式。如圖3.3.10所示，圖中間變量只供說(shuō)明用。細(xì)心的讀者可能已經(jīng)發(fā)現(xiàn)，DIF和DIT的流圖規(guī)律很相似。實(shí)際上，圖3.3.9和圖3.3.6是互易的，或稱(chēng)流圖轉(zhuǎn)置關(guān)系，把輸入名和輸出名對(duì)調(diào)，流向倒轉(zhuǎn)但增益和結(jié)構(gòu)不變，就可互相轉(zhuǎn)換。圖3.3.10和圖3.3.7也是互易的。圖3.3.10N=8時(shí)DIF的另一種流圖形式3.3.3N為組合數(shù)的FFT算法

以上討論的都是以2為基數(shù)的FFT算法，即N=2M。實(shí)際應(yīng)用時(shí)，有限長(zhǎng)序列的長(zhǎng)度N很大程度上由人為因素確定，因此多數(shù)場(chǎng)合可取N=2M，從而直接使用以2為基數(shù)的FFT算法。

如N不能人為確定，N的數(shù)值也不是以2為基數(shù)的整數(shù)次方，那么處理方法有以下兩種：

(1)補(bǔ)零：將x(n)補(bǔ)零，使N=2M。

例如N=30，補(bǔ)上x(chóng)(30)=x(31)=0兩點(diǎn)，使N=32=25，這樣可直接采用以2為基數(shù)M=5的FFT程序。有限長(zhǎng)序列補(bǔ)零后并不影響其頻譜X(ejw)，只是頻譜的采樣點(diǎn)數(shù)增加了，上例中由30點(diǎn)增加到32點(diǎn)，所以在許多場(chǎng)合這種處理是可接受的。

(2)如要求準(zhǔn)確的N點(diǎn)DFT值，可采用任意數(shù)為基數(shù)的FFT算法，其計(jì)算效率低于以2為基數(shù)的FFT算法。

如N為復(fù)合數(shù)，可分解為兩個(gè)整數(shù)P與Q的乘積，像前面以2為基數(shù)時(shí)一樣，F(xiàn)FT的基本思想是將DFT的運(yùn)算盡量減小。因此，在N=PQ的情況下，也希望將N點(diǎn)的DFT分解為P個(gè)Q點(diǎn)DFT或Q個(gè)P點(diǎn)DFT，以減少計(jì)算量。算法步驟是，先將n，k寫(xiě)成：式中：n0,k1分別為0,1,…,Q－1;n1,k0分別為0,1,…,P－1。

N點(diǎn)DFT可以重新寫(xiě)為(3.3.8)可以寫(xiě)成(3.3.9)再考慮到WNk0n1Q=WPk0n1，代入式(3.3.9)，得(3.3.10)令則有(3.3.11)由式(3.3.10)、式(3.3.11)可見(jiàn)，N點(diǎn)DFT變成Q個(gè)P點(diǎn)DFT、P個(gè)Q點(diǎn)DFT兩級(jí)運(yùn)算。下面以P=3，Q=4，N=12為例，說(shuō)明其算法處理過(guò)程，如圖3.3.11所示。

(1)先將x(n)通過(guò)

x(n1Q+n0)改寫(xiě)成x(n1，n0)。因?yàn)镼=4,n1=0、1、2，n0=0、1、2、3，故輸入是按自然順序的，即：

x(0,0)=x(0)x(0,1)=x(1)x(0,2)=x(2)x(0,3)=x(3)

x(1,0)=x(4)x(1,1)=x(5)x(1,2)=x(6)x(1,3)=x(7)

x(2,0)=x(8)x(2,1)=x(9)x(2,2)=x(10)x(2,3)=x(11)

(2)求Ｑ個(gè)Ｐ點(diǎn)的DFT。

(３)X1(k0，n0)乘以WＮk0n0得到X1′(k0，n0)。

(４)求P個(gè)Ｑ點(diǎn)的DFT，參變量是k0，

(５)將X2(k0，k1)通過(guò)X(k0+k1P)恢復(fù)為X(k)。圖3.3.11N=12時(shí)的FFT算法

N為組合數(shù)時(shí)的FFT運(yùn)算量分析：

(１)求Ｑ個(gè)P點(diǎn)DFT需要ＱP2次復(fù)數(shù)乘法和Ｑ·P·(P－1)次復(fù)數(shù)加法。

(２)乘Ｎ個(gè)W因子需要Ｎ次復(fù)數(shù)乘法。

(３)求P個(gè)Ｑ點(diǎn)DFT需要PQ2次復(fù)數(shù)乘法和P·Ｑ(Ｑ－1)次復(fù)數(shù)加法。

總的復(fù)數(shù)乘法量：QP2+N+PQ2=N(P+Q+1)

總的復(fù)數(shù)加法量：Q·P(P－1)+P·Q·(Q－1)=N(P+Q－2)3.3.4Chirp-Z變換

采用FFT可以算出全部N點(diǎn)DFT值，即Z變換X(z)在z平面單位圓上的等間隔取樣值。問(wèn)題的提出：

(1)不需要計(jì)算整個(gè)單位圓上Z變換的取樣，如對(duì)于窄帶信號(hào)，只需要對(duì)信號(hào)所在的一段頻帶進(jìn)行分析。這時(shí)，希望頻譜的采樣集中在這一頻帶內(nèi)，以獲得較高的分辨率，而頻帶以外的部分可不考慮。

(2)對(duì)其他圍線上的Z變換取樣感興趣，例如語(yǔ)音信號(hào)處理中，需要知道Z變換的極點(diǎn)所在頻率，如極點(diǎn)位置離單位圓較遠(yuǎn)，則其單位圓上的頻譜就很平滑；如果采樣不是沿單位圓而是沿一條接近這些極點(diǎn)的弧線進(jìn)行，則極點(diǎn)所在頻率上將出現(xiàn)明顯的尖峰，由此可較準(zhǔn)確地測(cè)定極點(diǎn)頻率。

(3)要求能有效地計(jì)算當(dāng)N是素?cái)?shù)時(shí)，序列的DFT。

算法原理：

已知x(n)，0≤n≤N－1，令zk=AW－k，k=0，…，M－1，M為采樣點(diǎn)數(shù)，A、W為任意復(fù)數(shù)，(3.3.12)式(3.3.12)中：A0表示起始取樣點(diǎn)的半徑長(zhǎng)度，通常A0≤1；θ0表示起始取樣點(diǎn)z0的相角；φ0表示兩相鄰點(diǎn)之間的等分角；W0為螺旋線的伸展率，W0<1則線外伸，W0>1則線內(nèi)縮(反時(shí)針)，W0=1則表示半徑為A0的一段圓弧，若A0=1，則這段圓弧是單位圓的一部分，如圖3.3.12所示。圖3.3.12Z平面上的頻率采樣形式計(jì)算Z變換在采樣點(diǎn)zk的值：(3.3.13)顯然，按照以上公式計(jì)算出全部M點(diǎn)采樣值需要NM

次復(fù)乘和(N－1)M次復(fù)加，當(dāng)N及M較大時(shí)，計(jì)算量迅速增加，以上運(yùn)算可轉(zhuǎn)換為卷積形式，從而可采用FFT進(jìn)行，這樣可大大提高計(jì)算速度，nk可以用以下表示式來(lái)替換：則(3.3.14)(3.3.15)

(3.3.16)又由于式中轉(zhuǎn)角意味著頻率?？梢?jiàn)，系統(tǒng)的單位脈沖響應(yīng)的相角隨時(shí)間呈線性增加，與線性調(diào)頻信號(hào)相似，因此稱(chēng)為Chirp-Z變換。

3.4離散傅立葉變換的實(shí)際應(yīng)用問(wèn)題

3.4.1頻譜泄漏(leakage)

實(shí)際上，頻譜泄漏是由于時(shí)域截?cái)嘣斐傻模魏螁我活l率的周期信號(hào)被截取成有限長(zhǎng)而成為非周期信號(hào)，它就無(wú)法保持單一頻率成分，頻譜線必然會(huì)向周?chē)孤┒蔀轭l帶，其微小幅度(能量)能延伸擴(kuò)散到很高頻區(qū)，如圖3.4.1(d)所示。圖3.4.1信號(hào)截?cái)嘁鸬念l譜泄漏對(duì)于序列x(n)的截?cái)嗖僮?，我們常用N點(diǎn)的矩形窗RN(n)去乘x(n)，截?cái)嗪蟮男盘?hào)頻譜是原有信號(hào)的頻譜X(ejw)與矩形窗譜WR(ejω)的卷積結(jié)果。這個(gè)WR(ejw)的特性顯然對(duì)信號(hào)頻譜的改變有著重大影響。圖3.4.2是N=8矩形窗的時(shí)域和頻域圖。RN(n)的DTFT為(3.4.1)當(dāng)|ω|→π時(shí)，到達(dá)折疊頻率，此時(shí)若N為偶數(shù)，則WR(ejπ)=

0；若N為奇數(shù)，則|WR(ejπ)|=1；而當(dāng)ω=2πi/N時(shí)(i為不等于零的整數(shù))，則WR(ejw)=0，即是幅度過(guò)零點(diǎn)，而直流點(diǎn)ω=0處，WR(ejω)=WR(1)=N。圖3.4.2N=8矩形窗的時(shí)域和頻域圖3.4.2分辨率及補(bǔ)零方法

頻譜分析的分辨率包括兩個(gè)方面的含義：頻率分辨率是指在分析頻帶范圍內(nèi)，頻率樣點(diǎn)間的間隔，即頻點(diǎn)密集度；幅度分辨力也稱(chēng)為分析精度，是指對(duì)于真正的頻譜幅度的逼近能力，幅度分辨力越高，越接近真實(shí)頻譜值，精度越好，誤差越小。前者由fs和N共同決定，后者由數(shù)據(jù)信息量大小決定，通常是原始信號(hào)的采集長(zhǎng)度。做DFT時(shí)，對(duì)序列填補(bǔ)零值可以改變其對(duì)DTFT的采樣密度。因此人們常常有一種誤解，認(rèn)為通過(guò)補(bǔ)零可以提高DFT的頻譜分辨率。我們知道，在數(shù)據(jù)x(n)后面填補(bǔ)若干個(gè)0，根據(jù)DTFT的定義，是不會(huì)對(duì)其DTFT結(jié)果X(ejω)帶來(lái)任何影響的。因此，盡管DFT的點(diǎn)數(shù)增大了，也只是改變了DFT計(jì)算頻譜點(diǎn)的位置和密度，是對(duì)同一個(gè)DTFT采用了不同的頻率采樣密度而已，頻譜分析的精度(幅度和相位的準(zhǔn)確程度)沒(méi)有絲毫改變，它完全取決于有效數(shù)據(jù)的長(zhǎng)度。不同長(zhǎng)度的實(shí)際x(n)其DTFT的結(jié)果X(ejω)是不一樣的，有效序列越長(zhǎng)，分析精度就越高。為此，我們特地把補(bǔ)0后得到的DFT頻譜稱(chēng)為高密度頻譜(Highdensityspectrum)，而將增加有效數(shù)據(jù)長(zhǎng)度獲得的DFT頻譜稱(chēng)為高精度頻譜(Highresolutionspectrum)。不過(guò)要注意的是，應(yīng)該在數(shù)據(jù)的尾部補(bǔ)零，按照DTFT式子，如果從數(shù)據(jù)的頭部補(bǔ)零，相當(dāng)于把有用的數(shù)據(jù)推后，從而帶來(lái)附加的滯后相位。也不能在數(shù)據(jù)之間插0，那會(huì)完全改變?cè)蛄?，屬于插值處理，是另外一種信號(hào)處理方法。

DFT還有一種現(xiàn)象，就是所謂的柵欄效應(yīng)，凡是用離散的值來(lái)近似連續(xù)的量都會(huì)有這種現(xiàn)象。我們知道DTFT的頻譜是連續(xù)且周期的，我們?cè)谄?π主周期內(nèi)將頻率均勻采樣成為DFT，如果頻率采樣點(diǎn)不夠密，就很有可能在兩個(gè)頻率點(diǎn)之間漏掉一個(gè)又窄又高的頻譜成分，我們沒(méi)有觀察到，這就好比從柵欄后面看景物似的，有些物體被欄板寬度所遮而看不到。幸好這個(gè)問(wèn)題不嚴(yán)重也容易解決，增加點(diǎn)數(shù)使得頻點(diǎn)間距變小，或者在序列后面添加少量0，使得所有頻率點(diǎn)都稍微移動(dòng)位置，這樣就能讓原本被遺漏的譜峰落在分析頻率點(diǎn)上，獲得計(jì)算。由于當(dāng)前硬件技術(shù)的發(fā)展，采樣點(diǎn)數(shù)可以做得很大，故柵欄效應(yīng)問(wèn)題已經(jīng)不突出了。

【例3.4.1】

對(duì)于有限長(zhǎng)指數(shù)信號(hào)x(n)=0.9n[U(n)－U(n－N)]和y(n)=0.9n[U(n)－U(n－2N)]，請(qǐng)分析其頻譜X(k)和Y(k)。

解我們對(duì)x(n)后面補(bǔ)N個(gè)0，使得其跟y(n)一樣長(zhǎng)；同時(shí)還可以調(diào)整N，以提高精度。分析如下：

x(t)=e－at≡0.9t，t=nT，T=1

MATLAB程序：

N=16；n=0:N-1；m=0:2*N-1;

x=0.9.^n;y=0.9.^m;

x0=[x,zeros(1,N)]；%序列后補(bǔ)16個(gè)0；

figure(1);

subplot(2,1,1);stem(m,x0);xlabel(′n′);ylabel(′x0(n)′);

subplot(2,1,2);stem(m,y);xlabel(′n′);ylabel(′y(n)′);

X0=fft(x0);Y=fft(y)；%求出高密度頻譜X0，高精度頻譜Y

figure(2);

subplot(2,1,1);stem(m,abs(X0));xlabel(′k′);ylabel(′|X0(k)|′);

subplot(2,1,2);stem(m,abs(Y));xlabel(′k′);ylabel(′|Y(k)|′);圖3.4.3實(shí)指數(shù)序列的高密度頻譜和高精度頻譜圖3.4.4是沒(méi)有補(bǔ)零的x(n)及其DFT頻譜X(k)，只有16點(diǎn)，它恰恰就是X0(k)的偶數(shù)序號(hào)的內(nèi)容，它們二者的頻譜包絡(luò)線是一樣的，說(shuō)明添0后的X0(k)對(duì)真正的頻譜的逼近程度沒(méi)有提高。這是顯然的，因?yàn)樘?并未給序列增加額外有效信息。因此，實(shí)際應(yīng)用中，在允許的情況下，可充分利用硬件資源，以足夠高的采樣率采集足夠長(zhǎng)的信息，盡量提高分辨率。圖3.4.4實(shí)指數(shù)序列的16點(diǎn)采樣值及其頻譜3.4.3DFT的處理增益

如果一個(gè)1V的直流電壓信號(hào)被采樣，然后運(yùn)用N點(diǎn)的DFT計(jì)算，將會(huì)在k=0處得到N×1V的數(shù)值，這就是DFT的處理增益，也叫運(yùn)算增益。對(duì)其他頻率點(diǎn)±k處的情況也是一樣，只不過(guò)每個(gè)k上的增益量為0.5N，考慮到對(duì)稱(chēng)的負(fù)頻率部分的0.5N，可以認(rèn)為DFT對(duì)所有各個(gè)信號(hào)分量都放大N倍，這從能量的角度也很容易理解。我們可以利用N點(diǎn)DFT的內(nèi)在相關(guān)增益來(lái)檢測(cè)隱藏于噪聲中的信號(hào)能量，從而把信號(hào)從噪聲背景中提取出來(lái)。一個(gè)N點(diǎn)DFT對(duì)信號(hào)的加工過(guò)程，還可以看成是N個(gè)窄帶的帶通濾波器并聯(lián)工作，它輸出N個(gè)頻率分量。每個(gè)窄帶濾波器就是DFT的一個(gè)輸出頻率單元，中心頻率就在分析頻點(diǎn)處，即kfs/N。隨著N的增大，濾波器的增益變大，而帶寬變窄。這在能量檢測(cè)方面十分有效，因?yàn)榻档蛶挸丝梢詾V除通帶內(nèi)的背景噪聲外，還可以提高頻率分辨率。我們用一個(gè)例子來(lái)說(shuō)明。

【例3.4.2】

在一個(gè)疊加有隨機(jī)噪聲的單一頻率的正弦波信號(hào)中，提取正弦頻譜。隨機(jī)噪聲服從0均值和方差1的N(0,1)正態(tài)分布(normaldistribution)，即高斯分布(Gaussiandistribution)。

N=64；n=0:N-1；%數(shù)據(jù)長(zhǎng)度64點(diǎn)

Noise=random(′N(xiāo)ormal′,0,1,1,N);

%(-2~2)之間的N(0,1)隨機(jī)

數(shù)生成

x=sin((2*pi*20/64).*n)；%歸一化fs=1，分64頻率點(diǎn)。第20點(diǎn)的

信號(hào)頻率為20/64=0.3125

xN=0.5*x+2*Noise；%將幅度為0.5的信號(hào)淹沒(méi)在幅度為其4

倍的噪聲中

X=fft(x);XN=fft(xN);

Xa=abs(X);PX=Xa.^2；%原始信號(hào)的頻譜幅度和功率(幅度平方)

XNa=abs(XN);PXN=XNa.^2；%加有噪聲信號(hào)的頻譜幅度和功率(幅度平方)

figure(1);

subplot(2,1,1);plot(n,20*log10(PX));xlabel(′k′);ylabel(′Px(k)′);

subplot(2,1,2);

plot(n,20*log10(PXN));

xlabel(′k′);

ylabel(′Pxs(k)′);

程序運(yùn)行后的結(jié)果如圖3.4.5所示，(a)圖為單一頻率信號(hào)的功率譜Px(k)，(b)圖是攜帶噪聲的信號(hào)的功率譜Pxs(k),后者已看不出來(lái)原來(lái)信號(hào)究竟在哪里了。我們?cè)黾佑^察時(shí)間,把采樣數(shù)據(jù)增加到256，情況就會(huì)好一些，如圖3.4.6所示。圖3.4.564點(diǎn)DFT計(jì)算的信號(hào)功率頻譜圖3.4.6256點(diǎn)DFT的信號(hào)功率頻譜繼續(xù)增加數(shù)據(jù)長(zhǎng)度到N=1024，DFT的處理增益就顯示威力了，提高了SNR，如圖3.4.7所示，正弦頻譜的位置在Pxs(k)中明顯從噪聲里展露出來(lái)。

我們有一個(gè)結(jié)論，DFT頻率點(diǎn)輸出噪聲的標(biāo)準(zhǔn)差RMS(均方根)與成比例，而包含信號(hào)諧波的頻率點(diǎn)的DFT輸出幅度與N成正比，二者增長(zhǎng)率不同，因此，輸入樣點(diǎn)N越大，DFT計(jì)算出來(lái)的信號(hào)幅度和噪聲幅度之比(即SNR)就越大。每增加一倍的點(diǎn)數(shù)，信噪比SNR增大3dB左右，即SNR2N=SNRN+20lg。圖3.4.71024點(diǎn)DFT的信號(hào)功率頻譜

3.5快速傅立葉變換FFT典型用法

3.5.1IDFT的快速算法

FFT可以用來(lái)計(jì)算IDFT，我們知道DFT的正反變換只有兩點(diǎn)不同：一個(gè)是變換指數(shù)上的正負(fù)號(hào)，另一個(gè)就是比例因子N。因?yàn)?/p>

這說(shuō)明把頻譜X(k)先求共軛，然后做FFT得到中間結(jié)果，將它求共軛后除以N就是x(n)。3.5.2實(shí)數(shù)序列的FFT

任何實(shí)數(shù)都可看成虛部為零的復(fù)數(shù)。例如：求某實(shí)信號(hào)x(n)的復(fù)譜，可認(rèn)為是將實(shí)信號(hào)加上數(shù)值為零的虛部變成復(fù)信號(hào)(x(n)+j0)，再用復(fù)數(shù)形式FFT求其離散傅立葉變換。這種做法很不經(jīng)濟(jì)，因?yàn)榘褜?shí)序列變成復(fù)序列，存儲(chǔ)器要增加一倍，且計(jì)算機(jī)運(yùn)行時(shí)，即使虛部為零，也要進(jìn)行涉及虛部的運(yùn)算，浪費(fèi)了運(yùn)算量。

合理的解決方法是利用復(fù)數(shù)形式FFT對(duì)實(shí)數(shù)據(jù)進(jìn)行有效計(jì)算，下面介紹兩種方法。

(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)可通過(guò)一次FFT運(yùn)算同時(shí)獲得。首先將x(n)、y(n)分別當(dāng)作一復(fù)序列g(shù)(n)的實(shí)部及虛部，即令g(n)=x(n)+jy(n)，經(jīng)過(guò)FFT運(yùn)算可獲得g(n)的DFT值:

G(k)=X(k)+jY(k)

(3.5.2)

利用離散傅立葉變換的共軛對(duì)稱(chēng)性：(3.5.3)通過(guò)g(n)的FFT運(yùn)算結(jié)果G(k)，由上式可得到X(k)的值。同理，通過(guò)G(k)，由上式也可得到Y(jié)(k)的值，(3.5.4)

(2)用一個(gè)N點(diǎn)的FFT運(yùn)算獲得一個(gè)長(zhǎng)度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－1

x2(n)=x(2n+1)n=0，1，…，N－1

然后將x1(n)及x2(n)組成一個(gè)復(fù)序列：

y(n)=x1(n)+jx2(n)

通過(guò)N點(diǎn)FFT運(yùn)算可得到：

Y(k)=X1(k)+jX2(k)

(3.5.5)

根據(jù)前面的討論，可得到

(3.5.6)現(xiàn)在，為求2N

點(diǎn)x(n)所對(duì)應(yīng)

X(k)，需求出

X(k)與X1(k)、X2(k)的關(guān)系。由定義有：(3.5.7)令代入式(3.5.7)可得：(3.5.9)(3.5.10)(3.5.8)3.5.3線性卷積的FFT計(jì)算

我們已經(jīng)知道，DFT的卷積定理是對(duì)應(yīng)于循環(huán)卷積的，雖然與線性卷積不同，但可以通過(guò)在序列后面補(bǔ)若干個(gè)0，就能使用循環(huán)卷積來(lái)計(jì)算線性卷積。也就是把FFT用在線性卷積上，從而發(fā)揮出它的快速優(yōu)勢(shì)，避免因直接計(jì)算線性卷積而付出大量的運(yùn)算時(shí)間。

設(shè)N點(diǎn)x(n)與M點(diǎn)h(n)，做線性卷積時(shí)得到輸出y(n)，共L=N+M－1點(diǎn)。

用L點(diǎn)的FFT完成卷積的過(guò)程是：

(1)x(n)補(bǔ)M－1個(gè)0到L點(diǎn)后，用FFT求出X(k)。

(2)h(n)補(bǔ)N－1個(gè)0到L點(diǎn)后，用FFT求出H(k)。

(3)將X(k)和H(k)對(duì)應(yīng)相乘得到Y(jié)(k)，并求其L點(diǎn)的IFFT，即可獲得線性卷積結(jié)果y(n)。

1.重疊相加法

假設(shè)無(wú)限長(zhǎng)序列x(n)，我們將其切成小段，每段N點(diǎn)。有如下表示：(3.5.11)那么離散線性系統(tǒng)h(n)的輸出就是(3.5.12)如圖3.5.1和圖3.5.2是長(zhǎng)的輸入信號(hào)x(n)切成N＝10點(diǎn)的片段和短的脈沖響應(yīng)h(n)，M＝7的情況。圖3.5.3是各小段輸出結(jié)果的重疊相加過(guò)程。圖3.5.1長(zhǎng)的輸入信號(hào)x(n)和短的脈沖響應(yīng)h(n)圖3.5.2輸入信號(hào)x(n)切成各小段xi(n)長(zhǎng)度N=10圖3.5.3輸出各小段yi(n)長(zhǎng)度N+M－1=16，頭尾疊加M－1=6點(diǎn)的情況

2.重復(fù)保留法

重復(fù)保留法與上面這個(gè)方法稍有不同，在分段序列xi(n)的后面不用補(bǔ)零，而是將前小段用過(guò)的數(shù)據(jù)尾部保留M－1點(diǎn)下來(lái)，再添上本小段新的數(shù)據(jù)N點(diǎn)，形成N+M－1點(diǎn)后參與卷積運(yùn)算，接下來(lái)，本段的尾部也要留下M－1給下一段用。因此保留的數(shù)據(jù)被用了2遍，在片段結(jié)果yi(n)中就應(yīng)該將其拋棄，它位于每段結(jié)果yi(n)的前部M－1點(diǎn)。這樣將剔除后的yi(n)與前面的yi－1(n)尾部直接連接即可。這個(gè)方法的名稱(chēng)指的是對(duì)輸入數(shù)據(jù)片段的形成特點(diǎn)。它跟前者不一樣，少了數(shù)據(jù)相加的環(huán)節(jié)，比較省事。要注意的是，在開(kāi)始的第一段，因?yàn)槭亲钋懊?，并沒(méi)有數(shù)據(jù)留下來(lái)，只好用M－1個(gè)0來(lái)填充，這個(gè)0是補(bǔ)在序列前面的！輸入序列小段如圖3.5.4所示。輸出構(gòu)成如圖3.5.5所示。圖3.5.4輸入各小段xi(n)