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

2、沖激函數或白噪聲選擇一個好的模型，在輸入是沖激函數或白噪聲的情況下，使其輸出等于所研究的信號，至少也的情況下，使其輸出等于所研究的信號，至少也是對該信號的一個良好近似。是對該信號的一個良好近似。(2) 利用已知的自相關函數或數據求模型的參數。利用已知的自相關函數或數據求模型的參數。(3) 利用求出的模型參數或數據估計該信號的功率譜。利用求出的模型參數或數據估計該信號的功率譜。 (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 若能確定模型中各參數若能確定模型中各參數和和就可以求得功率譜就可以求得功率譜Px(W W)描述p階AR模型的差分方程為 )(1knkynakyppn將上式兩端同乘以ykm再求數學期望 )(1kmkyEnkymkyEnakymkyEppn )(1mRnmRnamRyyppnymRynmRy mRy )(1mRnmRnamRyyppny 由于輸入是白噪聲信號 R n= 2

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

5、數獲得功率譜Py(W W)的估計值。的估計值。 一階一階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) 計算自相關函數的估計值計算自相關函數的估計值(2) 求解一階模型參數關函數的估計值求解一階模型參數關函數的估計值0 1 ) 1 (1yyRRa)1 (1(02121aRy L-D算法估計功率譜的步驟算法估計功率譜的步驟(4) 求出功率譜估計求出功率譜估計21j2ARe)(1)(pnnpnaPWW(3)由遞推算法求解由遞推算法求解p階模型參數階模型參數 L-D算法估計功率譜的步

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

8、序列的功率譜。算法估計該序列的功率譜。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例例: :前向線性預測前向線性預測(p=2) 2)2( 1) 1 ( 22kyakyaky f2kykyke2)2( 1) 1 (22kyakyaky2)2( 1) 1 (22f2kyakyakykemin2f2f2keEE 若若 2, 1, 0f2nnkykeE則估計的均方誤差達到最小。則估計的均方誤差達到最小。 且且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例例: :后向線性預測后向線性預測(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比較前向比較前向預測系統和后預測系統和后向向預測系統的預測系統的方程后可得方程后可得 系數系數a2(1),a2(2)滿足的方程相同滿足的方程相同 2b2f2 EE Ab(z)=z 2A(z) 在計算相關函數估計時，對在計算相關函數估計時，對N個觀測數據以個觀測數據以外的數據作零的假設，故譜估計誤差較大。外的數據作零的假設，故譜估計誤差較大。 直接從觀測的數據利用線性預測器的前

12、向和直接從觀測的數據利用線性預測器的前向和后向預測的總均方誤差之和為最小的準則來估計后向預測的總均方誤差之和為最小的準則來估計反射系數，進而通過反射系數，進而通過L-D算法的遞推公式求出算法的遞推公式求出AR模型的優化參數。模型的優化參數。2階前向預測誤差 2) 2( 1) 1 (22f2 kyakyakyke1階后向預測誤差 ) 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階預測器的反射系數階預測器的反射系數 一般地 1b1f1f keKkekepppp同理可得后向預測誤差的遞推公式 1f1b1bkeKkekepppp Kp=ap(p)為p階預測器的反射系數。 1b1f1f keKkekepppp 1f1b1bkeKkekepppp 前向和后向均方預測誤差的總和為前向和后向均方預測誤差的總和為 2b2f1kekeEppNpkp 由由0221pbpbppfpfpNpkppKkekeKkekeKE可得可得 1 1221211111kekekekeKbpfpNpkbpfpNpkp(1) 確定初始條件確定初始條件b0f0k

14、ykeke 121020kyNNk (2) 從從p=1開始迭代計算開始迭代計算: 計算計算AR模型參數模型參數1122b12f11b1f11 kekekekeKppNpkppNpkp遞推遞推p階均方誤差階均方誤差2122)1 (pppK(3) 遞推高一階前、后向預測誤差遞推高一階前、后向預測誤差 111keKkekebppfpfp 111keKkekefppbpbp(4) 若階數小于若階數小于p，則階數加，則階數加1，回到步驟，回到步驟(2)進行下進行下一次迭代，直到達到預定階數一次迭代，直到達到預定階數p。(5) 估計功率譜估計功率譜21j2ARe)(1)( pnnpnaPWWx：進行功率譜估計的輸入有限長序列；p: 模型的階數NFFT：DFT的點數；Fs ：繪制功率譜曲線的抽樣頻率，默認值為1；Pxx：功率譜估計值； f：Pxx值所對應的頻率點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);