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

　　

【正文】 ubplot(334)。plot(f4,10*log10(Pxx4))。title(39。hamming窗 39。)。 基于 MATLAB 的譜估計實現(xiàn) 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窗 后，功率譜的主瓣寬度比較窄，頻率泄露的情況相對不太明顯（也就是除主瓣以外的譜 線的幅度較主瓣幅度相比較?。?。 在上述程序中包含兩個頻率相差較小的正弦信號，但是運用不同窗函數(shù)，估計的分辨率也不同。對圖 36的主瓣進行放大后如圖 37所示 . 基于 MATLAB 的譜估計實現(xiàn) 16 圖 37 從圖 37可以看出，只有運用矩形窗和凱撒窗才能分辨出這兩個頻率相差較小的信號。所以運用矩形窗和凱撒窗進行功率譜估計的頻率分辨率較高。但是在實際應(yīng)用中，幅度相差 100倍甚至更高的正弦分量同時存在于一個信號中是可能的，這種情況下，由于矩形窗的旁瓣幅度較大，難以區(qū)分這樣的兩個正弦分量，故估計的效果較差，故還應(yīng)考慮旁瓣幅度較小的窗函數(shù)。 所以， 在實際應(yīng)用中，具體選擇哪種窗函數(shù)，應(yīng)該從旁瓣幅度和分辨率兩個角度綜合考慮到底選用哪一個窗函數(shù)進行功率譜估計的效果較好。 在采用周期圖法進行功率譜估計時，可以將平均周期圖法和加窗周期圖法這兩種方法結(jié)合在一起，效果會更好。對于此方法，在進行計算機仿真時，只需將平均周期圖法中， window參數(shù)改成相應(yīng)的窗函數(shù)即可。 BT 法 BT 法 BT 法首先計算隨機信號 []Xn的自相關(guān)函數(shù) [ ] [ ( ) ( ) ]XR m E X n X n m??，若[]Xm Rm????? ??? ，可由 []XRm的傅里葉變換的得到它的功率譜，即： ( ) [ ] jmXXmG R m e ???? ????? ? BT 法的估計步驟如下： 基于 MATLAB 的譜估計實現(xiàn) 17 1，先估計隨機信號 []Xn的自相關(guān)函數(shù)。自相關(guān)函數(shù)由下式估計： 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 個有限值得到自相關(guān)函數(shù)的估計 ?[]XRm，因此它與信號本身的自相關(guān)函數(shù) []XRm會有一定程度的差別。下面討論 ?[]XRm接近 []XRm的程度。 首先求 ?[]XRm的均值，有 101?[ [ ] ] { [ ] [ ] }NmXnE R m E X n X x mNm?????? ? 101 ? [ ] [ ]Nm XXn R m R mNm?????? ? 即自相關(guān)函數(shù)的估計值 ?[]XRm的數(shù)學期望等于序列 []Xn的自相關(guān)函數(shù)的真實值。因此?[]XRm是自相關(guān)函數(shù) []XRm的無偏估計。 其次求 ?[]XRm的方差，根據(jù)方差的定義 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 的譜估計實現(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??此。依此類推，對不同的 l 的情況，出現(xiàn)的次數(shù)將 為 ()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 的譜估計實現(xiàn) 19 式 (36)估計自相關(guān)函數(shù)的方法，雖然當 N 較大且 Nm?? 時能得到一致估計，但當 m 接近 N 時， ?[]XRm的估計方差變大，因而不能得到有用的估計。因此可以按下式估計自相關(guān)函數(shù) 139。01? ?[ ] [ ] [ ] [ ] , 1NmXXnNmR m X n X n m R m m NNN????? ? ? ? ?? 其均值 39。?[ [ ] ] [ ]XXNmE R m R mN?? 這相當于真值 []XRm用 三 角 窗函 數(shù) 加 權(quán)。 上 式 表明 ： 當 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作為自相關(guān)函數(shù)的估計，而較少使用 ?[]XRm。所以 ?[]XRm可表示為： 101? [ ] [ ] [ ]NmXnR m X n X n mN?????? BT 法的仿真 仿真程序如下： Fs=500。%采樣頻率 n=0:1/Fs:1。 %產(chǎn)生含有噪聲的序列 xn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+randn(size(n))。 nfft=512。 cxn=xcorr(xn)。 %計算序列的自相關(guān)函數(shù) CXk=fft(cxn,nfft) Pxx=abs(CXk)。%CXk估計出來以后是復數(shù)，此處是對它取模 index=0:round(nfft/2 1)。 k=index*Fs/nfft。%k相當于頻率 基于 MATLAB 的譜估計實現(xiàn) 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（采樣點數(shù)為 500） 圖 39(采樣點數(shù)為 5000) 從前述討論可知， BT法可通過增加采樣點數(shù) N 來使估計的方差減小，從而使性能得到改進。如圖 39 所示，當采樣點數(shù)為 5000 時，主瓣寬度明顯變窄，旁瓣的起伏程度也較小。 同周期圖法一樣， BT法也存在譜分辨率的問題。解決的 方法也可以在時域采用不同的窗函數(shù)，從而使得分辨率有所提高。 基于 MATLAB 的譜估計實現(xiàn) 21 4 AR 模型法 AR 模型法的基本理論 任何具有功率譜密度的隨機信號都可以看成由一白噪聲 []n? 激勵一物理網(wǎng)絡(luò)所形成。如果能根據(jù)已觀察到的數(shù)據(jù)估計出這一物理網(wǎng)絡(luò)的模型參數(shù)，就不必認為 N 個以外的數(shù)據(jù)全為零，這就有可能克服經(jīng)典譜估計的缺點。 實際應(yīng)用中所遇到的隨機過程常?？梢钥闯捎砂自肼?[]n? 經(jīng)意線性系統(tǒng)形成的。這個線性系統(tǒng)的傳遞函數(shù)為 00()()()qkk