1、1功率譜估計功率譜估計 - -參數估計方法參數估計方法2周期圖法的不足周期圖法的不足v估計方法的方差性能差估計方法的方差性能差 在功率譜密度計算中沒有實現求均值的運算v分辨率低分辨率低 樣本數據x(n)是有限長的，相當于在無限長樣本數據中加載了窗函數(矩形窗、Hanning等)3參數模型功率譜估計參數模型功率譜估計vMA模型模型vAR模型模型vARMA模型模型4平穩隨機信號的參數模型平穩隨機信號的參數模型v如果一個寬平穩隨機信號如果一個寬平穩隨機信號x(n)通過一個線性時通過一個線性時不變系統不變系統(LSI)h(n)，則系統輸出，則系統輸出y(n)也是寬平也是寬平穩隨機過程，并且穩隨機過程，

2、并且y(n)的功率譜密度和的功率譜密度和x(n)的功的功率譜密度滿足下式：率譜密度滿足下式：v其中其中Pyy、Pxx分別為系統輸出、輸入的功率譜密分別為系統輸出、輸入的功率譜密度，而度，而H(w)為系統脈沖響應的傅立葉變換。為系統脈沖響應的傅立葉變換。2( )( ) |( )|yyxxhPwPwHw5平穩隨機信號的參數模型平穩隨機信號的參數模型v如果系統輸入為白噪聲信號如果系統輸入為白噪聲信號u(n)，其功率譜密，其功率譜密度為常數度為常數2，則輸出信號功率譜密度，則輸出信號功率譜密度Pxx(w)完全完全由系統傳遞函數由系統傳遞函數|H(w)|2決定，因此我們通過對決定，因此我們通過對H(w)

3、進行建模，從而得到輸出信號的功率譜密進行建模，從而得到輸出信號的功率譜密度。度。H(z)u(n)x(n)6平穩隨機信號的參數模型平穩隨機信號的參數模型v在上圖中，輸入在上圖中，輸入u(n)為白噪聲信號，其方差為為白噪聲信號，其方差為2 ，則系統輸出，則系統輸出x(n)的功率譜密度的功率譜密度Pxx(w)為：為：22( )|( )|xxPwH w7平穩隨機信號的參數模型平穩隨機信號的參數模型v因此我們利用確定性系統傳遞函數因此我們利用確定性系統傳遞函數H(z)的特性的特性去表征隨機信號去表征隨機信號x(n)的功率譜密度，稱為參數的功率譜密度，稱為參數模型功率譜估計。模型功率譜估計。v參數模型功率

4、譜估計的步驟：參數模型功率譜估計的步驟： 對H(z)選擇合適的模型：MA模型、AR模型、ARMA模型 根據已知樣本數據x(n)，或者x(n)的自相關函數，確定H(z)的參數 利用H(z)估計x(n)的功率譜。8平穩隨機信號的參數模型平穩隨機信號的參數模型vH(z)的模型：的模型： AR模型：auto-Regressive 此模型只有極點，沒有零點，對應其幅度譜結構存在譜峰11( )1piiiH za z9平穩隨機信號的參數模型平穩隨機信號的參數模型 MA模型：Moving-Average 此模型只有零點，沒有極點，對應幅度譜結構中存在譜谷點。1( )1qiiiH zb z 10平穩隨機信號的參

5、數模型平穩隨機信號的參數模型 ARMA模型： 此模型同時有零點、極點，對應幅度譜結構中存在譜峰、譜谷111( )1qiiipiiib zH za z11系統模型系統模型v對于一階全極點傳遞函數對于一階全極點傳遞函數 傳遞函數所對應的幅度響應實際上是：11( )1H zaz|1|( )|zH zzaza12v當當a013v當當ap，因此我們利用，因此我們利用p個估計的自相關函數，個估計的自相關函數，可以對可以對mp所有的自相關函數所有的自相關函數rxx (m)進行延拓，進行延拓，從而提高了自相關函數窗的長度，增加了功率從而提高了自相關函數窗的長度，增加了功率譜估計的頻域分辨率。譜估計的頻域分辨率

6、。1( )()pxxk xxkrma rmk 44AR模型階數模型階數p的選擇的選擇v如果模型的階數過小，則會增加對功率譜的平如果模型的階數過小，則會增加對功率譜的平滑作用，降低譜的分辨率滑作用，降低譜的分辨率v但如果階數太高，雖然會降低預測誤差的方差，但如果階數太高，雖然會降低預測誤差的方差，但會導致譜峰的分裂，增加估計誤差。但會導致譜峰的分裂，增加估計誤差。v這是由于階數實際上對應于譜結構中的譜峰情這是由于階數實際上對應于譜結構中的譜峰情況。況。45AR模型階數模型階數p的選擇的選擇46AR模型階數模型階數p的選擇的選擇47AR模型階數模型階數p的選擇的選擇v在進行在進行AR譜估計時，首先

7、需要確定階數譜估計時，首先需要確定階數p。p的的選擇可以基于以下三種準則進行。選擇可以基于以下三種準則進行。v最終預測誤差準則最終預測誤差準則(FPE)v其中其中k為階數，為階數，N為樣本數據為樣本數據x(n)的長度，而的長度，而k表示表示k階階AR模型得到的白噪聲方差。模型得到的白噪聲方差。v上式最小值對應的階數為最終選擇的階數。上式最小值對應的階數為最終選擇的階數。( )kNkFPE kNk48AR模型階數模型階數p的選擇的選擇v阿凱克信息論準則阿凱克信息論準則(AIC)v同樣選擇使上式最小的同樣選擇使上式最小的k值作為模型的階數。值作為模型的階數。vAIC準測和準測和FPE準則在樣本數據

8、準則在樣本數據x(n)長度較長時，長度較長時，估計得到的模型階數相似。對于較短的樣本數估計得到的模型階數相似。對于較短的樣本數據，建議使用據，建議使用AIC準則。準則。( )ln2kAIC kNk49AR模型階數模型階數p的選擇的選擇v自回歸傳遞函數準則自回歸傳遞函數準則(CAT)v同樣使得上式最小的同樣使得上式最小的k為模型階數。為模型階數。111( )kiikNiCAT kNN50AR模型參數的求解模型參數的求解v自相關法自相關法v利用利用Yule-Walker方程得到方程得到AR模型參數模型參數ai：221( )|1|ARpjwkkkPwa e12(0)(1)(1)(1)(1)(0)(2

9、)(2)(1)(2)(0)( )xxxxxxxxxxxxxxxxpxxxxxxxxarrrprarrrprarprprrp 51AR模型參數的求解模型參數的求解vYule-Walker方程中的自相關函數方程中的自相關函數rxx(m)為有偏為有偏估計值：估計值：21(0)( )pxxk xxkra rk1 | |01( )( ) ()Nmxxnrmx n x nmN 52AR模型估計功率譜密度模型估計功率譜密度v系數矩陣不僅僅是對稱的，而且沿著和主對角系數矩陣不僅僅是對稱的，而且沿著和主對角線平行的任意一條對角線上的元素都相等，這線平行的任意一條對角線上的元素都相等，這樣的矩陣稱為樣的矩陣稱為T

10、oeplitz矩陣，可以利用矩陣，可以利用Levinson-Durbin遞推算法得到遞推算法得到p個參數個參數ai以及方差以及方差2。53AR模型估計功率譜密度模型估計功率譜密度v如果采用有偏估計得到自相關函數，就可以利如果采用有偏估計得到自相關函數，就可以利用用Levinson-Durbin高效的求解高效的求解AR模型參數，模型參數，并且可以保證求解的系數并且可以保證求解的系數ai在單位圓內，即保證在單位圓內，即保證AR模型的穩定性。這種方法稱為自相關法模型的穩定性。這種方法稱為自相關法v同時自相關法計算的白噪聲信號功率會隨著階同時自相關法計算的白噪聲信號功率會隨著階數的增加而減小或者保持不

11、變。數的增加而減小或者保持不變。54AR模型估計功率譜密度模型估計功率譜密度v但自相關法也存在一定的問題，由于在求解自但自相關法也存在一定的問題，由于在求解自相關函數的時候，進行了矩形加窗處理，降低相關函數的時候，進行了矩形加窗處理，降低了分辨率。了分辨率。v同時當樣本數據長度較短時，估計誤差會比較同時當樣本數據長度較短時，估計誤差會比較大，出現譜峰偏移和譜線分裂。大，出現譜峰偏移和譜線分裂。55AR模型估計功率譜密度模型估計功率譜密度v協方差法協方差法(1)(1,1)(2,1)( ,1)(0,1)(2)(1,2)(2,2)( ,2)(0,2)(1, )(2, )( , )(0, )( )px

12、xxxpxxxxxxxxparrr prarrr prrprpr p prpap 56AR模型估計功率譜密度模型估計功率譜密度v其中的自相關函數為：其中的自相關函數為：v同時白噪聲的方差為：同時白噪聲的方差為：1*1( , )()()Nxn pr k lx nl x nkNp21(0,0)(0, )pxxk xxkra rk57AR模型估計功率譜密度模型估計功率譜密度v計算自相關函數時，樣本數據的取值范圍與自計算自相關函數時，樣本數據的取值范圍與自相關法不同，這樣保證了不對樣本數據進行矩相關法不同，這樣保證了不對樣本數據進行矩形窗的截斷，因此如果樣本函數的長度較短時，形窗的截斷，因此如果樣本函

13、數的長度較短時，可以獲得比自相關法更好的譜分辨率。如果樣可以獲得比自相關法更好的譜分辨率。如果樣本函數的長度遠遠大于階數本函數的長度遠遠大于階數p時，自相關法和協時，自相關法和協方差法的性能是差不多的。方差法的性能是差不多的。v同時協方差法求解的是非同時協方差法求解的是非Toeplitz陣，不能用迭陣，不能用迭代的方法計算，因此運算復雜度較大。同時也代的方法計算，因此運算復雜度較大。同時也不能像自相關法一樣保證不能像自相關法一樣保證AR模型的穩定性。模型的穩定性。58AR模型估計功率譜密度模型估計功率譜密度v修正的協方差法修正的協方差法v與協方差法類似，自相關函數的求解修正為：與協方差法類似，

14、自相關函數的求解修正為：v同時估計的白噪聲方差為：同時估計的白噪聲方差為：1*1*01( , )()()2()() ()Nxn pNpnr k lx nl x nkNpx nk x nl 21(0,0)(0, )pxxk xxkra rk59AR模型估計功率譜密度模型估計功率譜密度v修正的協方差法從線性預測的角度分析，實際修正的協方差法從線性預測的角度分析，實際上是同時進行前向、后向預測，因此其估計譜上是同時進行前向、后向預測，因此其估計譜的分辨率比較高，譜峰的偏移也比較小。的分辨率比較高，譜峰的偏移也比較小。v但缺點同樣是需要求解非但缺點同樣是需要求解非Toeplitz陣，計算比較陣，計算比

15、較復雜。復雜。60AR模型估計功率譜密度模型估計功率譜密度vBurg遞推法遞推法以上提到的自相關法、協方差法和修正的協方差以上提到的自相關法、協方差法和修正的協方差都要估計樣本數據的自相關函數，如果能免去都要估計樣本數據的自相關函數，如果能免去自相關函數的求解，從而直接根據樣本函數得自相關函數的求解，從而直接根據樣本函數得到到AR模型參數模型參數ai，從而可以減少中間步驟，提，從而可以減少中間步驟，提高譜估計的性能。同時采用前向、后向線性預高譜估計的性能。同時采用前向、后向線性預測。測。這種算法對短數據的功率譜估計比自相關函數法這種算法對短數據的功率譜估計比自相關函數法要準確。要準確。61應用

16、應用v針對含噪正弦信號針對含噪正弦信號v數據數據64點，采用分段平均周期圖法和點，采用分段平均周期圖法和AR模型模型1122( )5*sin()5*sin()( )x nwnwnv n62應用應用vAR模型，模型，周期圖法周期圖法0123456-10-505101520253063應用應用v數據長度為數據長度為64，采用分段的周期圖法和，采用分段的周期圖法和AR模型模型( )( )(2)x nw nw n( )22cos2xxPww64應用應用vAR模型，模型，周期圖法周期圖法0123456-15-10-5051065應用應用v數據長度為數據長度為64點，取不同階數點，取不同階數(4、20)的

17、的AR模型模型對譜估計的影響。對譜估計的影響。( )( )0.18* (2)x nv nx n66應用應用01234567-2-1.5-1-0.500.511.522.567MA模型估計功率譜密度模型估計功率譜密度v與與AR模型一樣，首先推導模型一樣，首先推導MA參數參數bi與樣本數與樣本數據據x(n)的正則方程。的正則方程。v首先，首先，MA模型參數為：模型參數為：1( )1qiiiH zb z 68MA模型估計功率譜密度模型估計功率譜密度v輸入白噪聲信號輸入白噪聲信號u(n)、MA模型以及輸出模型以及輸出x(n)之之間為線性卷積的關系：間為線性卷積的關系：v與與AR模型進行相同的分析，得到

18、：模型進行相同的分析，得到：1( )( )()qkkx nu nu nk b200,1,., ( )0q mkm kkxxb bmqrmmq69MA模型估計功率譜密度模型估計功率譜密度vMA模型只有模型只有q個零點，并且注意個零點，并且注意MA的正則方的正則方程，程，MA模型計算的自相關函數的取值范圍為模型計算的自相關函數的取值范圍為-q q，并且類似于，并且類似于MA模型參數模型參數bi的自相關函數的自相關函數v這樣估計的功率譜密度為：這樣估計的功率譜密度為：221( )|1|qjwkkkP wb e70MA模型估計功率譜密度模型估計功率譜密度v注意到，根據計算的自相關函數注意到，根據計算的

19、自相關函數rxx(m)得到的功得到的功率譜為：率譜為：v因此從譜估計的角度，因此從譜估計的角度，MA模型譜估計等效于經模型譜估計等效于經典譜估計中的自相關法，譜估計的分辨率低。典譜估計中的自相關法，譜估計的分辨率低。( )( )qjwmBTxxmqPwrm e220|qjwkkkb e71ARMA模型估計功率譜密度模型估計功率譜密度vARMA模型實際上模型實際上AR模型和模型和MA模型的綜合，模型的綜合，其正則方程為：其正則方程為：v其中其中h(k)為為ai和和bi的函數，因此該方程為非線性的函數，因此該方程為非線性方程，求解較為復雜。方程，求解較為復雜。2101()( )0,1.,( )()

20、pq mk xxkkkxxpk xxka rmkh k bmqrma rmkmq72最大熵譜估計方法最大熵譜估計方法vrxx(k)的最大熵外推法的最大熵外推法v經典譜估計中是零值外推，對于窄帶信號是很經典譜估計中是零值外推，對于窄帶信號是很不精確的，如何對不精確的，如何對rxx(k)進行外推？進行外推？v這里這里re(k)表示自相關函數的外推值表示自相關函數的外推值1(1)| | (1)( )( )( )NjwkjwkxxxxekNkNPwrk er k e73最大熵譜估計方法最大熵譜估計方法v對對re(k)的約束條件是什么？的約束條件是什么？ 保證得到的功率譜密度是實數，并且是非負的。 使隨

21、機信號x(n)的熵最大，等價為使得x(n)盡可能的白化，對功率譜而言，使得估計的功率譜Pxx(w)盡可能平坦。74最大熵譜估計方法最大熵譜估計方法v對于能量有限的信號，具有高斯分布的隨機信對于能量有限的信號，具有高斯分布的隨機信號號x(n)具有最大的熵率，并且具有最大的熵率，并且x(n)是高斯是高斯AR過過程，即程，即x(n)的功率譜是全極點形式的譜結構。的功率譜是全極點形式的譜結構。75最大熵譜估計方法最大熵譜估計方法v根據以上要求外推的自相關函數根據以上要求外推的自相關函數rxx(k)：v而而ap為自相關正則方程的解：為自相關正則方程的解：1( )()0pxxk xxkrka rkll 1

22、1(0)(1)( )1(1)(0)(1)00( )(1)(0)xxxxxxppxxxrrrparrrparprpr 76最大熵譜估計方法最大熵譜估計方法v最大熵譜估計的解釋是：根據給定隨機信號最大熵譜估計的解釋是：根據給定隨機信號x(n)，利用對，利用對x(n)AR模型的限制，對自相關函模型的限制，對自相關函數進行外推，并且假設數進行外推，并且假設x(n)為高斯分布。為高斯分布。v實際上最大熵譜估計與實際上最大熵譜估計與Yule-Walker方法估計功方法估計功率譜是等價的。率譜是等價的。77最大熵譜估計方法最大熵譜估計方法v這里對最大熵譜估計方法的描述，用于解釋該這里對最大熵譜估計方法的描述

23、，用于解釋該方法譜估計的實質問題。方法譜估計的實質問題。v最大熵估計外推對數據強加了一個全極點模型，最大熵估計外推對數據強加了一個全極點模型，因此因此MEM估計是否優于傳統方法，取決于所分估計是否優于傳統方法，取決于所分析的信號類型，以及信號模型逼近析的信號類型，以及信號模型逼近AR過程的程過程的程度度78特征分解法譜估計特征分解法譜估計v對于帶有白噪聲的正弦波組合，由于正弦波之對于帶有白噪聲的正弦波組合，由于正弦波之間是非諧波的關系間是非諧波的關系 ，因此不能用基于傅立葉變，因此不能用基于傅立葉變換的周期圖法進行分析。而特征分解法可以得換的周期圖法進行分析。而特征分解法可以得到比到比AR模型

24、更高的分辨率，特別是信噪比比較模型更高的分辨率，特別是信噪比比較低的時候，譜估計的效果比較理想。低的時候，譜估計的效果比較理想。79特征分解法譜估計特征分解法譜估計v一階諧波過程：一階諧波過程：v其中復指數其中復指數A1=|A1|ej?1 , ?1是均勻分布的隨機變是均勻分布的隨機變量，量，w(n)是方差為是方差為w2 的白噪聲的白噪聲11( )( )jnwx nAew n80特征分解法譜估計特征分解法譜估計v復數隨機信號的自相關矩陣定義為：復數隨機信號的自相關矩陣定義為：*(0)(1)(1)(1)(0)(2)(1)(2)(0)xxxxxxxxxxxxxxxxxxxxrrrMrrrMRrMrM

25、r81特征分解法譜估計特征分解法譜估計vRss的秩為的秩為1.)(|)(2211keAkrwjkwxxxxsswwRRR111111(1)(2)21(1)(2)11|1jwj Mwjwj Mwssj Mwj MweeeeRAee82特征分解法譜估計特征分解法譜估計v定義定義v則則Rss可以用可以用e1表示：表示：v因此因此Rss的非零特征值為的非零特征值為M|A1|21111,.,jwjMw TeeeHsseeAR1121| 22111 1111|HssReAee eM Ae83特征分解法譜估計特征分解法譜估計vRxx和和Rss的特征向量是一致的。的特征向量是一致的。vRxx的特征根是的特征根

26、是Rss的特征根和噪聲方差之和的特征根和噪聲方差之和222()()()xxisswwisswissiw isswiRvRRvRI vR vvv22max1|wM A84特征分解法譜估計特征分解法譜估計v根據根據Rx 的特征值和特征矢量獲得關于的特征值和特征矢量獲得關于x(n)的參的參數：數： Rx進行特征值分解，最大的特征值為M|A1|2+w2 ，其它的特征值均為w2 利用Rx 的特征值求信號功率|A1|2和噪聲方差：2min2maxmin1|wAM85特征分解法譜估計特征分解法譜估計 最大特征值對應的特征矢量為e1，則e1第二個系數為 ejw1, 其頻率即為w186特征分解法譜估計特征分解法

27、譜估計v基于信號自相關矩陣分解的頻率估計算法：將基于信號自相關矩陣分解的頻率估計算法：將樣本空間分為信號子空間和噪聲子空間，然后樣本空間分為信號子空間和噪聲子空間，然后用頻率估計函數估計頻率值。用頻率估計函數估計頻率值。v假設隨機信號假設隨機信號x(n)由由p個復指數信號和白噪聲信個復指數信號和白噪聲信號組成：號組成：v其中其中s(n)為正弦信號，為正弦信號，v(n)為白噪聲為白噪聲1( )( )( )( )ipjnwiix ns nv nAev n87特征分解法譜估計特征分解法譜估計vX(n)信號的自相關函數：信號的自相關函數：v其中其中Pi是功率：是功率： Pi =|Ai|221( )(

28、)ipjkwxxiwirkPek 88特征分解法譜估計特征分解法譜估計v則則x(n)的自相關矩陣為：的自相關矩陣為：(1)2(2)2(1)(2)2iiiiiijwj Mwiwiiiiijwj Mwiiwiiiixxj Mwj MwiiiwiiiPPePePePPeRPePeP89特征分解法譜估計特征分解法譜估計v其中其中ei為：為：2Hxxsswwi i iwiRRRPee2(1)1,.,iiijwj wj MwTieeee90特征分解法譜估計特征分解法譜估計v設設vi是是Rss的特征矢量：的特征矢量：vvi也是也是Rxx的特征矢量，并且的特征矢量，并且Rxx的特征根是的特征根是Rss特特征根

29、和噪聲方差之和：征根和噪聲方差之和：222()()xx isswwissiwissiiw iiwiR vRRvR vIvvvv2xsiiw91特征分解法譜估計特征分解法譜估計vRss的秩為的秩為p，因此，因此Rss有有p個非零特征根，因此個非零特征根，因此Rxx特征根分為特征根分為p個大于個大于w2 的特征根，和的特征根，和M-p個個為為w2 的特征根的特征根。v對應于特征根的分類，特征矢量也分為兩類。對應于特征根的分類，特征矢量也分為兩類。實際上，大于實際上，大于w2的特征根和特征矢量對應信號的特征根和特征矢量對應信號子空間，而等于子空間，而等于w2的特征根和特征矢量對應噪的特征根和特征矢量對應噪聲子空間。聲子空間。v這樣可以根據自相關函數這樣可以根據自相關函數Rxx的特征根求解，最的特征根求解，最小的特征根就是白噪聲信號的方差。小的特征根就是白噪聲信號的方差。92特征分解法譜估計特征分解法譜估計v如果自相關矩陣如果自相關矩陣Rxx的維數的維數M=p+1，則，則p個大于個大于w2 的特征根，的特征根，1個等于個等于w2 的特征根，因此信號的特征根，因此信號子空間的特征矢量為子空間的特征矢量為p個，噪聲子空間的特征個，噪聲子空間的特征矢量為矢量為1個，定