第3章離散傅里葉變換3.1信號頻域分析的幾種形式3.2離散傅里葉變換的定義3.3DFT與DTFT、Z變換的關(guān)系3.4離散傅里葉變換的性質(zhì)3.5用DFT計(jì)算數(shù)字頻譜的誤差及解決方法3.6用MATLAB實(shí)現(xiàn)離散傅里葉變換習(xí)題

3.1信號頻域分析的幾種形式

1.連續(xù)時(shí)間非周期信號的傅里葉變換(FT)

連續(xù)時(shí)間信號x(t)的傅里葉變換X(Ω)定義為反變換由下式確定

如圖3-1所示為連續(xù)信號x(t)及其傅里葉變換X(Ω)。由該圖可以看出,連續(xù)信號的傅里葉變換也是連續(xù)的,這不便于計(jì)算機(jī)處理。圖3-1連續(xù)信號及其傅里葉變換

2.連續(xù)時(shí)間周期信號的傅里葉級數(shù)表示(FS)

連續(xù)時(shí)間周期信號x(t)的傅里葉級數(shù)表示把x(t)展開為復(fù)指數(shù)信號ejkΩ0t的加權(quán)和,即

設(shè)周期信號x(t)的周期為Ｔ,稱Ω0=2π/T為基波頻率,簡稱為基頻；對k>1的一般情況,稱kΩ0為k次諧波頻率。k次諧波的系數(shù)X(kΩ0)由下式確定

圖3-2所示為連續(xù)周期信號x(t)及其傅里葉級數(shù)表示X(kΩ0)。由該圖可以看出,連續(xù)時(shí)間周期信號的傅里葉級數(shù)把頻域變成離散的了,但時(shí)域仍然是連續(xù)的,這也不便于計(jì)算機(jī)處理。圖3-2連續(xù)時(shí)間周期信號及其傅里葉級數(shù)

3.序列的離散時(shí)間傅里葉變換(DTFT)

對連續(xù)時(shí)間信號進(jìn)行采樣,得到在時(shí)間上離散的信號即序列。非周期序列x(n)可以通過離散時(shí)間傅里葉變換得到其頻譜,離散時(shí)間傅里葉變換X(ejω)定義如下：

離散時(shí)間傅里葉反變換定義為

非周期序列x(n)的傅里葉變換X(ejω)是ω的連續(xù)函數(shù)且隱含著周期性,X(ejω)的周期為2π,這個(gè)周期性質(zhì)表明為了分析目的僅僅需要X(ejω)的一個(gè)周期即可。如圖3-3所示

為離散序列及其離散時(shí)間傅里葉變換,從圖中可以看出頻域的周期性是由時(shí)域的離散化所致。但是從圖中可以看出頻域還是連續(xù)的,同樣不便于計(jì)算機(jī)處理。

圖3-3離散序列及其離散時(shí)間傅里葉變換

圖3-4離散周期序列及其離散傅里葉級數(shù)

至此,我們已經(jīng)介紹了信號頻域分析的四處形式,總結(jié)如下：

?連續(xù)時(shí)間非周期信號的傅里葉變換其時(shí)域和頻域都是連續(xù)的；

?連續(xù)時(shí)間周期信號的傅里葉級數(shù)表示其時(shí)域是連續(xù)的,但頻域是離散的；

?非周期序列的離散時(shí)間傅里葉變換其時(shí)域是離散的,但頻域是連續(xù)的；

?周期序列的離散時(shí)間傅里葉級數(shù)表示其時(shí)域和頻域都是離散的。

信號在時(shí)域和頻域都具有連續(xù)與離散、周期與非周期兩種特征,其對應(yīng)規(guī)律如下：如果信號在頻域是離散的,則該信號在時(shí)域必然表現(xiàn)為周期性；反之,如果信號在時(shí)域是離散的,則該信號在頻域必然是周期的。所以,周期的離散信號(序列)經(jīng)過DFS表示,其頻域一定既是周期又是離散的。存在這種對應(yīng)規(guī)律的原因是：一個(gè)域的離散化導(dǎo)致另一個(gè)域的周期化；一個(gè)域的周期化導(dǎo)致另一個(gè)域的離散化。

3.2離散傅里葉變換的定義3.2.1DFT與DFS的關(guān)系DFS將信號的時(shí)域和頻域都離散化了,但兩者都是周期性的。實(shí)際中遇到的序列往往都是有限長的序列,一般并不滿足周期性。但是我們知道,一個(gè)周期序列和其周期的離散頻譜所攜帶的信息都體現(xiàn)在一個(gè)周期,即周期序列和有限長序列在攜帶信息上并無本質(zhì)區(qū)別。這樣通過將有限長序列人為地延拓成為周期序列,然后對其離散傅里葉級數(shù)取它的主值區(qū)間,就可以得到有限長序列的離散頻譜。所以,離散傅里葉變換為有限長序列提供了一個(gè)有效的頻率分析方法。

圖3-5

DFS與DFT的關(guān)系示意圖

3.2.2

DFT的定義

對于有限長序列x(n),0≤n≤N-1,其N點(diǎn)離散傅里葉變換和離散傅里葉反變換的定義分別為

DFT的定義可以寫成矩陣形式

X=Wx

其中,x為輸入序列組成的向量,即x=[

x(0),x1),…,x(N-1)]T；X為DFT序列組成的向量,即X=[X(0),X(1),…,X(N-1)]T；W為DFT矩陣,即

3.2.3

DFT的物理意義

式(3-13)看似比較復(fù)雜,它是數(shù)字信號處理中數(shù)字譜分析最有用的工具。為了理解式(3-13),通過歐拉公式e-jθ

=cosθ-jsinθ,可以將其轉(zhuǎn)化為

把式(3-13)中復(fù)雜的指數(shù)轉(zhuǎn)化為實(shí)數(shù)部分與虛數(shù)部分,其中X(k)為DFT第k個(gè)輸出部分,k為DFT在頻域的輸出序號,k=0,1,2,3,…,N-1。x(n)則為輸入樣本序列,n為時(shí)域輸入序號,n=0,1,2,3,…,N-1。數(shù)值N是一個(gè)很重要的參數(shù),因?yàn)樗鼪Q定了所需要輸入樣本的大小頻域輸出結(jié)果的分辨率和計(jì)算N點(diǎn)DFT函數(shù)所要的時(shí)間。

由此可見,DFT的每一個(gè)輸出值X(k)都由輸入序列值與不同頻率的正弦和余弦乘積之和來確定。不同的頻率取決于原始信號的采集頻率fs和DFT的樣本數(shù)量N。比如,給定一個(gè)速率為500每秒的樣本連續(xù)信號,然后對其執(zhí)行16點(diǎn)DFT數(shù)據(jù)抽樣,則正弦函數(shù)的頻率為fs/N=500/16或者31.25Hz。其他X(k)的分析頻率則是基頻的整數(shù)倍,例如：

圖3-6模擬頻率與數(shù)字頻率之間的定標(biāo)關(guān)系

例3-1已知序列x(n)=δ(n),求它的N點(diǎn)DFT。

δ(n)的X(k)如圖3７所示。這是一個(gè)很特殊的例子,它表明對序列δ(n)來說,不論對它進(jìn)行多少點(diǎn)的DFT,所得結(jié)果都是一個(gè)離散矩形序列,這是因?yàn)棣?n)的頻譜是均勻譜。

圖3-7序列δ(ｎ)及其離散傅里葉變換

圖3-8有限長序列及其DFT

3.3

DFT與DTFT、Z變換的關(guān)系

前面已討論了DFT與DFS的關(guān)系,現(xiàn)在討論DFT與DTFT、Z變換的關(guān)系。

1.DFT與Z變換的關(guān)系若x(n)是一個(gè)長度為N的有限長序列,對x(n)進(jìn)行Z變換,有

圖3-9DFT與DTFT、Z變換的關(guān)系

2.DFT與DTFT變換的關(guān)系

由于DTFT,即序列的傅里葉變換X(ejω)是單位圓上的Z變換,根據(jù)式(3-26),DFT與DTFT的關(guān)系為

上式說明X(k)也可以視為序列x(n)的傅里葉變換X(ejω)在區(qū)間[0,2π]上的N點(diǎn)等間隔采樣,其采樣間隔為ωN=2π/N,圖3-9(b)給出了X(ejω)和X(k)的關(guān)系。DFT的變換區(qū)間長度N不同,對X(ejω)在區(qū)間[0,2π]上的采樣間隔和采樣點(diǎn)數(shù)不同,所以DFT的變換結(jié)果也不同。

該序列的N=8點(diǎn)DFT為

圖3-10給出了矩形序列R6(n)的DTFT以及N=8、N=16、N=32、N=64不同點(diǎn)數(shù)的DFT的幅度值。從中可以看出,DFT就是對序列頻譜的等間隔采樣。

圖3-10矩形序列R6(n)的DTFT以及不同點(diǎn)數(shù)的DFT

3.4離散傅里葉變換的性質(zhì)

1．線性特性若w(n)=ax(n)+by(n),則

DFT[w(n)]=W(k)=aX(k)+bY(k)(3-36)

如果x(n)和y(n)的序列長度N1、N2不同,則選擇N=ｍax［N1,N2］作為變換長度,而短序列通過補(bǔ)0達(dá)到N點(diǎn)。

2.DFT隱含的周期性

離散傅里葉變換定義了時(shí)域N個(gè)離散點(diǎn)到頻域N個(gè)頻率點(diǎn)的一一映射,但是定義式本身隱含著周期性,即

x(n)和X(k)都是周期為N的序列。

離散傅里葉變換定義中的正變換和反變換分別限定0≤n≤N-1和0≤k≤N-1。用計(jì)算機(jī)或數(shù)字處理器對信號進(jìn)行頻譜分析時(shí),一般都期望信號在時(shí)域和頻域都是離散和有限的,因而這種限定是自然而然和合情合理的。離散傅里葉變換潛在的周期性直接導(dǎo)出了它具有許多優(yōu)異的特性,利用這些特性就得到了離散傅里葉變換的快速算法,它使得計(jì)算DFT高效而方便。我們可以把x(n)和Ｘ(k)看成是其中的任意一個(gè)周期,只要進(jìn)行DFT計(jì)算時(shí)取的點(diǎn)數(shù)N不小于對信號的采樣點(diǎn)數(shù)M,DFT和IDFT都是精確的。

3.循環(huán)(圓周)移位的概念

在定義序列x(n)的N點(diǎn)離散傅里葉變換時(shí),把序號n約束為0≤n≤N-1。期望移位后新序列的序號n也滿足前述約束,顯然使用模運(yùn)算符定義循環(huán)移位即可滿足這個(gè)約束條件。定義循環(huán)移位運(yùn)算如下：

圖3-11給出了循環(huán)移位運(yùn)算的實(shí)例。其中圖3-11(a)為原序列x(n),圖3-11(b)為x1(n)=x(〈n-2〉6),圖3-11(c)為x2(n)=x(〈n+2〉6)。圖中的n表示序列的序號,序號按順時(shí)針依次排列。x(〈n-2〉6)使得x(n)的所有取值依順時(shí)針移位2個(gè)點(diǎn)；而x(〈n+2〉6)使得x(n)的所有取值依逆時(shí)針移位2個(gè)點(diǎn)。這樣移動(dòng)的結(jié)果是所有取值的相對位置不變,而對應(yīng)的新序列的序號保持不動(dòng),把序號與取值對應(yīng)起來就構(gòu)成了新序列,如圖3-11(b)中,新的序列為x1(0)=x(4),x1(1)=x(5),x1(2)=x(0),x1(3)=x(1),x1(4)=x(2),x1(5)=x(3)。

圖3-11循環(huán)(圓周)移位運(yùn)算示意圖

循環(huán)移位可以想象將序列x(n)順序排列在N等分的圓周上,x(n-m)是將序列各元素在圓周上按順時(shí)針旋轉(zhuǎn)m個(gè)位置,而x(n+m)是將序列各元素在圓周上按反時(shí)針方向旋轉(zhuǎn)m個(gè)位置,所以,循環(huán)移位又叫圓周移位。

為簡便起見,記

同理

4.時(shí)域循環(huán)移位性質(zhì)

若x(n)?X(k),則

或者簡記為

序列在時(shí)域循環(huán)移位m個(gè)單位,其DFT的相位頻譜各次諧波相位平移,但幅度頻譜不變。

5.頻域循環(huán)移位性質(zhì)

若x(n)?X(k),則

或者簡記為

序列在時(shí)域被調(diào)制相當(dāng)于頻譜的平移,包括幅度頻譜和相位頻譜。

6.循環(huán)卷積定理

1)兩個(gè)有限長序列的循環(huán)卷積

2)時(shí)域循環(huán)卷積定理

7.共軛對稱性

設(shè)x?(n)為x(n)實(shí)共軛復(fù)序列,若x(n)?X(k),則有

同理

在討論DTFT的對稱性時(shí),也提到過共軛對稱序列和共軛范對稱序列,那里的對稱性是指關(guān)于縱坐標(biāo)的對稱性。DFT的對稱性有所不同,序列x(n)和其離散傅里葉變換X(k)都是有限長序列,其定義的區(qū)間為0≤n≤N-1,所以,DFT的對稱性是指關(guān)于N/2的對稱性。

8.帕薩瓦爾定理

設(shè)x(n)和y(n)的N點(diǎn)離散傅里葉變換分別為X(k)和Y(k),則有如下帕薩瓦爾定理

特別地,當(dāng)x(n)=y(n)時(shí),帕薩瓦爾定理變?yōu)?/p>

式(3-57)表明,序列在時(shí)域和頻域的能量是不變的,故也稱為能量守恒定理。

3.5用DFT計(jì)算數(shù)字頻譜的誤差及解決方法

首先通過一個(gè)例題來了解利用DFT計(jì)算數(shù)字頻譜的好處。

圖3-12混合頻率信號的時(shí)域、頻域分析

在實(shí)際應(yīng)用中,原始信號可能是無限長的連續(xù)時(shí)間信號,為了利用離散傅里葉變換計(jì)算這樣信號的離散頻譜,還必須進(jìn)行一些必要的處理,而這些處理過程可能對原始信號的真實(shí)頻譜產(chǎn)生誤差,下面討論誤差產(chǎn)生的原因和解決的辦法。

如果原始信號是連續(xù)的模擬信號xa(t),第一步是對其進(jìn)行等間隔采樣得到采樣信號,則由時(shí)域采樣定理知采樣信號的頻譜是xa(t)頻譜的周期延拓。為了避免頻譜混疊,要求采樣頻率不小于信號最高頻率的兩倍。一般來說,原始信號是有限長的,所以其頻譜寬度是無限的,因此,在對有限長信號進(jìn)行采樣前先進(jìn)行抗混疊濾波,一方面使得采樣頻率無需非常高,另外一方面避免了頻譜混疊。通過對連續(xù)信號xa(t)進(jìn)行采樣就得到了離散序列x(n)=xa(nT)。

第二步是對x(n)進(jìn)行截短處理,這是因?yàn)閤(n)太長則不利于處理。截短處理在時(shí)域進(jìn)行,它通過把x(n)乘以一個(gè)長度為N的窗函數(shù)w(n)得到有限長序列x(n)w(n)。x(n)w(n)的離散時(shí)間傅里葉變換(DTFT)是x(n)的離散時(shí)間傅里葉變換X(ω)與w(n)的離散時(shí)間傅里葉變換W(n)之卷積。因?yàn)橛邢揲L序列w(n)的頻譜寬度是無限的,所以乘積x(n)w(n)的頻譜也是無限的,即頻譜“擴(kuò)散”(拖尾,展寬)了,或者說頻譜“泄露”了。為了減小頻譜泄露,可以采取兩個(gè)方法：其一,增加窗函數(shù)的寬度,即取更多的數(shù)據(jù)；其二,采用性能更優(yōu)的窗函數(shù),這些窗函數(shù)具有比矩形窗旁瓣小主瓣窄的優(yōu)點(diǎn),這必將減小頻譜泄露,在FIR數(shù)字濾波器設(shè)計(jì)一章中會(huì)詳細(xì)講解這個(gè)問題。

作為實(shí)例,下面考慮用矩形窗函數(shù)RL(n)對離散余弦序列cos(ω0n)進(jìn)行截短處理后的頻譜,截短得到的序列為cos(ω0n)RL(n)。圖3-13(a)給出了cos(nπ/3)的頻譜,圖3-13(b)為cos(nπ/3)截短處理后的頻譜,由圖可以清楚地看出頻譜發(fā)生了泄漏。

圖3-13

cos(nπ/3)的DTFT及截短處理后的DTFT

由于x(n)w(n)的頻譜依然是連續(xù)的,所以要進(jìn)行第三步處理：對x(n)w(n)的頻譜進(jìn)行離散化,即等間隔采樣得到M個(gè)離散點(diǎn)。與時(shí)域的等間隔采樣導(dǎo)致頻域的周期延拓相對偶,這種頻域的等間隔采樣導(dǎo)致了時(shí)域的周期延拓。在講解離散傅里葉變換的物理意義時(shí),已經(jīng)指出只要一個(gè)周期內(nèi)的采樣點(diǎn)數(shù)M≥N,則時(shí)域信號不會(huì)發(fā)生混疊。

以下介紹兩個(gè)有關(guān)分辨率的概念?！邦l率分辨率”是指所用的算法將所分析的時(shí)域信號頻譜中兩個(gè)靠得很近的譜峰分辨開來的能力,“頻率分辨率”通常也稱為“物理頻率分辨

率”,以便與“計(jì)算頻率分辨率”相區(qū)別?！坝?jì)算頻率分辨率”是指所分析的信號的離散頻譜中相鄰點(diǎn)間的頻率間隔。

將x(t)用間隔Ts=采樣得到x(n),采樣頻率為fs=1/Ts,則能得到總的采樣點(diǎn)數(shù)為M=Tu/Ts=Tufs。我們知道在對x(n)進(jìn)行N點(diǎn)離散傅里葉變換時(shí),只要N≥M。這里的計(jì)算頻率分辨率為

物理頻率分辨率為

將M=Tufs代入式(3-60)得到的物理頻率分辨率與式(3-58)給出的一致,這說明不能光靠增加采樣點(diǎn)數(shù)M來提高物理頻率分辨率,這是因?yàn)樵谛盘柕某掷m(xù)時(shí)間Tu保持不變的前提下,采樣點(diǎn)數(shù)M增加,則采樣間隔Ts=減小,但M與Ts=兩者的乘積保持為Tu不變。

在進(jìn)行離散傅里葉變換時(shí),通過在有效數(shù)據(jù)后面補(bǔ)充一些零來達(dá)到改善頻譜的目的,但這并不能提高算法的物理頻率分辨率,這是因?yàn)槲锢眍l率分辨率由有效的數(shù)據(jù)點(diǎn)數(shù)M決定。從根本上說,補(bǔ)零并沒有提供新的信息,不能提高分辨率就理所當(dāng)然了,但這樣做能提高算法的計(jì)算頻率分辨率。此外補(bǔ)零可以使得總的點(diǎn)數(shù)N是2的整數(shù)次冪,便于用下一章要講解的快速傅里葉變換來高效地計(jì)算離散傅里葉變換。

我們知道x(n)的N點(diǎn)離散傅里葉變換X(k)與x(n)的離散時(shí)間傅里葉變換在頻率點(diǎn)2πk/N的取值相等。通過這N個(gè)離散的頻率點(diǎn)處的離散時(shí)間傅里葉變換來觀察序列的頻譜,就像通過一個(gè)“柵欄”觀賞風(fēng)景一樣,只能在這些離散的頻率點(diǎn)處看到真實(shí)的景象,稱這種現(xiàn)象為“柵欄效應(yīng)”。顯然如果N足夠大(如果有效數(shù)據(jù)長度M不變,可以通過補(bǔ)零增大N),這些頻率點(diǎn)分布得足夠稠密,則可減小柵欄效應(yīng)。補(bǔ)零可以把原來由于頻率點(diǎn)錯(cuò)位而攔住的有效頻率成分顯現(xiàn),但并不能提高算法的物理頻率分辨率,或者說原來分不開的兩個(gè)頻峰,補(bǔ)零后依然不能分開。

圖3-14為序列補(bǔ)零對DTFT與DFT的影響。圖3-14(a)為在有效長度為Ｍ的序列尾部補(bǔ)充了N1-M個(gè)零點(diǎn),圖3-14(b)為對應(yīng)的DTFT。由式(3-62)知,補(bǔ)零后的DTFT與沒有補(bǔ)零的DTFT一致,對圖3-14(b)所示的DTFT在一個(gè)周期內(nèi)進(jìn)行間隔為2π／N1等間隔采樣,采樣點(diǎn)對應(yīng)的值就是N1點(diǎn)離散傅里葉變換X1(k),如圖3-14(c)所示。對原序列尾部補(bǔ)充了N2-M(其中N2＞N1)個(gè)零點(diǎn)得到圖3-14(d),圖3-14(e)為對應(yīng)的DTFT,與圖31４(b)一致。對圖3-14(e)所示的DTFT在一個(gè)周期內(nèi)進(jìn)行間隔為2π／N2等間隔采樣,采樣點(diǎn)對應(yīng)的值就是N2點(diǎn)離散傅里葉變換X2(k),如圖3-14(f)所示。

圖3-14序列補(bǔ)零對DTFT與DFT的影響

在計(jì)算DFT時(shí),需要選擇合適的采樣點(diǎn)數(shù)M,下面給出具體的方法。設(shè)X(ω)的截止頻率為fc,這個(gè)可以事先確定,因?yàn)橐阎獂(t)就可以得到Ｘ(ω)及其截止頻率fc。為了避免混疊,選擇采樣頻率fs滿足2.5fc≤fs≤3.0fc。若給定頻率分辨率Δfp,可得需要的最少采樣點(diǎn)數(shù)Ｍ為

3.6用MATLAB實(shí)現(xiàn)離散傅里葉變換

本節(jié)通過MATLAB實(shí)現(xiàn)對連續(xù)信號進(jìn)行DFT的分析。

例3-9設(shè)待分析的信號是三個(gè)正弦波的加權(quán)和,即

圖3-15依次為采樣得到的離散序列的40、100、150和200點(diǎn)DFT。由圖3-15可以看出DFT的點(diǎn)數(shù)比所需的采樣點(diǎn)數(shù)少時(shí),算法的分辨率不足以分辨出所有的頻率分量。

圖3-15例3-9進(jìn)行不同點(diǎn)數(shù)的DFT分析結(jié)果

圖3-15例3-9進(jìn)行不同點(diǎn)數(shù)的DFT分析結(jié)果

例3-10設(shè)x(n)＝R4(n),計(jì)算該序列在變換區(qū)間分別為N=8、N=16、N=32、N=64點(diǎn)的DFT,畫出其幅度頻譜和相位頻譜,并進(jìn)一步討論序列補(bǔ)零對頻譜分辨率的影響。

需要注意的是,當(dāng)幅度采樣是零時(shí),相應(yīng)的相位并不是零,這是由于用MATLAB計(jì)算相位部分所用的特定算法所致。圖3-16所示的DFT圖,圖中也用虛線示出X(ejω)的圖以供比較。從圖3-16中可見，X4(k)正確給出了X(ejω)

的４個(gè)樣本值,但僅有一個(gè)是非零樣本,即一個(gè)常數(shù)(或直流DC)信號,這就是由于DFTX4(k)所預(yù)期的在k=0(或ω=0)時(shí)有一個(gè)非零樣本,而在其他頻率沒有值的緣故。

圖3-16例3-10在N=4點(diǎn)時(shí)的DFT結(jié)果

那么如何得到DTFTX(ejω)

的其他樣本？很明顯應(yīng)該增加采樣的密度,即應(yīng)增大Ｎ。設(shè)想現(xiàn)用兩倍的點(diǎn)數(shù),這可以通過將x(n)補(bǔ)上4個(gè)零值而構(gòu)成一個(gè)8點(diǎn)的序列來完成

這是一種非常重要的運(yùn)算,稱為補(bǔ)零運(yùn)算。在實(shí)際中為了得到信號更密的譜,這一運(yùn)算是必須的。

圖3-17例3-10在N=8點(diǎn)時(shí)的DFT結(jié)果

繼續(xù)下去,若將x(n)補(bǔ)12個(gè)零值而作為一個(gè)16點(diǎn)序列,那么頻率分辨率是ω1=2π/16＝π/8和W16=e-jπ/8,因此得到一個(gè)更為密集的譜,譜樣本相距π/8。X16(k)如圖3-18所示。

圖3-18例3-10在N=16點(diǎn)時(shí)的DFT結(jié)果

通過以上例題,可以得出如下幾點(diǎn)結(jié)論。

(1)補(bǔ)零就是將更多的零值補(bǔ)在原序列的后面,所得到的更長一些的序列對離散時(shí)間傅里葉變換提供了更為密集的樣本。在MATLAB中,補(bǔ)零使用zeros函數(shù)實(shí)現(xiàn)。

(2)補(bǔ)零給出了一種高密度的譜,并對畫圖提供了一種更好的展現(xiàn)形式,但是它并沒有給出高分辨率的譜,因?yàn)闆]有任何新的信息附加到這個(gè)信號上,而僅在數(shù)據(jù)中添加了額外的零值。

(3)為了得到更高分辨率的譜,就必須要從實(shí)驗(yàn)或觀察中獲取更多的數(shù)據(jù),也存在一些高級的方法是利用附加的邊緣信息或非線性技術(shù)來獲得的。

習(xí)題

第4章快速傅里葉變換4.1直接計(jì)算DFT的運(yùn)算量及改進(jìn)的途徑4.2基2時(shí)間抽取的快速傅里葉變換4.3基2頻率抽取的快速傅里葉變換4.4離散傅里葉反變換的快速算法4.5用FFT計(jì)算序列的線性卷積4.6用MATLAB實(shí)現(xiàn)快速傅里葉變換習(xí)題

4.1直接計(jì)算DFT的運(yùn)算量及改進(jìn)的途徑

由x(n)的N點(diǎn)離散傅里葉變換定義式

可以看出，為了完成每一個(gè)值X,k)的計(jì)算，需要N次復(fù)數(shù)乘法x,n)WnkN，n=0，1，…，N-1，然后把這N個(gè)乘積相加，這需要N-1次復(fù)數(shù)加法，因此，計(jì)算全部N點(diǎn)離散傅里葉變換Ｘ,k)總共需要N2次復(fù)數(shù)乘法和N,N-1)次復(fù)數(shù)加法，當(dāng)N較大時(shí)，運(yùn)算量是很大的。例如當(dāng)N=1024時(shí)，需要完成1048756次復(fù)數(shù)乘法，即100多萬次復(fù)數(shù)乘法。可見離散傅里葉變換的計(jì)算量隨N2變化，當(dāng)N較大時(shí)，計(jì)算量相當(dāng)驚人，所以要把較長點(diǎn)數(shù)的離散傅里葉變換分解為兩個(gè)較短點(diǎn)數(shù)的離散傅里葉變換，進(jìn)而得到快速傅里葉變換算法。那么，如何減少DFT的計(jì)算量呢？觀察離散傅里葉變換定義式，可利用WnkN的固有性質(zhì)來實(shí)現(xiàn)。

利用WnkN的對稱性、周期性和可約性，可以將長序列的DFT分解為短序列的DFT，從而減小計(jì)算量。本章只講解在時(shí)間和頻率進(jìn)行基2的分解算法，即參與運(yùn)算的點(diǎn)數(shù)經(jīng)過一次分解變?yōu)樵瓉淼囊话?，所以?jì)算量大約也可以減半，經(jīng)過連續(xù)的分解，最后只進(jìn)行2點(diǎn)離散傅里葉變換。這種分解既可以在時(shí)域進(jìn)行，也可以在頻域進(jìn)行。

4.2基2時(shí)間抽取的快速傅里葉變換

設(shè)N為2的整數(shù)次冪，即N=2M。在計(jì)算N點(diǎn)離散傅里葉變換時(shí)，按序列x(n)序號的奇偶，把x(n)分成兩個(gè)長度均為N/2的序列x1(n)和x2(n)，即顯然x1(n)和x2(n)剛好取遍序列x(n)的N個(gè)序列值。

由離散傅里葉變換的定義得

式(4-17)和式(4-18)的計(jì)算可用如圖4-1所示的蝶形運(yùn)算符號來示意。顯然，圖4-1(b)和圖4-1(a)等價(jià)，而圖4-1(b)的運(yùn)算結(jié)構(gòu)像蝴蝶，所以把這種運(yùn)算稱為“蝶形運(yùn)算”。

圖4-1基２時(shí)間抽取FFT的蝶形運(yùn)算

圖4-2為采用蝶形運(yùn)算對N=8點(diǎn)序列進(jìn)行一次時(shí)間抽取計(jì)算DFT的示意圖。圖4-2進(jìn)行一次時(shí)間抽取后的信號流圖(N=8)

在第一次抽取(分解)中，蝶形運(yùn)算中的復(fù)數(shù)乘法

總共需要兩個(gè)N/2次蝶形運(yùn)算，所以總的運(yùn)算量為N/2次復(fù)數(shù)乘法。為計(jì)算N/2點(diǎn)離散傅里葉變換X1(k)和X2(k)需要的運(yùn)算量為次復(fù)數(shù)乘法。與直接進(jìn)行N點(diǎn)離散傅里葉變換的運(yùn)算量相比，經(jīng)過一次分解后的運(yùn)算量降低一半。

為了計(jì)算X1(k)，同樣按序號的奇偶把x1(n)分解成兩個(gè)長度均為N/4的序列，即

先計(jì)算x3(n)和x4(n)的N/4點(diǎn)離散傅里葉變換X3(k)和X4(k)，然后用蝶形算法即可完成X1(k)的計(jì)算，有

用同樣的方法可以完成X2(k)的計(jì)算。需要注意的是，每次分解都是把上一次分解得到的序列按序號的奇偶分解成兩個(gè)序列，然后對每一對序列進(jìn)行蝶形運(yùn)算。圖4-3為采用

蝶形運(yùn)算進(jìn)行兩次時(shí)間抽取完成N=8點(diǎn)DFT的信號流圖。

圖4-3進(jìn)行二次時(shí)間抽取后的信號流圖(N=8)

最后，將兩點(diǎn)的DFT用蝶形運(yùn)算表示，圖4-4為完全分解三次得到的8點(diǎn)FFT信號流圖。

圖4-4時(shí)域抽取完全分解得到的8點(diǎn)FFT

設(shè)序列長度滿足N=2M，這樣徹底分解為兩點(diǎn)DFT需要進(jìn)行M級次分解。由圖4-4可以看出，每一級都包含N/2個(gè)蝶形，每個(gè)蝶形需要１次復(fù)數(shù)乘法，復(fù)數(shù)乘法次數(shù)為NM/2=N1gN/2。直接計(jì)算N=2M點(diǎn)的ＤＦＴ需要N2次復(fù)數(shù)乘法，所以FFT相對于DFT的復(fù)數(shù)乘法運(yùn)算效率為

當(dāng)N=210時(shí)，上式約等于205；N=211時(shí)，上式約等于372。顯然N越大，F(xiàn)FT的運(yùn)算效率越高。

4.3基2頻率抽取的快速傅里葉變換

基2時(shí)間抽取的FFT算法將輸入序列x(n)，根據(jù)序號的奇偶逐次分解成較短的子序列完成X(k)的計(jì)算，而基2頻率抽取的FFT算法，根據(jù)輸出序列X(k)序號的奇偶逐次分解成較短的子序列。下面先介紹基2頻率抽取的第一次分解過程。

顯然在區(qū)間0≤n≤N/2-1內(nèi)，x(n)對應(yīng)原序列的前半部分序列值，而xn+N/2)對應(yīng)原序列的后半部分序列值。為了得到x(n)的N點(diǎn)離散傅里葉變換X(k)，只需先求得N/2點(diǎn)離散傅里葉變換X1(k)和X2(k)，它們分別與X(k)奇數(shù)序號和偶數(shù)序號的N/2個(gè)序列值對應(yīng)。由x(n)和

得到x1(n)和x2(n)也可以由“蝶形運(yùn)算”完成，如圖4-5所示。

圖4-5基2頻率抽取FFT的蝶形運(yùn)算

圖4-6為采用蝶形運(yùn)算對原序列進(jìn)行一次頻率抽取計(jì)算DFT的示意圖。圖4-6進(jìn)行一次頻率抽取計(jì)算N=8點(diǎn)DFT

這種分解一直進(jìn)行下去。圖4-7為采用蝶形運(yùn)算進(jìn)行兩次頻率抽取計(jì)算N=8點(diǎn)DFT的信號流圖。與時(shí)間抽取算法完全一致的是，當(dāng)N=2M時(shí)這樣的分解可以進(jìn)行M次。每次

分解都是把上一次分解得到的序列前后兩半序列分解成兩個(gè)序列，對得到的兩個(gè)序列再進(jìn)行蝶形運(yùn)算。全部分解完成后對得到的所有2點(diǎn)序列進(jìn)行2點(diǎn)DFT變換即可。圖4-8為采用蝶形運(yùn)算進(jìn)行頻率抽取完全分解后得到的N=8點(diǎn)FFT的信號流圖。對照圖4-4和圖4-8可以看出，頻率抽取的運(yùn)算量與時(shí)間抽取相同。

圖4-7進(jìn)行二次頻率抽取后的信號流圖(N=8)

圖4-8

頻率抽取完全分解得到的8點(diǎn)FFT

4.4離散傅里葉反變換的快速算法

利用離散傅里葉正、反變換定義的對稱性，并借助于FFT可以很方便地實(shí)現(xiàn)離散傅里葉反變換的快速算法———IFFT?；仡欕x散傅里葉反變換的定義式

可以看出，式(4-35)右邊即為X?(k)的DFT(這里時(shí)域序號為k，頻域序號為n)。這就得到了第一種實(shí)現(xiàn)IFFT的方法，其步驟為：

(1)對輸入序列X(k)取共軛得到X?(k)；

(2)對X?(k)利用FFT計(jì)算N

點(diǎn)DFT；

(3)將前述結(jié)果除以N

并取共軛即可得X(k)的離散傅里葉反變換x(n)。

對離散傅里葉反變換式(4-33)右邊取兩次共軛得

式中，右邊大括號內(nèi)的部分即為X?(k)的DFT(這里時(shí)域序號為k，頻域序號為n)。這就得到了第二種實(shí)現(xiàn)IFFT的方法，其步驟為：

(1)對輸入序列X(k)取共軛得到X?(k)；

(2)對X?(k)利用FFT計(jì)算Ｎ點(diǎn)DFT；

(3)將對前述結(jié)果取共軛并除以Ｎ即可得X(k)的離散傅里葉反變換x(n)。

另外，也可以這樣理解離散傅里葉反變換式(4-33)：

4.5用FFT計(jì)算序列的線性卷積

線性時(shí)不變系統(tǒng)對任意輸入的響應(yīng)是輸入序列和系統(tǒng)沖激響應(yīng)序列的線性卷積，如果在一定的條件下線性卷積與循環(huán)卷積結(jié)果一致，則可以借助于FFT完成線性卷積的計(jì)算。分析表明，只要循環(huán)卷積的點(diǎn)數(shù)不比有限長序列卷積結(jié)果的長度小，則循環(huán)卷積與線性卷積相等。

首先討論兩個(gè)有限長序列的線性卷積。設(shè)有限長序列x1(n)的長度為M1，x2(n)的長度為M2，其線性卷積的長度為L=M1+M2-1。這時(shí)x1(n)和x2(n)的線性卷積yL(n)為

回顧x1(n)和x2(n)兩者N點(diǎn)循環(huán)卷積yc(n)的定義，有

式中，N≥max(M1，M2)?？紤]到循環(huán)卷積結(jié)果yc(n)的序號取值范圍為0≤n≤N-1，可以把式(4-4