【正文】 主瓣寬度比較窄，頻率泄露的情況相對不太明顯（也就是除主瓣以外的譜線的幅度較主瓣幅度相比較?。?。triang窗39。figure(8)plot(f7,10*log10(Pxx7))。kaiser窗39。figure(7)plot(f6,10*log10(Pxx6))。hann窗39。figure(6)plot(f5,10*log10(Pxx5))。hamming窗39。figure(5)plot(f4,10*log10(Pxx4))。blackman窗39。figure(4)plot(f3,10*log10(Pxx3))。bartlett窗39。figure(3)plot(f2,10*log10(Pxx2))。boxcar窗39。figure(2)plot(f1,10*log10(Pxx1))。)。title(39。subplot(337)。kaiser窗39。plot(f6,10*log10(Pxx6))。)。title(39。subplot(335)。hamming窗39。plot(f4,10*log10(Pxx4))。)。title(39?；?MATLAB 的譜估計實現(xiàn)15subplot(333)。bartlett窗39。plot(f2,10*log10(Pxx2))。)。title(39。figure(1)subplot(331)。[Pxx6,f6]=periodogram(xn,window6,nfft,Fs)。[Pxx4,f4]=periodogram(xn,window4,nfft,Fs)。%直接法[Pxx2,f2]=periodogram(xn,window2,nfft,Fs)。window7=triang(length(xn))。window5=hann(length(xn))。window3=blackman(length(xn))。window1=boxcar(length(xn))。%產(chǎn)生含有噪聲的序列xn=cos(2*pi*40*n)+cos(2*pi*41*n)+3*cos(2*pi*90*n)+*randn(size(n))。 加窗周期圖的仿真采用加窗周期圖方法進(jìn)行功率譜估計的仿真程序如下：Fs=600。figure(1)plot(k,plot_Pxx) 數(shù)據(jù)重疊方法的仿真圖如圖 35 所示。k=index*Fs/nfft。[Pxx,f]=pwelch(xn,window,noverlap,nfft)。window=boxcar(100)。xn=cos(2*pi*40*n)+5*cos(2*pi*90*n)+*randn(size(n))。程序如下：Fs=600。end。 sx=sx+p。for k=1:K [p,f]=periodogram(xn(n1:(n1+L1)))。f=0。 function [sx,f]=bart(xn,K)L=floor(length(xn)/K)。%k相當(dāng)于頻率,256是PERIODOGRAM（在bart函數(shù)中調(diào)用的）函數(shù)中默認(rèn)的nfft 的值plot_Pxx=10*log10(Pxx(index+1))。index=0:round(256/2 1)。K=10。 %采樣頻率n=0:1/Fs:1。如圖 41，真實的譜線應(yīng)該是在頻率為 40Hz 和 90Hz 處為兩個沖激，而周期圖法估計出的功率譜在這兩個基于 MATLAB 的譜估計實現(xiàn)13頻率出的主瓣非常寬，失真較明顯。0 50 100 15060504030201001020圖33由前述討論可知，在周期圖方法中，無論怎樣增大數(shù)據(jù)采樣點數(shù) ，功率譜估N計的方差一直不為零。[Pxx,f]=periodogram(xn,window,nfft,Fs)。 xn=cos(2*pi*40*n)+5*cos(2*pi*90*n)+*randn(size(n))。仿真程序如下：Fs=600。采用加窗周期圖方法的數(shù)學(xué)表達(dá)式為 10[]NjwnnXe?????? ??其中， 為窗函數(shù)。幾種常用窗函數(shù)的數(shù)學(xué)表達(dá)式如下：矩形窗： ()1,0nN????三角窗（bartlett）: 2/(1),0()/2nnN??????漢寧窗： 1()cos,21nN?????????????????哈明窗： 2(),0nN???????布萊克曼窗： 24().,011nnN?????????????????????基于 MATLAB 的譜估計實現(xiàn)11在MATLAB中還給出了幾種常用的窗函數(shù)，他們的時域和頻域特性分別如下：50 0 5001boxcar域域域域域域50 0 5001bartlett域域域域域域50 0 5001blackman域域域域域域50 0 5001hamming域域域域域域50 0 5001hann域域域域域域50 0 5001kaiser域域域域域域50 0 5001triang域域域域域域圖31500 0 5000204060 boxcar域域域域域域500 0 5000102030 bartlett域域域域域域500 0 5000102030blackman域域域域域域500 0 5000102030hamming域域域域域域500 0 5000102030 hann域域域域域域500 0 50002040 kaiser域域域域域域500 0 5000102030 triang域域域域域域圖32從這幾種窗函數(shù)的頻域特性可以看出，矩形窗的主瓣寬度最小，因此運(yùn)用矩形窗函數(shù)以后，譜分辨率是最高的，但旁瓣幅度較大，功率泄露情況也比較嚴(yán)重。這就是說，若信號中包含兩個以上正弦信號且頻率非常接近時，兩個分量的頻譜疊加的結(jié)果，兩個正弦信號不可能分辨，使信號模糊，降低了分辨率。而 函數(shù)的主瓣不是無限窄，主瓣寬度反比于數(shù)據(jù)記錄長度 ，而在實際中又不可Sa N能獲得很長的數(shù)據(jù)長度。時域上乘以一矩形窗相當(dāng)于頻域與一個 函數(shù)進(jìn)行卷積。為了在平均周期圖的基礎(chǔ)上更進(jìn)一步降低用周期圖法進(jìn)行功率譜估計的方差，可以采用的方法是使每一個分段的數(shù)據(jù)有一些重疊，使得觀測的數(shù)據(jù)段的數(shù)目增多，以此來達(dá)到進(jìn)一步減小方差的目的。 平均周期圖由于周期圖估計的方差總是不趨于零，所以需要對其進(jìn)行改進(jìn)，其中一種方法就是平均周期圖法。而且，無論 取何值，??4r[()]0XI???N總與 在同一數(shù)量級。,。1??根據(jù)式(32)可得 122122[()]()()jjNNNEIEXe?????????? 2[][]NNnkpqRnkpRq?基于 MATLAB 的譜估計實現(xiàn)9 (35) 12()()[][]jnkjpqEXnkpqe????其中 為矩形窗函數(shù)。按方差的定義應(yīng)有2X? 22var[()][()][()]NNNIEII????由于假設(shè) 零均值，方差為 的實高斯白噪聲序列，根據(jù)式(34)可知，當(dāng)[]nX?時， 。下面討論譜估計的方差。 211[()][()][()]j jjNNNEIXeEXe?????1100[()()]jnjknk????? [][[]]jnjknkEdXedXeN??????????? ()1[[][]]jknnk ??????令 代入上式得mkn??1[()][][]]jmNmnEIdmEXne???? ???? ???? ?? ?? （33）[]jNXmbRe??????式中 1[][]]Nnbdm???? ,10N?????其 他由于 為兩個矩形窗函數(shù)的卷積，因此它是一個三角窗函數(shù)，其傅里葉變[]Nbm換為基于 MATLAB 的譜估計實現(xiàn)8 22sin11()[]jnNn NBde??????????????????????式(33)為兩個函數(shù)乘積的傅里葉變換，等于他們各自傅里葉變換的卷積，即 [()]()NXEIBG???由上式可知，除非 是一個沖激函數(shù)，否則 將不等于 ，故B[()]NEI()XG?周期圖作為功率譜的估計是有偏差的。由于序列的傅里葉變化可利用 計算而提高運(yùn)算效率，這是周期圖FT法的主要優(yōu)點。當(dāng) 為實()jNXe?[]Xn()jN[]n序列時， 也為實序列，根據(jù)傅里葉變換的性質(zhì)， 的傅里葉變換為 。N?[]定義矩形窗函數(shù)如下 1,01[]nNd??????其 他則式(31)可寫成 1?[][][][]XnRmdXnmN??????令 ，上式可簡寫為[][]dnXY? 11?[][]XnRYY???????即 可看成是 與 的卷積。下面逐一闡述。(1)k?Xk相合性是對一個估計量的基本要求，若估計量不具有相合性，那么不論將樣本容量 取得多么大，都不能將 估計得足夠準(zhǔn)確，這樣的估計量是不可取的。12(,)n ?即，若對于任意 都滿足：對于任意 ，有???0???lim{}1nP??????則稱 是 的相合估計量。這樣，對估計量又有下述相合性的要求。?1??2 相合性無偏性和有效性都是在樣本容量 固定的前提下提出的。??可以證明，樣本均值（也就是總體的一階矩）是總體均值的無偏估計量，樣本方差是總體方差的無偏估計量。 是包含在總體 的分布中的待估參12,nX…X???X基于 MATLAB 的譜估計實現(xiàn)5數(shù)，這里 是 的取值范圍。它給出了平穩(wěn)隨機(jī)過程的時域特性與頻域特性之間的聯(lián)系，它是分析隨機(jī)信號的一個最重要，最基本的公式。?4，功率譜密度可積。 2，功率譜密度是 的實函數(shù)。)XnXREXxfndx????? 離散隨機(jī)過程的功率譜密度 離散隨機(jī)過程功率譜密度定義設(shè) 為寬平穩(wěn)離散時間隨機(jī)過程，或簡稱為寬平穩(wěn)隨機(jī)序列，具有零均值，[]Xn其自相關(guān)函數(shù)為基于 MATLAB 的譜估計實現(xiàn)4(){[]]}XRmEnTXm??或簡寫為 (21)()[]]X當(dāng) 滿足 的條件時，我們定義 的功率譜密度為 的()XRm()XR??????[]Xn()XRm離散傅里葉變換，并記為 ()G?()jmTXXmRe??????式中， 是隨機(jī)序列相鄰各值的時間間隔。[][]nmnmEXYE?nXmY離散時間隨機(jī)過程的均方值定義為 222[](。[][]nmnmXYY?（2） ，即 乘以一個常數(shù) 的均值等于 的均值乘以此常Ea?XanX基于 MATLAB 的譜估計實現(xiàn)3數(shù)。)nxnXmExfnd?????若 是單值函數(shù)，則 構(gòu)成一個新的離散時間隨機(jī)過程，其均值可定義()gA{}g為 []()。本篇論文主要涉及的信號類型屬于離散隨機(jī)信號，它一般是對一個連續(xù)信號進(jìn)行抽樣，對抽樣所得的序列進(jìn)行功率譜估計，故它屬于連續(xù)隨機(jī)序列。0,12,kX…， ks（4） 離散型隨機(jī)序列：相應(yīng)于時間和狀態(tài)都是離散的情況。（3） 連續(xù)隨機(jī)序列，隨機(jī)過程 在任一離散時刻的狀態(tài)是連續(xù)型隨機(jī)變()t量，也就是時間離散，狀態(tài)連續(xù)的情況?；?MATLAB 的譜估計實現(xiàn)22 隨機(jī)過程及估計理論基本知識 隨機(jī)過程的基本知識 隨機(jī)過程的分類按照時間和狀態(tài)（一般稱隨機(jī)過程 在 的可能取值為它的狀態(tài)）是連續(xù)()iXtit?的還是離散的來分類，可以分成以下四類：（1） 連續(xù)型隨機(jī)過程： 對于任意的 ， 都是連續(xù)型隨機(jī)變量，()ttT?()Xt也就是時間和狀態(tài)都是連續(xù)的情況。在現(xiàn)代譜估計的方法中，不再認(rèn)為在觀測到的 個數(shù)據(jù)以外的數(shù)據(jù)全為零，因此克服了經(jīng)典法的缺點，并且提高的譜估計的分辨率，能給出較好的估計。另一種是直接法，通過對觀測到的數(shù)據(jù)樣本進(jìn)行傅里葉變換，然后將所得到結(jié)果的幅值平方后得到功率譜估計，這種方法又稱為周期圖法。經(jīng)典譜估計法通常又分為兩種，一種是間接法，它是由布萊克曼（）和圖基（）提出，又稱為 BT 法。一類是非參數(shù)化方法，另一類是參數(shù)化方法。例如，根據(jù)信號、干擾與噪聲的功率譜，可以設(shè)計適當(dāng)?shù)臑V波器，以盡量不失真地重現(xiàn)信號，而最大幅度地抑制干擾與噪聲。在電子戰(zhàn)中，譜分析可用來對目標(biāo)進(jìn)行分類識別等。IIABSTRACTPower Spectral Estimation is a fundamental topic in the signal processing field. Actually