北京交通大學(xué)(數(shù)字信號(hào)處理研究生課程)ch6_3現(xiàn)代功率譜估計(jì)

1、近代數(shù)字信號(hào)處理近代數(shù)字信號(hào)處理(Advanced Digital Signal Processing) 電子信息工程學(xué)院電子信息工程學(xué)院信號(hào)與圖像處理研究室1. 方差性能不好，不是方差性能不好，不是Px(W W)的一致估計(jì)的一致估計(jì)2. 平滑周期圖和平均周期圖改善了周期圖的方差平滑周期圖和平均周期圖改善了周期圖的方差性能，但卻降低了譜分辨率和增大了偏差。性能，但卻降低了譜分辨率和增大了偏差。3. 可能使短序列的功率譜估計(jì)出現(xiàn)錯(cuò)誤的結(jié)果可能使短序列的功率譜估計(jì)出現(xiàn)錯(cuò)誤的結(jié)果將觀測(cè)數(shù)據(jù)以外的數(shù)據(jù)一律視為零，與實(shí)際不符。將觀測(cè)數(shù)據(jù)以外的數(shù)據(jù)一律視為零，與實(shí)際不符。(1) 選擇一個(gè)好的模型，在輸入是

2、沖激函數(shù)或白噪聲選擇一個(gè)好的模型，在輸入是沖激函數(shù)或白噪聲的情況下，使其輸出等于所研究的信號(hào)，至少也的情況下，使其輸出等于所研究的信號(hào)，至少也是對(duì)該信號(hào)的一個(gè)良好近似。是對(duì)該信號(hào)的一個(gè)良好近似。(2) 利用已知的自相關(guān)函數(shù)或數(shù)據(jù)求模型的參數(shù)。利用已知的自相關(guān)函數(shù)或數(shù)據(jù)求模型的參數(shù)。(3) 利用求出的模型參數(shù)或數(shù)據(jù)估計(jì)該信號(hào)的功率譜。利用求出的模型參數(shù)或數(shù)據(jù)估計(jì)該信號(hào)的功率譜。 (AutoRegressive model )(111)(0zAzazHpnnn (Moving Average model) qlllzbzH11)(AutoRegressive- Moving Average mod

3、el )()(1)(00zAzBzazbzHpnnnqlll 若輸入白噪聲的功率譜若輸入白噪聲的功率譜2)( W W P則輸出序列的功率譜為則輸出序列的功率譜為2j22j)e()()e()(W W W W W WW WHPHPx 若能確定模型中各參數(shù)若能確定模型中各參數(shù)和和就可以求得功率譜就可以求得功率譜Px(W W)描述p階AR模型的差分方程為 )(1knkynakyppn將上式兩端同乘以ykm再求數(shù)學(xué)期望 )(1kmkyEnkymkyEnakymkyEppn )(1mRnmRnamRyyppnymRynmRy mRy )(1mRnmRnamRyyppny 由于輸入是白噪聲信號(hào) R n= 2

4、d dn mRmhmRy 對(duì)于因果系統(tǒng), p階AR模型的自相關(guān)函數(shù)與模型參數(shù)的關(guān)系 pmmnmRnamRyppny, 1 , 0)(21 d d Yule-Walker(Y-W) 方程2mh Y-W 方程的矩陣表示pmmnmRnamRyppny, 1 , 0)(21 d d 例: p=3 時(shí)的Y-W 方程 000)3()2()1(101231012210132102333 aaaRRRRRRRRRRRRRRRRyyyyyyyyyyyyyyyy 若已知若已知Ryn ，由，由Y-W方程解出各參數(shù)方程解出各參數(shù)a3(1), a3 (2), a3 (3), 2。則可由則可由AR模型參數(shù)獲得功率譜模型參

5、數(shù)獲得功率譜Py(W W)的估計(jì)值。的估計(jì)值。 一階一階Y-W方程的解方程的解0) 1 (10 1 1 0211aRRRRyyyy解此方程得 0 1 ) 1 (1yyRRa) 1 ( 1 0121aRRyy) 1 (1(021aRy 二階二階Y-W方程的解方程的解00)2() 1 (10 1 2 1 0 1 2 1 02222aaRRRRRRRRRyyyyyyyyy 1 0 1 20)2(2222yyyyyRRRRRa211 1 ) 1 (2yyRaR 1 02 1 1 0) 1 (222yyyyyyRRRRRRa) 1 ()2() 1 (121aaa)2(2) 1 ( 1 02222aRaR

6、Ryyy2122)2(1 a p階階Y-W方程的遞推解方程的遞推解21111)()(pyppnypnpRnapRpa) 1, 2 , 1()()()()(11pnnpapananapppp2122)(1 ppppa(1) 計(jì)算自相關(guān)函數(shù)的估計(jì)值計(jì)算自相關(guān)函數(shù)的估計(jì)值(2) 求解一階模型參數(shù)關(guān)函數(shù)的估計(jì)值求解一階模型參數(shù)關(guān)函數(shù)的估計(jì)值0 1 ) 1 (1yyRRa)1 (1(02121aRy L-D算法估計(jì)功率譜的步驟算法估計(jì)功率譜的步驟(4) 求出功率譜估計(jì)求出功率譜估計(jì)21j2ARe)(1)(pnnpnaPWW(3)由遞推算法求解由遞推算法求解p階模型參數(shù)階模型參數(shù) L-D算法估計(jì)功率譜的步

7、驟算法估計(jì)功率譜的步驟21111)()(pyppnypnpRnapRpa) 1, 2 , 1()()()()(11pnnpapananapppp2122)(1 ppppax：進(jìn)行功率譜估計(jì)的輸入有限長(zhǎng)序列；p: 模型的階數(shù)NFFT：DFT的點(diǎn)數(shù)；Fs ：繪制功率譜曲線的抽樣頻率，默認(rèn)值為1；Pxx：功率譜估計(jì)值； f：Pxx值所對(duì)應(yīng)的頻率點(diǎn)Pxx,f = pyulear(x,p,NFFT,Fs)例：例：利用利用L-D算法進(jìn)行譜估計(jì)算法進(jìn)行譜估計(jì)一序列含有白噪聲和兩個(gè)頻率的余弦信號(hào)，一序列含有白噪聲和兩個(gè)頻率的余弦信號(hào)，)4 . 0cos()3 . 0cos(kkkkx 利用利用L-D算法估計(jì)該

8、序列的功率譜。算法估計(jì)該序列的功率譜。N=128;p=40;NFFT=2048;Fs=2; n=0:N-1;randn(state,0);x=cos(0.3*pi*n)+cos(0.4*pi*n)+randn(size(n);P,f=periodogram(x,NFFT,2) ;Py,fy=pyulear(x,p,NFFT,2);subplot(211);plot(f,10*log(P);grid;title( Periodogram );axis(0 1 -60 60);subplot(212);plot(fy,10*log(Py);grid;title( Yule);axis(0 1 -6

9、0 60);00.10.20.30.40.50.60.70.80.91-60-40-200204060 Periodogram 00.10.20.30.40.50.60.70.80.91-60-40-200204060 Yule2212)2() 1 (1zaza例例: :前向線性預(yù)測(cè)前向線性預(yù)測(cè)(p=2) 2)2( 1) 1 ( 22kyakyaky f2kykyke2)2( 1) 1 (22kyakyaky2)2( 1) 1 (22f2kyakyakykemin2f2f2keEE 若若 2, 1, 0f2nnkykeE則估計(jì)的均方誤差達(dá)到最小。則估計(jì)的均方誤差達(dá)到最小。 且且minf22f2

10、f2kykeEkeEE2)2( 1) 1 (22f2kyakyakyke 1f2kykeE0 1 )2(0) 1 ( 1 22yyyRaRaR 2f2kykeE00)2( 1 ) 1 (222yyyRaRaRf2f2kykeEE 2)2( 1 ) 1 (022yyyRaRaR00)2() 1 (10 1 2 1 0 1 2 1 0f222EaaRRRRRRRRRyyyyyyyyy例例: :后向線性預(yù)測(cè)后向線性預(yù)測(cè)(p=2) yk12 kyyk)2( 1) 1 (2 22kyakyaky2 2b2kykyke)2( 1) 1 (222kyakyaky)()()(b2bzYzEzA2122) 1

11、()2(zzaa)(12zAz由由1 , 0, 0b2nnkykeE2minb22b2b2kykeEkeEE可得可得00)2() 1 (10 1 2 1 0 1 2 1 0b222EaaRRRRRRRRRyyyyyyyyy比較前向比較前向預(yù)測(cè)系統(tǒng)和后預(yù)測(cè)系統(tǒng)和后向向預(yù)測(cè)系統(tǒng)的預(yù)測(cè)系統(tǒng)的方程后可得方程后可得 系數(shù)系數(shù)a2(1),a2(2)滿足的方程相同滿足的方程相同 2b2f2 EE Ab(z)=z 2A(z) 在計(jì)算相關(guān)函數(shù)估計(jì)時(shí)，對(duì)在計(jì)算相關(guān)函數(shù)估計(jì)時(shí)，對(duì)N個(gè)觀測(cè)數(shù)據(jù)以個(gè)觀測(cè)數(shù)據(jù)以外的數(shù)據(jù)作零的假設(shè)，故譜估計(jì)誤差較大。外的數(shù)據(jù)作零的假設(shè)，故譜估計(jì)誤差較大。 直接從觀測(cè)的數(shù)據(jù)利用線性預(yù)測(cè)器的前

12、向和直接從觀測(cè)的數(shù)據(jù)利用線性預(yù)測(cè)器的前向和后向預(yù)測(cè)的總均方誤差之和為最小的準(zhǔn)則來(lái)估計(jì)后向預(yù)測(cè)的總均方誤差之和為最小的準(zhǔn)則來(lái)估計(jì)反射系數(shù)，進(jìn)而通過(guò)反射系數(shù)，進(jìn)而通過(guò)L-D算法的遞推公式求出算法的遞推公式求出AR模型的優(yōu)化參數(shù)。模型的優(yōu)化參數(shù)。2階前向預(yù)測(cè)誤差 2) 2( 1) 1 (22f2 kyakyakyke1階后向預(yù)測(cè)誤差 ) 1 ( 11b1kyakyke L-D算法的遞推公式 ) 1 () 2() 1 () 1 (1212aaaa 2) 2( 1) 1 () 2() 1 (2121f2 kyakyaaakyke 1) 1 ( 2) 2( 1) 1 (121 kyakyakyaky 1b

13、12f1f2 keKkekeK2=a2(2) 2階預(yù)測(cè)器的反射系數(shù)階預(yù)測(cè)器的反射系數(shù) 一般地 1b1f1f keKkekepppp同理可得后向預(yù)測(cè)誤差的遞推公式 1f1b1bkeKkekepppp Kp=ap(p)為p階預(yù)測(cè)器的反射系數(shù)。 1b1f1f keKkekepppp 1f1b1bkeKkekepppp 前向和后向均方預(yù)測(cè)誤差的總和為前向和后向均方預(yù)測(cè)誤差的總和為 2b2f1kekeEppNpkp 由由0221pbpbppfpfpNpkppKkekeKkekeKE可得可得 1 1221211111kekekekeKbpfpNpkbpfpNpkp(1) 確定初始條件確定初始條件b0f0k

14、ykeke 121020kyNNk (2) 從從p=1開始迭代計(jì)算開始迭代計(jì)算: 計(jì)算計(jì)算AR模型參數(shù)模型參數(shù)1122b12f11b1f11 kekekekeKppNpkppNpkp遞推遞推p階均方誤差階均方誤差2122)1 (pppK(3) 遞推高一階前、后向預(yù)測(cè)誤差遞推高一階前、后向預(yù)測(cè)誤差 111keKkekebppfpfp 111keKkekefppbpbp(4) 若階數(shù)小于若階數(shù)小于p，則階數(shù)加，則階數(shù)加1，回到步驟，回到步驟(2)進(jìn)行下進(jìn)行下一次迭代，直到達(dá)到預(yù)定階數(shù)一次迭代，直到達(dá)到預(yù)定階數(shù)p。(5) 估計(jì)功率譜估計(jì)功率譜21j2ARe)(1)( pnnpnaPWWx：進(jìn)行功率譜估計(jì)的輸入有限長(zhǎng)序列；p: 模型的階數(shù)NFFT：DFT的點(diǎn)數(shù)；Fs ：繪制功率譜曲線的抽樣頻率，默認(rèn)值為1；Pxx：功率譜估計(jì)值； f：Pxx值所對(duì)應(yīng)的頻率點(diǎn)N=512;NFFT=1024;Fs=2;p=40;n=0:N-1;randn(state,0);x=cos(0.3*pi*n)+cos(0.32*pi*n)+randn(size(n);P,f=pyulear(x,p,NFFT,2);Pw, f2=pburg(x,p,NFFT,2);subplot(211);plot(f,10*log(P);grid;title(L-D);axis(0 1 -30 60);