第六章 數字信號分析(I)-DFT與FFT

第六章數字信號分析（Ⅰ）——DFT與FFT●

計算機技術的發(fā)展是數字信號分析的基礎?！?/p>

數字濾波優(yōu)于模擬濾波，（1）速度快，例如，采用數字信號分析技術對于1024采樣點進行A/D轉換，僅需4～15μs，進行FFT運算須250ms，較快的只需數毫秒；一個蝶形FFT硬件運算只需2μs。（2）分辨力高，在高頻段（50kHz)可達25Hz；在超低頻段可達0.0025Hz?！駭底中盘柗治鲆话惆ǎ侯l譜分析與數字濾波等主要內容。前者又包含相關與統(tǒng)計分析等。第一節(jié)模擬信號離散化●

本章內容：主要介紹離散Fourier變換與快速Fourier變換的基本原理及應用。一、A/D與D/A轉換1、A/D轉換過程（1）采樣（或稱抽樣）連續(xù)時間信號x(t)→離散采樣信號x(n△t)。

量化誤差呈等概率均勻分布，其概率密度函數p(R)=1/R，最大量化誤差為±0.5R，其均方差為（2）關于量化誤差問題或量化增量R取決于A/D轉換器位數。例如，8位A/D，R=Vref/28。（3）編碼——將離散幅值經過量化以后變?yōu)槎M制數字D。信號幅值2、D/A轉換過程譯碼——把數字信號RD恢復為有限幅值A的過程，即

D/A轉換過程包括：譯碼與波形復原。波形復原——把離散幅值恢復為連續(xù)波形的過程，由保持電路實現。例如，零階保持與一階多角保持等。前者是在兩個采樣值之間，令輸出保持上一個采樣值的值；后者是在兩個采樣值之間，令輸出為兩個采樣值的線性插值。經過保持變換構成的信號存在不連續(xù)點，用模擬低通濾波器消除輸出波形的不連續(xù)點。

采樣過程——是通過采樣脈沖序列與連續(xù)時間信號x(t)相乘來完成的。根據采樣脈沖序列的形狀，分為理想脈沖采樣與矩形脈沖采樣。二、采樣信號的Fourier變換1、時域采樣（1）理想脈沖采樣采樣脈沖序列：采樣信號：采樣信號頻譜沿頻率軸每隔一個采樣頻率ωs，重復出現一次，即頻譜產生周期延拓。幅值被Cn所加權，故頻譜形狀不變。（2）矩形脈沖采樣周期矩形脈沖序列的傅立葉變換：故有：可見：Xs(ω)是X(ω)在以ωs為周期的重復過程中，其幅值按sinc(nωsτ/2)規(guī)律變化的函數。2、頻域采樣采樣脈沖序列采樣間隔則有可見：若頻譜X(ω)被間隔為ω1的脈沖序列在頻域中采樣，則在時域中等效于x(t)以T1(=2π/ω1)為周期而重復。就是說，周期信號的頻譜是離散的。根據因此，

上述分析證明：信號的時域與頻域呈采樣（離散）與重復（周期)關系。采樣定理：ωs≥2ωm或fs≥2fm。因為時域采樣間隔決定于fs，故又稱為時域采樣定理。三、采樣定理1、頻混現象F[xs(t)]為周期譜，其周期ωs=2π/Tsωs>2ωm,周期譜圖相互分離ωs<2ωm,周期譜圖相互重疊2、采樣頻率3、信號復原為了從Xs(ω)中無失真地選出X(ω)，用頻域矩形窗函數H(ω)與Xs(ω)相乘，得

實現方法：將采樣信號xs(t)通過傳遞函數為H(ω)的理想低通濾波器，則在濾波器輸出端可得到頻譜為X(ω)的連續(xù)信號x(t)。理想濾波器的H(ω)為：所以有：●連續(xù)信號可以展成正交采樣函數（sinc(t)型函數）的無窮級數，級數的系數等于采樣值x(nTs)。也可以說，若在采樣信號xs(t)的每個采樣值上畫一個峰值為x(nTs)的sinc(t)型函數波形，則合成波形就是x(t)。所以，若xs(t)通過理想低通濾波器時，每個采樣值產生一個脈沖響應，這些響應進行疊加就得到x(t)，從而達到由xs(t)恢復x(t)的目的?！?/p>

為什么稱sinc(t)為內插函數？所謂內插是指從已知離散點的值，求在離散點之間另外一些點的值時，所采用的數學插值法。理論上講，對不在取樣時刻任意點的數值應該是無限加權樣值的總和，但由于這里的內插函數是衰減的，因此，實際上可由該點附近的一組有限值之和而得到良好的逼近?！癞敠豷=2ωm，ωc=ωm，各個采樣的沖激響應零點恰好落在采樣時刻上。就采樣點的數值而言，在這種情況下，各個沖激響應互相不產生“串擾”?！裆鲜隼脼V波器由采樣信號恢復原信號的方法，又稱為惠特克波形重構法或理想內插法。若x(t)是時域有限信號，并集中在-tm～tm的時間范圍內。若在頻域中以不大于1/2tm的頻率間隔對X(ω)進行采樣，則采樣信號X1(ω)可唯一地表示原信號?！駯艡谛l域采樣后，只能獲得采樣點的頻率成分，其余的頻率成分一概被舍去。這猶如透過柵欄觀賞光景，只能看到一部分，就可能使一部分有用的頻率成分被漏掉，而丟掉了部分有用信息，此現象稱為柵欄效應。4、頻域采樣定理第二節(jié)離散Fourier變換（DiscreteFourierTransform,DFT)●

DFT并非泛指對任意離散信號取傅里葉積分或傅里葉級數，而是專指適用于計算機計算的FT。這是因為：●

DFT表達式的導出方法：

（1）從離散時間序列的Z變換基礎上導出，即有限長序列的離散傅立葉變換可解釋為它的Z變換在單位圓上的采樣；原因：對x(t)進行FT或IFT運算時，無論在時域或在頻域，都需進行包括(-∞,∞)區(qū)間的積分運算。而要實現DFT，則必須做到：1）把連續(xù)信號（時域或頻域）改造為離散數據；2）把計算范圍收縮到一個有限區(qū)間；3）實現正、逆傅立葉變換運算。

（2）把DFT作為連續(xù)信號傅立葉變換的一種特殊情況來導出。其物理概念較前者更清楚。一、離散Fourier變換關系式1、時域采樣2、時域截斷對連續(xù)時間信號進行DFT,一般包括時域采樣,時域截斷,頻域采樣三個步驟.

皺紋問題：因矩形函數有突變階躍點，在時域截斷后，反映在頻域上產生皺紋。此即由于Gibbs現象產生的能量泄漏效應。

★采樣信號經過截斷處理后，雖然在時域為有限長的離散樣本，但頻域內仍為連續(xù)函數，若實現逆變換，還必須改造頻域函數為有限離散值。3、頻域采樣采樣脈沖序列：令頻域采樣脈沖序列為δ1(f)，1）根據頻域采樣定理（f0≤1/2tm)，選取采樣間隔f0=1/T0(T0為時域截斷信號分布區(qū)間，它相當于2tm)。2）又根據FT的對稱性，對應的時域函數為：采樣信號：被δ1(f)采樣后的頻域采樣信號為：逆傅立葉變換表明：是周期為T0的離散函數，每個周期內有N個離散點。δ函數的卷積)特性積分是在一個周期T0內進行，即r=0δ函數的篩選(積分)特性由于是周期函數，所以其傅立葉變換也是等間隔脈沖序列，即T0=NTs這表明：是

的傅立葉變換，在頻域內（k區(qū)間(-∞,∞))是一個被Ck所加權的周期性離散脈沖序列，每一個周期內有k=N個樣本點。

Ck表示頻譜中的一個周期內N個采樣點的復數值，以下用X(kf0)來表示同理可證明X(kf0)的逆傅立葉變換：通常將上兩式寫成：

這就是DFT，它通過連續(xù)FT，將N個時域采樣點與N個頻率采樣點聯(lián)系起來，是連續(xù)FT的一種特殊情況.

●此式物理意義是十分明確的，在變換過程中，都僅僅涉及到了處理區(qū)間(0≤n≤N-1;0≤k≤N-1)的N個x(n)和N個X(k)值。實質上是表明了兩個N維矢量的相互線性變換（映射）。其中，x(n)可以分解為N個諧波復指數序列，每個諧波分量的頻率為f0或kf0，復振幅為(1/N)X(k)。同樣，X(k)被分解為N個復指數序列之和，每個分量的復振幅是x(n)，頻率是nf0?！駥嶋H信號x(n)在非處理區(qū)間，可能是零(時域有限信號)，可能是周期的(周期信號)，可能是非零、非周期的(隨機信號)，甚至是不定的，但只要x(n)在所關心的處理區(qū)間是確定的，就可沿用上述DFT與IDFT關系。二、DFT與FT的關系DFT與FT之間是一個近似，因為DFT需要采樣與截斷。而近似的程度是被分析波形的函數。1、頻域有限的周期信號，時域截斷長度等于周期時●時域截斷波形的FT為X(f)*δ0(f)*U(f)，與X(f)相比，產生嚴重畸變?！癞斶@一畸變圖形被δ1(f)采樣后[圖(g)]，畸變卻避免了。這是因為頻域采樣脈沖的間隔為1/T0，在這些采樣頻率坐標點1/T0,2/T0，…處，圖(e)中的實線除了在±1/T0點有數值外，其余點都是零。

●●x(t)經采樣、截斷后的幅值關系設x(t)的時域幅值為A，★直接作FT變換以后頻域幅值為A/2。x(t)經時域采樣、截斷再卷積以后([x(t)δ0(t)u(t)]*δ1(t))其幅值為AT0，不是原來的A。x(t)在頻域內的X(f)經過卷積及頻域采樣后([X(f)*δ0(f)*U(f)]δ1(f))其幅值為AT0/2Ts，也不是原來X(f)的A/2。因此，如果希望用DFT來計算FT，就必須對DFT的變換偶對間的常數因子作一些調整，即將離散時間函數乘以因子Ts/T0，這樣就得到頻率函數所要求的A/2的常數因子，所以:●以上分析表明，對頻域有限的周期信號，當截斷長度等于其周期（或周期的整數倍）時，DFT與FT之間的差別僅僅是一個比例因子T0/Ts。因此，欲使DFT與FT等價，就必須：（1）時間函數x(t)是周期性的；（2）x(t)是頻域有限信號；（3）采樣頻率至少是x(t)的上限頻率的2倍；（4）截斷函數u(t)必須正好在x(t)的一個周期內(或整數倍周期)是非零的.★更正2、頻域有限的周期信號，時域截斷長度不等于周期時

●時域內：信號周期為T’0，截斷長度為T0，且T’0≠T0。頻域內：X(f)*δ0(f)*U(f)，注意這個卷積圖形是頻率間隔為1/T’0的脈沖的卷積，而sinc(t)型函數的零點是1/T0,2/T0，…。卷積以后的疊加圖形在1/T0,2/T0，…等點處并不是零。這表明：由于截斷后的波形不是周期的整數倍，故其積分平均值不為零。當這一波形被頻率間隔為1/T0的脈沖作頻率采樣以后，零點不再恰好與每個樣本點（脈沖采樣點）相重合，故引起DFT與FT之間的差異。產生這一現象的原因解釋：從時域看，不按周期的整數倍對x(t)進行采樣和截斷，會產生具有間斷點的周期函數(圖中(g))，這些劇烈變化，將在頻域中產生附加的頻率成分。從頻域看，時域截斷等效于sinc(t)型函數與X(f)的單個脈沖的卷積，結果頻率函數不再是一個脈沖，而是頻率的連續(xù)函數，在原來脈沖的位置上，這個連續(xù)函數具有局部最大值，還有一系列稱為旁瓣的峰值。這些旁瓣在頻域抽樣后就造成附加的頻率成分，這就是所謂泄漏效應，因為時域中截斷是必須的，所以泄漏效應是DFT所固有的。3、時域有限而頻域無限信號

●這類信號在時域采樣后必然會產生頻混現象，因此采樣時間間隔的選擇必須使頻混現象減小到允許的限度以內。

●如果對有限長的波形進行采樣，選擇的采樣點數N正好等于它的樣本點總數，則不必在時域進行截斷。截斷被省略，經時域采樣后的函數作傅立葉變換以后[圖(c)]，與頻域采樣脈沖δ1(f)相乘，這個乘積在時域等效于圖(c)和(d)所示的時間函數的卷積。最后所產生的波形是周期的，周期由原函數的N個樣本確定，所以它是原信號x(t)的復制品。這個周期函數的傅立葉變換即為圖(e)所示的頻域采樣后的波形。

●對于這類函數，如果選擇N等于時域有限信號的樣本點數，則誤差僅由頻混造成。如果選擇采樣間隔足夠小，就可以減少由頻混引起的誤差。因此，在這種情況下DFT與FT可很好的一致起來。(a)(b)(c)(d)(e)4、一般周期信號

●對于時域無限頻域無限的一般周期信號，采樣后存在著頻混效應，如果選取截斷函數精確地等于周期的整數倍，那么將不產生泄漏效應。這種情況下變換的誤差源主要是頻混效應，若時域截斷不等于周期，則產生泄漏效應.

●實際情況中所遇到的信號，往往既不是有限時間，也不是有限帶寬或周期性的信號。這類任意信號在作離散傅立葉變換中，存在著頻混和泄漏效應。減小采樣間隔可以減小頻混效應；改善截斷函數，即選擇合適的窗函數，可以減小或抑制泄漏現象。5、任意信號★可見，如果處理恰當，在許多應用中都可以用DFT來得到本質上和FT等價的結果。其中，值得記住的一個重要概念是，DFT意味著在時域、頻域兩方面都周期化，而時域函數的N個樣本點則表示周期化以后所形成的新周期函數的一個周期?！镒詈髴赋?，DFT實際上是建立了函數x(t)的N個時間樣本點與N個頻率樣本點之間的互換關系，利用這一關系可從x(n)計算X(k)，也可以從X(k)計算x(n)。這種相互關系，在數據處理的數字量分析法中經常遇到，例如，從數據的相關函數計算功率譜密度函數，也可以從功率譜密度函數計算相關函數。三、DFT的性質（與連續(xù)傅立葉變換的性質類同）1、線性2、時移特性有限長序列x(n)的移位x(n)x(n)有限長序列的圓移位(或循環(huán)移位)

●時移特性——若將x(t)沿時間軸位移t0，則其FT要乘以因子e-j2πft0，即與此類同，如果DFT[x(n)]=X(k)，則3、頻移特性DFT：若DFT[x(n)]=X(k)

則DFT[x(n)ej2πl(wèi)n/N]=X(k-l)

或

DFT[x(n)W-ln]=X(k-l)

FT：若F[x(t)]=X(f)

則F[x(t)ej2πf0t]=X[(f-f0)]

時間函數x(n)乘以指數項ej2πl(wèi)n/N，則DFT向右圓移l單位。4、離散卷積

依據運算方式不同，可分為線卷積與圓卷積。（1）時域線卷積或已知x(n)和h(n)兩有限序列，卷積如圖示x(n)1234h(n)12344812163691224681234302011420114y(n)nnn（2）時域圓卷積用圓卷積方法計算上例，公式如下：y(n)=x(n)

*

h(n)圓卷積符號

可見，圓卷積與線卷積所得結果不同。這是因為線卷積過程中，經反折再向右平移，在左端將依次留出空位。而圓卷積過程中，經反折的圓移序列向右移去的樣值又從左端循環(huán)出現，這使得兩種情況下相乘疊加而得之數值不同。01234

為解決圓卷積與線卷積結果不同的問題，將x(n)和h(n)都適當地補一些零值，以擴展其長度。那么在作圓卷積時，向右移去的零值，從左端出現仍取零值，這樣就與線卷積的情況相同.補零擴展的條件為：

x(n)的樣本點數h(n)的樣本點數432121000438300000432110004324000000432132100041262000043210432100012620004321004321000830004321432100016941000432100043210004000x(m)x(m)x(m)x(m)x(m)x(m)x(m)h(0-m)h(1-m)h(2-m)h(3-m)h(4-m)h(5-m)h(6-m)411411202030●若選取的L不夠長，圓卷積將首尾交疊混淆，其結果與線卷積不一致（這可看做是一種混疊現象）。取L≥N1+N2-1，可避免這一現象?！?/p>

圓卷積可利用快速傅立葉變換技術，實現快速卷積。因此，對于有限長序列求線卷積的問題，可轉化為圓卷積來求解，以便利用快速傅立葉變換技術。（3）離散時域卷積定理運用這一定理，可對兩個時域周期序列x(n)與h(n)分別計算離散傅立葉變換，再將結果相乘，然后計算乘積的離散傅立葉逆變換，即可得兩個時域周期序列的卷積。這一定理為用快速傅立葉變換計算時域卷積提供了依據。（4）離散頻域卷積定理即兩個周期為N的時域周期采樣函數，它們的乘積的離散傅立葉變換等于它們的離散傅立葉變換的卷積。5、離散相關定理變換對稱為離散相關定理，即兩個周期為N的時域周期序列，它們的時域離散相關的離散傅立葉變換等于它們的離散傅立葉變換的乘積。運用這個定理，可以等效地在頻域中確定相關性。連續(xù)函數時域卷積定理：離散時域卷積定理：6、巴什瓦定理對于離散信號，時域功率和頻域功率之間的關系，由下式給出：若x(n)為實序列第三節(jié)快速Fourier變換（FastFourierTransform,FFT)

●

FFT是一種減少DFT計算時間的算法。在FFT出現之前，雖然DFT為離散信號的分析從理論上提供了變換工具，但是很難實現，因為計算時間很長。例如，對采樣點N=1000，DFT算法運算量約需200萬次，而FFT僅約需1.5萬次。

●

FFT方法于1965年由美國庫利-圖基(J.W.Cooley-J.W.Tukey)首先提出，曾被認為是信號分析技術的劃時代的進步。一、FFT算法的基本原理1、DFT的計算量N×N對稱方陣，Wnk=[Wnk]TN×N對稱方陣，W-nk=[W-nk]T計算量：每計算一個X(k)值，需進行N次復數相乘和(N-1)次復數相加。當計算X(0)，X(1)，…共N個X(k)值時，需要N2次復數相乘，N(N-1)次復數相加。矩陣[W]與[x(n)]相乘過程中存在不必要的重復運算，這是簡化運算的關鍵.以N=4為例討論如下：分析：（1）不必要的計算①W0=1，②

WN/2=[e-j2π/N]N/2=-1①

Wnk的周期性，即Wnk=Wn(k+N)=Wk(n+N)，當N=4，有W2=W6，W1=W9等（2）可利用的特性②

Wnk的對稱性，即W(nk+N/2)=-Wnk，當N=4，有W3=-W1，W2=-W0等原計算式周期性簡化對稱性簡化這就是庫利-圖基FFT算法的基本思想。

2、減小運算工作量的途徑

FFT算法有多種變型，其算法很多，但每種變型的建立，多是考慮了被分析數據的特性，或者利用計算機特性，或者利用專用計算機FFT硬件特性等。本節(jié)以基2FFT算法作為討論的起點，因為它包含了FFT算法的基本要素，運算過程比較單純，適于人們學習。3、FFT計算方法

●基2算法設一個點序列x(n)，要求采樣點數N=2M，M為正整數。

●基2算法的出發(fā)點把N點DFT運算分解為兩組N/2點的DFT運算，即把x(n)按n為偶數和n為奇數分解為兩部分。將x(n)的DFT計算分為奇、偶兩部分以2r表示偶數n，2r+1表示奇數n，r=0,1,2,…,(N/2-1)一個N點的DFT被分解為兩個N/2點的DFT

●必須注意

G(k)和H(k)只有N/2個點，k=0,1,2,…,N/2-1。而X(k)卻需要N個點，k=0,1,2,…,N-1。如果以G(k)和H(k)表達全部X(k)，應利用G(k)和H(k)的兩個重復周期，由周期性可知●加權系數WkN為：●計算X(k)的全部關系式：G(k)和H(k)可分別看成是序列x(2r)與x(2r+1)的N/2點DFT

此式表明，一個N點的DFT可分解成兩個N/2點的DFT

●采用蝶形流程圖，以N=4為例，說明上式的計算過程。由G(0),G(1),H(0),H(1)計算X(0),X(1),X(2),X(3)蝴蝶結

一個蝴蝶結包括兩次復數乘法，兩次復數加法，但其中有重復。H(0)與W04相乘以及與-W04相乘可以改成只與W04相乘，再分別加減，這樣就使運算量減少至只有一次復數乘法和兩次復數加(減)法。同理，用第二個蝴蝶結計算X(1)和X(3)，也只有一次乘法，兩次加(減)法。這樣，由G(k)和H(k)計算X(k)的過程中，包含N/2個蝴蝶結運算，共N/2次復數乘法，和N次復數加(減)法?！裼蓌(0),x(1),x(2),x(3)計算G(0),G(1),H(0),H(1)時的計算量矩陣表示W12=e-j2π/2=-1=-W02實現N=4的FFT計算全過程，所需乘法運算次數為2×N/2(4次)，加法運算次數為2×N(8次)。顯然，這比用DFT直接計算所需(乘法N2=16，加法N(N-1)=12）運算工作量大為減小。2次乘法，4次加減法●N=23=8的FFT計算量分三級蝶形運算，每級需乘法N/2次，加法N次。全部運算3×N/2=12次乘法，3×N=24次加法。而直接DFT運算量N2=64次乘，N(N-1)=56次加?！馧=2M的FFT計算量運算分解為M級蝶形圖，每級包含N/2次乘，N次加，故算法的工作量為復數乘法復數加法●直接DFT與FFT算法所需乘法次數的比較。當N=256，二者比值為64；當N=1024，其比值為204.8二、FFT算法的應用

FFT是實施DFT的一種快速算法，FFT的應用實質上是DFT的應用。FFT可直接用來處理離散信號數據,也可用于對連續(xù)時間信號分析的逼近?！緦嵗?】單邊指數函數的DFT函數傅立葉變換1、對連續(xù)時間變換的逼近FT圖DFT圖令β=1，a=1令β=1，a=1，利用DFT方法，按下列步驟得到離散變換圖：①時域采樣，設N=32，Ts=0.25，得到每點樣值x(nTs)；②計算

可見，DFT是FT的逼近，其實部是偶函數，在頻域k=N/2點對稱，在k>N/2時，代表了負頻率點處理的結果。而虛部為奇函數，k>N/2處是負頻率處理結果?！緦嵗?】方波的諧波分析將DFT用于方波的諧波分析，需要計算以得到各次諧波系數值。對方波在時域的采樣點數N=32。在k=N/2點處對稱?？梢钥闯?，低次諧波比較逼近，而高次諧波有誤差，這是由于頻率混疊效應所致。雖然可以通過提高采樣頻率來減少這一現象，但不可能完全避免，因為周期方波為時域無限、頻域無限信號。2、卷積運算（1）快速卷積運算過程長度為N1的序列x(n)和長度為N2的序列h(n)卷積，其結果y(n)長度為N1+N2-1●卷積運算中，每個x(n)的樣值必須與每個h(n)的樣值相乘，因此，共需要N1×N2次乘法運算。

●如果把線卷積改為求圓卷積，并借助FFT技術，可減少運算量?！窨焖倬矸e運算過程實現快速卷積算法中，由于利用了DFT分析，即時域或頻域都是周期性的離散數據，當對他們作卷積運算時，將出現一種周期數據之間的疊帶求和現象，給計算結果帶來一種所謂的環(huán)繞誤差。下面分析其產生原因和避免方法。第一步：利用FFT算法計算兩信號的DFT第二步：在各頻率點處兩信號的變換相乘第三步：運用IFFT算法，計算變換式乘積的反變換

實現這一過程共需兩次FFT和一次IFFT運算（相當于三次FFT運算），此外，完成X(k)與H(k)兩序列相乘，需作N次乘法。在一般的有限沖激響應（FIR）數字濾波器中，由h(n)求H(k)這一步是預先設計好的，數據已置于存貯器中，故實際只需兩次FFT的運算量。如果假定N1=N2=N，則全部復數乘法運算次數為。可見，隨N值增大，計算量顯著減少，故圓卷積的方案可以快速完成卷積運算。（2）圓卷積的環(huán)繞誤差●環(huán)繞誤差產生原因

圖①為兩個非周期離散序列的卷積。采用直接線卷積或補零圓卷積方法很容易求得其計算結果[圖中(g)]。兩個非周期離散序列的卷積①兩個周期離散序列的圓卷積②

原因：當采用DFT分析方法，上述兩個非周期離散信號被改造為時域、頻域相對應的周期離散信號，導致卷積結果與圖①不同（見圖2）。主要區(qū)別：當h(n-m)向右移動時，h(n-m)另一周期的一部分進入到求和區(qū)域，導致錯誤的計算結果，被稱為“環(huán)繞誤差”或“疊帶效應”?！癖苊猸h(huán)繞誤差的方法

方法：對x(n)與h(n)分別在尾部填補N1(x(n)的樣點數)與N2(h(n)的樣點數）零值點，即使其周期加倍。如果x(n)與h(n)的長度相等，則都加長N點。采用補點方法后所得計算結果如圖。

補點后的副作用：對x(n)與h(n)進行補點后，避免了環(huán)繞誤差的同時，使的DFT(或FFT)算法所需容量加倍，在各DFT表示式中，必須用2N來代替N，在各函數的尾部補填足夠的零值使有效周期為2N，這就能夠對兩個含有N點的函數進行正常的卷積運算。兩個含有N點的非周期性函數的離散卷積給出一個具有2N-1點新函數，當在原函數上填加N個零點后，所得圓卷積的周期為2N，比原來的非周期信號的卷積多一個點，每周內多一個附加零點。如果這兩函數相當靠近，但長度不等，則首先將短函數的尾部補零使與長函數的長度相等，然后再補零到2N-1，這也即對較短信號補充的零點數總共超過了N個?？偨Y以上各點，快速卷積過程可按如下步驟進行：（a）用補零法修正x(n)和h(n)，以避免環(huán)繞誤差的出現。（b）用FFT算法計算兩個修正后的函數的DFT，得到X(k)與H(k)（c）將X(k)與H(k)相乘，得到（d）利用FFT算法，計算出Y(k)的IDFT，即3、相關運算相關函數的數字計算方法有時域直接計算與FFT快速算法。直接計算方法是依據下述定義進行的，即互相關函數（在數字信號分析中，一般用符號r）自相關函數：這種計算方法與卷積運算相類同（卷積多一個時間反折），也是一個乘、加序列，所需計算量很大?！锔锔?/p>

●相關函數的FFT算法，依據的是維納-辛欽關系，即自相關函數或互相關函數可以由功率譜密度或互譜密度函數來求得。

①這種方法是一種迂回的方法，但它比直接時域計算方法快5～100倍。

②當運用FFT方法計算相關函數時，也必須注意到環(huán)繞誤差的影響，它類似于圓卷積中的誤差。解決的方法也是對時間序列x(n)與y(n)補零擴展。第一步：對x(n)和y(n)作FFT分析，得到復頻譜X(k)與Y(k)第二步：對X(k)與Y(k)作共軛乘積，得到或當x(n)與y(n)相同時，得到自功率譜密度:

第三步：作IFFT分析，從功率譜密度獲得相關函數，即

★本節(jié)內容：隨機信號的譜分析與譜估計技術，即對功率譜密度的傳統(tǒng)估計方法。

★DFT與FFT是信號處理的重要工具，尤其是DFT的基本概念與算法，為信號頻譜分析的各種應用鋪平了道路。

★傳統(tǒng)譜分析方法，是基于Fourier變換的譜分析方法，包括相關函數法與周期圖法。第四節(jié)譜分析與譜估計

★

譜分析與譜估計在生產實踐與科學研究中獲得了日益廣泛的應用?！纠?】在聲納系統(tǒng)中，通過對噪聲信號進行譜分析，以判斷水面艦艇或潛艇的運動速度、方向、位置、大小等；【例2】對飛機、輪船、汽車、汽輪機、電機、機床等主體或部件進行實際運動的譜分析，可以提供設計數據和檢驗設計效果，或者尋找振源和診斷故障。

1）相關函數法(又稱BT法)1958年由布萊克曼-圖基(Blackman-Tukey)提出。它是通過統(tǒng)計分析，從時域上先求信號的自相關函數，再作Fourier變換，求得功率譜估計值。2）周期圖法它是直接將數據進行Fourier變換，再取其幅度平方，得到信號的功率譜密度。此法是在1965年由Cooley-Tukey提出的FFT方法問世以后，被用于譜估計?！锾貏e注意一、周期圖法作功率譜估計

★自相關函數作為時移的函數是最能較完整地表征隨機信號的特定統(tǒng)計平均量值的。而隨機信號的功率譜密度，正是自相關函數的Fourier變換。

★對于隨機信號而言，其自身的Fourier變換是不存在的，只能用功率譜密度來表征它的統(tǒng)計平均譜特性。根據相關定理與維納-辛欽關系式（參見圖2-5）易于證明隨機信號序列x(n)的功率譜密度：1)“泄漏”問題以上兩種傳統(tǒng)方法本質上是一樣的，都認為“有限長的數據段，可以看作是無限長的取樣序列給予開窗截斷后的結果”。不論是數據開窗，還是自相關函數開窗，在頻率域內都會發(fā)生“泄漏”現象，即功率譜主瓣內的能量泄漏到旁瓣內。這樣，弱信號的主瓣很容易被強信號的旁瓣所淹沒或歪曲，造成譜的模糊與失真。

2）所有旁瓣抑制技術，都是以損失譜分辨率為代價。3）在譜分析應用中，頻率分辨率與低旁瓣一樣是個重要指標，有時甚至更重要。因此，解決高分辨率與低旁瓣的矛盾是譜分析中的一個重點問題。

●

為什么叫周期圖法？

由于序列x(n)的離散傅里葉變換X(k)具有周期函數的性質，故把它稱為長度為N的實平穩(wěn)隨機信號序列x(n)的周期圖。

●

周期圖法作譜估計時，存在的兩個主要問題：功率譜密度的統(tǒng)計變異性和能量泄漏。前者是統(tǒng)計誤差，是由于在功率譜測量中收集到的數據數量有限，而引起的不確定度。后者是譜分析中所固有的，它將造成估計的偏度誤差。因此，實際應用中對周期圖法進行修改，以盡量減小估計誤差。其方法是：（1）采取平均化處理，以減少統(tǒng)計變異性；（2）采用窗處理，以減少泄漏。其估計值周期圖法計算功率譜密度流程

統(tǒng)計變異性產生原因——對于有限長數據的處理，由于數據的概率性質，而產生了統(tǒng)計性誤差。二、平均化處理（以減少統(tǒng)計變異性）一般用變異系數（Coefficientofvariation，或稱為標準化標準偏差）來表征譜估計的變異性質，其定義為譜估計的均方差均值當用周期圖法作譜估計時，周期圖是一個復數，故而有可以證明，實部XR(k)與虛部XI(k)是等方差和零均值的兩個不相關的隨機變量。由于Fourier變換是線性運算，如果被分析數據x(n)是高斯分布，則XR(k)與XI(k)也是高斯隨機變量。|X(k)|2=X2R(k)+X2I(k)表明，譜估計值是兩個獨立高斯變量的平方和，它相當于具有兩個自由度的卡埃平方(Chi-Square)分布.從概率統(tǒng)計學可知，卡埃平方分布自由度為n的卡埃平方變量n為獨立變量數，即自由度

zi服從正態(tài)分布的獨立隨機變量1、譜估計的變異性變量χ2的均值和方差：變異系數

顯然，用式作譜估計時，相當于具有自由度n=2（實部和虛部），其估計的變異系數εr=1。這表明估計的相對誤差達到100%，即估計的變異性和被估計量一樣大，這樣的估計是不合用的。為此可采用平均化處理方法來提高估計精度.平均周期圖的方法——是將序列x(n)分段，求各段周期圖，再進行平均。2、平均化處理方法設序列x(n)[或x(t)]總體長度為N（或T），將其分為q段，每段長度為Nq（或Te），對每段數據Nq作譜估計，得到

，此時有如果各段頻率分量的實部XRj(k)與虛部XIj(k)是互為獨立的隨機變量，那么有此時χ2分布的自由度n=2q。將各段估計譜在對應的頻率上作q個估計量的平均，得

N此時變異系數為

●分析：對于連續(xù)隨機過程，樣本總體長度為T(T=N△t)，△t為采樣間隔，分段長度Te(T/q)，那么分析帶寬Be=1/Te，故有。這一關系式與模擬分析方法中所得到的結論是一致的，欲提高譜估計精度，必須同時考慮到樣本總體長度T與頻率分析帶寬Be。三、窗口函數因不可能對無限長信號進行分析，故必須截斷。假定截斷區(qū)間為(-T,T)，因為對|τ|>T時的Rx(τ)值假定為零，所得到的估計譜為近似譜，即正弦信號的真譜與估計譜之間的關系1、真譜與估計譜根據維納-辛欽定理，理論譜密度的定義為：兩者間總是有誤差例如，正弦信號的真譜與估計譜之間的關系。估計譜是真譜與窗譜的卷積，即理想窗譜W(ω)應為δ函數。此時估計譜與真譜完全相同。但要得到δ函數形窗譜，其時域窗口必然是無限寬，這不可能。實際窗口為有限寬，窗譜W(ω)為sinc(t)型函數,故有主瓣能量泄漏到旁瓣.

2、窗函數★窗函數設計的目的——改善窗譜形狀。

★

窗函數的基本要求——窗譜的主瓣要窄且高，以提高分辨率；旁瓣應小，正負交替接近相等，以減小泄漏或負譜現象。

★加窗的其他作用——可抑制噪聲，提高頻率辨識能力。例如，1）沖擊測量或脈沖激振時，實際有用信號延續(xù)時間很短。如果采樣時間較長，可利用窗口控制，避免在余留樣本時間內，無實際信號輸入時而有噪聲混入。2）在系統(tǒng)辨識中，小阻尼系統(tǒng)的響應衰減緩慢。如在較短時間內截斷，會丟掉有用信息；若加指數窗口，使其增加衰減速度，既保留了有用信息，又可防止噪聲信號混入?！锔纳拼白V形狀的基本思想——改善截斷處的不連續(xù)狀態(tài)。Gibbs現象的研究已經表明：時域內的間斷，反映到頻域，必然發(fā)生振蕩現象；反之，頻域內的間斷，反映到時域，也同樣發(fā)生振蕩現象。

1）冪窗：采用時間變量某種冪次的函數，如矩形、三角形、梯形或其他時間t的高次冪2）三角函數窗：應用正弦或余弦函數等組合成復合函數，如漢寧窗、海明窗等3）指數窗：采用指數時間函數，如e-at，如高斯窗等★窗函數主要類型（1）矩形窗（時間變量的零次冪窗）

優(yōu)點是主瓣比較集中；缺點是旁瓣較高，并有負旁瓣，導致變換中帶進了高頻干擾和泄漏，甚至出現負譜現象。（2）三角窗（也稱費杰（Fejer）窗，是冪窗的一次方形式）特點：主瓣寬約為矩形窗的兩倍，旁瓣小，且無負旁瓣（3）漢寧(Hanning)窗（又稱升余弦窗）漢寧窗w(t)●漢寧窗的譜窗為3個矩形窗的頻譜之和，即3個sinc(t)型函數之和。括號中的兩項相對于第一個譜窗向左、右各移動了π/T，從而使旁瓣互相抵消，消去高頻干擾和漏能.

●漢寧窗與矩形窗的譜圖對比圖(a)為W(ω)－ω關系；圖(b)為相對幅度(相對于主瓣衰減)－logω關系?？梢钥闯?，漢寧窗主瓣加寬（第一個零點在2π/T處）并降低，旁瓣則顯著減小。第一個旁瓣衰減-32dB，而矩形窗第一個旁瓣衰減-13dB。此外，漢寧窗的旁瓣衰減速度也較快，約為60dB/10oct，而矩形窗為20dB/10oct。由以上比較可知，從減小泄漏觀點出發(fā)，漢寧窗優(yōu)于矩形窗。但漢寧窗主瓣加寬，相當于分析帶寬加寬，頻率分辨力下降。octave——倍頻程（4）海明(Hamming)窗（也是余弦窗的一種，為改進的升余弦窗）

海明窗與漢寧窗都是余弦窗，只是加權系數不同。海明窗的加權系數使旁瓣更小。分析表明，海明窗的第一旁瓣衰減為-42dB。海明窗的頻譜也是由3個矩形時窗的頻譜合成，但其旁瓣衰減速度為20dB/10oct，比漢寧窗衰減速度(約為60dB/10oct)慢。（5）高斯窗（是一種指數窗）a為常數，決定了函數曲線衰減的快慢。適當選取a值，可使截斷點（T為有限值）處的函數值比較小，則截斷造成的影響就比較小。高斯窗譜無負的旁瓣，主瓣較寬，故頻率分辨力低。第一旁瓣衰減達-55dB。高斯窗函數常用來截斷一些非周期信號，如指數衰減信號等。

●除以上常用窗函數外，還有多種窗函數，如帕仁(Parzen)窗、布萊克曼(Blackman)窗、凱塞(Kaiser)窗等?！裎宸N典型窗函數的性能特點

【例】矩形窗、三角窗、漢寧窗分析同一信號數據。(a)為被分析信號的真譜；(b)和(c)是用兩種時寬(T/2與T)的矩形窗分析的結果?？梢?，時域窗口窄，分辨率低，相鄰的兩譜線不能分辨。(d)和(e)是分別用三角窗與漢寧窗分析的結果，與圖(c)相比，三角窗與漢寧窗分辨率低,不能分辨相鄰譜線，但由于旁瓣衰減快,譜的分布區(qū)域窄而邊沿清晰。窗函數類型-3dB帶寬等效噪聲帶寬旁瓣幅度(dB)旁瓣衰減速度(dB/10oct)矩形三角形漢寧海明高斯0.89B1.28B1.20B1.30B1.55BB1.33B1.23B1.36B1.64B-13-27-32-42-55-20-60-60-20-20(a)(b)(c)(d)(e)

●窗函數的選擇應考慮被分析信號的性質與處理要求?！魞H要求精確讀出主瓣頻率，而不考慮幅值精度，可選主瓣寬度較窄的矩形窗，例如測量物體的自振頻率等；▲如果分析窄帶信號，且有較強的干擾噪聲，則選用旁瓣幅度小的窗函數，如漢寧窗、三角窗等；▲對于隨時間按指數衰減的函數，可采用指數窗來提高信噪比。一、最大熵譜估計的基本原理

●

本節(jié)將依據已經證明的最大熵定理（3-4）闡明這一方法的計算原理。

●最大熵譜分析方法把信息熵的概念引入信號處理中，有時又稱為現代時序譜分析方法。這是一種把自相關函數外推的方法。在分析過程中，沒有固定的窗口函數。在每一步外推自相關函數中，使估計的相關函數包含過程的信息最多，即要求在過程的熵達到最大的條件下，確定未知的自相關函數值，借以達到譜估計的逼真和穩(wěn)定度最好的目的。換句話說，就是采用譜熵為最大的準則來估計功率譜。前面已經證明，N維高斯隨機序列信源的相對熵為【參見式（3-15）】

：|R|(或記為det[R])表示矩陣[R]的行列式。自相關矩陣是