中科院研究生院-現(xiàn)代數(shù)字信號處理課件Chapter-4

1、第四章 功 率 譜 估 計 4.1 引言 4.2 經(jīng)典譜估計 4.3 現(xiàn)代譜估計中的參數(shù)建模 4.4 AR模型譜估計方法 4.5 其他譜估計方法4.1 引 言功率譜定義估計質(zhì)量評價功率譜估計的方法功率譜估計的應(yīng)用本章討論的主要內(nèi)容1、功率譜的定義信號的功率譜和其自相關(guān)函數(shù)服從一對傅里葉變換關(guān)系 對于平穩(wěn)隨機(jī)信號，服從各態(tài)歷經(jīng)定理，集合平均可以用時間平均代替令l=n+m, 那么 2、估計質(zhì)量評價無偏性：一致性：3、功率譜估計的方法經(jīng)典譜估計方法間接方法：BT法直接方法：周期圖法現(xiàn)代譜估計方法參數(shù)法：ARMA模型法AR模型、MA模型、ARMA模型非參數(shù)法：諧波分解法、多分量法1、經(jīng)典譜估計方法BT

2、法：先按照有限個觀測數(shù)據(jù)估計自相關(guān)函數(shù)，再計算功率譜；周期圖法：直接對觀測數(shù)據(jù)進(jìn)行處理，計算功率譜。經(jīng)典譜估計方法的特點(diǎn)：都采用傅立葉變換方法，物理概念比較清楚；頻率分辨率低；估計量的方差和分辨率是一對矛盾。2、現(xiàn)代譜估計方法以信號模型為根底，估計功率譜的問題轉(zhuǎn)化成由觀測數(shù)據(jù)估計信號模型參數(shù)的問題?，F(xiàn)代譜估計方法的特點(diǎn)：頻率分辨率較經(jīng)典法高；缺乏如何選擇信號模型的理論指導(dǎo)。4、功率譜估計的應(yīng)用在信號處理的許多場所，要求預(yù)先知道信號的功率譜密度(或自相關(guān)函數(shù))；常常利用功率譜估計來得到線性系統(tǒng)的參數(shù)估計；從寬帶噪聲中檢測窄帶信號。 5、本章討論的主要內(nèi)容主要內(nèi)容：BT法、周期圖法、改進(jìn)的周期圖法

3、、AR模型法、最大熵譜估計法、基于矩陣特征分解的方法。分析方法：介紹各種估計方法的原理，根據(jù)估計質(zhì)量評價準(zhǔn)那么，分析討論其估計性能。4.2 經(jīng) 典 譜 估 計BT法周期圖法改進(jìn)的周期圖法4.2.1 BT法BT法是先估計自相關(guān)函數(shù), 然后進(jìn)行傅里葉變換得到功率譜。有偏自相關(guān)函數(shù)估計的誤差相對較小，是一種漸近一致估計： BT法的加權(quán)協(xié)方差譜估計式中 -(M-1)m(M-1) 其它 , MN 窗函數(shù)w(m)的傅里葉變換必須是非負(fù)的。4.2.2 周期圖法 周期圖法的定義如下: 1. 周期圖與BT法的等價關(guān)系 令 m=k-n, 即k=m+n，那么 利用有偏自相關(guān)函數(shù)的BT法和周期圖法是等價的。 2. 周

4、期圖法譜估計質(zhì)量分析 1) 周期圖的偏移 式中 上式在頻域表示為： 式中 周期圖的統(tǒng)計平均值等于它的真值卷積三角譜窗函數(shù)，因此周期圖是有偏估計，但當(dāng)N時，wB(m)1，三角譜窗函數(shù)趨近于函數(shù)，周期圖的統(tǒng)計平均值趨于它的真值，因此周期圖屬于漸近無偏估計。 2) 周期圖的方差 為分析簡單起見， 假設(shè)x(n)是實(shí)的零均值的正態(tài)白噪聲信號，方差是x2，即功率譜是常數(shù)x2 ，其周期圖用IN()表示，N表示觀測數(shù)據(jù)的長度。 周期圖的均值式中 這里由于對信號作了實(shí)白噪聲的假設(shè)，才有無偏估計的結(jié)果。利用正態(tài)白噪聲、多元正態(tài)隨機(jī)變量的多階矩公式，有 周期圖的均方值將上式代入周期圖的均方值公式中， 得到 將=1=

5、2代入上式，得到 信號的功率譜真值是2x,說明周期圖的方差很大，周期圖的均方誤差也是非常大。 用這種方法估計的功率譜在2x附近起伏很大，故周期圖是非一致估計，是一種很差的功率譜估計方法。 圖 4.2.2 白噪聲的周期圖 4.2.3 經(jīng)典譜估計方法改進(jìn)Bartlett平均周期圖法窗口處理法平均周期圖Welch法修正的周期圖求平均法存在問題：BT法和周期圖法估計功率譜都不是一致估計，頻率分辨率低。解決方法：對周期圖進(jìn)行修正，使其滿足一致估計條件。可以采用平滑處理的方法，使其方差減小。1. Bartlett平均周期圖法 主要思想：對序列x(n)進(jìn)行L次獨(dú)立觀測或?qū)⑵浞殖蒐段，計算每組觀測數(shù)據(jù)的周期圖

6、，再將L個周期圖加和后求平均。 假設(shè)隨機(jī)信號x(n)的觀測數(shù)據(jù)區(qū)間為：0nM-1，共進(jìn)行了L次獨(dú)立觀測，得到L組記錄數(shù)據(jù)，每一組記錄數(shù)據(jù)用xi(n), i=1, 2, 3, ,L表示； 或?qū)﹂L為N的數(shù)據(jù)x(n)分成L段，每段有M個數(shù)據(jù)，N=LM，第i段數(shù)據(jù)表示為xi(n)= x(n+iM-M)。 第i組的周期圖用下式表示： 估計方法： 將得到的L個周期圖進(jìn)行平均，作為信號x(n)的功率譜估計， 公式如下： 偏移分析： 估計效果分析： 平均周期圖仍然是有偏估計，偏移和每一段的數(shù)據(jù)個數(shù)M有關(guān)； 偏移的大小反映分辨率的上下。 方差分析： 平均周期圖的估計方差是周期圖的方差的1/L，L越大方差越小，功

7、率譜越平滑；相應(yīng)的，M越小，偏移越大，分辨率越低；估計的均方誤差也減少； 以分辨率的降低換取了估計方差的減少，估計量的方差和分辨率是一對矛盾。圖 4.2.3 平均周期圖法 2、窗口處理法平均周期圖主要思想：用一適當(dāng)?shù)墓β首V窗函數(shù)W(ej)與周期圖進(jìn)行卷積，來到達(dá)使周期圖平滑的目的的。 式中 -(M-1)nM-1 估計方法：那么 周期圖的窗函數(shù)法就是前面提到的BT法的加權(quán)協(xié)方差譜估計。又 偏移分析： 估計效果分析：可得 周期圖的窗函數(shù)法仍然是有偏估計, 其偏移和wB(m)、w(m)兩個窗函數(shù)有關(guān)。 如果w(m)窗的寬度比較窄，M比N小得多，這樣|m|p 4.4 AR譜估計的方法AR譜估計方法可歸

8、結(jié)為求解AR模型系數(shù)或線性預(yù)測器系數(shù)的問題。AR模型參數(shù)估計方法：信號預(yù)測誤差最小原那么或預(yù)測誤差功率最小自相關(guān)法Levison遞推法Burg法協(xié)方差法修正協(xié)方差法前后向線性預(yù)測最小二乘法最大熵原那么最大熵譜估計方法1、 自相關(guān)法列文森Levinson遞推 估計方法：自相關(guān)法的出發(fā)點(diǎn)是選擇AR模型的參數(shù)使預(yù)測誤差功率最??；采用Levison-Durbin遞推方法求解Yule-Walker方程得到AR模型參數(shù)。 預(yù)測誤差功率為 假設(shè)信號x(n)的數(shù)據(jù)區(qū)間在0nN-1范圍，有P個預(yù)測系數(shù)，N個數(shù)據(jù)經(jīng)過沖激響應(yīng)為api(i=0,1, 2, , P)的濾波器， 輸出預(yù)測誤差e(n)的長度為N+P， 因

9、此應(yīng)用下式計算: 預(yù)測誤差功率最小，得到 采用Levinson-Durbin遞推法求解Yule-Walker方程： 由k=1開始遞推，遞推到k=p，依次得到a11,21，a21,a22, 22，ap1,ap2,app,2p。 AR模型的各個系數(shù)以及模型輸入白噪聲方差求出后， 信號功率譜用下式計算： 圖 4.5.1 利用列文森遞推法計算功率譜的流程圖 性能分析：該方法需要基于有限的觀測數(shù)據(jù)估計自相關(guān)序列，當(dāng)數(shù)據(jù)長度較短時，估計誤差會比較大，AR參數(shù)的計算就會引入很大的誤差。從而導(dǎo)致功率譜估計出現(xiàn)譜線分裂與譜峰頻率偏移等現(xiàn)象。2、 伯格Burg遞推法 估計方法：直接由時間序列計算AR模型參數(shù)的方法

10、，求前、后向預(yù)測誤差平均功率最小時的反射系數(shù)kp，進(jìn)而求AR模型參數(shù)ak和2w。 設(shè)信號x(n)觀測數(shù)據(jù)區(qū)間為：0nN-1，前向、后向預(yù)測誤差功率分別用p,e和p,b表示，預(yù)測誤差平均功率用p為 其中，前向、后向預(yù)測誤差公式分別為 求預(yù)測誤差平均功率p最小時的反射系數(shù)kp，令 基于反射系數(shù)kp，由Levinson-Durbin遞推關(guān)系求AR模型參數(shù)ak和2w，進(jìn)而求得功率譜Pxx圖 4.5.2 伯格遞推法流程圖 性能分析：該方法防止了采用有限數(shù)據(jù)估計自相關(guān)函數(shù)的計算，適合短序列參數(shù)估計，克服了L-D遞推中的某些缺點(diǎn)，計算量小。但對正弦信號的譜估計，仍存在某些譜線分裂與頻率偏移現(xiàn)象。3、 協(xié)方差

11、法與修正協(xié)方差法1. 協(xié)方差法 估計方法：利用使預(yù)測誤差功率最小的方法求模型參數(shù) 該公式中使用的觀測數(shù)據(jù)均已得到，不需要在數(shù)據(jù)兩端補(bǔ)充零點(diǎn)， 因此比較自相關(guān)法去掉了加窗處理的不合理假設(shè)。式中 即通過求解以下方程組求apk 性能分析：適用于非平穩(wěn)信號；一些實(shí)驗結(jié)果說明它的分辨率優(yōu)于自相關(guān)法，另外對于純粹弦信號數(shù)據(jù)，可以有效地估計正弦信號的頻率。 2. 修正協(xié)方差法前后向線性預(yù)測最小二乘法： 估計方法：修正協(xié)方差法使用前向和后向預(yù)測誤差平均值最小的方法， 估計AR模型的參數(shù)，進(jìn)而估計信號的功率譜。 前向和后向預(yù)測誤差功率pe、pb分別用下式表示： 預(yù)測誤差平均功率最小經(jīng)過簡化，得到 令 將上式寫成

12、矩陣形式求apk ： 性能分析：該方法去掉了Burg法所用的Levinson的約束條件，估計得到的譜在譜線分裂和頻率偏移時較Burg法有較大改善；該方法也適用于非平穩(wěn)信號。 幾種方法的比較：自相關(guān)法可以用Levinson遞推算法，運(yùn)算量小，但分辨率受窗長度的限制；協(xié)方差法，去除了自相關(guān)法加窗處理的不合理假設(shè)，分辨率高，運(yùn)算量較大；修正協(xié)方差法，分辨率高，在譜線分裂和偏移上較Burg法有較大改善，運(yùn)算量大；Burg算法，可用改進(jìn)的Levinson遞推算法，分辨率高，但對正弦信號存在譜線分裂和偏移現(xiàn)象。 例 信號的四個觀察數(shù)據(jù)為x(n)=x(0), x(1), x(2), x(3)=2, 4, 1

13、, 3， 分別用自相關(guān)法和協(xié)方差法估計AR1模型參數(shù)。 解1 自相關(guān)法： (2) 協(xié)方差法: 圖 4.5.3 AR模型階次太小時的平滑作用4、 關(guān)于AR模型階次的選擇 如果是純P階AR信號，應(yīng)選擇模型階次k P。 如果選擇模型階次kP時，將產(chǎn)生對譜的平滑作用，降低譜的分辨率。如果選擇kP，且假定觀測的數(shù)據(jù)沒有誤差沒有干擾，估計的參數(shù)應(yīng)是： 對于白噪中的AR信號，其階次的選擇應(yīng)折衷考慮。如選擇AR模型，其階次應(yīng)加大，較低的階次會使譜估計產(chǎn)生偏移， 降低分辨率。信噪比愈低，平滑作用愈嚴(yán)重，愈需要高的階次， 因此信噪比低應(yīng)選高的階次。階次愈高，分辨率愈高；但階次太高，會使估計誤差加大，譜峰分裂。最終

14、預(yù)測誤差FPE準(zhǔn)那么 阿凱克信息論準(zhǔn)那么自回歸傳遞函數(shù)準(zhǔn)那么CAT 4.5 其他譜估計方法估計思想：采用最大熵原那么，外推自相關(guān)函數(shù)方法估計信號功率譜。它基于將的有限長度自相關(guān)序列以外的數(shù)據(jù)用外推的方法求得，而不是把它們當(dāng)作是零。 4.5.1 最大熵譜估計 1. 利用最大熵的原那么外推自相關(guān)函數(shù) 按照Shannon對熵的定義， 當(dāng)隨機(jī)變量X取離散值時，熵的定義為 式中pi是出現(xiàn)狀態(tài)i的概率。當(dāng)X取連續(xù)值時，熵的定義為 式中, p(x)是X的概率密度函數(shù)， 假設(shè)x(n)是零均值正態(tài)分布的平穩(wěn)隨機(jī)序列，它的N維高斯概率密度函數(shù)為 式中 先討論一維高斯分布的信號的熵，然后推廣到N維。 同理可求得N維

15、高斯分布信號的熵為式中det(Rxx(N)表示矩陣Rxx(N)的行列式，由上式說明為使熵最大，要求det(Rxx(N)最大。 用最大熵方法外推rxx(N+1)： 設(shè)rxx(N+1)是信號自相關(guān)函數(shù)的第N+2個值，根據(jù)自相關(guān)函數(shù)的性質(zhì)，由N+2個自相關(guān)函數(shù)組成的矩陣為 為選擇rxx(N+1)使det(Rxx(N+1)最大， 解以下方程： 用數(shù)學(xué)歸納法，得到 上式是rxx(N+1)的一次函數(shù)，可以解出rxx(N+1)。以此類推，可推出任意多個其它自相關(guān)函數(shù)值，而不必假設(shè)它們?yōu)榱悖?這就是最大熵譜估計的根本思想。2. 最大熵譜估計與AR模型譜估計的等價性 AR模型信號自相關(guān)函數(shù)與模型參數(shù)服從Yule

16、-Walker方程m1 m=0 將m1的情況寫成： 如果從N個線性方程中解得的N個AR參數(shù)a1, a2, aN值，代入上式并將其整理成行列式的形式，即可得 這就證明了信號為高斯分布時，AR模型功率譜估計和最大熵譜估計的等價性。 通過解Yule-Walker方程，解出模型參數(shù)，最大熵譜估計用下式計算信號功率譜： 4.5.2 特征分解法譜估計該方法特別適用于對多正弦加白噪聲序列進(jìn)行譜分析，可以得到比AR模型法更高的分辨率和更準(zhǔn)確的頻率估計，尤其在信噪比低時更有效。估計思想：將白噪聲加正弦波作為一特殊的ARMA模型，用特征方程求它的參數(shù)，計算出正弦波的頻率。1 正弦波用退化AR模型表示 設(shè)P個實(shí)正弦

17、波組成的信號用下式表示： 式中，初相位i是在區(qū)間-，均勻分布的隨機(jī)變量，令x(n)=sin(n+), 那么上式變?yōu)?這樣上式的特征多項式為 由下面三角恒等式：- Z變換 兩個根分別是:z1=ej,z2=e-j,它們共軛成對，且模為1，由這兩個根可以確定正弦波的頻率。 比照AR模型的系統(tǒng)函數(shù)， 可以把正弦波信號用一個特殊的AR2模型表示。該AR模型的鼓勵白噪聲方差趨于0，極點(diǎn)趨于單位圓。通常稱為退化的AR模型。x(n)+a1x(n-1)+a2x(n-2)=w(n) 對于P個實(shí)正弦波, 特征多項式是 上式是z-1的2P階多項式，可以表示為 考慮到根共軛成對，也可表示為 P個正弦波組合的模型用2P階差分方程描述 2 白噪聲中正弦波組合用一特殊