7-醫(yī)學信號處理現(xiàn)代譜估計.ppt

1、1,第七章功率譜估計的現(xiàn)代方法 現(xiàn) 代 譜 估 計,2,經(jīng)典譜估計以傅立葉變換為基礎，具有計算效率高的優(yōu)點，但是由于將未觀測數(shù)據(jù)認為0和數(shù)據(jù)加窗，而具有頻率分辨率低、旁瓣泄漏等嚴重的缺陷。 現(xiàn)代譜估計與經(jīng)典譜估計不同，它以參數(shù)模型為基礎，能夠得到小方差和高分辨率，特別是數(shù)據(jù)長度很短的情況，更具優(yōu)勢。,7.1 概述,3,現(xiàn)代譜估計法的基本思想：,處理步驟： 1 確定或選擇一個合適的模型依賴于對所研究隨機過程進行理論分析和實驗研究； 2 根據(jù)觀測數(shù)據(jù)估計模型參數(shù)涉及各種算法的研究； 3 由模型參數(shù)計算功率譜。 關鍵 1、模型選擇問題（AR， MA ，ARMA） 2、參數(shù)確定方法（導致產(chǎn)生了各種算法

2、）,4,7.2 自回歸模型(AR)譜估計,數(shù)字系統(tǒng)的數(shù)學模型:有理分式傳遞函數(shù)的模型如下圖：,式中ak為自回歸系數(shù)，稱為AR系數(shù)；bk為滑動平均系數(shù)，稱為MA系數(shù)。,模型傳遞函數(shù)為：,5,有理分式傳遞函數(shù)的模型的差分方程為：,令a0=1有：,w(n)為高斯白噪聲，,6,求功率譜的實質(zhì)變?yōu)榇_定系統(tǒng)參數(shù)的問題,模型的功率譜密度：即系統(tǒng)輸出功率譜和輸入功率譜之間關系為（假定h(n)為實序列）：,7,如果除b0外其它的MA系數(shù)都等于0，即,AR（p）模型,全極點模型,8,自回歸模型,9,如果除a0外其它的AR系數(shù)都等于0，即,全零點模型,MA（q）模型,10,如果除a0=1和b0=1外其它的AR系數(shù)和

3、MA系數(shù)都不全等于0，即,稱為ARMA（p,q）模型，即極點-零點模型。,11,到底選擇什么模型？三種模型之間關系如何？Word分解定理,Wold分解定理：任何一個有限方差的平穩(wěn)ARMA過程可以分為完全隨機的部分和確定的部分。 推論：任何有限方差的ARMA或MA平穩(wěn)過程可以用一個無限階的AR模型表示；同樣，任何ARMA或AR模型可以用一個無限階的MA模型表示。因 此，如果在這三個模型中選了一個與信號不匹配的模型，利用高的階數(shù)仍然可以得到對信號的好的逼近。,12,結論： 由于對AR模型參數(shù)的估計，得到的是線性方程。故AR模型比ARMA以及MA模型有計算上的優(yōu)點，即只需解一組線性方程，而ARMA或

4、MA模型一般需要解一組非線性方程。同時，實際的物理系統(tǒng)往往是全極點系統(tǒng)。AR模型得到了深入的研究和廣泛的應用。,13,已知：自相關函數(shù),要求： AR模型的階數(shù)p，以及p個AR 參數(shù)a(i),激勵源方差,Yule-Walker方程,7.3 AR模型的Yule-Walker方程,14,7.3.1 Yule-Walker方程的推導,1. 對 進行求逆z變換 2. 直接由模型差分方程推導，把模型的差分方程代入x(n)的自相關函數(shù),15,如何根據(jù)自相關函數(shù)確定系統(tǒng)參數(shù),16,可見，AR模型輸出信號的自相關函數(shù)具有遞推性質(zhì)，即：,Yule-Walker方程（Y-W方程）,17,選擇m0的前P個方程并寫成單

5、一正規(guī)矩陣的形式為：,以上利用了自相關函數(shù)的偶對稱性。Y-W方程表明：只要已知輸出平穩(wěn)隨機信號的自相關函數(shù)，就能求出AR模型中的參數(shù)ak，并且需要的觀測數(shù)據(jù)較少。,18,AR模型譜估計,N個樣值 x(0),x(1)x(N),自相關函數(shù) R(0),R(1).R(N),AR模型參數(shù)和 a1,a2，ap 激勵源方差,功率譜密度,Y-W方程,19,Yule-Walker方程的求解,1、采用高斯消元法，解線性方程組常用方法，運算量數(shù)量級為p的三次方。 2、用Levinson-Durbin算法，Y-W方程的高效解法，即按階次進行遞推運算量數(shù)量級為p的二次方。,7.3.2 Levinson-Durbin算法

6、,20,Levinson-Durbin遞推算法： 算法的關鍵就是要推導出由第K階AR模型的參數(shù)計算第k+1階AR模型AR(k+1)參數(shù)的迭代計算公式。,首先以AR(0)和AR(1)模型參數(shù)作為初始條件，計算AR(2)模型參數(shù)，然后根據(jù)這些參數(shù)計算AR(3)模型參數(shù)，等等，一直到計算出AR(p)模型參數(shù)為止。,21,22,23,遞推公式為：,其中akk稱為反射系數(shù),將所估計的模型參數(shù)代入即可計算功率譜估計值：,24,AR模型參數(shù)和 a1,a2，ap 激勵源方差,功率譜密度,AR模型譜估計,25,給定初始值和AR模型的階數(shù)p，可按照L-D算法流程進行估計，流程終止規(guī)則為 或,MATLAB里有專門實

7、現(xiàn)L-D算法的函數(shù)可估計AR模型參數(shù)： a E=aryule(x,p),a為模型參數(shù)，E為噪聲方差。 分析： AR模型的穩(wěn)定性；L-D算法的收斂性。,26,AR模型譜估計的L-D算法流程,、給定N個觀察數(shù)據(jù)xN(n),n=0,1,N-1; 、由xN(n)估計自相關函數(shù)值 ，m=0,1,p; 、利用L-D遞推算法，根據(jù) 計算AR模型參數(shù)的估計值。首先令p=1，按下式計算a11和,然后，使p=p+1,按下式計算app,api, 、重復以上遞推過程，直到滿足p=m或者 。 、代入 計算公式估計功率譜。,27,例7-1、已知實數(shù)據(jù)序列,的自相關為：,用Levinson-Durbin遞推算法求AR模型的

8、參量：,解：,28,29,一、AR模型的穩(wěn)定性具有下面性質(zhì)：,H(z)全部的極點在單位圓內(nèi) 自相關矩陣正定 激勵信號方差隨階次增加而遞減,7.3.3 確定AR模型的階,30,階太低，功率譜平滑的太厲害，平滑后的譜分辨不出真實譜中的兩個峰； 階太高，可以提高譜估計的分辨率，但會出現(xiàn)許多虛假譜峰。,真實譜,虛假譜峰,二、有關AR模型的階的問題：,31,所以，估計一個AR(p）過程，選取AR(k）階數(shù)要求 kp , 但k不能太大。如果估計精確的話， kp時，AR(p)模型參數(shù)估計為：,AR模型譜估計方法，既要估計模型參數(shù)，又要估計模型的階，在這樣復雜的情況下，如何評價各種譜估計的性能，目前尚無定論。

9、,32,三、確定AR模型的階的方法 一般的觀察方法，簡單而直觀,不斷增加階數(shù)，觀察預測誤差功率，下降到最小時，對應的階選為模型的階； 不斷增加階數(shù)，觀察各階模型預測誤差序列的周期圖，最接近平坦時對應于最佳的階；,33,1、FPE(最終預測誤差）,N為觀測數(shù)據(jù)長度， 為擬合殘差方差，隨階增加而減小，F(xiàn)PE的最小值對應的階數(shù)為最后確定的階。,四、確定AR模型的階的方法根據(jù)誤差準則確定,34,2、Akaike（AIC)信息準則,AIC(i)= 最小所對應的階。,i為模型的階， 為模型誤差，隨著階的增加而減小，而式中第二項隨階次增加而增加。 AIC定義式有一個最小值。 適用于AR模型。,35,此外，還

10、有CAT等準則 。,通過實驗發(fā)現(xiàn)：在將這些準則用于估計AR模型的階，對于實際數(shù)據(jù)，所得到的譜估計結果常常無太大區(qū)別。 對于短數(shù)據(jù)，以上準則都不理想。 在實際應用中，應該參照實驗結果對模型的階加以適當調(diào)整。,36,7.4 線性預測譜估計,假設x(n)是一個N階AR過程，現(xiàn)在時刻x(n)的值 可以由過去N個時刻的取樣值的加權來預測，加權系數(shù)為-ak，那么,N階線性預測器：可看作用序列x(n-N),x(n-N-1) , ,x(n-1)激勵一個沖擊響應為-ak的線性時不變系統(tǒng)的輸出值。,x(n-N),x(n-N-1) , ,x(n-1),-ak,預測誤差為：,預測誤差功率為：,38,確定系數(shù)ak的一個

11、原則是使預測誤差功率最小。根據(jù)這一原則推導出的預測器系數(shù)-ak與x(n)的自相關序列Rxx(m)之間的關系為：,將兩個關系式寫成矩陣展開式分別為：,40,將(1)和(2)兩個關系式合并為一個式子：,41,將(3)寫成矩陣展開形式為：,可以看出：N階線性預測器的系數(shù)ak與AR模型中的AR系數(shù)相等，預測誤差概率最小值Pmin與AR模型中的輸入噪聲方差 相等。所以，線性預測譜估計與AR譜估計是等效的。,42,熵是信息量的一種量度，也是不確定性的一種量度。信息量與事件發(fā)生概率之間有類似于反比例的關系，信息量與概率之間存在對數(shù)關系。復合事件的信息量等于各獨立事件信息量之和。對于事件A有：,7.5 最大熵

12、譜估計（MESE） Maximum Entropy Spectral Estimation）,7.5.1 按最大熵譜外推自相關函數(shù),43,N個符號組成信號系統(tǒng)傳遞消息，每個符號出現(xiàn)的概率為pi，接收到第i個符號的信息量為I(i)，消息中總的平均信息量為：,這個平均信息量稱為具有符號i和概率pi的信號系統(tǒng)的熵。對于隨機過程，應該用聯(lián)合概率密度函數(shù)來定義熵： x0,x1,xN為隨機過程的N+1個取樣值。,44,對于零均值的高斯平穩(wěn)隨機過程則有：,其中X=x0,x1,xN為由N+1個取樣值構成行矩陣，XH為X的共軛轉(zhuǎn)置矩陣。,45,detR(N)是行列式的值。于是有下列式子：,46,均值為0的高斯平

13、穩(wěn)隨機過程的熵的表達式，它是R(N)的函數(shù)。,47,最大熵譜估計：為了使得H取得最大值，應當使detR(N)取最大值。根據(jù)外推或預測方法，求出使detR(N+1) 取最大值的Rxx(N+1)：,由,得到：,48,結論：,上述方程為Rxx(N+1)的一元一次方程，可由已知或估計的N+1個自相關值Rxx(0)、 Rxx(1)、 Rxx(N)求出Rxx(N+1)。以此類推。 可以證明這種按最大熵外推自相關函數(shù)的結果與AR模型是等價的。所以，上式實質(zhì)為Yulerwalker方程。,49,實際上，假定在線性預測譜估計中，用外推的方法得到了Rxx(N+1)，即：,7.5.2 MESE與AR譜估計等效,50

14、,使聯(lián)立方程有非零解的充分必要條件是系數(shù)行列式等于0：,結論 因此：最大熵譜估計與AR模型譜估計和線性預測譜估計是等效的 。 此外，還可以證明：AR譜估計等效于最佳白化處理。,51,7.6.1 預測誤差格型濾波器,7.6 預測誤差格型濾波器及伯格(Burg)遞推算法,已知n個觀測數(shù)據(jù)x(1),x(2),x(n-1)，利用p階線性預測濾波器估計x(n)為 估計誤差為：,52,代入有：,53,上式中,因此可得：,54,格型前向預測誤差濾波器,55,格型預測誤差濾波器傳遞函數(shù)為,同時有：,56,也就是相當于：,57,后向預測誤差,58,格型后向預測誤差濾波器傳遞函數(shù)為,又因：,因此：,59,7.6.

15、2 Burg遞推算法Kp的確定,根據(jù)信號的有限個取樣值估計AR模型參數(shù)的方法自相關法、協(xié)方差法和Burg遞推法。 自相關法和協(xié)方差法都是直接估計AR參數(shù)，而Burg法是先估計反射系數(shù)，然后利用Levinson-Durbin遞推算法由反射系數(shù)求得AR參數(shù)。 Burg法可以保證：,60,算法準則是前向均方誤差和后向均方誤差之和最小。,如果用前向預測方法以均方誤差最小為準則確定Kp，用 表示為：,令,61,如果用后向預測方法以均方誤差最小為準則確定Kp，用 表示為：,令,62,Burg算法準則是前向均方誤差和后向均方誤差之和最小，令,63,于是得到：,對于平穩(wěn)隨機序列，集合平均可用時間均值來代替，故

16、上式為：,通過觀察發(fā)現(xiàn)，下式成立：,64,對于長度為N的有限長序列，有下面一系列遞推關系式。當p=1時：,65,于是有：,66,于是,其中,67,繼續(xù)代入求得K3，b3(n)，e3(n)，以此類推。便可以由x(n)求得各階Kp以及前向與后向誤差及其各個akk。,68,Burg法估計AR(p)模型參數(shù)的具體步驟為： 1、確定初始條件： 2、按照公式 計算Kp。 3、按照公式 計算ep(n)和bp(n)。 4、計算均方誤差： 5、p=p+1。 6、重復第25步，直至滿足條件為止。,69,例5-2、設N=5的數(shù)據(jù)記錄為x(0)=1,x(1)=2,x(2)=3,x(3)=4,x(4)=5,AR模型的階

17、次p=3,試用相關函數(shù)法確定AR參量及預測值 .,解：先由數(shù)據(jù)求自相關函數(shù)式：,70,用Levinson-Durbin遞推算法求AR模型的參量分別是：,71,根據(jù)所得的AR(3)參量，預測值：,若使用的是二階線性預測器，有例51所得的結果，則,可分別由前向與后向預測得到如下：,72,例5-3、設仍利用例52中的記錄數(shù)據(jù)，試用伯格法求AR(2)的參量。 解：用上述遞推公式，i=1時：,e1(n)和b1(n),73,p=2時：,若使用此二階線性預測，可得：,74,算法比較 Levinson-Durbin Burg算法 真正樣值,25.9090 3.0650,25.5403 0.17415,1.27

18、00 0.8549 1.0000,2.8983 4.5825 5.0000,顯然，伯格算法要比萊文森德賓算法優(yōu)越得多。短數(shù)據(jù)！,75,比較Welch方法和Burg方法在噪聲信號的功率譜估計中的效果。 為高斯型白噪聲。,例子：現(xiàn)代譜估計和經(jīng)典譜估計方法的比較。,利用MATLAB中的pburg和pwelch函數(shù)分別用Burg方法和Welch方法對上述噪聲信號進行功率譜估計并比較結果。,76,function =burgwelchpsd() fs=1000; t=0:1/fs:1; xn=sin(2*pi*80*t)+2*sin(2*pi*140*t)+randn(size(t); plot(xn)

19、; pw,f=pwelch(xn,fs,twosided); pb1,f=pburg(xn,17,fs,twosided); pb2,f=pburg(xn,13,fs,twosided); figure w=10*log10(pw) 10*log10(pb1) 10*log10(pb2); plot(f,w); grid xlabel(frequency(Hz); ylabel(amplitude(dB); axis(0 200 -50 0); legend(welch方法,Burg方法高階,Burg方法低階);,77,信號：,78,很明顯，Burg方法比Welch方法更光滑。但是，當AR模型

20、階數(shù)降低時，譜峰的頻移越來越明顯，頻率分辨率降低。,79,7.7 AR模型譜估計存在的問題,7.7.1 譜線分裂,由正弦信號疊加噪聲構成的隨機信號，在下列四種情況下容易出現(xiàn)譜線分裂的現(xiàn)象，即譜線頻率偏移或出現(xiàn)兩個靠得很近的譜峰。 高信噪比；正弦信號分量的初始相位是/4的奇數(shù)倍； 數(shù)據(jù)長度為正弦分量的1/4周期的奇數(shù)倍； AR模型參數(shù)的數(shù)目與數(shù)據(jù)的個數(shù)相比的百分比較大，即二者大小可比擬。,80,對于Burg算法，譜線分裂是由于第一個反射系數(shù)K的計算誤差引起的，K1的估計并沒有使預測誤差功率最小。 改善措施： 1、用解析信號代替實值信號，克服信號相位的影響； 2、調(diào)整修正反射系數(shù)，以使預測誤差功率達到最小。,81,7.7.2 附加噪聲使分