【導讀】率譜估計的各種原理方法，比較分析功率譜估計的各自特點。數選取，最后通過實例加以說明。第5周熟悉MATLAB的相關基本知識。第6周閱讀相關參考文獻。第8周進行經典譜估計系統(tǒng)仿真。第10周研究現代經典譜估計的方法。第13周總結以上討論結果，寫出論文初稿。第14周對畢業(yè)論文進行修改，并進行后期檢查。第16周完成畢業(yè)設計答辯。并指出其中的重點難點，檢查上次布置的任務。本計劃為開題之初所定，后續(xù)會根據具體情況隨時調整。內容之一，在軍事，生物醫(yī)學，通信等領域得到了較為廣泛的應用。前者的主要方法有BTPSD估計法和周期圖法;后者的主要方法有最大熵。其中周期圖法和AR模型法是用得較多且最。供用于功率譜估計的函數。的了解，在大學期間，對隨機信號分析，MATLAB這兩門課程已有初步了解，開始的第1周周五之前獨立撰寫完成，并交指導教師審閱。晰的認識，故應復習隨機信號基本理論和MATLAB語言的使用。

　　

【正文】 ,window6,nfft,Fs)。 [Pxx7,f7]=periodogram(xn,window7,nfft,Fs)。 figure(1) subplot(331)。plot(f1,10*log10(Pxx1))。title(39。boxcar窗 39。)。 subplot(332)。plot(f2,10*log10(Pxx2))。title(39。bartlett窗 39。)。 subplot(333)。plot(f3,10*log10(Pxx3))。title(39。blackman窗 39。)。 subplot(334)。plot(f4,10*log10(Pxx4))。title(39。hamming窗 39。)。 基于 MATLAB 的譜估計實現 15 subplot(335)。plot(f5,10*log10(Pxx5))。title(39。hann窗 39。)。 subplot(336)。plot(f6,10*log10(Pxx6))。title(39。kaiser窗 39。)。 subplot(337)。plot(f7,10*log10(Pxx7))。title(39。triang窗 39。)。 %plot(f,10*log10(Pxx))。 figure(2) plot(f1,10*log10(Pxx1))。title(39。boxcar窗 39。)。 figure(3) plot(f2,10*log10(Pxx2))。title(39。bartlett窗 39。)。 figure(4) plot(f3,10*log10(Pxx3))。title(39。blackman窗 39。)。 figure(5) plot(f4,10*log10(Pxx4))。title(39。hamming窗 39。)。 figure(6) plot(f5,10*log10(Pxx5))。title(39。hann窗 39。)。 figure(7) plot(f6,10*log10(Pxx6))。title(39。kaiser窗 39。)。 figure(8) plot(f7,10*log10(Pxx7))。title(39。triang窗 39。)。 仿真圖如下： 圖 36 從圖 36可以看出，采用 blackman窗和 hamming窗 后，功率譜的主瓣寬度比較窄，頻率泄露的情況相對不太明顯（也就是除主瓣以外的譜線的幅度較主瓣幅度相比較?。? 在上述程序中包含兩個頻率相差較小的正弦信號，但是運用不同窗函數，估計 的分辨率也不同。對圖 36的主瓣進行放大后如圖 37所示 . 基于 MATLAB 的譜估計實現 16 圖 37 從圖 37可以看出，只有運用矩形窗和凱撒窗才能分辨出這兩個頻率相差較小的信號。所以運用矩形窗和凱撒窗進行功率譜估計的頻率分辨率較高。但是在實際應用中，幅度相差 100倍甚至更高的正弦分量同時存在于一個信號中是可能的，這種情況下，由于矩形窗的旁瓣幅度較大，難以區(qū)分這樣的兩個正弦分量，故估計的效果較差，故還應考慮旁瓣幅度較小的窗函數。 所以，在實際應用中，具體選擇哪種窗函數，應該從旁瓣幅度和分辨率兩個角度綜合考慮到底選用哪一個窗函數進行功率 譜估計的效果較好。 在采用周期圖法進行功率譜估計時，可以將平均周期圖法和加窗周期圖法這兩種方法結合在一起，效果會更好。對于此方法，在進行計算機仿真時，只需將平均周期圖法中， window參數改成相應的窗函數即可。 BT 法 BT 法 BT 法首先計算隨機信號 []Xn 的自相關函數 [ ] [ ( ) ( ) ]XR m E X n X n m??，若[]Xm Rm????? ??? ，可由 []XRm的傅里葉變換的得到它的功率譜，即： ( ) [ ] jmXXmG R m e ???? ?? ??? ? BT 法的估計步驟如下： 基于 MATLAB 的譜估計實現 17 1，先估計隨機信號 []Xn的自相關函數。自相關函數由下式估計： 101? [ ] [ ] [ ] , 1NmXnR m X n X n m m NNm???? ? ? ?? ? (36) 式中， ? []XRm的長度為 21N? 。 2，求 ? []XRm的離散傅里葉變換，得到 []Xn功率譜的估值： ( ) [ ] , 1M jmXXmMG R m e M N??? ???? ? ?? 由于式 (36)僅利用了 []Xn的 N 個有限值得到自相關函數的估計 ? []XRm，因 此它與信號本身的自相關函數 []XRm會有一定程度的差別。下面討論 ? []XRm接近 []XRm的程度。 首先求 ? []XRm的均值，有 101?[ [ ] ] { [ ] [ ] }NmXnE R m E X n X x mNm?????? ? 101 ? [ ] [ ]Nm XXn R m R mNm?????? ? 即自相關函 數的估計值 ? []XRm的數學期望等于序列 []Xn的自相關函數的真實值。因此? []XRm是自相關函數 []XRm的無偏估計。 其次求 ? []XRm的方差，根據方差的定義 22? ? ?v a r [ [ ] ] [ [ ] ] [ [ ] ]X X XR m E R m E R m?? (37) 由式（ 36）可知，上式等號右邊第一項可表示為 1122 001?[ [ ] ] { [ ] [ ] [ ] [ ] }()N m N mX nkE R m E X n X n m X k X k mNm? ? ? ???? ? ?? ?? (38) 當 []Xn是零均值正態(tài)序列時，它的各高階矩都可以用其一階和二階矩來表示。對于零均值正態(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 的譜估計實現 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? ）的情況將出現 Nm? 次， 1l? 的情況將出現 1Nm??此。依此類推，對不同的 l 的情況，出現的次數將為 ()N m l??，于是上式可寫成 1 22 ( 1 )1?v a r [ [ ] ] { ) [ [ ] [ ] [ ] }()NmX X X Xl N mR m N m l R l R l m R l mNm??? ? ? ?? ? ? ? ? ?? ? 1 22 ( 1 )1 1 [ [ ] [ ] [ ] ]()NmX X Xl N m ml R l R l m R l mN m N??? ? ? ?? ? ?? ? ? ? ??????? 1 22 ( 1 ) [ [ ] [ ] [ ] ]()NmX X Xl N mN R l R l m R l mNm??? ? ? ?? ? ? ?? ? 當 Nm?? 時，上式以 1N 趨于零。即 ?lim {v a r[ [ ]]} 0XN Rm?? ? 故 ? []XRm滿足一致估計的條件。 基于 MATLAB 的譜估計實現 19 式 (36)估計自相關函數的方法，雖然當 N 較大且 Nm?? 時能得到一致估計，但當 m 接近 N 時， ? []XRm的估計方差變大，因而不能得到有用的估計。因此可以按下式估計自相關函數 139。01? ?[ ] [ ] [ ] [ ] , 1NmXXnNmR m X n X n m R m m NNN????? ? ? ? ?? 其均值 39。?[ [ ] ] [ ]XXNmE R m R mN?? 這相當于真值 []XRm用 三 角 窗函 數 加 權。 上 式 表明 ： 當 0m? 時，由于39。?[ [ ] ] [ ]XXE R m R m? ， 39。? []XRm是 []XRm的有偏估計；若 m 取有限值，則當 N?? 時，39。?[]XRm的偏差為零，即 39。?lim [ [ ] ] [ ] 0XXN E R m R m?? ??， 39。? []XRm是 []XRm的漸近無偏估 計，其估計方差 239。?v a r [ [ ] ] v a r [ [ ] ] v a r [ [ ] ]X X XNmR m R m R mN? ? ??????? 由于 39。? []XRm是 []XRm的漸近無偏一致估計，且其估計方差小于 ? []XRm的估計方差，所以一般用 39。? []XRm作為自相關函數的估計，而較少使用 ? []XRm。所以 ? []XRm可表示為： 101? [ ] [ ] [ ]NmXnR m X n X n mN?????? BT 法的仿真 仿真程序如下： Fs=500。%采樣頻率 n=0:1/Fs:1。 %產生含有噪聲的序列 xn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+randn(size(n))。 nfft=512。 cxn=xcorr(xn)。 %計算序列的自相關函數 CXk=fft(cxn,nfft) Pxx=abs(CXk)。%CXk估計出來以后是復數，此處是對它取模 index=0:round(nfft/2 1)。 k=index*Fs/nfft。%k相當于頻率 基于 MATLAB 的譜估計實現 20 plot_Pxx=10*log10(Pxx(index+1))。 figure(1)。 plot(k,plot_Pxx)。 仿真圖如圖 38 所示 0 50 100 150 200 25020253035404550550 50 100 150 200 25025303540455055 圖 38（采樣點數為 500） 圖 39(采樣點數為 5000) 從前述討論可知， BT 法可通過增加采樣點數 N來使估計的方差減小，從而使性能得到改進。如圖 39所示，當采樣點數為 5000 時，主瓣寬度明顯變窄，旁瓣的起伏程度也較小。 同周期圖法一樣， BT 法也存在譜分辨率的問題。解決的方法也可以在時域采用不同的窗函數，從而使得分辨率有所提高。