【正文】 圖表整潔，布局合理，文字注釋必須使用工程字書寫，不準(zhǔn)用徒手畫3）畢業(yè)論文須用 A4 單面打印，論文 50 頁以上的雙面打印4）圖表應(yīng)繪制于無格子的頁面上5）軟件工程類課題應(yīng)有程序清單，并提供電子文檔1）設(shè)計(jì)（論文）2）附件：按照任務(wù)書、開題報(bào)告、外文譯文、譯文原文（復(fù)印件）次序裝訂3）其它。：任務(wù)書、開題報(bào)告、外文譯文、譯文原文（復(fù)印件） 。保密的論文（設(shè)計(jì)）在解密后適用本規(guī)定。有 權(quán)將論文（設(shè)計(jì)）用于非贏利目的的少量復(fù)制并允許論文（設(shè)計(jì)）進(jìn)入學(xué)校圖書館被查閱。對本論文（設(shè)計(jì)）的研究做出重要貢獻(xiàn)的個(gè)人和集體，均已在文中作了明確說明并表示謝意。北京:化學(xué)工業(yè)出版社，2022.(3)丁玉美等，數(shù)字信號處理——時(shí)域離散信號處理[M].西安:西安電子科技大學(xué)出版社，2022.(4)盛驟等，概率論與數(shù)理統(tǒng)計(jì)[M].北京:高等教育出版社，2022.(5)姚武川，姚天任，經(jīng)典譜估計(jì)方法的 MATLAB 分析[J] .華中理工大學(xué)學(xué)報(bào)，2022，28(4).(6)黃志宇等，隨機(jī)信號的功率譜估計(jì)及 MATLAB 的實(shí)現(xiàn)[J] .現(xiàn)代電子技術(shù)，2022，3(134).基于 MATLAB 的譜估計(jì)實(shí)現(xiàn)33畢業(yè)論文（設(shè)計(jì)）原創(chuàng)性聲明本人所呈交的畢業(yè)論文（設(shè)計(jì)）是我在導(dǎo)師的指導(dǎo)下進(jìn)行的研究工作及取得的研究成果。基于 MATLAB 的譜估計(jì)實(shí)現(xiàn)31致謝我要感謝邢老師和我的同學(xué)們，他們在一些數(shù)學(xué)推導(dǎo)以及 MATLAB 的仿真上給我提供了有力幫助。AR 模型法有效地避免了經(jīng)典譜估計(jì)的缺點(diǎn)，且在分辨率等性能上能夠獲得較好效果。周期圖法直接對觀測數(shù)據(jù)進(jìn)行傅里葉變換，從而求得功率譜估計(jì)，雖然效率較高，但它不是功率譜密度的一致估計(jì)，且估計(jì)性能較差，所以需要進(jìn)行平均修正。從本篇論文可以看出，BT 法先利用觀測數(shù)據(jù)求出信號自相關(guān)函數(shù)的估計(jì)量，然后對其進(jìn)行傅里葉變換，得出功率譜估計(jì)。因此，在缺乏先驗(yàn)信息的情況下，可以試著用不同方法確定模型的階，并對最終的譜估計(jì)結(jié)果進(jìn)行分析和選擇。但遺憾的是這個(gè)函數(shù)往往不能給出理想的 AR 模型的階數(shù)。如圖 47，估計(jì)選擇的階數(shù)是 80，雖然頻率為 40Hz 和 41Hz 的兩個(gè)譜峰較為明顯，但是其他頻率譜峰的起伏性也明顯增加，因此，增加了虛假的頻率成分。如圖 45，估計(jì)選擇的是階數(shù)為 50，估計(jì)的信號接近真實(shí)信號，可以分辨出頻率為 40Hz 和 41Hz 的兩個(gè)譜峰。當(dāng)選擇的模型階次較大時(shí)會出現(xiàn)虛假的頻率成分。一般模型的最好選擇是先p驗(yàn)未知的，實(shí)際中需預(yù)先選定模型的階次。)%half表示單邊譜基于 MATLAB 的譜估計(jì)實(shí)現(xiàn)280 10 20 30 40 50 60 70 80 90 1004035302520151050Frequency (Hz)Power/frequency (dB/Hz)Power Spectral Density Estimate via Burg圖 44從圖 44 可以看出，burg 遞推法估計(jì)出的譜線在寬度和幅度上面的性能都要優(yōu)于 YuleWalker 方法。%order1就是AR模型數(shù)學(xué)公式中的Pfigure(1)pburg(xn,order1,nfft,fs,39。nfft=512。xn=cos(2*pi*40*n)+cos(2*pi*41*n)+3*cos(2*pi*90*n)+*randn(size(n))。burg遞推法的仿真程序如下:fs=200。所以這兩種方法適合于估計(jì)正弦波的頻率，但產(chǎn)生的虛假譜峰也比較多。%AR模型協(xié)方差算法subplot(212)。subplot(211)。nfft=512。n=0:1/fs:1。 協(xié)方差方法與修正協(xié)方差方法協(xié)方差方法和修正協(xié)方差方法克服了 YuleWalker 方法中存在的頻率偏移的問題且在分辨率方面有良好的性能，但虛假譜峰的問題依然沒有較好的解決。YuleWalker方法的仿真圖如下：0 10 20 30 40 50 60 70 80 90 1003025201510505Frequency (Hz)Power/frequency (dB/Hz)Power Spectral Density Estimate via YuleWalker 0 10 20 30 40 50 60 70 80 90 1002520151050510Frequency (Hz)Power/frequency (dB/Hz)Power Spectral Density Estimate via YuleWalker 圖 41（采樣點(diǎn)為 50） 圖 42（采樣點(diǎn)為 200）從圖 41 可以看出，當(dāng)采樣點(diǎn)為 50 時(shí)，在頻率為 40Hz 處的譜峰明顯向左偏移，且在頻率為 6070Hz 之間出現(xiàn)了虛假譜峰。nfft=512。n=0:1/fs:1。1pN??為了求出 最小時(shí)的模型參數(shù)，可先求出參數(shù) ，再根據(jù)解 YuleWalke 方程p?a中導(dǎo)出的遞推公式，進(jìn)一步解出其他參數(shù)，從而可根據(jù)觀測數(shù)據(jù)直接計(jì)算 AR 模型的參數(shù)。在這兩個(gè)公式中，信號項(xiàng)的自變量最大的是 ，最小的是 。伯格遞推法則由信號的觀測數(shù)據(jù)直接計(jì)算 AR 模型的參數(shù)。再根據(jù)最小均方誤差準(zhǔn)則可求得修正協(xié)方差方法中的模型參數(shù)。其預(yù)測誤差的平均功率為基于 MATLAB 的譜估計(jì)實(shí)現(xiàn)25()ppfb????式中， 和 分表表示前向和后向預(yù)測誤差功率。由于使用的 個(gè)數(shù)()pN?據(jù)都是已知的，對 以外的數(shù)據(jù)不需要做任何假設(shè)，相當(dāng)于避免了加窗概念的引入，p所以得到的預(yù)測誤差只存在于 之間，所以，預(yù)測誤差功率的公式如下式：~1npN?? 21121[][][]pN knpnpeXan?????????運(yùn)用最小均方誤差準(zhǔn)則可求得模型參數(shù)。在 YuleWalker 中，沒有用到前、后向預(yù)測的概念，所以由觀測數(shù)據(jù)求預(yù)測誤差功率的公式為： 22111[][][]pinnieXanN????????????假設(shè)信號 的數(shù)據(jù)區(qū)在 范圍，有 個(gè)預(yù)測系數(shù)， 個(gè)數(shù)據(jù)經(jīng)過沖激[]Xn0?N相應(yīng)為 的濾波器，輸出預(yù)測誤差 的長度為 ,應(yīng)用下式計(jì)算：(0,12,)piap?…[]eP? 2112022[][][]pNPNPinnieXan?????????顯然， 的長度長于數(shù)據(jù)的長度，上式中數(shù)據(jù) 的兩端需補(bǔ)充零點(diǎn)，這相當(dāng)于[]en []無窮長的信號經(jīng)過加窗處理，得到長度為 的數(shù)據(jù)。,39。[][]39。[]pkXnan????這樣，前向預(yù)測是由 預(yù)測 ，后向[],[],[2],[1]Xn?…?[]Xn預(yù)測是由 預(yù)測 。[]pka???一般，前向，后向預(yù)測用同一數(shù)據(jù)進(jìn)行，即利用 [],1],[2],Xnn?…進(jìn)行預(yù)測，為此，將上式改為[1]Xnp?? ,1?39。[] ,1?39。用均方誤差最小準(zhǔn)則求前向預(yù)測誤差濾[1][2],[]Xnn??…波器的最佳系數(shù) ,pk 2,{[]}0,1,fpkEepa??…從而求得估計(jì)系數(shù) 。前向預(yù)測就是指利用 預(yù)測 ，其估計(jì)值 和[1],[2],[]XnXnp??…[]n?[]Xn預(yù)測誤差 可用下式表示[]fpen,1?[][]pknan??? (47),1?[][][]pfp keXaXn???其中 就是指預(yù)測系數(shù)與階數(shù) 有關(guān)。對于實(shí)序列，由于 ，因此只[]XRm?要已知或估計(jì)出 個(gè)自相關(guān)函數(shù)的值，可由該方程解出 個(gè)模型參數(shù)? 1?，根據(jù)這些參數(shù)即可得到隨機(jī)信號的功率譜估計(jì)。2??ka[]Xn AR 模型的 YuleWalker 方法階 AR 模型的系統(tǒng)函數(shù)和差分方程分別如式(41),(42) 所示，模型輸出的功p率譜密度如式(43)所示。其系統(tǒng)函數(shù)只有極點(diǎn)，沒有零點(diǎn)。2??kab若 ， ，除 外其余的 都為 0，此時(shí)系統(tǒng)函數(shù)就變?yōu)?1b?00k (41)1()pkkBzHAaz????其差分方程為 (42)1[][][pkXnan????其中， 為白噪聲通過物理網(wǎng)絡(luò) 后的輸出信號，也就是待估計(jì)的信號。N實(shí)際應(yīng)用中所遇到的隨機(jī)過程常?？梢钥闯捎砂自肼?經(jīng)意線性系統(tǒng)形成的?；?MATLAB 的譜估計(jì)實(shí)現(xiàn)214 AR 模型法 AR 模型法的基本理論任何具有功率譜密度的隨機(jī)信號都可以看成由一白噪聲 激勵一物理網(wǎng)絡(luò)所[]n?形成。同周期圖法一樣，BT 法也存在譜分辨率的問題。仿真圖如圖 38 所示0 50 100 150 200 25020253035404550550 50 100 150 200 25025303540455055 圖 38（采樣點(diǎn)數(shù)為 500） 圖 39(采樣點(diǎn)數(shù)為 5000)從前述討論可知，BT 法可通過增加采樣點(diǎn)數(shù) N 來使估計(jì)的方差減小，從而使性能得到改進(jìn)。figure(1)?；?MATLAB 的譜估計(jì)實(shí)現(xiàn)20k=index*Fs/nfft。 %計(jì)算序列的自相關(guān)函數(shù)CXk=fft(cxn,nfft)Pxx=abs(CXk)。nfft=512。%采樣頻率n=0:1/Fs:1。所以 可表示為：39。?var[]var[]var[]XXX????????由于 是 的漸近無偏一致估計(jì)，且其估計(jì)方差小于 的估計(jì)方差，39。?li[]0NERm????39。?[]X N??的偏差為零，即 ， 是 的漸近無偏估39。上式表明：當(dāng) 時(shí)，由于[]Xm0?， 是 的有偏估計(jì)；若 取有限值，則當(dāng) 時(shí)，39。 0? ?[][]][],1NX Xn RN????????其均值 39。?[]XR基于 MATLAB 的譜估計(jì)實(shí)現(xiàn)19式(36)估計(jì)自相關(guān)函數(shù)的方法，雖然當(dāng) 較大且 時(shí)能得到一致估計(jì)，Nm?但當(dāng) 接近 時(shí)， 的估計(jì)方差變大，因而不能得到有用的估計(jì)。依此類推，對不同?1l的 的情況，出現(xiàn)的次數(shù)將為 ，于是上式可寫成l Nm122()?var[]{)[][][]}()X XXlNRmlRlmRl?? ??? 122()[][][]()mXXlNlll???????? 122()[][][]()mXXlNRllRlm?????當(dāng) 時(shí)，上式以 趨于零。?[]XRm其次求 的方差，根據(jù)方差的定義?[]X (37)22???var[][][]XXXmERm??由式（36）可知，上式等號右邊第一項(xiàng)可表示為 (38)1220?[]{[]][]}()NXnkERmnk????當(dāng) 是零均值正態(tài)序列時(shí)，它的各高階矩都可以用其一階和二階矩來表示。首先求 的均值，有?[]XRm 10?[]{[]]}NmXnEEXx???? 10?[]NmXnR??即自相關(guān)函數(shù)的估計(jì)值 的數(shù)學(xué)期望等于序列 的自相關(guān)函數(shù)的真實(shí)值。?[]Xm22，求 的離散傅里葉變換，得到 功率譜的估值：[] (),1MjmXXmGReN???????由于式(36)僅利用了 的 個(gè)有限值得到自相關(guān)函數(shù)的估計(jì) ，因此它[]nN?[]XRm與信號本身的自相關(guān)函數(shù) 會有一定程度的差別。 BT 法 BT 法BT 法首先計(jì)算隨機(jī)信號 的自相關(guān)函數(shù) ，若[]Xn[]())]XRmEnm??，可由 的傅里葉變換的得到它的功率譜，即：[]XmR??????Rm ()[]jmXXmGRe????????BT 法的估計(jì)步驟如下：基于 MATLAB 的譜估計(jì)實(shí)現(xiàn)171，先估計(jì)隨機(jī)信號 的自相關(guān)函數(shù)。在采用周期圖法進(jìn)行功率譜估計(jì)時(shí)，可以將平均周期圖法和加窗周期圖法這兩種方法結(jié)合在一起，效果會更好。但是在實(shí)際應(yīng)用中，幅度相差100倍甚至更高的正弦分量同時(shí)存在于一個(gè)信號中是可能的，這種情況下，由于矩形窗的旁瓣幅度較大，難以區(qū)分這樣的兩個(gè)正弦分量，故估計(jì)的效果較差，故還應(yīng)考慮旁瓣幅度較小的窗函數(shù)。對圖36的主瓣進(jìn)行放大后如圖37所示.基于 MATLAB 的譜估計(jì)實(shí)現(xiàn)16圖37從圖37可以看出，只有運(yùn)用矩形窗和凱撒窗才能分辨出這兩個(gè)頻率相差較小的信號。仿真圖如下：圖36 從圖36可以看出，采用blackman窗和hamming窗后，功率譜的