ch11_1經(jīng)典功率譜估計(jì)

1、biaE)(var2EE的一致估計(jì)為則稱若, 0varlim, 0biaN)(DTFTxxPnR (1) 由隨機(jī)序列一個(gè)樣本的由隨機(jī)序列一個(gè)樣本的N個(gè)觀測(cè)值計(jì)算自相關(guān)函個(gè)觀測(cè)值計(jì)算自相關(guān)函數(shù)的估計(jì)數(shù)的估計(jì)nRx(2) 對(duì)對(duì) 進(jìn)行進(jìn)行DTFT即得該隨機(jī)序列的功率譜估即得該隨機(jī)序列的功率譜估計(jì)計(jì))( xPnRxXk是寬平穩(wěn)各態(tài)遍歷隨機(jī)信號(hào)，是寬平穩(wěn)各態(tài)遍歷隨機(jī)信號(hào)，xk是其一個(gè)樣本是其一個(gè)樣本NNkNxnkxkxNnR121lim已知已知xk的的N個(gè)觀測(cè)值個(gè)觀測(cè)值x0,x1,xN 1，則自相關(guān)函數(shù)的估計(jì)為則自相關(guān)函數(shù)的估計(jì)為110nkxkxNnRNkx*1nxnxN1) 1(NnNnRx的計(jì)算過(guò)程

2、的計(jì)算過(guò)程10Nn11NnkxnkxkxNnR 110nkxkxNnRnNkx0) 1(nN1110NnnkxkxNnRnNkxnNlnlxlxN101例例已知平穩(wěn)各態(tài)遍歷的實(shí)隨機(jī)序列已知平穩(wěn)各態(tài)遍歷的實(shí)隨機(jī)序列Xk的單一樣的單一樣本的本的N個(gè)觀測(cè)值為個(gè)觀測(cè)值為xk=1, 0， 1，試計(jì)算該隨機(jī)，試計(jì)算該隨機(jī)序列的自相關(guān)函數(shù)估計(jì)。序列的自相關(guān)函數(shù)估計(jì)。解:1, 0, 2, 0, 131nRx110nkxkxNnRNkx231nxnx1, 0, 1nx , 1, 0 , 12nx例例已知平穩(wěn)各態(tài)遍歷的實(shí)隨機(jī)序列已知平穩(wěn)各態(tài)遍歷的實(shí)隨機(jī)序列Xk的單一樣本的的單一樣本的N個(gè)觀測(cè)個(gè)觀測(cè)值為值為xk=1

3、, 0， 1，試計(jì)算該隨機(jī)序列的自相關(guān)函數(shù)估計(jì)。，試計(jì)算該隨機(jī)序列的自相關(guān)函數(shù)估計(jì)。解：解： format rat x=1 0 -1; Rx,n=xcorr(x,biased)Rx = -1/3 0 2/3 0 -1/3 n = -2 -1 0 1 2 1) 對(duì)xk補(bǔ)零形成L點(diǎn)序列) 12(NLkxLDFTkxmXLLIDFT12mXNnRLx2)3)2j)e(1DTFT XNnRx 自相關(guān)函數(shù)估計(jì)的計(jì)算自相關(guān)函數(shù)估計(jì)的計(jì)算已知平穩(wěn)各態(tài)遍歷的實(shí)隨機(jī)序列已知平穩(wěn)各態(tài)遍歷的實(shí)隨機(jī)序列Xk的單一樣本的的單一樣本的N個(gè)觀個(gè)觀測(cè)值為測(cè)值為xk=1, 0， 1，試計(jì)算該隨機(jī)序列的自相關(guān)函數(shù)，試計(jì)算該隨機(jī)序

4、列的自相關(guān)函數(shù)估計(jì)。估計(jì)。x=1 0 -1;X=fft(x,8);Rx1=ifft(abs(X).2)/3;Rx=fftshift(Rx1);Rx1 = 0.6667 0 -0.3333 0 0 0 -0.3333 0Rx =0.0000 0 -0.3333 0 0.6667 0 -0.3333 0),(DTFT,DFTmPPnRkxxxx 估計(jì)估計(jì) 101NkxnkxkxNnRnLLnxxnRPje )(e Re2 0 j11 nNnnRR 相關(guān)法相關(guān)法計(jì)算功率譜估計(jì)計(jì)算功率譜估計(jì)已知實(shí)平穩(wěn)隨機(jī)序列已知實(shí)平穩(wěn)隨機(jī)序列Xk單一樣本的單一樣本的N個(gè)觀測(cè)值為個(gè)觀測(cè)值為xk=1, 0， 1，試?yán)茫?/p>

5、試?yán)孟嚓P(guān)法相關(guān)法估計(jì)其功率譜。估計(jì)其功率譜。解：Xk的自相關(guān)函數(shù)估計(jì)值為1, 0, 2, 0, 131nRxDTFT)(nRPxx eRe22312j )2cos(132 對(duì) 進(jìn)行傅里葉變換得Xk的功率譜估計(jì)nRx功率譜估計(jì)的質(zhì)量與自相關(guān)函數(shù)估計(jì)的質(zhì)量密切相關(guān)功率譜估計(jì)的質(zhì)量與自相關(guān)函數(shù)估計(jì)的質(zhì)量密切相關(guān))(1var2nrRnrRrRNnRrbianRNnnRxxN，偏差偏差、方差趨于零，是一致估計(jì)。方差趨于零，是一致估計(jì)。N 固定時(shí)，固定時(shí)，n N，偏差較大，偏差較大其它 01, 1 , 0NkkxkxN方法基礎(chǔ)方法基礎(chǔ):*1nxnxNnRNNx由由*DTFT1)(nxnxNPNNx )e

6、()e(1jj NNXXN2j)e(1 NXN 已知已知稱稱 用用IN( )表示表示10jje DTFT)e (NkkNNNkxkxX2jjDTFT)e (1)()e (NxNNXNPXkx 功率譜估計(jì)其中其中2DFTL1mXNmPmXkxNxNN 功功率率譜譜估估計(jì)計(jì)點(diǎn)點(diǎn) 102jeDFTLkmkLNNNkxkxmX 周期圖法計(jì)算功率譜估計(jì)周期圖法計(jì)算功率譜估計(jì)已知實(shí)平穩(wěn)隨機(jī)序列已知實(shí)平穩(wěn)隨機(jī)序列Xk單一樣本的單一樣本的N個(gè)觀測(cè)值為個(gè)觀測(cè)值為xk=1, 0， 1，試?yán)弥芷趫D法估計(jì)其功率譜。，試?yán)弥芷趫D法估計(jì)其功率譜。對(duì)對(duì)xk進(jìn)行離散時(shí)間傅里葉變換進(jìn)行離散時(shí)間傅里葉變換(DTFT)2 j1

7、0jje1e )e (NkkNkxX)e ()e (1)e (1)(j*j2j NNNNXXNXNI )e1)(e1(312j2j )2cos(132功率譜估計(jì)為：功率譜估計(jì)為：njNNnNenRNnNIE)(1) 1(N，EIN( )= Px( )，漸進(jìn)無(wú)偏估計(jì)漸進(jìn)無(wú)偏估計(jì)sin)sin(1)(var24NNINN增加，增加， 方差不減小，不是一致估計(jì)方差不減小，不是一致估計(jì)利用利用MATLAB計(jì)算計(jì)算周期圖x：進(jìn)行功率譜估計(jì)的輸入有限長(zhǎng)序列；NFFT：DFT的點(diǎn)數(shù)；Fs ：繪制功率譜曲線的抽樣頻率，默認(rèn)值為1；Pxx：功率譜估計(jì)值； F：Pxx值所對(duì)應(yīng)的頻率點(diǎn)Pxx,F = periodo

8、gram(x,NFFT,Fs)利用周期圖法估計(jì)利用周期圖法估計(jì)零均值方差為零均值方差為1平穩(wěn)高斯白噪平穩(wěn)高斯白噪聲的功率譜。取序列的長(zhǎng)度分別為聲的功率譜。取序列的長(zhǎng)度分別為N=64，128，256，512。并分析譜估計(jì)質(zhì)量。并分析譜估計(jì)質(zhì)量。N=512;Nfft=2048;randn(state,0);x=randn(1,N);P,F=periodogram(x,Nfft,2) ;plot(F,10*log10(P);A=axis;axis(A(1:2) -40 20);grid;title(N=,num2str(N);xlabel(Frequency);ylabel(Power Spectr

9、al(dB);平穩(wěn)高斯白噪聲功率譜估計(jì)結(jié)果(周期圖法)00.20.40.60.81-40-20020N=64FrequencyPower Spectral(dB)00.20.40.60.81-40-20020N=128FrequencyPower Spectral(dB)00.20.40.60.81-40-20020N=256FrequencyPower Spectral(dB)00.20.40.60.81-40-20020N=512FrequencyPower Spectral(dB) 波動(dòng)的大小不隨數(shù)據(jù)長(zhǎng)度波動(dòng)的大小不隨數(shù)據(jù)長(zhǎng)度N的增加而減小，即周期圖法譜估的增加而減小，即周期圖法譜估計(jì)的

10、方差較大，且不隨計(jì)的方差較大，且不隨N的增加而減小。的增加而減小。 周期圖法進(jìn)行功率譜估計(jì)周期圖法進(jìn)行功率譜估計(jì), 方差不隨方差不隨N的增加減小。的增加減小。如何提高譜估計(jì)質(zhì)量如何提高譜估計(jì)質(zhì)量?1) 對(duì)自相關(guān)函數(shù)估計(jì)值加窗對(duì)自相關(guān)函數(shù)估計(jì)值加窗2) 將將N個(gè)觀測(cè)值分段個(gè)觀測(cè)值分段, 計(jì)算各段的周期圖計(jì)算各段的周期圖, 再取平均再取平均對(duì)自相關(guān)函數(shù)估計(jì)值加窗對(duì)自相關(guān)函數(shù)估計(jì)值加窗, 將誤差較大的估計(jì)值截去將誤差較大的估計(jì)值截去nxNNnMnRnwPj1)1(e )(PM( )波動(dòng)比波動(dòng)比IN ( )小，可證是一致估計(jì)小，可證是一致估計(jì)降低了頻率分辨率降低了頻率分辨率1, 1 , 0; 1, 1

11、 , 0;MkAikiMxkxi第第i段序列的周期圖為段序列的周期圖為2j)e (1)( MiMXMI 平均周期圖平均周期圖)(1)(10 iMAiAMIAP 將隨機(jī)序列將隨機(jī)序列Xk的的N個(gè)觀測(cè)值分成個(gè)觀測(cè)值分成A段段平均周期圖法平均周期圖法：減小方差：減小方差：增加估計(jì)的偏差，：增加估計(jì)的偏差，降低了譜的分辨率降低了譜的分辨率：分段即加窗，段越多，窗越短，主瓣寬度越大：分段即加窗，段越多，窗越短，主瓣寬度越大：將各段數(shù)據(jù)有一定程度的重疊：將各段數(shù)據(jù)有一定程度的重疊 X:進(jìn)行功率譜估計(jì)的輸入有限長(zhǎng)序列；進(jìn)行功率譜估計(jì)的輸入有限長(zhǎng)序列； WINDOW: (1)給定長(zhǎng)度給定長(zhǎng)度 L，則用指定長(zhǎng)度

12、的，則用指定長(zhǎng)度的Hamming窗。窗。 (2)給定序列給定序列 ，則用該序列為窗。，則用該序列為窗。 NOVERLAP指定分段重疊的樣本數(shù)。缺省時(shí)指定分段重疊的樣本數(shù)。缺省時(shí)NOVERLAP=L/2，L為窗的長(zhǎng)度。為窗的長(zhǎng)度。：DFT的點(diǎn)數(shù)，的點(diǎn)數(shù)，NFFTL。 ：繪制功率譜曲線的抽樣頻率，默認(rèn)值為：繪制功率譜曲線的抽樣頻率，默認(rèn)值為1；：功率譜估計(jì)值；：功率譜估計(jì)值；F：Pxx值所對(duì)應(yīng)的頻率點(diǎn)值所對(duì)應(yīng)的頻率點(diǎn)：一序列含有白噪聲和兩個(gè)頻率的余弦信號(hào)，設(shè)一序列含有白噪聲和兩個(gè)頻率的余弦信號(hào)，設(shè))34. 0cos()3 . 0cos(kwkkkx分別采用周期圖法和分別采用周期圖法和Welch法估

13、計(jì)該序列的功率譜。法估計(jì)該序列的功率譜。N=1024;NFFT=1024;Fs=2*pi;L=128;n=0:N-1;randn(state,0);x=cos(0.3*pi*n)+cos(0.4*pi*n)+randn(size(n);window=boxcar(L);P,f=periodogram(x,NFFT,2) ;Pw, f2=pwelch(x,window,NFFT,2);subplot(211);plot(f2,10*log(Pw);grid;title(Welch);axis(0 1 -60 60);subplot(212);plot(f,10*log(P);grid;title(Periodogram);axis(0 1 -60 60);00.20.40.60.81-50050Welch00.20.40.60.81-50050Periodogram 周期圖法譜估計(jì)曲線的波動(dòng)很大，即估計(jì)的方差較大。周期圖法譜估計(jì)曲線的波動(dòng)很大，即估計(jì)的方差較大。 WelchWelch法譜估計(jì)曲線較為平滑，方差減小，但分辨率降低。法譜估計(jì)曲線較為平滑，方差減小，但分辨率降低。 對(duì)對(duì)WelchWelch法，當(dāng)數(shù)據(jù)分段數(shù)增加，各段數(shù)據(jù)長(zhǎng)度較短時(shí)，譜法，當(dāng)數(shù)據(jù)分段數(shù)增加，各段數(shù)據(jù)長(zhǎng)度較短時(shí)，譜的分辨率明顯下降，而譜估計(jì)曲線較為平滑，方差較??；的分辨率明顯下降，而譜估計(jì)曲