【正文】 傅里葉變換對(duì)，即 ()()jXXGRed????????以及 1()()2jXXRed????????? 這一關(guān)系就是著名的維納——辛欽定理。??無(wú)偏性的描述如下若估計(jì)量 的數(shù)學(xué)期望 存在，且對(duì)于任意 有12?(,)nX?…?[]E?????則稱 是 的無(wú)偏估計(jì)量。我們自然希望隨著樣n本容量的增大，一個(gè)估計(jì)量的值穩(wěn)定于待估參數(shù)的真值。??可以證明，樣本 階矩是總體 的 階矩的相合估計(jì)量?；?MATLAB 的譜估計(jì)實(shí)現(xiàn)6 周期圖 周期圖的定義對(duì)于一個(gè)離散隨機(jī)過(guò)程，由公式(21)，以及樣本的一階矩就是對(duì)總體期望的無(wú)偏估計(jì)可得： (31)10?[][]]NmXnRX????其中 是信號(hào)的采樣點(diǎn)數(shù)， 是對(duì)此離散時(shí)間隨機(jī)過(guò)程自相關(guān)函數(shù)的估計(jì)量。[]n []Y?()jNe??根據(jù)功率譜密度與自相關(guān)函數(shù)的關(guān)系，功率譜估計(jì)為自相關(guān)函數(shù)估計(jì)的傅里葉變換，即 211?()()()jj jXNNGeXe?????直接將 模的平方除以 求得的功率譜估計(jì)的方法稱為周期圖法，其結(jié)jNe?果用 表示，即()I基于 MATLAB 的譜估計(jì)實(shí)現(xiàn)7 (32)21?()()jXNNGIXe???周期圖法是利用數(shù)據(jù)的傅里葉變換直接求得的，而不再計(jì)算自相關(guān)函數(shù)，所以又稱直接法。當(dāng) 時(shí)N?? lim[]li1Nmb????????????故 li()NB??此時(shí) (34)li[]()XEIG???因此，周期圖作為功率譜估計(jì)，當(dāng) 時(shí)是無(wú)偏的，即漸近無(wú)偏。N??2()()NXEIG?為了求 ，先求 在兩個(gè)頻率 和 處的協(xié)方差，最后令2[]NI12。,npqkq??其 他將上式代入式(36)中，得 12 12411())()21200[()]NNjnkjnkXNnknkEI ee???????????? ??? ?? ?? 2211422()()sisinnXNN?? ?????? ????????? ??????? ?將 代入上式，可得12?? 224sin()[()]NXEI??????????????????? 24si()var[()]1nNXNI????????由上式可見(jiàn)，當(dāng) 時(shí)， 。平均周期圖法就是將長(zhǎng)度為 的觀測(cè)數(shù)據(jù)截?cái)喑?段，每一段NK長(zhǎng)度為 ,分別計(jì)算每一段的周期圖，最后進(jìn)行平均得到平均周期圖的估計(jì)，其數(shù)學(xué)L表達(dá)式如下： 1() 20?|()|LmjnpermnGXe??????基于 MATLAB 的譜估計(jì)實(shí)現(xiàn)101()0??()KmavperperG????由此可見(jiàn)，估計(jì)的方差很明顯多了一項(xiàng) ，這樣方差便隨著分段數(shù)目的增加而下降。在理想的情況下應(yīng)Sa該是原信號(hào)乘以幅值始終為1的信號(hào)，此信號(hào)在頻域表現(xiàn)為一沖激函數(shù)，主瓣無(wú)限窄。對(duì)其進(jìn)行改進(jìn)的辦法是將矩形窗換為其它的窗函數(shù)。[]n? 周期圖的仿真在 MATLAB 中用周期圖方法進(jìn)行功率譜估計(jì)的函數(shù)為 periodogram。nfft=1024。所以周期圖不是功率譜的一致估計(jì)。%產(chǎn)生含有噪聲的序列xn=cos(2*pi*40*n)+ 5*cos(2*pi*90*n)+*randn(size(n))。k=index*Fs/256。sx=0。 p=p/K。0 50 100 150 200 250 30030201001020300 50 100 150 200 250 3003020100102030圖 34 圖 35從圖 34 和圖 33 的對(duì)比中很明顯可以看出，圖 34 頻率為 40Hz 和 90Hz 處的兩個(gè)基于 MATLAB 的譜估計(jì)實(shí)現(xiàn)14主瓣的寬度較圖 33 而言明顯變窄，更接近于原信號(hào)。nfft=512。index=0:round(nfft/2 1)。從此圖與圖 34 的對(duì)比中可以看出，采用數(shù)據(jù)重疊方法后，主瓣寬度變得更窄，更接近于原信號(hào)。nfft=1024。window4=hamming(length(xn))。[Pxx1,f1]=periodogram(xn,window1,nfft,Fs)。[Pxx5,f5]=periodogram(xn,window5,nfft,Fs)。plot(f1,10*log10(Pxx1))。subplot(332)。)。blackman窗39。title(39。plot(f5,10*log10(Pxx5))。subplot(336)。 )。triang窗39。title(39。title(39。title(39。title(39。title(39。title(39。title(39。在上述程序中包含兩個(gè)頻率相差較小的正弦信號(hào)，但是運(yùn)用不同窗函數(shù)，估計(jì)的分辨率也不同。所以，在實(shí)際應(yīng)用中，具體選擇哪種窗函數(shù)，應(yīng)該從旁瓣幅度和分辨率兩個(gè)角度綜合考慮到底選用哪一個(gè)窗函數(shù)進(jìn)行功率譜估計(jì)的效果較好。自相關(guān)函數(shù)由下式估計(jì)：[]Xn (36)10?[]],1NmnRXN??????式中， 的長(zhǎng)度為 。因?[]XR[]n此 是自相關(guān)函數(shù) 的無(wú)偏估計(jì)。即Nm??li{var[]}0XN??故 滿足一致估計(jì)的條件。?[][]XXNmERR??這相當(dāng)于真值 用三角窗函數(shù)加權(quán)。 39。?[]XRm[] ?[]XRm所以一般用 作為自相關(guān)函數(shù)的估計(jì)，而較少使用 。%產(chǎn)生含有噪聲的序列xn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+randn(size(n))。%CXk估計(jì)出來(lái)以后是復(fù)數(shù)，此處是對(duì)它取模index=0:round(nfft/2 1)。plot(k,plot_Pxx)。解決的方法也可以在時(shí)域采用不同的窗函數(shù)，從而使得分辨率有所提高。[]n這個(gè)線性系統(tǒng)的傳遞函數(shù)為 0()qkkpbzBHzAa???當(dāng)輸入白噪聲的功率譜密度 時(shí)，輸出的功率譜密度為2()G?? 222()()jjXBeHeA???如果能確定 與 ， 的值，通過(guò)上式就可得到所需信號(hào)的功率譜密度。模型輸出功率譜密p度為 221())XjpjkkGAeae????????基于 MATLAB 的譜估計(jì)實(shí)現(xiàn)22(43)只要能求得 及所有 的值，就可以得到隨機(jī)信號(hào) 的功率譜。212{,}pa??…基于 MATLAB 的譜估計(jì)實(shí)現(xiàn)23 AR 模型的協(xié)方差方法與修正協(xié)方差方法介紹兩種方法之前首先介紹一下前、后向線性預(yù)測(cè)誤差濾波器的基本知識(shí)。,pk如果利用 數(shù)據(jù)預(yù)測(cè) ，則成為后向預(yù)測(cè)，其估計(jì)[1],[2],[]XnXnp??…[]Xn值用 表示?39。[]39。[]pbp kenXnpXaXn????????同樣，利用最小均方誤差準(zhǔn)則可求得 。而在協(xié)方差算法中運(yùn)用了前向預(yù)測(cè)算法求模型參數(shù)。pf?b由數(shù)字濾波器的相關(guān)理論可知，前向線性預(yù)測(cè)誤差濾波器的系數(shù)和后向線性預(yù)測(cè)誤差濾波器的系數(shù)互為共軛。設(shè)信號(hào) 觀測(cè)數(shù)據(jù)區(qū)間為： ，前向，后向預(yù)測(cè)誤差功率分別[]Xn01nN??用 和 表示，公式分別為：pf?b (49) 12[]Nfpf pneP???? (410)12[]bpbpn?? ()ppfb????我們已經(jīng)知道，前向，后向預(yù)測(cè)誤差公式分別為式(47)和式(48)。 AR 模型法的仿真 YuleWalker 方法YuleWalker 方法的仿真程序如下：fs=200。pyulear(xn,10,nfft,fs)如前所述，對(duì)于YuleWalker方法，當(dāng)數(shù)據(jù)很長(zhǎng)，估計(jì)較準(zhǔn)確，但是當(dāng)估計(jì)數(shù)據(jù)較短時(shí)，估計(jì)誤差較大，會(huì)出現(xiàn)譜峰頻率偏移與虛假譜線的問(wèn)題。仿真程序如下：fs=200。order=50。pmcov(xn,order,nfft,fs)%AR模型改進(jìn)協(xié)方差算法基于 MATLAB 的譜估計(jì)實(shí)現(xiàn)27仿真圖如下：0 10 20 30 40 50 60 70 80 90 10040302010010Frequency (Hz)Power/frequency (dB/Hz)Power Spectral Density Estimate via Covariance0 10 20 30 40 50 60 70 80 90 10040302010010Frequency (Hz)Power/frequency (dB/Hz)Power Spectral Density Estimate via Modified Covariance圖 43從上圖可以看出，協(xié)方差方法和修正協(xié)方差方法在分辨率方面的性能甚至比加窗周期圖的分辨率還要好。n=0:1/fs:1。order1=50。 AR 模型階的選擇在 AR 模型中，模型階次 的選擇是一個(gè)關(guān)鍵問(wèn)題。這一點(diǎn)可由下圖說(shuō)明。0 10 20 30 40 50 60 70 80 90 100302520151050510Frequency (Hz)Power/frequency (dB/Hz)Power Spectral Density Estimate via YuleWalker圖 45（階數(shù)為 50）基于 MATLAB 的譜估計(jì)實(shí)現(xiàn)290 10 20 30 40 50 60 70 80 90 10035302520151050510Frequency (Hz)Power/frequency (dB/Hz)Power Spectral Density Estimate via YuleWalker0 10 20 30 40 50 60 70 80 90 10035302520151050510Frequency (Hz)Power/frequency (dB/Hz)Power Spectral Density Estimate via YuleWalker 圖 46（階數(shù)為 30） 圖 47（階數(shù)為 80） 在 MATLAB 的 hosa 工具箱中給出了 AR 模型的定階函數(shù)：arorder?；?MATLAB 的譜估計(jì)實(shí)現(xiàn)305 結(jié)論本篇論文介紹了功率譜估計(jì)的幾種常用方法，包括周期圖及其修正，BT 法，AR 模型法，給出了詳細(xì)的估計(jì)步驟和估計(jì)性能，對(duì)每一種方法及相應(yīng)的改進(jìn)方法進(jìn)行了基于 MATLAB 的仿真，并結(jié)合仿真圖進(jìn)一步討論了每種方法的估計(jì)質(zhì)量。此外，為了提高譜分辨率，還進(jìn)行了加窗修正?；?MATLAB 的譜估計(jì)實(shí)現(xiàn)32參考文獻(xiàn)(1)張旭東等，離散隨機(jī)信號(hào)處理[M].北京:清華大學(xué)出版社，2022.(2)景占榮等，信號(hào)檢測(cè)與估計(jì)[M]。 作者簽名： 日期： 畢業(yè)論文（設(shè)計(jì)）授權(quán)使用說(shuō)明本論文（設(shè)計(jì)）作者完全了解**學(xué)院有關(guān)保留、使用畢業(yè)論文（設(shè)計(jì)）的規(guī)定，學(xué)校有權(quán)保留論文（設(shè)計(jì)）并向相關(guān)部門送交論文（設(shè)計(jì)）的電子版和紙質(zhì)版。 作者簽名： 指導(dǎo)教師簽名： 日期： 日期： 基于 MATLAB 的譜估計(jì)實(shí)現(xiàn)34注 意 事 項(xiàng)（論文）的內(nèi)容包括：1）封面（按教務(wù)處制定的標(biāo)準(zhǔn)封面格式制作）2）原創(chuàng)性聲明3）中文摘要（300 字左右） 、關(guān)鍵詞4）外文摘要、關(guān)鍵詞 5）目次頁(yè)（附件不統(tǒng)一編入）6）論文主體部分：引言（或緒論） 、正文、結(jié)論7）參考文獻(xiàn)8）致謝9）附錄（對(duì)論文支持必要時(shí)）：理工類設(shè)計(jì)（論文）正文字?jǐn)?shù)不少于 1 萬(wàn)字（不包括圖紙、程序清單等） ，文科類論文正文字?jǐn)?shù)不少于 萬(wàn)