【正文】 學(xué)生姓名： 鈕金鑫 學(xué)號： 05052112 指導(dǎo)教師： 邢務(wù)強(qiáng) 報告日期： 2021 年 3月 16日 1． 本課題所涉及的問題及應(yīng)用現(xiàn)狀綜述 功率譜估計 (PSD) 是用有限長的數(shù)據(jù)來估計信號的功率譜 , 它對于認(rèn)識一個隨機(jī)信號或其他應(yīng)用方面來講都是重要的 , 是數(shù)字信號處理的重要研究內(nèi)容之一，在軍事，生物醫(yī)學(xué)，通信等領(lǐng)域得到了較為廣泛的應(yīng)用。 第 15 周 總結(jié)、 完成畢業(yè)設(shè)計論文。 第 11 周 進(jìn)行現(xiàn)代譜估計系統(tǒng)仿真。 第 7 周 研究經(jīng)典譜估計的方法。 西 安 郵 電 學(xué) 院 畢 業(yè) 設(shè) 計（論 文） 題 目： 基于 MATLAB 的譜估計實現(xiàn) 系 別： 通信工程 專 業(yè)： 通信工程 班 級： 通工 0509 班 學(xué)生姓名： 鈕金鑫 導(dǎo)師姓名： 邢務(wù)強(qiáng) 職稱： 講師 起 止時間： 2021 年 3 月 10 日至 2021 年 6 月 11 日 西 安 郵 電 學(xué) 院 畢業(yè)設(shè)計 (論文 )任務(wù)書 學(xué)生姓名 鈕金鑫 指導(dǎo)教師 邢務(wù)強(qiáng) 職稱 講師 系別 通信工程系 專業(yè) 通信工程 題目 基于 Matlab 的譜估計實現(xiàn) 任務(wù)與要求 掌握隨機(jī)信號功率譜的相關(guān)理論及其性質(zhì)，在此基礎(chǔ)上討論功率譜估計的各種原理方法，比較分析功率譜估計的各自特點。 第 6 周 閱讀相關(guān)參考文獻(xiàn) 。 第 10 周 研究現(xiàn)代經(jīng)典譜估計 的方法。 第 14 周 對畢業(yè)論文進(jìn)行修改，并進(jìn)行后期檢查。 對計劃的說明 本計劃為開題之初所定，后續(xù)會根據(jù)具體情況隨時調(diào)整。 后者的主要方法有最大熵譜分析法 (AR 模型法 )、 Pisarenko 諧波分解法、 Prony 提取極點法、 Prony 譜線分解法以及 Capon 最大似然法。 （ 3）， 學(xué)習(xí) MATLAB 的相關(guān)知識。 說明： 本報告必須由承擔(dān)畢業(yè)論文 (設(shè)計 )課題任務(wù)的學(xué)生在畢業(yè)論文 (設(shè)計 ) 正式 開始的第 1 周周五之前獨立撰寫完成，并交指導(dǎo)教師審閱。 第 9～ 10周 學(xué)習(xí) MATLAB 的相關(guān)知識。 第 15 周 總結(jié)、 完成畢業(yè)設(shè)計論文。 本論文研究了功率譜估計的幾種常用的方法，包括經(jīng)典法和參數(shù)法中的 AR 模型法，對每種方法的估計質(zhì)量做了數(shù)學(xué)推導(dǎo)，給出仿真程序及仿真圖，在仿真圖的基礎(chǔ)上對每種方法的性能進(jìn)行了討論。 關(guān)鍵字：功率譜估計；周期圖； BT 法； AR模型法。功率譜估計在各種隨機(jī)信號處理中得到了十分 廣泛的應(yīng)用。非參數(shù)化方法又叫經(jīng)典譜估計法，它實質(zhì)上仍依賴于傳統(tǒng)的傅里葉變換法。 經(jīng)典法的缺點表現(xiàn)為，除了得到的 N 個數(shù)據(jù)以外，序列的其他值均被認(rèn)為是零（或者等效，序列的自相關(guān)函數(shù)值除了能估計出的有限個值之外，其他的值被當(dāng)做零），但序列或其自相關(guān)函數(shù)的那些未能觀測到或未估計出來的值，實際上并不全是零。 （ 2） 離散型隨機(jī)過程： ()Xt 對于任意的 tT? ， ()Xt 都是離散型隨機(jī)變量，也就是時間連續(xù)，狀態(tài)離散的情況。例如，對連續(xù)型隨機(jī)序列再進(jìn)行量化，即得到離散隨機(jī)序列。 )nXE g X g x f x n d x????? ? 均值有下列性質(zhì)： （ 1） [ ] [ ] [ ]n m n mE X Y E X E Y? ? ?，即和的均值等于均值的和。 )nX n XE X x f x n d x??????? ? 離散時間隨機(jī)過程的方差定義為 22[ ] [ ( ) ]nnX n n XD X E X m? ? ? ? 由于和的均值等于均值的和，所以容易證明上式可寫成 2 2 2 2 2[ ] [ ]n n nX n n X XE X E X m??? ? ? ? 2nX? 為非負(fù)函數(shù)，其平方根稱作離散時間隨機(jī)過程的標(biāo)準(zhǔn)差或均方差，即 2 []nnX X nDX???? 一般來說，均值，均方值， 方差都是 n 的函數(shù)，但對平穩(wěn)離散隨機(jī)過程來說，它們與 n 無關(guān)，都是常數(shù)，即 []Xnm E X? 22[]XnEX? ? 22[( ) ]X n XE X m? ?? 離散時間隨機(jī)過程的自相關(guān)函數(shù)定義為 121 2 1 2 1 2 1 2 1 2( , ) [ ] ( , 。 3， 對于實隨機(jī)過程來說，功率譜密度是 ? 的偶函數(shù)。 估計量的評價標(biāo)準(zhǔn) 無偏性 設(shè) 12, , , nX X X… 是總體 X 的一個樣本。 有效性 設(shè) 1 1 1 2? ? ( , , , )nX X X??? …與 2 2 1 2? ? ( , , , )nX X X??? …都是 ? 的無偏估計量，若對于任意??? ，有 12? ?[ ] [ ]DD??? 且至少對于某一個 ??? 上式中的不等號成立，則稱 1?? 較 2?? 有效。 設(shè) 12? ?( , , , )nX X X??? … 為參 數(shù) ? 的估計量，若對于任意 ??? ，當(dāng) n?? 時12?( , , , )nX X X? … 依概率收斂于 ? ，則稱 ?? 是 ? 的相合估計量。 3 經(jīng)典譜估計及其仿真 經(jīng)典譜估計方法分為兩種： BT 法和周期圖法。而 []Ym的傅里葉變換 10[ ] [ ] [ ] [ ] ( )Nj m j m j m jNm m mY m e d m X m e X m e X e? ? ? ?? ? ? ? ?? ? ?? ? ? ? ? ? ?? ? ?? ? ? ()jNXe? 是有限長序列 []Xn的傅里葉變換，顯然 ()jNXe? 是周期性的。 周期圖的性能 為了了解周期圖法的譜估計效果，我們來討論它的估計均值和方差。 為分析簡單起見，通常假設(shè) []Xn零均值，方差為 2X? 的實高斯白噪聲序列，即功率譜密度為常數(shù) 2X? 。利用正態(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 的譜估計實現(xiàn) 9 [ [ ] [ ] ] [ [ ] [ ] ]E X n X q E X k X p? 40X????? , 。顯然，周期圖不是功率譜估計的一致估計，所以這種估計不是最好的估計方法。 加窗周期圖 周期圖法相當(dāng)于只取 了原信號的一段進(jìn)行估計，也就是說對原信號加了一個矩形窗。如果原來真是的功率譜是窄的，那么與主瓣卷積后會使功率向附近頻域擴(kuò) 散，卷積以后使得加窗的信號的頻譜變成 Sa 函數(shù)頻譜的形狀，即有一定寬度，存在若干旁瓣的光滑脈沖那樣。 旁瓣較低的窗是 blackman窗和 hamming窗 ， 但這 兩種窗的主瓣寬度比矩形窗的主瓣寬度大，因此分辨率沒有矩形窗的分辨率高。 n=0:1/Fs:1。 figure(1) plot(f,10*log10(Pxx))。 基于 MATLAB 的譜估計實現(xiàn) 13 平均周期圖的仿真 仿真程序如下： Fs=600。 Pxx=bart(xn,10)。 figure(1) plot(k,plot_Pxx)。 n1=1。 n1=n1+L。 n=0:1/Fs:1。 noverlap=50。 plot_Pxx=10*log10(Pxx(index+1))。 %采樣頻率 n=0:1/Fs:1。 window2=bartlett(length(xn))。 window6=kaiser(length(xn),)。 [Pxx3,f3]=periodogram(xn,window3,nfft,Fs)。 [Pxx7,f7]=periodogram(xn,window7,nfft,Fs)。boxcar窗 39。title(39。plot(f3,10*log10(Pxx3))。 subplot(334)。)。hann窗 39。title(39。plot(f7,10*log10(Pxx7))。 %plot(f,10*log10(Pxx))。)。)。)。)。)。)。)。所以運用矩形窗和凱撒窗進(jìn)行功率譜估計的頻率分辨率較高。對于此方法，在進(jìn)行計算機(jī)仿真時，只需將平均周期圖法中， window參數(shù)改成相應(yīng)的窗函數(shù)即可。下面討論 ? []XRm接近 []XRm的程度。對于零均值正態(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??此。因此可以按下式估計自相關(guān)函數(shù) 139。?[ [ ] ] [ ]XXE R m R m? ， 39。? []XRm是 []XRm的漸近無偏估 計，其估計方差 239。所以 ? []XRm可表示為： 101? [ ] [ ] [ ]NmXnR m X n X n mN?????? BT 法的仿真 仿真程序如下： Fs=500。 cxn=xcorr(xn)。%k相當(dāng)于頻率 基于 MATLAB 的譜估計實現(xiàn) 20 plot_Pxx=10*log10(Pxx(index+1))。如圖 39所示，當(dāng)采樣點數(shù)為 5000 時，主瓣寬度明顯變窄，旁瓣的起伏程度也較小