ch64現(xiàn)代功率譜估計(jì)

數(shù)字信號(hào)處理

(DigitalSignalProcessing)

信號(hào)與系統(tǒng)系列課程組國(guó)家電工電子教學(xué)基地?cái)?shù)字信號(hào)處理

(DigitalSignalProcess1離散隨機(jī)序列的特征描述平穩(wěn)隨機(jī)序列通過(guò)LTI系統(tǒng)經(jīng)典功率譜估計(jì)

現(xiàn)代功率譜估計(jì)隨機(jī)信號(hào)的功率譜估計(jì)離散隨機(jī)序列的特征描述隨機(jī)信號(hào)的功率譜估計(jì)2現(xiàn)代譜估計(jì)簡(jiǎn)介問(wèn)題提出平穩(wěn)隨機(jī)信號(hào)的參數(shù)模型AR模型參數(shù)與自相關(guān)函數(shù)的關(guān)系Y-W方程的L-D遞推算法AR模型參數(shù)與線性預(yù)測(cè)濾波器的關(guān)系伯格(Burg)遞推算法現(xiàn)代譜估計(jì)簡(jiǎn)介問(wèn)題提出3問(wèn)題提出經(jīng)典法存在問(wèn)題:1.方差性能不好，不是Px(W)的一致估計(jì)2.平滑周期圖和平均周期圖改善了周期圖的方差性能，但卻降低了譜分辨率和增大了偏差。3.可能使短序列的功率譜估計(jì)出現(xiàn)錯(cuò)誤的結(jié)果出現(xiàn)問(wèn)題的原因:將觀測(cè)數(shù)據(jù)以外的數(shù)據(jù)一律視為零，與實(shí)際不符。問(wèn)題提出經(jīng)典法存在問(wèn)題:1.方差性能不好，不是Px(4參數(shù)模型法的基本思想根據(jù)所研究信號(hào)的先驗(yàn)知識(shí)，對(duì)觀測(cè)數(shù)據(jù)以外的數(shù)據(jù)作出某種比較合理的假設(shè)。假設(shè)信號(hào)是白噪聲通過(guò)LTI系統(tǒng)產(chǎn)生的。由觀測(cè)數(shù)據(jù)估計(jì)LTI系統(tǒng)模型的參數(shù)。最后由LTI系統(tǒng)模型的參數(shù)得出功率譜。hh[k]y[k]輸入白噪聲的自相關(guān)函數(shù)參數(shù)模型法的基本思想根據(jù)所研究信號(hào)的先驗(yàn)知識(shí)5平穩(wěn)隨機(jī)信號(hào)的參數(shù)模型AR模型(AutoRegressivemodel)MA模型(MovingAveragemodel)ARMA模型(AutoRegressive-MovingAveragemodel)平穩(wěn)隨機(jī)信號(hào)的參數(shù)模型AR模型MA模型ARMA模6若輸入白噪聲的功率譜則輸出序列的功率譜為若能確定模型中各參數(shù){an,bn}和s2就可以求得功率譜Px()平穩(wěn)隨機(jī)信號(hào)的參數(shù)模型若輸入白噪聲的功率譜則輸出序列的功率譜為若能確定模型中各參數(shù)7AR模型參數(shù)與自相關(guān)函數(shù)的關(guān)系描述p階AR模型的差分方程為

將上式兩端同乘以y[k-m]再求數(shù)學(xué)期望

AR模型參數(shù)與自相關(guān)函數(shù)的關(guān)系描述p階AR模型的差分方程為8AR模型參數(shù)與自相關(guān)函數(shù)的關(guān)系由于輸入是白噪聲信號(hào)Rh[n]=s2d[n]對(duì)于因果系統(tǒng),p階AR模型的自相關(guān)函數(shù)與模型參數(shù)的關(guān)系

Yule-Walker(Y-W)方程AR模型參數(shù)與自相關(guān)函數(shù)的關(guān)系由于輸入是白噪聲信號(hào)Rh9AR模型參數(shù)與自相關(guān)函數(shù)的關(guān)系Y-W方程的矩陣表示例:p=3時(shí)的Y-W方程若已知Ry[n]，由Y-W方程解出各參數(shù)a3(1),a3(2),

a3(3),s2。則可由AR模型參數(shù)獲得功率譜Py()的估計(jì)值。AR模型參數(shù)與自相關(guān)函數(shù)的關(guān)系Y-W方程的矩陣表示例:p10Y-W方程的L-D遞推算法一階Y-W方程的解解此方程得

Y-W方程的L-D遞推算法一階Y-W方程的解解此方程得11Y-W方程的L-D遞推算法二階Y-W方程的解Y-W方程的L-D遞推算法二階Y-W方程的解12Y-W方程的L-D遞推算法

p階Y-W方程的遞推解Y-W方程的L-D遞推算法p階Y-W方程的遞推解13Y-W方程的L-D遞推算法(1)計(jì)算自相關(guān)函數(shù)的估計(jì)值(2)求解一階模型參數(shù)關(guān)函數(shù)的估計(jì)值

L-D算法估計(jì)功率譜的步驟Y-W方程的L-D遞推算法(1)計(jì)算自相關(guān)函數(shù)的估計(jì)值(14(4)求出功率譜估計(jì)Y-W方程的L-D遞推算法(3)由遞推算法求解p階模型參數(shù)

L-D算法估計(jì)功率譜的步驟(4)求出功率譜估計(jì)Y-W方程的L-D遞推算法(3)由遞15利用L-D算法估計(jì)頻譜的MATLAB函數(shù)x：進(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ì)頻譜的MATLAB函數(shù)x：進(jìn)行功率譜估計(jì)的16例：利用L-D算法進(jìn)行譜估計(jì)一序列含有白噪聲和兩個(gè)頻率的余弦信號(hào)，利用L-D算法估計(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([01-6060]);subplot(212);plot(fy,10*log(Py));grid;title('Yule');axis([01-6060]);例：利用L-D算法進(jìn)行譜估計(jì)一序列含有白噪聲和兩個(gè)頻率的余弦17譜估計(jì)結(jié)果——p=40,N=12800.10.20.30.40.50.60.70.80.91-60-40-200204060Periodogram00.10.20.30.40.50.60.70.80.91-60-40-200204060Yule譜估計(jì)結(jié)果——p=40,N=12800.10.20.30.418AR模型參數(shù)與前向線性預(yù)測(cè)濾波器的關(guān)系前向線性預(yù)測(cè)濾波器前向預(yù)測(cè)誤差前向預(yù)測(cè)誤差濾波器系統(tǒng)函數(shù)例:前向線性預(yù)測(cè)(p=2)

AR模型參數(shù)與前向線性預(yù)測(cè)濾波器的關(guān)系前向線性預(yù)測(cè)濾波器19AR模型參數(shù)與前向線性預(yù)測(cè)濾波器的關(guān)系均方誤差最小的前向預(yù)測(cè)誤差濾波器正交準(zhǔn)則若

則估計(jì)的均方誤差達(dá)到最小。且AR模型參數(shù)與前向線性預(yù)測(cè)濾波器的關(guān)系均方誤差最小的前向預(yù)20預(yù)測(cè)系數(shù)及均方誤差的確定

AR模型參數(shù)與前向線性預(yù)測(cè)濾波器的關(guān)系預(yù)測(cè)系數(shù)及均方誤差的確定AR模型參數(shù)與前向線性預(yù)測(cè)濾波器的21后向線性預(yù)測(cè)濾波器后向預(yù)測(cè)誤差后向預(yù)測(cè)誤差濾波器系統(tǒng)函數(shù)AR模型參數(shù)與后向線性預(yù)測(cè)濾波器的關(guān)系例:后向線性預(yù)測(cè)(p=2)

后向線性預(yù)測(cè)濾波器后向預(yù)測(cè)誤差后向預(yù)測(cè)誤差濾波器系統(tǒng)22預(yù)測(cè)系數(shù)及均方誤差的確定AR模型參數(shù)與后向線性預(yù)測(cè)濾波器的關(guān)系由可得比較前向預(yù)測(cè)系統(tǒng)和后向預(yù)測(cè)系統(tǒng)的方程后可得

系數(shù){a2(1),a2(2)}滿(mǎn)足的方程相同

Ab(z)=z-2A(z)預(yù)測(cè)系數(shù)及均方誤差的確定AR模型參數(shù)與后向線性預(yù)測(cè)濾波器的關(guān)23伯格(Burg)遞推算法L-D算法缺點(diǎn)：在計(jì)算相關(guān)函數(shù)估計(jì)時(shí)，對(duì)N個(gè)觀測(cè)數(shù)據(jù)以外的數(shù)據(jù)作零的假設(shè)，故譜估計(jì)誤差較大。伯格(Burg)遞推算法基本思想：直接從觀測(cè)的數(shù)據(jù)利用線性預(yù)測(cè)器的前向和后向預(yù)測(cè)的總均方誤差之和為最小的準(zhǔn)則來(lái)估計(jì)反射系數(shù)，進(jìn)而通過(guò)L-D算法的遞推公式求出AR模型的優(yōu)化參數(shù)。伯格(Burg)遞推算法L-D算法缺點(diǎn)：24伯格(Burg)遞推算法前向預(yù)測(cè)誤差的遞推公式2階前向預(yù)測(cè)誤差

1階后向預(yù)測(cè)誤差

L-D算法的遞推公式

K2=a2(2)2階預(yù)測(cè)器的反射系數(shù)

伯格(Burg)遞推算法前向預(yù)測(cè)誤差的遞推公式2階前向預(yù)測(cè)25伯格(Burg)遞推算法預(yù)測(cè)誤差的遞推公式一般地同理可得后向預(yù)測(cè)誤差的遞推公式

Kp=ap(p)為p階預(yù)測(cè)器的反射系數(shù)。

伯格(Burg)遞推算法預(yù)測(cè)誤差的遞推公式一般地同理可得后26伯格(Burg)遞推算法預(yù)測(cè)誤差濾波器的格形結(jié)構(gòu)

伯格(Burg)遞推算法預(yù)測(cè)誤差濾波器的格形結(jié)構(gòu)27伯格(Burg)遞推算法反射系數(shù)Kp

的確定前向和后向均方預(yù)測(cè)誤差的總和為

由可得伯格(Burg)遞推算法反射系數(shù)Kp的確定前向和后向均方預(yù)28伯格(Burg)遞推算法步驟(1)確定初始條件(2)從p=1開(kāi)始迭代計(jì)算:計(jì)算AR模型參數(shù)遞推p階均方誤差伯格(Burg)遞推算法伯格(Burg)遞推算法步驟(1)確定初始條件(2)29(3)遞推高一階前、后向預(yù)測(cè)誤差(4)若階數(shù)小于p，則階數(shù)加1，回到步驟(2)進(jìn)行下一次迭代，直到達(dá)到預(yù)定階數(shù)p。(5)估計(jì)功率譜伯格(Burg)遞推算法步驟伯格(Burg)遞推算法(3)遞推高一階前、后向預(yù)測(cè)誤差(4)若階數(shù)小于p，30[Pxx,f]=pburg(x,p,NFFT,Fs)Burg算法估計(jì)頻譜的MATLAB函數(shù)x：進(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]=pburg(x,p,NFFT,Fs)B31利用Burg法進(jìn)行譜估計(jì)程序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([01-3060]);subplot(212);plot(f2,10*log(Pw));grid;title('Burg');axis([01-3060]);利用Burg法進(jìn)行譜估計(jì)程序N=512;NFFT=1024;32AR模型階數(shù)p=50的譜估計(jì)結(jié)果00.10.20.30.40.50.60.70.80.91-200204060L-D00.10.20.30.40.50.60.70.80.91-200204060BurgAR模型階數(shù)p=50的譜估計(jì)結(jié)果00.10.20.30.433AR模型階數(shù)p=80的譜估計(jì)結(jié)果00.10.20.30.40.50.60.70.80.91-200204060L-D00.10.20.30.40.50.60.70.80.91-200204060BurgAR模型階數(shù)p=80的譜估計(jì)結(jié)果00.10.20.30.434數(shù)字信號(hào)處理

(DigitalSignalProcessing)

信號(hào)與系統(tǒng)系列課程組國(guó)家電工電子教學(xué)基地?cái)?shù)字信號(hào)處理

(DigitalSignalProcess35離散隨機(jī)序列的特征描述平穩(wěn)隨機(jī)序列通過(guò)LTI系統(tǒng)經(jīng)典功率譜估計(jì)

現(xiàn)代功率譜估計(jì)隨機(jī)信號(hào)的功率譜估計(jì)離散隨機(jī)序列的特征描述隨機(jī)信號(hào)的功率譜估計(jì)36現(xiàn)代譜估計(jì)簡(jiǎn)介問(wèn)題提出平穩(wěn)隨機(jī)信號(hào)的參數(shù)模型AR模型參數(shù)與自相關(guān)函數(shù)的關(guān)系Y-W方程的L-D遞推算法AR模型參數(shù)與線性預(yù)測(cè)濾波器的關(guān)系伯格(Burg)遞推算法現(xiàn)代譜估計(jì)簡(jiǎn)介問(wèn)題提出37問(wèn)題提出經(jīng)典法存在問(wèn)題:1.方差性能不好，不是Px(W)的一致估計(jì)2.平滑周期圖和平均周期圖改善了周期圖的方差性能，但卻降低了譜分辨率和增大了偏差。3.可能使短序列的功率譜估計(jì)出現(xiàn)錯(cuò)誤的結(jié)果出現(xiàn)問(wèn)題的原因:將觀測(cè)數(shù)據(jù)以外的數(shù)據(jù)一律視為零，與實(shí)際不符。問(wèn)題提出經(jīng)典法存在問(wèn)題:1.方差性能不好，不是Px(38參數(shù)模型法的基本思想根據(jù)所研究信號(hào)的先驗(yàn)知識(shí)，對(duì)觀測(cè)數(shù)據(jù)以外的數(shù)據(jù)作出某種比較合理的假設(shè)。假設(shè)信號(hào)是白噪聲通過(guò)LTI系統(tǒng)產(chǎn)生的。由觀測(cè)數(shù)據(jù)估計(jì)LTI系統(tǒng)模型的參數(shù)。最后由LTI系統(tǒng)模型的參數(shù)得出功率譜。hh[k]y[k]輸入白噪聲的自相關(guān)函數(shù)參數(shù)模型法的基本思想根據(jù)所研究信號(hào)的先驗(yàn)知識(shí)39平穩(wěn)隨機(jī)信號(hào)的參數(shù)模型AR模型(AutoRegressivemodel)MA模型(MovingAveragemodel)ARMA模型(AutoRegressive-MovingAveragemodel)平穩(wěn)隨機(jī)信號(hào)的參數(shù)模型AR模型MA模型ARMA模40若輸入白噪聲的功率譜則輸出序列的功率譜為若能確定模型中各參數(shù){an,bn}和s2就可以求得功率譜Px()平穩(wěn)隨機(jī)信號(hào)的參數(shù)模型若輸入白噪聲的功率譜則輸出序列的功率譜為若能確定模型中各參數(shù)41AR模型參數(shù)與自相關(guān)函數(shù)的關(guān)系描述p階AR模型的差分方程為

將上式兩端同乘以y[k-m]再求數(shù)學(xué)期望

AR模型參數(shù)與自相關(guān)函數(shù)的關(guān)系描述p階AR模型的差分方程為42AR模型參數(shù)與自相關(guān)函數(shù)的關(guān)系由于輸入是白噪聲信號(hào)Rh[n]=s2d[n]對(duì)于因果系統(tǒng),p階AR模型的自相關(guān)函數(shù)與模型參數(shù)的關(guān)系

Yule-Walker(Y-W)方程AR模型參數(shù)與自相關(guān)函數(shù)的關(guān)系由于輸入是白噪聲信號(hào)Rh43AR模型參數(shù)與自相關(guān)函數(shù)的關(guān)系Y-W方程的矩陣表示例:p=3時(shí)的Y-W方程若已知Ry[n]，由Y-W方程解出各參數(shù)a3(1),a3(2),

a3(3),s2。則可由AR模型參數(shù)獲得功率譜Py()的估計(jì)值。AR模型參數(shù)與自相關(guān)函數(shù)的關(guān)系Y-W方程的矩陣表示例:p44Y-W方程的L-D遞推算法一階Y-W方程的解解此方程得

Y-W方程的L-D遞推算法一階Y-W方程的解解此方程得45Y-W方程的L-D遞推算法二階Y-W方程的解Y-W方程的L-D遞推算法二階Y-W方程的解46Y-W方程的L-D遞推算法

p階Y-W方程的遞推解Y-W方程的L-D遞推算法p階Y-W方程的遞推解47Y-W方程的L-D遞推算法(1)計(jì)算自相關(guān)函數(shù)的估計(jì)值(2)求解一階模型參數(shù)關(guān)函數(shù)的估計(jì)值

L-D算法估計(jì)功率譜的步驟Y-W方程的L-D遞推算法(1)計(jì)算自相關(guān)函數(shù)的估計(jì)值(48(4)求出功率譜估計(jì)Y-W方程的L-D遞推算法(3)由遞推算法求解p階模型參數(shù)

L-D算法估計(jì)功率譜的步驟(4)求出功率譜估計(jì)Y-W方程的L-D遞推算法(3)由遞49利用L-D算法估計(jì)頻譜的MATLAB函數(shù)x：進(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ì)頻譜的MATLAB函數(shù)x：進(jìn)行功率譜估計(jì)的50例：利用L-D算法進(jìn)行譜估計(jì)一序列含有白噪聲和兩個(gè)頻率的余弦信號(hào)，利用L-D算法估計(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([01-6060]);subplot(212);plot(fy,10*log(Py));grid;title('Yule');axis([01-6060]);例：利用L-D算法進(jìn)行譜估計(jì)一序列含有白噪聲和兩個(gè)頻率的余弦51譜估計(jì)結(jié)果——p=40,N=12800.10.20.30.40.50.60.70.80.91-60-40-200204060Periodogram00.10.20.30.40.50.60.70.80.91-60-40-200204060Yule譜估計(jì)結(jié)果——p=40,N=12800.10.20.30.452AR模型參數(shù)與前向線性預(yù)測(cè)濾波器的關(guān)系前向線性預(yù)測(cè)濾波器前向預(yù)測(cè)誤差前向預(yù)測(cè)誤差濾波器系統(tǒng)函數(shù)例:前向線性預(yù)測(cè)(p=2)

AR模型參數(shù)與前向線性預(yù)測(cè)濾波器的關(guān)系前向線性預(yù)測(cè)濾波器53AR模型參數(shù)與前向線性預(yù)測(cè)濾波器的關(guān)系均方誤差最小的前向預(yù)測(cè)誤差濾波器正交準(zhǔn)則若

則估計(jì)的均方誤差達(dá)到最小。且AR模型參數(shù)與前向線性預(yù)測(cè)濾波器的關(guān)系均方誤差最小的前向預(yù)54預(yù)測(cè)系數(shù)及均方誤差的確定

AR模型參數(shù)與前向線性預(yù)測(cè)濾波器的關(guān)系預(yù)測(cè)系數(shù)及均方誤差的確定AR模型參數(shù)與前向線性預(yù)測(cè)濾波器的55后向線性預(yù)測(cè)濾波器后向預(yù)測(cè)誤差后向預(yù)測(cè)誤差濾波器系統(tǒng)函數(shù)AR模型參數(shù)與后向線性預(yù)測(cè)濾波器的關(guān)系例:后向線性預(yù)測(cè)(p=2)

后向線性預(yù)測(cè)濾波器后向預(yù)測(cè)誤差后向預(yù)測(cè)誤差濾波器系統(tǒng)56預(yù)測(cè)系數(shù)及均方誤差的確定AR模型參數(shù)與后向線性預(yù)測(cè)濾波器的關(guān)系由可得比較前向預(yù)測(cè)系統(tǒng)和后向預(yù)測(cè)系統(tǒng)的方程后可得

系數(shù){a2(1),a2(2)}滿(mǎn)足的方程相同

Ab(z)=z-2A(z)預(yù)測(cè)系數(shù)及均方誤差的確定AR模型參數(shù)與后向線性預(yù)測(cè)濾波器的關(guān)57伯格(Burg)遞推算法L-D算法缺點(diǎn)：在計(jì)算相關(guān)函數(shù)估計(jì)時(shí)，對(duì)N個(gè)觀測(cè)數(shù)據(jù)以外的數(shù)據(jù)作零的假設(shè)，故譜估計(jì)誤差較大。伯格(Burg)遞推算法基本思想：直接從觀測(cè)的數(shù)據(jù)利用線性預(yù)測(cè)器的前向和后向預(yù)測(cè)的總均方誤差之和為最小的準(zhǔn)則來(lái)估計(jì)反射系數(shù)，進(jìn)而通過(guò)L-D算法的遞推公式求出AR模型的優(yōu)化參數(shù)。伯格(Burg)遞推算法L-D算法缺點(diǎn)：58伯格(Burg)遞推算法前向預(yù)測(cè)誤差的遞推公式2階前向預(yù)測(cè)誤差

1階后向預(yù)測(cè)誤差

L-D算法的遞推公式

K2=a2(2)2階預(yù)測(cè)器的反射系數(shù)

伯格(Burg)遞推算法前向預(yù)測(cè)誤差的遞推公式2階前向預(yù)測(cè)59伯格(Burg)遞推算法預(yù)測(cè)誤差的遞推公式一般地同理可得后向預(yù)測(cè)誤差的遞推公式

Kp=ap(p)為p階預(yù)測(cè)器的反射系數(shù)。

伯格(Burg)遞推算法預(yù)測(cè)誤差的遞推公式一般地同理可得后60伯格(Burg)遞推算法預(yù)測(cè)誤差濾波器的格形結(jié)構(gòu)

伯格(Burg)遞推算法預(yù)測(cè)誤差濾波器的格形結(jié)構(gòu)61伯格(Burg)遞推算法反射系數(shù)Kp

的確定前向和后向均方預(yù)測(cè)誤差的總和為

由可得伯格(Burg)遞推算法反射系數(shù)Kp的確定前向和后向均方預(yù)62伯格(Burg)遞推算法步驟(1)確定初始條件(2)從p=1開(kāi)始迭代計(jì)算:計(jì)算AR模型參數(shù)遞推p階均方誤差伯格(Burg)遞推算法伯格(Burg)遞推算法步驟(1)確定初始條件(2)63(3)遞推高一階前、后向預(yù)測(cè)誤差(4)若階數(shù)小于p，則階數(shù)加1，回到步驟(2)進(jìn)行