【正文】 模的平方除以 N 求得的功率譜估計(jì)的方法稱為周期圖法，其結(jié)果用 ()NI ? 表示，即 21? ( ) ( ) ( )jX N NG I X eN ????? (32) 周期圖法是利用數(shù)據(jù)的傅里葉變換直接求得的，而不再計(jì)算自相關(guān)函數(shù)，所以又稱 直接法。 周期圖的性能 為了了解周期圖法的譜估計(jì)效果，我們來討論它的估計(jì)均值和方差。 當(dāng) N?? 時(shí) lim [ ] lim 1 1NNN mbm N? ? ? ? ??? ? ????? 故 lim ( ) ( )NN B ? ? ??? ? 此時(shí) l i m [ ( ) ] ( )NXN E I G???? ? (34) 因此，周期圖作為功率譜估計(jì)，當(dāng) N?? 時(shí)是無偏的，即漸近無偏。 為分析簡單起見，通常假設(shè) []Xn零均值，方差為 2X? 的實(shí)高斯白噪聲序列，即功率譜密度為常數(shù) 2X? 。 為了求 2[ ( )]NEI? ，先求 ()NI ? 在兩個(gè)頻率 1? 和 2? 處的協(xié)方差，最后令 12? ? ???。利用正態(tài)白噪聲、多元正態(tài)隨機(jī)變量的多階矩公式，有 [ [ ] [ ] [ ] [ ] ] [ [ ] [ ] ] [ [ ] [ ] ] [ [ ] [ ] ] [ [ ] [ ] ]E X n X k X p X q E X n X k E X p X q E X n X p E X k X q?? 基于 MATLAB 的譜估計(jì)實(shí)現(xiàn) 9 [ [ ] [ ] ] [ [ ] [ ] ]E X n X q E X k X p? 40X????? , 。 ,k n p q p n q k q n p k? ? ? ? ? ?其 他 將上式代入式 (36)中，得 1 2 1 24 1 1 1 1( ) ( ) ( ) ( )212 2 0 0 0 0[ ( ) ( ) ]N N N Nj n k j n kXNN n k n kE I I N e eN ? ? ? ????? ? ? ?? ? ? ? ? ?? ? ? ???? ? ?????? ? ? ? 221 2 1 241 2 1 2( ) ( )sin sin221( ) ( )sin sinXNNNN? ? ? ?? ? ? ? ?????? ? ? ???? ? ? ?? ? ???? ? ? ?? ? ? ?? ? ? ??? 將 12? ? ???代入上式，可得 224 s in ( )[ ( ) ] 2 s in ( )NX NEI N ??? ????????????????? 24 s in ( )v a r [ ( ) ] 1s in ( )NX NI N ??? ????????????????? 由上式可見，當(dāng) N?? 時(shí)， 4var[ ( )] 0NXI ????。顯然，周期圖不是功率譜估計(jì)的一致估計(jì)，所以這種估計(jì)不是最好的估計(jì)方法。平均周期圖法就是將長度為 N 的觀測數(shù)據(jù)截?cái)喑?K 段，每一段長度為 L ,分別計(jì)算每一段的周期圖，最后進(jìn)行平均得到平均周期圖的估計(jì)，其數(shù)學(xué)表達(dá)式如下： 1( ) 201? ( ) | ( ) |Lm j np e r mnG X n eL?? ? ??? ? 1 ()01? ?( ) ( )K ma v p e r p e rmGGK????? ? 基于 MATLAB 的譜估計(jì)實(shí)現(xiàn) 10 由此可見，估計(jì)的方差很明顯多了一項(xiàng) 1K ，這樣方差便隨著分段數(shù)目的增加而下降。 加窗周期圖 周期圖法相當(dāng)于只取了原信號的一段進(jìn)行估計(jì)，也就是說對原信號加了一個(gè)矩形窗。在理想的情況下 應(yīng)該是原信號乘以幅值始終為 1的信號，此信號在頻域表現(xiàn)為一沖激函數(shù)，主瓣無限窄。如果原來真是的功率譜是窄的，那么與主瓣卷積后會使功率向附近頻域擴(kuò)散，卷積以后使得加窗的信號的頻譜變成 Sa 函數(shù)頻譜的形狀，即有一定寬度，存在若干旁瓣的光滑脈沖那樣。 對其進(jìn)行改進(jìn)的辦法是將矩形窗換為其它的窗函數(shù)。 旁瓣較低的窗是 blackman窗和 hamming窗 ， 但這兩種窗的主瓣寬度比矩形窗的主瓣寬度大，因此分辨率沒有矩形窗的分辨率高。 周期圖的仿真 在 MATLAB 中用周期圖方法進(jìn)行功率譜估計(jì)的函數(shù)為 periodogram。 n=0:1/Fs:1。 nfft=1024。 figure(1) plot(f,10*log10(Pxx))。所以周期圖不是功率譜的一致估計(jì)。 基于 MATLAB 的譜估計(jì)實(shí)現(xiàn) 13 平均周期圖的仿真 仿真程序如下： Fs=600。%產(chǎn)生含有噪聲的序列 xn=cos(2*pi*40*n)+ 5*cos(2*pi*90*n)+*randn(size(n))。 Pxx=bart(xn,10)。 k=index*Fs/256。 figure(1) plot(k,plot_Pxx)。 sx=0。 n1=1。 p=p/K。 n1=n1+L。 0 50 100 150 200 250 300 3 0 2 0 1 001020300 50 100 150 200 250 300 3 0 2 0 1 00102030 圖 34 圖 35 從圖 34 和圖 33 的對比中很明顯可以看出，圖 34 頻率為 40Hz 和 90Hz 處的兩個(gè)主瓣的寬度較圖 33 而言明顯變窄，更接近于原信號。 n=0:1/Fs:1。 nfft=512。 noverlap=50。 index=0:round(nfft/2 1)。 plot_Pxx=10*log10(Pxx(index+1))。從此圖與圖 34的對比中可以看出，采用數(shù)據(jù)重疊方法后，主瓣寬度變得更窄，更接近于原信號。 %采樣頻率 n=0:1/Fs:1。 nfft=1024。 window2=bartlett(length(xn))。 window4=hamming(length(xn))。 window6=kaiser(length(xn),)。 [Pxx1,f1]=periodogram(xn,window1,nfft,Fs)。 [Pxx3,f3]=periodogram(xn,window3,nfft,Fs)。 [Pxx5,f5]=periodogram(xn,window5,nfft,Fs)。 [Pxx7,f7]=periodogram(xn,window7,nfft,Fs)。plot(f1,10*log10(Pxx1))。boxcar窗 39。 subplot(332)。title(39。)。plot(f3,10*log10(Pxx3))。blackman窗 39。 subplot(334)。title(39。)。plot(f5,10*log10(Pxx5))。hann窗 39。 subplot(336)。title(39。)。plot(f7,10*log10(Pxx7))。triang窗 39。 %plot(f,10*log10(Pxx))。title(39。)。title(39。)。title(39。)。title(39。)。title(39。)。title(39。)。title(39。)。 在上述程序中包含兩個(gè)頻率相差較小的正弦信號，但是運(yùn)用不同窗函數(shù)，估計(jì)的分辨率也不同。所以運(yùn)用矩形窗和凱撒窗進(jìn)行功率譜估計(jì)的頻率分辨率較高。 所以， 在實(shí)際應(yīng)用中，具體選擇哪種窗函數(shù)，應(yīng)該從旁瓣幅度和分辨率兩個(gè)角度綜合考慮到底選用哪一個(gè)窗函數(shù)進(jìn)行功率譜估計(jì)的效果較好。對于此方法，在進(jìn)行計(jì)算機(jī)仿真時(shí)，只需將平均周期圖法中， window參數(shù)改成相應(yīng)的窗函數(shù)即可。自相關(guān)函數(shù)由下式估計(jì)： 101? [ ] [ ] [ ] , 1NmXnR m X n X n m m NNm???? ? ? ?? ? (36) 式中， ?[]XRm的長度為 21N? 。下面討論 ?[]XRm接近 []XRm的程度。因此?[]XRm是自相關(guān)函數(shù) []XRm的無偏估計(jì)。對于零均值正態(tài)序列，有 1 2 3 4 1 2 3 4 1 3 2 4 1 4 2 3[ ] [ ] [ ] [ ] [ ] [ ] [ ]E X X X X E X X E X X E X X E X X E X X E X X? ? ? 所以 基于 MATLAB 的譜估計(jì)實(shí)現(xiàn) 18 { [ ] [ ] [ ] [ ] } [ [ ] [ ] ] [ [ ] [ ] ]E X n X n m X k X k m E X n X n m E X k X k m? ? ? ? ? ? [ [ ] [ ] ] [ [ ] [ ] ]E X n X k E X n m X k m? ? ? [ [ ] [ ] ] [ [ ] [ ] ]E X n X k m E X n m X k? ? ? 22[ ] [ ] [ ] [ ]X X X XR m R k n R k m n R k m n? ? ? ? ? ? ? ? (39) 將式 (39)代入式 (38)得 112 2 22 001?[ [ ] ] [ ] [ [ ] [ ] [ ] ]()N m N mX X X X XnkE R m R m R k n R k m n R k m nNm? ? ? ???? ? ? ? ? ? ? ?? ?? (310) 將式 (310)代入式 (37)，并注意到 22?[ [ ]] [ ]XXE R m R m? ，得 11 22 001?v a r [ [ ] ] [ [ ] [ ] [ ] ]()N m N mX X X XnkR m R k n R k m n R k m nNm? ? ? ???? ? ? ? ? ? ?? ?? 令 l k n??，顯然 l 的最小值為 ( 1)Nm? ? ? ，最大值為 ( 1)Nm??，且 0l ? （即 kn? ）的情況將出現(xiàn) Nm? 次， 1l? 的情況將出現(xiàn) 1Nm??此。即 ?lim { v a r[ [ ]]} 0XN Rm?? ? 故 ?[]XRm滿足一致估計(jì)的條件。因此可以按下式估計(jì)自相關(guān)函數(shù) 139。?[ [ ] ] [ ]XXNmE R m R mN?? 這相當(dāng)于真值 []XRm用 三 角 窗函 數(shù) 加 權(quán)。?[ [ ]] [ ]XXE R m R m? ， 39。?[]XRm的偏差為零，即 39。?[]XRm是 []XRm的漸近無偏估 計(jì)，其估計(jì)方差 239。?[]XRm是 []XRm的漸近無偏一致估計(jì)，且其估計(jì)方差小于 ?[]XRm的估計(jì)方差，所以一般用 39。所以 ?[]XRm可表示為： 101? [ ] [ ] [ ]NmXnR m X n X n mN?????? BT 法的仿真 仿真程序如下： Fs=500。 %產(chǎn)生含有噪聲的序列 xn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+randn(size(n))。 cxn=xcorr(xn)。%CXk估計(jì)出來以后是復(fù)數(shù)，此處是對它取模 index=0:round(nfft/2 1)。%k相當(dāng)于頻率 基于 MATLAB 的譜估計(jì)實(shí)現(xiàn) 20 plot_Pxx=10*log10(Pxx(index+1))。 plot(k,plot_Pxx)。如圖 39 所示，當(dāng)采樣點(diǎn)數(shù)為 5000 時(shí)，主瓣寬度明顯變窄，旁瓣的起伏程度也較小。解決的 方法也可以在時(shí)域采用不同的窗函數(shù)，從而使得分辨率有所提高。如果能根據(jù)已觀察到的數(shù)據(jù)估計(jì)出這一物理網(wǎng)絡(luò)的模型參數(shù)，就不必認(rèn)為 N 個(gè)以外的數(shù)據(jù)全為零，這就有可能克服經(jīng)典譜估計(jì)的缺點(diǎn)。這個(gè)線性系統(tǒng)的傳遞函數(shù)為 00()()()qkkk