非參數(shù)譜估計算法與實現(xiàn)

非參數(shù)譜估計的算法與實現(xiàn)內(nèi)容摘要：功率譜估計是隨機信號分析中的一個重要內(nèi)容，主要研究信號在頻域中的各種特征，目的是根據(jù)有限觀測數(shù)據(jù)在頻域內(nèi)提取被噪聲污染的有用信號。本設(shè)計針對非參數(shù)譜估計的算法與實現(xiàn)展開研究。首先討論了一種最常見的譜估計法一一周期圖法。由于周期圖法是有偏估計，其方差過大，分辨率較低，因此，需要研究基于周期圖的各種改進方法。主要的改進途徑包括改善窗口形狀、對數(shù)據(jù)進行平均和平滑等。為此，討論了基于周期圖的其他幾種非參數(shù)譜估計算法（即BT法、平均周期圖法，包括Bartlett法和Welch法）及其通過MATLAB的仿真實現(xiàn)。此外，信號來波方向估計是功率譜估計算法在空間領(lǐng)域的推廣應(yīng)用，具體就是以多元天線陣結(jié)合現(xiàn)代數(shù)字信號處理為基礎(chǔ)的新型測向技術(shù)。為此，本設(shè)計還討論比較了采用Capon法、MUSIC算法實現(xiàn)波達方向估計的實現(xiàn)。關(guān)鍵詞：周期圖法平均周期圖法BT譜估計Capon法Music算法波達方向AlgorithmsandImplementationsofNon-ParameterSpectralEstimationsAbstract：Powerspectrumestimateisanimportantpartoftherandomsignalanalysis,whichprimarilystudiestherandomsignalcharacteristicsinthefrequencydomainandaimsatrecoveringtheusefulsignalinthefrequencydomaincontaminatedwithnoisebasedonthelimitedobserveddata.Thedesignfocusesonthealgorithmandimplementationofthenon-parameterspectralestimations.First,themostgeneralspectralestimation-periodogrammethodisdiscussed.Astheperiodogrammethodisabiasedestimatewiththelargevarianeeandlowresolution,theimprovedmethodsbasedontheperiodgrammethodneedtobestudied.Themainimprovementsincludechangingtheshapeofwindows,averagingandsmoothingthedata,etc.Asaresult,severalothernon-parameterspectralestimationalgorithmsbasedontheperiodgrammethod(suchastheBTmethod,theaveragingperiodgrammethod,includingtheBartlettmethodandWelchmethod)andtheirrealizationsbythesoftwareofMATLABaregiven.Moreover,thedirectionofarrival(DOA)estimationistoapplythepowerspectrumestimationalgorithmsintothespacedomain,i.e.,DOAisthenewdirectionfindingtechniquebasedonthemultiplesensorarraycombinedwiththemoderndigitalsignalprocessing.Accordingly,therealizationsofDOAestimationalgorithmsbasedonthemethodofCaponandMUSICarediscussedandcomparedinthisdesign.Keywords:PeriodgramAveragingPeriodgramBTSpectrumCaponmethodMusicAlgorithmDirectionofarrival(DOA)TOC\o"1-5"\h\z1緒論 11.1引言 1非參數(shù)譜估計 1空間譜估計 11.2本文的主要工作 2\o"CurrentDocument"2幾種常用的非參數(shù)譜估計算法與實現(xiàn) 22.1引言 22.2周期圖法 32.2.1周期圖譜估計 32.2.2周期圖法的性質(zhì) 32.2.3周期圖法在MATLAB中的仿真實現(xiàn) 42.3BT法 5BT法譜估計 5BT譜估計的性質(zhì) 6BT譜估計在MATLAB中的仿真實現(xiàn) 72.4平均周期圖法 7Bartlett方法及在MATLAB中的仿真實現(xiàn) 8Welch方法及在MATLAB中的仿真實現(xiàn) 102.5小結(jié) 13\o"CurrentDocument"3空間譜估計的算法與實現(xiàn) 143.1引言 143.2全向天線的波束下傾 143.2.1全向天線的波束下傾的基本原理 143.2.2全向天線的波束下傾的仿真實現(xiàn) 153.3天線的波達方向估計 153.3.1天線的波達方向估計的基本原理 16求波達方向的方法 173.4天線陣的波束形成 203.5小結(jié) 23\o"CurrentDocument"4結(jié)論 23致謝 錯誤！未定義書簽。參考文獻： 23附錄 24非參數(shù)譜估計的算法與實現(xiàn)1.1引言1.1.1非參數(shù)譜估計信號的頻譜分析是研究信號特性的重要手段之一。對于確定性信號，可以用傅里葉變換來考察其頻譜性質(zhì)。在通信系統(tǒng)中，最常見的信號往往不是確定信號，而是具有某種統(tǒng)計特性的隨機信號。對于隨機信號，由于它一般不是周期的，持續(xù)時間無限長，具有無限長能量的功率信號，不滿足傅里葉變換條件，隨機信號也不存在解析表達式，因此，對隨機信號，我們不能像確定信號那樣進行頻譜分析。然而，雖然隨機信號的頻譜不存在，但如果隨機信號是平穩(wěn)的，其相關(guān)函數(shù)就是確定的，那么相關(guān)函數(shù)的傅里葉變換就是它的功率譜密度函數(shù)，簡稱功率譜。通?？梢酝ㄟ^求隨機信號的功率譜來對這類信號進行頻譜分析。功率譜是頻率的函數(shù)，其物理單位是W/Hz,反映了隨機信號各頻率成份功率能量的分布情況，可以揭示信號中隱含的周期性及相鄰譜峰等有用信息。功率譜估計的應(yīng)用極其廣泛，例如，在語音信號識別、雷達雜波分析、地震勘探信號處理、水聲信號處理、系統(tǒng)辨識中非線性系統(tǒng)識別、物理光學(xué)中透鏡干涉、流體力學(xué)的內(nèi)波分析、太陽黑子活動周期研究等許多領(lǐng)域，發(fā)揮了重要作用。然而，實際應(yīng)用中的平穩(wěn)隨機信號通常是有限長的，只能根據(jù)有限長信號估計原信號的真實功率譜，這就是功率譜估計問題。功率譜估計分為非參數(shù)模型譜估計和參數(shù)模型譜估計，本文主要研究非參數(shù)模型譜估計.非參數(shù)模型譜估計又稱經(jīng)典譜估計，主要方法有：周期圖法，平均周期圖法，BT法及改進的周期圖法，如Bartiett法和Welch法。1.1.2空間譜估計電磁波到達方向(DOA)估計在無線電通信、雷達、聲納和導(dǎo)航等領(lǐng)域有著廣泛的應(yīng)用。由于空間信號存在多源可能，加上電磁波傳播的多徑效應(yīng)，將會給波束法測向帶來較大測量誤差，并受干擾信號的影響嚴重。利用陣列天線進行測向，其空間分辨率因陣列孔徑尺寸的選擇受經(jīng)典瑞利限的限制，在瑞利限以內(nèi)的空間目標則不能分辨。為了突破瑞利限約束，提高角度分辨力，新的測向技術(shù)不斷涌現(xiàn)。在波束形成技術(shù)、零點技術(shù)和時域譜估計技術(shù)的基礎(chǔ)上發(fā)展起來的空間譜估計技術(shù)是基于陣列信號處理技術(shù)的一種超分辨測向技術(shù)，用于提高在處理帶寬內(nèi)空問信號角度的估計精度和分辨力?？臻g譜分析采用了類似時域譜估計中的非線性處理，從而產(chǎn)生了一些特殊的算法，如多重信號分類法(MUSIC)、旋轉(zhuǎn)不變技術(shù)的參數(shù)估計法(ESPRIT)、最小內(nèi)積法JVW)、投影矩陣法和矩陣分解法等。這些方法都巧妙地利用接收信號的相關(guān)矩陣的特征結(jié)構(gòu)。其中，MUSIC算法將陣列輸出數(shù)據(jù)的自相關(guān)矩陣進行特征分解，然后利用這兩個子空間的正交性來估計信號的參數(shù)。這一算法的提出開創(chuàng)了空間譜估計算法研究的新時代，成為空間譜估計理論體系中的標志性算法。這種算法也是本文會研究到的?？臻g譜估計通過陣列信號處理，估計出陣列信號的空間譜，即信號源在空間的分布情況。采用該技術(shù)不僅可以提高測向精度和分辨力，而且加快了信號處理的運算速度。1.2本文的主要工作第一章，扼要介紹了譜估計，非參數(shù)譜估計以及空間譜估計等概念。第二章，首先研究周期圖法的算法與實現(xiàn)，用MATLAB對其仿真，分析仿真結(jié)果，總結(jié)周期圖法的優(yōu)缺點，再引出其他幾種對周期圖法的改進方法，如，BT法，平均周期圖法(Bartiett法以及Welch法)。從各種估計法的算法著手，結(jié)合MATLAB仿真，分析仿真結(jié)果，與周期圖法作比較，總結(jié)各種算法對周期圖法作出了哪些改進。第三章，首先引入什么是空間譜估計，簡介其主要應(yīng)用。采用Capon算法和MUSIC算法仿真實現(xiàn)了陣列側(cè)向以及波束形成。第四章，總結(jié)了本課題研究過程中的主要工作、理論結(jié)論。2幾種常用的非參數(shù)譜估計算法與實現(xiàn)2.1引言譜估計的非參數(shù)方法完全依賴于PSD的定義式來實現(xiàn)功率譜的估計，這些方法構(gòu)成了PSD估計的“經(jīng)典手段”本章主要討論幾種非參數(shù)譜估計方法、這些方法的性質(zhì)以及這些方法在MATLAB中的仿真。這里首先介紹一種最常用的譜估計算法，即由(2)式直接導(dǎo)出的周期圖法。對于足夠長的數(shù)據(jù)長度，周期圖法能夠獲得滿意的分辨率，但是不是可靠的譜估計算法，這是因為它的方差很大，而且并不隨著數(shù)據(jù)長度的增加而降低。由于周期圖法的估計方差很大，促使人們不斷地研究新的改進方法，試圖以降低分辨率為代價來獲得較小的估計方差。目前已經(jīng)提出了一些改進的方法，我們僅討論其中最常見的幾種?？梢宰C明，在大量數(shù)據(jù)樣本的條件下，所有這些方法在性質(zhì)和性能方面或多或少是等價的。

2.2周期圖法2.2.1周期圖譜估計一種信號功率譜密度估計方法。它的特點是：為得到功率譜估值，先取信號序列的離散傅里葉變換，然后取其幅頻特性的平方并除以序列長度 ，由于序列的離散傅里葉變換具有周期性，因而這種功率譜也具有周期性，常稱為周期圖。Schuster于1899年首先提出周期圖法，也稱直接法，取平穩(wěn)隨機信號Y(t)的有限個觀察值y(O), y(1)， ，y(T-1)對功率譜?6)進行估計，周期圖方法依賴于PSD的定義(2)式。略去只有信號樣本信息｛yOh的條件下不能i=1實現(xiàn)的期望和極限運算，我們得到?6)=—蘭y(t丄h(周期圖) ⑶pNi=1用來確定時間序列中可能存在“隱周期”，是周期圖譜估計子 (3)式最早的應(yīng)用之一，這也許可以認為這種方法的由來。2.2.2周期圖法的性質(zhì)分析?卩6)的統(tǒng)計特性的重要性在于，證明了周期圖作為PSD估計的不可靠性，此外:還提供了一些如何改進周期圖以便獲得更好的譜估計的內(nèi)部機理。我們用兩個部分來分析周期圖方法的性質(zhì)：偏差分析和方差分析。通常用偏差和方差這兩種基本度量來刻劃估計子的性能，其主要原因是一個估計的總平方誤差就是它的偏差的平方與方差之和。為了說明這個問題，設(shè)a為任意一個帶估計量，a是啊的一個估計。估計的均方差(MSE)為MSE=E,=EQa-MSE=E,=EQa-EQa\+EQa=EQa-EQaa+ReEla-EQa通過分別考慮MSE的偏差和方差分量，我們可以獲得誤差來源的深入了解，并得到減少誤差的方法。估計的均值：?6)=丄?6)*w6)其中WC)其中WC)是窗函數(shù)o(n)的傅里葉變換。當NTa時，E?6)-》?Co)是無偏估計，當N為有限值時，是有偏估計，偏差為：bias$6)=丄?(b^W6-九)d(x)—?6)

2兀bias估計的方差:var?Co)=―^F?(X)D6_X)D6+九加+E?6) ⑺_ 」2兀N_K 0 0 _其中，化(o)是矩形窗dQ(n)=1，|n<\N_1|的傅里葉變換，可見，$6)不是?(o)的一致估計。隨著N增大，譜估計的起伏增大，NTa時，i?(o)A?2(o)o估計的分辨率：除了通過偏差和方差來分析周期圖法的性質(zhì)，我們還可以通過分辨率來分析。數(shù)據(jù)窗為長度為N的矩形窗時，Res{?(o)}=0.89蘭，N增大時，分辨率N提高，但會使?6)的起伏加劇，可見方差與分辨率是一對矛盾。y周期圖法應(yīng)用比較廣泛，主要是由于它與序列的頻譜有對應(yīng)關(guān)系，可以采用FFT快速算法來計算。但是，這種方法需要對無限長的平穩(wěn)序列進行截斷，相當于對其加矩形窗，使之成為有限長數(shù)據(jù)。同時，這也意味著對自相關(guān)函數(shù)加三角窗，使功率譜與窗函數(shù)卷積，從而產(chǎn)生頻譜泄漏，容易使弱信號的主瓣被強信號的旁瓣所淹沒，造成頻譜的模糊和失真，使得譜分辨率較低。2.2.3周期圖法在MATLAB中的仿真實現(xiàn)仿真實例1用周期圖法對信號x(n)二、：20sin(2兀0.2n)+^2sin(2兀0.213n)+w(n)進行功率譜估計，其中w(n)為白噪聲。對信號x(n)的采樣率為f=1Hz，滿足取樣定理。程序代碼見附錄I。sn=(0,N-1)，當N=256和N=1024時，MATLAB仿真得到的功率譜估計圖分別為：

1ijAm'.A'.-A.1Jp-rV.LhJ'1,,'ijrT11II■■r■uVI1III111111I11200-20-400 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1頻率（Hz）1024頻率（Hz）1024點功率譜圖1周期圖法得出的功率譜估計)Bd(譜率功隨著采樣點數(shù)增加，該估計是漸進無偏的。從圖1中可以看出，采用周期圖法估計得到的功率譜很不平滑，也就是估計的協(xié)方差比較大，而且采用增加采樣點的方法也不能使周期圖變得更加平滑，這是周期圖法的缺點。周期圖法的優(yōu)點是能應(yīng)用離散傅里葉變換的快速算法來進行估值。這種方法適用于長信號序列的情況，在有足夠的序列長度時，應(yīng)用改進的周期圖法，可以得到較好的功率譜估值，因而應(yīng)用很廣。BT法BT法即Blackman-Turkey法，這種方法是Blackman與Turkey兩人于1958年提出的，故以他們的名字命名。這一節(jié)，我們研究BT法，并與周期圖法作比較。BT法譜估計正如我們已經(jīng)看到的，即使有很長的樣本長度，方差仍然很大時周期圖是這

種譜估計子存在的主要問題。周期圖譜估計較差的統(tǒng)計質(zhì)量可以直觀地解釋為是

由這樣兩個方面的原因同時引起的，即e6）中末端滯后?沁n）的血）估計精度

c低，以及大量（即便是少量）的協(xié)方差估計誤差在e6）中累積相加。ce6）e6）的定義公式ce6）=匸r（k）e-iOkck=-（N-1）（相關(guān)圖）是基于相關(guān)的psd定義（1）直接導(dǎo)出的相關(guān)圖譜估計算法。通過在e6）的定義式c?6)=￡w(k)r(k)e-iOkBT⑼>M),（8）式中截短求和的方法可以減少這兩方面的影響。這個思想就導(dǎo)出了BT估計子,其中M<N,t（?6)=￡w(k)r(k)e-iOkBT⑼>M),同時w（k）隨k平滑衰減到零。由于w（k）在式（9）中對樣本協(xié)方差序列的滯后進行加權(quán)，故被稱為滯后窗。因為由這種方法求出的功率譜是通過自相關(guān)函數(shù)間接得到的，所以稱為間接法，又稱為自相關(guān)法或BT法。當M較小時，式（9）的計算量不是很大，因此，該方法是在FFT問世之前常用的譜估計方法。當M=N時，BT法與周期圖法估計的功率譜是一樣的；當M<N時，BT法的偏差大于周期圖法，在窗函數(shù)滿足一定條件下是漸進無偏估計；方差小于周期圖的方差；分辨率比周期圖法低，與窗函數(shù)選擇有關(guān)。令W6）為w（k）的DTFT,(10)W6)=w(k)e-曲=氏(k)e一血k(10)K=-(M-1)利用DTFT的性質(zhì)，我們可以寫成BT6)BT6)=￡①(k)r(k)e={..,0,0,w(-(M-1》...,w(M-1),0,0,...} 和...,0,0,r(-(N-1))...,r(N-1)0,0..J的積的DTFTDTFT*{DTFT*{DTFT(w(k))}(11)記DTFT{..,0,0,r(-(N—1))…,r(N—1)0,0...}=化6)，可得(12)e6）=06）*w6）=—G）*w（p）d^(12)BT P 2兀斤P—兀由于極大多數(shù)常用的窗函數(shù)的WC）在w=0處都有一個相對較窄的主峰。由（12）式可得BT譜估計子（9）式就是一個局部加權(quán)的平均周期圖。BT譜估計的性質(zhì)由于elw]>0，自然也要求06）>0。選擇適當?shù)拇昂瘮?shù)就能夠使估計普擁BT有這樣的性質(zhì)。

如果滯后窗｛w（k）｝式正半定的（即W（to）＞0），則加窗的協(xié)方差序列也是正半定的，這意味著對所有的?,e6）＞oBT即BT譜估計具有非負性。證明 既然（w（k）｝對于k=0是實對稱的,那么其DTFTWC）是一個實的偶函數(shù)。進一步w（k）是一個正半定序列，則對于所有?，W6）＞0。由于定義e o，由（12）式可知，w6）＞o即意味著e 0。P BT應(yīng)該注意的是，有些滯后窗，如三角窗，并不滿足（13）式中的假設(shè)，因此,這些窗函數(shù)的使用可能導(dǎo)致估計譜出現(xiàn)負值。BT譜估計在MATLAB中的仿真實現(xiàn)仿真實例2用BT法對信號x（n）^.20sin（2兀0.2n）+、2sin（2兀0.213n）+w（n）進行功率譜估計，其中w（n）為白噪聲。對信號x（n）的采樣率為f=1Hz，滿足取樣s定理，n=（0,N-1）。程序代碼見附錄II。當N=512時，MATLAB仿真得到的功率譜估計圖為：35512占功率1譜估計±30i1[251*LB譜20\\\15?'屮功105[10"50 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1頻率（Hz）圖2BT法得出的功率譜估計圖通過實驗仿真可以直觀地看出功率譜估計中BT法和周期圖法所得到的結(jié)果是一致的，其特點是離散性大，曲線粗糙，方差較大，但是分辨率較高。2.4平均周期圖法周期圖是信號功率譜的一個有偏估值；而且，當信號序列的長度增大到無窮時，估值的方差不趨于零。因此，隨著所取的信號序列長度的不同,所得到的周期圖也不同，這種現(xiàn)象稱為隨機起伏。由于隨機起伏大，使用周期圖不能得到比較穩(wěn)定的估值。一些學(xué)者對此作了改進。

為了減小隨機起伏，M.S.Bartiett特提出平均周期圖法，即先把信號序列分為若干段，對每段分別計算其周期圖，然后取各個周期圖的平均作為功率譜的估值。平均周期圖可以減小隨機起伏，但是，如果信號序列不是足夠長，由于每段序列長度變短，功率譜估值對不同頻率成分的分辨能力也隨之下降。另一種改進方法是將周期圖與一個適當?shù)念l域窗函數(shù)相褶積，從而對周期圖產(chǎn)生平滑作用，以減小隨機起伏。加窗處理的結(jié)果雖然可以使隨機起伏減小，但也會使周期圖的分辨能力下降。P.O.Welch提出一種把加窗處理與平均處理結(jié)合起來的方法。先把分段的數(shù)據(jù)乘以窗函數(shù)（進行加窗處理），分別計算其周期圖，然后進行平均。Welch方法是較常用的一種計算方法。為了得到較好的功率譜估值，加窗和平均處理均應(yīng)兼顧減小隨機起伏和保證有足夠的譜分辨率兩個方面。2.4.1Bartiett方法及在MATLAB中的仿真實現(xiàn)（一）Bartlett方法Bartie方法的基本思想很簡單：為了減少周期圖的起伏，把長度為N的觀測樣本分割為長度為M的L二N/M個子樣本，然后對從這些子樣本得到的周期圖對于每一個o取平均。將Bartlett方法表示成如下的數(shù)學(xué)形式。令y(t)=y(t)=y((j-1)M+1)jt=1,…，M(14)表示第j個子樣本的觀測值，并令0(￡D0(￡D)=—蘭y(t)e-iOtji=1(15)表示對應(yīng)的周期圖。Bartlett譜估計由下式給出0Co）=丄藝06） （16）BLjj=1由于Bartlett方法是在長度為M的數(shù)據(jù)段上計算的，故其分辨率應(yīng)該約為1/M。因此，與與原來周期圖方法比較，Bartlett方法的分辨率降低了L倍。作為對分辨率降低的回報，可以預(yù)料，Bartlett方法的方差減小了。事實上，可以發(fā)現(xiàn)Bartlett方法把周期圖的方差減小了L倍。在選擇M（或L）時，分辨率和方差之間的折中是很明顯的。（二） Bartlett方法的性質(zhì)正如我們所知道的，（16）式中的06）可以寫成j06)=j06)=j藝#(k)e-iojk=-(m-1)(17)其中bjC)｝是第j個子樣本的協(xié)方差序列。把(17)式代入(16)式得我們發(fā)現(xiàn)?6)在形式上與使用矩形窗的BT估計很相似。子樣本協(xié)方差rj(k)在j上的平均是ACSrG)的一個估計。但是，在(18)式中ACS的估計并沒有有效地利用可以得到的數(shù)據(jù)滯后的乘積yQy丄-k)，特別是當|k接近M-1時。實際上，對于k二M-1，只有大約可以滯后乘積的1/M被用來形成(18)式中ACS的估計。可以預(yù)料，這些滯后的方差要大于BT估計中對應(yīng)的r(k)滯后，類似地，?6)方差要比? 6)的方差大。此外，由于Bartlett方法使用了一個固定的矩形B BT滯后窗，所以其分辨率-泄漏折中的靈活性也不如BT法。由于上述原因，我們總結(jié)為：在形式上，(14)—(16)式定義的Bartlett估計與使用長度為M的固定矩形滯后窗的BT估計類似，但其方差略大于后者。由此和以前討論的BT譜估計的性質(zhì)可知，與基本的周期圖方法比較，Bartlett估計降低了分辨率，減少了方差(其倍數(shù)均為L二N/M)。矩形窗的主瓣要比極大多數(shù)其他滯后窗的主瓣窄。所以，在BT類譜估計中，Bartlett估計可望得到最少的譜峰兼并(即有最大的分辨率)，但有最明顯的泄漏。(三)Bartlett估計的MATLAB仿真仿真實例3用Bartlett估計對信號xC)=sin2k60t+2sin2兀110t+n(t)進行

功率譜估計，其中n(t)為高斯白噪聲。對信號xC)的采樣頻率為f=500Hz，滿足抽樣定理。N=500，源程序代碼見附錄III。得到的功率譜估計圖如下：Bartlett估計對周期圖法德改進的思想是將信號分段進行估計，并將這些估計結(jié)果進行平均而減少估計的方差，使得估計功率譜圖變得平滑。如上例中將501點的信號分為三段，分別作周期圖法估計，然后加以平均。功率譜估計結(jié)果如圖3所示，將信號分段得到的估計值顯然平滑了許多。-20-40i1|L圖3Bartlett法頻率譜估計圖-20-40i1|L圖3Bartlett法頻率譜估計圖-60o50 100 150 200 250 300 350 400 450 500頻率(Hz))Bd(譜率功Welch方法及在MATLAB中的仿真實現(xiàn)（一）Welch譜估計Welch方法是Bartlett方法經(jīng)過兩個方面的改進得到的。首先，Welch方法中的數(shù)據(jù)段允許重疊。其次，在計算在周期圖前，對每一個數(shù)據(jù)段加窗。為了用數(shù)學(xué)方式描述Welch方法，令(19)y（t）=y（（j—1）K+1），t=人…,M(19)j j=1,…,S表示第j個數(shù)據(jù)段。在（19）式中，（j-1）k是第j個序列起始點。若K=M，則序列間沒有重疊（但是都是鄰接的），我們得到用于Bartlett方法中的樣本分割（得到S=L=N/M個數(shù)據(jù)子樣本）。但是，在Welch方法中K的推薦值為K=M/2，在這種情況下，得到S沁2M/N個數(shù)據(jù)段（前后數(shù)據(jù)段之間重疊50%）。計算與yC）對應(yīng)的加窗周期圖j(20)弋v(t)y(t)e一血(20)ji=1其中P表示時間窗｛C）｝的功率:i=1通過對(20)式中的加窗周期圖取平均得到PSD的Welch估計:要解釋上述對Bartle方法的改進，從而導(dǎo)出Welch方法的理由很簡單。通過允許數(shù)據(jù)段之間的重疊，就有更多的周期圖參與(22)式中的平均，故可以降低PSD估計的方差。在周期圖計算中引入窗函數(shù)，就能夠更有效地控制PSD估計中的偏差/分辨率性質(zhì)。除此之外，時間窗還可用來對位于沒個子樣本末端的數(shù)據(jù)給予較小的加權(quán)，這樣，盡管子樣本間是重疊的，但可以使相鄰的子樣本間的相關(guān)性減弱。通過(22)式中平均。這種去相關(guān)的重要作用能夠更有效地降低方差。(二)Welch譜估計的MATLAB仿真仿真實例4用Welch譜估計對信號xtsin260t2sin2110tnt進行功率譜估計，其中nt為高斯白噪聲。對信號xt的采樣頻率為f=500Hz，滿s足抽樣定理。N=500,源程序代碼見附錄IV。得到的功率譜估計圖如下：周期圖法得到的500點功率譜OO8-11?1周期圖法得到的500點功率譜OO8-11?11'tI.1 1C1f\r.；?-150 100 150 200 250 300 350 400 450 500頻率(Hz)將500點信號分五段且每段之間有重疊得到的功率譜50 100 150 200 250 300 350 400 450 500頻率(Hz)圖4采用樣本混疊增加分段數(shù)得到的功率譜O-40O-6*譜率功-40-60-8000-800/|i-40-60-8000-800/|ihIInI'li|嘰、?><50 100 150 200 250 300 350 400 450 500頻率(Hz)采用welch平均修正周期圖法得出的功率譜估計-20-40-6050 100 150 200 250 300 350 400 450 500頻率（Hz）圖5Welch平均修正周期圖法得到的功率譜估計圖增加分段數(shù)可以進一步降低估計的方差，然而若每段中的數(shù)據(jù)點數(shù)太少，就會使得估計的頻率分辨率下降很多。在樣本信號總點數(shù)一定的條件下，我們可以采用使分段相互重疊的方法來增加分段數(shù)，而保持每段中信號點數(shù)不變，這樣就在保證頻率分辨率的前提下進一步降低估計的協(xié)方差。女f仿真實例4所示，將xn分為5段，而且每段之間有重疊，對每段分別求FFT后求平均值，得出功率估計值。功率譜估計結(jié)果如圖4所示，與圖3相比，顯然估計值又平滑了許多。另外，采用加窗法也可以降低估計的協(xié)方差，這也是窗函數(shù)應(yīng)用的一個方面。即在計算周期圖法之前，對數(shù)據(jù)分段并加非矩形窗(如海明窗、漢寧窗或凱瑟窗等)，然后再采用分段長度一半的混疊率，就能夠大大降低估計方差。這種方法稱為Welch平均修正周期圖法，簡稱Welch法。仿真實例4是采用海明窗的Welch法。譜估計結(jié)果如圖5所示，與圖3，圖4相比，顯然Welch法得出的功率譜估計顯然平滑性較好。不同窗函數(shù)的Welch譜估計的MATLAB仿真仿真實例5用Welch方法對信號x(n)=f20sin(2兀0.2n)+^2sin(2兀0.213n)+w(n)進行功率譜估計，其中w(n)為白噪聲。對信號x(n)的采樣率為f=1Hz，滿足取樣定理，n=(0,N-1)。程序代碼見附s錄V。當N=1000時，窗函數(shù)采用矩形窗、Hanning窗hamming窗，MATLAB仿真得到的功率譜估計圖分別為：

node譜率功0*—■—-f―七卄卜嚴*—?一^-F0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1)node譜率功頻率（Hznode譜率功0*—■—-f―七卄卜嚴*—?一^-F0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1)node譜率功頻率（Hz）hanning窗400 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1頻率（Hz）)DBU(譜率功y.V11k"，--*J40hamming窗200 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9頻率（Hz）10圖6Welch法加不同的窗函數(shù)得出的功率譜估計圖仿真結(jié)果分析：由圖6可以看出，使用不同的窗函數(shù)譜估計的質(zhì)量是不一樣的，矩形窗的主瓣較窄，分辨率較好，但方差較大，噪聲水平較高；而Haning窗和Hamming窗的主瓣較寬，分辨率較低，但方差較小，噪聲水平較低。因此，在進行譜分析時選擇何種窗函數(shù)，要視具體情況而定。如果強調(diào)高分辨率，能精確讀出主瓣頻率，而不關(guān)心幅度的精度，例如測量震動物體的自震頻率，可以選用主瓣寬度比較窄的矩形窗；對受到強干擾的窄帶信號，若干擾靠近信號，則可選用旁瓣幅度較小的窗函數(shù)，若離開通帶較遠，則可選用漸近線衰減速度比較快的窗函數(shù)。總之，要針對不同的信號和不同的處理目的來選擇合適的窗函數(shù)，這樣才能得到良好的效果。2.5小結(jié)本章先研究周期圖譜估計，從均值，方差，分辨率討論周期圖法的性能。估計的均值得出當NTa時，是無偏估計，當N為有限值時，是有偏估計。估計的方差得出$6）不是?6）的一致估計。隨著N增大，譜估計的起伏增大。N增大時，分辨率提高，但會使?6）的起伏加劇，可見方差與分辨率是一y對矛盾。為此，提出了改進周期圖法的幾種途徑，即改善窗口形狀、平均和平滑。BT法中涉及到兩次加窗：一次是用于截取數(shù)據(jù)，叫數(shù)據(jù)窗；一次是用于截取相關(guān)函數(shù)，叫延遲窗，延遲窗常用三角窗。相關(guān)法的兩次加窗將造成泄漏效應(yīng)，影響到譜峰分辨率和方差。平均就是將截取的數(shù)據(jù)段再分成L個小段，分別計算功率譜后取功率譜的平均，這種方法使估計的方差減少，但偏差加大，分辨率下降。平滑是用一個適當?shù)拇昂瘮?shù)與算出的功率譜進行卷積，使譜線平滑。這種方法得出的譜估計是無偏的，方差也小，但分辨率下降?，F(xiàn)在比較常用的改進方法是Welch法，又叫加權(quán)交疊平均法，簡記為WOSA法，這種方法以加窗（加權(quán)）求取平滑，以分段重疊求得平均，因此集平均與平滑的優(yōu)點于一體，同時也不可避免帶有兩者的缺點，因此歸根到底是一種折中。Welch法得出的是無偏譜估計，當段數(shù)L增大時方差減小，趨于一致估計，分段平均減小了由數(shù)據(jù)樣本本身的隨機性帶來的方差，段數(shù)越多方差越小，但分辨率下降。非參數(shù)譜估計法的基礎(chǔ)是FFT，它們出現(xiàn)較早，運算量較小，物理概念明確，便于工程實現(xiàn)，仍是目前常用的譜估計方法，有的已被固化到FFT儀和譜分析儀中。但有一些難以克服的缺點:譜的分辨率較低，它正比于 （N是數(shù)據(jù)長度）。加窗的壞影響不可避免。較寬的主瓣降低分辨力，較大的旁瓣有可能掩蓋真實譜中較弱的成分，或是產(chǎn)生假的峰值。而沒有一個窗函數(shù)能使譜估計在方差、偏差和分辨率各方面同時改善的，使用窗函數(shù)只是一種折中的技巧，不是解決問題的根本辦法。不是真實譜的一致估計，且N增大時譜曲線起伏加劇。基于上述原因，非參數(shù)譜估計還有繼續(xù)研究的必要性。3空間譜估計的算法與實現(xiàn)3.1引言本章將討論使用一個m元無源陣元組成的陣列對n個輻射源進行定位的問題。這里，從輻射源來的能量可能是聲音、電磁波等等，接收傳感器將該能量轉(zhuǎn)換為電信號。接收陣元可能為電磁天線、水聲傳感器或者地震檢波器。這類問題在雷達或聲納系統(tǒng)、通信、天體物理學(xué)、生物醫(yī)學(xué)、地震水下監(jiān)視以及其他許多領(lǐng)域具有廣泛的應(yīng)用。這一問題的基本內(nèi)容包括確定能量在空間（空氣、水下或地下）的分布，這里能量高度集中的點就是輻射源。這類問題稱為空間譜估計。之所以稱為譜估計還因為輻射源定位問題和前面討論的譜估計有著非常緊密的聯(lián)系。事實上我們將看到，前面所述的方法幾乎都可以用作推導(dǎo)輻射源定位問題的解。3.2全向天線的波束下傾3.2.1全向天線的波束下傾的基本原理廣播用的發(fā)射機天線建在天線塔上，希望覆蓋廣大的地域。在覆蓋區(qū)的邊沿，由于距離遠，地球表面彎曲，電磁波場強衰減很快，因此為了有效地覆蓋既定的區(qū)域，保證區(qū)域內(nèi)場強不低于特定值，通常采用波束下傾的方法，將軸向排列的半波振子天線，通過調(diào)節(jié)天線的軸向間距、饋電的相位，使得軸向天線陣列的方向圖實現(xiàn)波束下傾。軸向排列的半振子天線結(jié)構(gòu)，陣元之間的距離是d，垂直軸線與電磁波輻射方向的夾角是9，相鄰陣元饋入信號的相位差是a，相鄰陣元發(fā)出電磁波到達

+ejd+ej2$+ejd+ej2$+ +ej6-山丿E=E(23)(23)?=Bdco0-aP=竺H九3.2.2全向天線的波束下傾的仿真實現(xiàn)根據(jù)以上討論，編寫出繪制軸向排列天線陣列的方向圖程序。源程序代碼見附錄切。仿真結(jié)果如下圖：圖7軸向陣列波束下傾仿真圖3.3天線的波達方向估計天線陣列信號處理作為信號處理的一個重要分支，在通信、雷達、地震勘測、射電天文等領(lǐng)域獲得了廣泛應(yīng)用和迅速發(fā)展。它包括：信源定位——確定陣列到信源的仰角和方位角，甚至距離（若信源位于近場）信源分離——確定各個信源發(fā)射的信號波形。各個信源從不同方向到達陣列，即使它們在時域和頻域上是疊加的。信道估計——確定信源與陣列之間的傳輸信道的參數(shù)（多徑參數(shù)）。陣列信號處理時改善蜂窩和個人通信服務(wù)系統(tǒng)質(zhì)量、覆蓋范圍和容量的一個強有力的工具。建立在波達方向估計、波束形成基礎(chǔ)上的智能天線的應(yīng)用，抑制了干擾信號，改善了欲接收信號的信噪比，降低了數(shù)字通信系統(tǒng)的誤碼率。將接收天線陣列用于公眾網(wǎng)的反向連接（客戶到基站）時，多個接收天線能夠收集更多的信號能量，若天線在空間足夠分離或極化各異，則多個天線能夠提供很好的分集接收，并抑制多徑傳輸引起的衰落。這些好處可以擴大基站的覆蓋范圍，改善通信質(zhì)量。3.3.1天線的波達方向估計的基本原理全向天線不僅利用率不高，而且對各種信號不加區(qū)別地接收，降低了通信質(zhì)量。定點無線通信采用定向天線，大幅度地改善了通信質(zhì)量。面對眾多移動用戶的公眾通信網(wǎng)基站和專用移動通信網(wǎng)，采用天線指向即波束可變的天線（智能天線），可以使移動通信的質(zhì)量得到很大的改善。為使天線的波束指向可控，甚至形狀可控，采用陣列天線是合適的。在距離通信源足夠遠的空間里，可以將到達的電磁波視為平面波。對于等距離直線陣天線，由于調(diào)制在載波上的基帶信號碼元寬度與波速的乘積遠大于天線陣的尺寸，因此多個天線陣元上的信號的幅度可視為不變，而它們的載波的相位差則取決于其相互位置、尺寸、波長和到達方向。無線接收的無線信號中有許多成分，其中我們關(guān)心的是S信號。天線陣列各個陣元接收的電磁波信號因為陣元排列位置的不同帶來相位差。經(jīng)過特定參數(shù)的加權(quán)控制器w處理后，進一步改變了各個陣元輸出信號的相位和幅度。處理的目標是使得陣元輸出信號的和天線輸出Y中的S成分具有最大輸出。用S信號作為基準信號，反饋控制單元的功能就是將輸出信號Y與基準信號S的差值（即誤差信號￡），作為調(diào)節(jié)控制加權(quán)控制器w參數(shù)的依據(jù)。反饋控制是使8減小，Y中的S成分加大，也就是說，天線陣列接收方向圖指向了S信號的方向。以下通過兩個簡單的例子介紹智能天線中波達方向估計和波束形式的原理。等距離直線陣智能天線，我們把天線陣元沿x軸排列，從0到M-1。若有一平面波以9角（入射線與z軸的夾角）和0角（入射線與x軸的夾角）入射到陣列上，第K號陣元上產(chǎn)生的信號為x，它與0號陣元的相位差為：K2兀A=Kdsi9co0 （24）K九式中，九和d分別是入射波的波長和陣元的間距，A亦稱陣因子。計入陣因子的影響，第K號陣元的輸出時Ax，即u。為了使天線陣的輸出滿足需要，在Kk k每個陣元上，用加權(quán)因子w進行控制。這樣第K號陣元上輸出的信號為KWRARxK，即wu。若到達天線陣的信號是N個，則天線陣的輸出是N個信號在M個車元上的輸出的疊加。將問題簡化為xy平面的二維問題（sin9二1），并用解析式表達如下：X（n）=lx（n）x（n）,x（n）］ （25）1 2 N

(d)p-j2兀-co(S1其中，0為第一個信號的入射角。1A=[a,A，…,A]TOC\o"1-5"\h\z1 2 NW=W,W，…,W]1 2 My(t)二WhAXh=WhU(26)(27)(28)(29)其中，U=lu,u,(26)(27)(28)(29)01 M=1求多個信號到達的方向（波達方向）的方法有多種。下面討論其中兩種方法（Capon法和MUSIC算法）及其仿真驗證結(jié)果。3.3.2求波達方向的方法Capon法Capon法亦稱最小方差無畸變響應(yīng)MVDR（MinimumVarianceDistortion-LessResponse）。天線陣列中的陣元數(shù)決定了陣列方向圖設(shè)計中的自由度數(shù)。Capon法將陣列中可控的自由度用來形成期望的波束形狀，達到對有用信號進行提升和對無用信號進行抑制的目的S并將其優(yōu)化問題表達為mineK（h）2匸minwhrw （30）UUw w其約束條件為WhA（）=1?可以證明上式的解為-A-A-1皿A-a代入式（29），可以得到相應(yīng)到的功率為Pcap1Pcap1AhG)R-1A(|))UU(31)MUSIC法MUSIC法亦稱多重信號分類（MultipleSignalClassification）。當入射信號數(shù)N小于陣元數(shù)M時，MxM的R矩陣有與進入陣列的信號數(shù)目相等的UU非零特征值及M-N個為零的特征值。V v ?… v]（V是R相應(yīng)的噪聲特征矢量）n N+1 N+2 M n UU

間譜為因為A與V的正交性，分母很小，峰值很大，這樣可以得出MUSIC法的空間譜為(32)_AhG)aG)icAh口iA?)(32)其中，口i_VVh稱為噪聲子空間的正交投影估計。nn仿真以下是對七單元線性天線陣在四信號輸入情況下了仿真仿真實例6設(shè)信號1從兀/4方向入射，信號2從兀/3方向入射，信號3從兀/6方向入射，信號4從3兀/4方向入射。下面是用兩種方法，分別用兩種坐標（直角坐標和極坐標）求波達方向估計。程序代碼見附錄⑷。仿真結(jié)果如下：90 1499999999999999900000000000000090 14999999999999999000000000000000270圖8Capon法作出的波達方向估計（極坐標）0 20 406080100 120 140 160 180圖9MUSIC法作出的波達方向估計（直角坐標）90 1500000270圖10 MUSIC法作出的波達方向估計（極坐標）10°11I111111/-I111u1v1■■■-_-0 20 40 60 80 100 120 140 160 18010-5圖11Capon法作出的波達方向估計(直角坐標)從對以上四圖的分析可以得到如下結(jié)論：兩種方法估計得都比較準確。在兀/6,兀/4,兀13,3兀/4處有尖銳的方向圖線。MUSIC法方向圖像幅度更大。3.4天線陣的波束形成1.原理我們以等距離圓陣為例來討論天線陣的波束形成。我們?nèi)A陣的圓心為O,在陣列圓上任取一點B，把天線陣元順序定為0到M-1。若有一平面波以0角入射到陣列上，第K號陣元上產(chǎn)生的信號為x，它K與到達陣元中心的波前的相位差是A@)=exp〔-j2兀?cosG+0)]；九與r分別k I九k丿是入射波的波長與陣列圓的半徑，A?亦稱陣因子。為了使天線陣的輸出滿足需要，在每個陣元上加上加權(quán)因子W控制。這樣第K號陣元上輸出的信號為K?Ax。若到達天線陣的信號是N個，天線陣輸出的是N個信號的M個陣元KKK上的輸出的疊加。用解析表達式如下：TOC\o"1-5"\h\z\o"CurrentDocument"A@)=exp(-j2兀；cos(?+0)) (33)K1 九kA=LA,A，…,A]hK K1K2 MX(n)=lx(n),x(n)…,x(n)]H (34)0 1 N-1式中，0是K陣元以O(shè)B為基準順時針畫出的角度。KTOC\o"1-5"\h\z\o"CurrentDocument"A=Ia,A，…,A] （35）1 2 N\o"CurrentDocument"W=lw,w，…,w1H （36）0 1 M-1y（t）=WhAX=WhU （37）為了求得多個信號到達的方向（波達方向），可以采用上述的Capon、MUSIC兩種方法）。波束形成可以采用下面的方法：當有多（N）個信號輸入時，其中1個信號時我們關(guān)心的，N-1個信號是需要抑制的。方程組（38）描述了上述需求的約束條件（四個信號輸入，第一個信號是我們關(guān)心的，其余信號是需要抑制的）。TOC\o"1-5"\h\zWhU=1,0,0,。｝ （38）根據(jù)信號波達方向的知識U（U=AX）及約束條件求解方程組（38），可以得到 W=\w,w，…,w] （39）01 M-1代入式（37）可得到陣列輸出的方向特性。仿真下面是四個輸入信號八單元圓陣列天線陣的波達方向估計和波束形成。仿真實例7設(shè)信號1從2兀/3方向入射，信號2從兀/3方向入射，信號3從圖12Capon算法波達方向估計圖13MUSIC算法波達方向估計90 1.5270圖14波束形成3.5小結(jié)本章在對非參數(shù)譜估計的理解基礎(chǔ)上，討論了其在空間領(lǐng)域進行推廣，具體就是了采用Capon算法和MUSIC算法求天線陣的波達方向，討論了天線陣的波束形成，及其MATLAB的仿真實現(xiàn)。4結(jié)論本設(shè)計的主要任務(wù)是研究非參數(shù)譜估計的算法與實現(xiàn)，另外對功率譜估計算法在空間領(lǐng)域的推廣應(yīng)用尤其是陣列側(cè)向問題進行了討論。在經(jīng)典譜估計中，無論是周期圖法還是其改進的方法，都存在著頻率分辨率低、方差性能不好的問題，原因是譜估計時需要對數(shù)據(jù)加窗截斷，用有限個數(shù)據(jù)或其自相關(guān)函數(shù)來估計無限個數(shù)據(jù)的功率譜，這其實是假定了窗以外的數(shù)據(jù)或自相關(guān)函數(shù)全為零，這種假定是不符合實際的，正是由于這些不符合實際的假設(shè)造成了經(jīng)典譜估計分辨率較差。另外，經(jīng)典譜估計的功率譜定義中既無求均值運算又無求極限運算，因而使得譜估計的方差性能較差，當數(shù)據(jù)很短時，這個問題更為突出。如何選取最佳窗函數(shù)、提高頻譜分辨率，如何在短數(shù)據(jù)情況下提高信號譜估計質(zhì)量，還需要進一步研究。Capon法和MUSIC算法可以用于空間譜估計測向，充分利用了天線各個陣元從空間電磁場接收到的全部信息，而傳統(tǒng)測向方法僅僅利用了其中一部分信息，因而具有多方面的優(yōu)越性：1）可以實現(xiàn)對同一一信道中存在的多個信號同時測向；2）超分辨率測向，即對處于天線陣固有波束寬度以內(nèi)的兩個以上方向的來波同時測向；3）陣元位置可以隨意放置，各個陣元的方向特性也可以是任意的；4）具有優(yōu)良的測向靈敏度和準確度；5）可以較好地克服測向模糊問題。參考文獻：PetreStoica,RandolphL.Moses.Introductiontospectralanalysis,PrenticeHall,Inc.l997魏平，唐斌，陸明泉，唐海燕，談覓舟，肖先賜.譜估計原理[J].電子科技大學(xué)研究生教材，1999胡廣書.數(shù)字信號處理一理論、算法與實現(xiàn)[M].北京：清華大學(xué)出版社,2003.王曉峰，王炳和.周期圖及其改進方法中譜分析率的Matlab分析J].武警工程學(xué)院學(xué)報，2003（6）：75—77.姚武川，姚天任.經(jīng)典譜估計方法的Matlab分析J].華中理工大學(xué)學(xué)報，2000,28(4)：45—47.劉剛，呂新華，攸陽.陣列信號處理中基于1USIC算法的空間譜估計[J].微計算機信息，2006，22(4-3):302?303.附錄附錄］%直接法clearall;clc;%N=256Fs=l;N=256;n=0:N-l;xn=sqrt(20)*sin(2*pi*0.2*n)+sqrt(2)*sin(2*pi*0.213*n)+randn(size(n);nfft二256;Pxx=10*log10(abs(fft(xn,nfft).“2)/(N));f=(0:length(Pxx)T)*Fs/length(Pxx);subplot(211);plot(f,Pxx)xlabel('頻率(Hz)');ylabel('功率譜(dB)');title('256點功率譜')grid%N=1024Fs=1;N=1024;n=0:N-1;xn=sqrt(20)*sin(2*pi*0.2*n)+sqrt(2)*sin(2*pi*0.213*n)+randn(size(n);nfft=1024;Pxx=10*log10(abs(fft(xn,nfft).“2)/N);f=(0:length(Pxx)T)*Fs/length(Pxx);subplot(212);plot(f,Pxx)xlabel('頻率(Hz)');ylabel('功率譜(dB)');title('1024點功率譜')grid附錄II%間接法，512點功率譜估計clearall;clc;Fs=l;N=512;n=O:N-l;xn=sqrt(20)*sin(2*pi*0.2*n)+sqrt(2)*sin(2*pi*O.213*n)+randn(size(n);nfft=512;Pxx二10*log10(fft(xcorr(xn,'unbiased'),nfft));f=(0:length(Pxx)T)*Fs/length(Pxx);plot(f,Pxx)xlabel('頻率(Hz)');ylabel('功率譜(dB)');title('512點功率譜估計')grid附錄皿%Bartlettclearall;clc;%N=500fs=500;t=0:1/fs:1;F=0:1:fs;xn=sin(2*pi*60*t)+2*sin(2*pi*110*t)+randn(1,length(t));Pxx1=abs(1/fs*fft(xn))「2;Pxx2=(abs(fft(xn(1:167)))."2+abs(fft(xn(168:334)))「2+abs(fft(xn(335:501)))?遼)/3;subplot(211);plot(F,10*log10(Pxx1));xlabel('頻率(Hz)');ylabel('功率譜(dB)');title('周期圖法得到的500點功率譜')；subplot(212);plot(0:3:fs,10*log10(Pxx2));xlabel('頻率(Hz)');ylabel('功率譜(dB)');title('將500點信號分三段得到的功率譜')；%Welchclearall;clc;%N=500fs=500;t=O:l/fs:l;F=O:l:fs;xn=sin(2*pi*60*t)+2*sin(2*pi*110*t)+randn(l,length(t));Pxxl二abs(l/fs*fft(xn))「2;Pxx3=(abs(3/fs*fft(xn(1:167)))「2+abs(3/fs*fft(xn(83:249)))."2+abs(3/fs*fft(xn(168:334)))「2+abs(3/fs*fft(xn(250:416)))「2+abs(3/fs*fft(xn(335:501))).V)/5;subplot(211);plot(F,10*log10(Pxx1));xlabel('頻率(Hz)');ylabel('功率譜(dB)');title('周期圖法得到的500點功率譜')；subplot(212);plot(0:3:fs,10*log(Pxx3))xlabel('頻率(Hz)');ylabel('功率譜(dB)');title('將500點信號分五段且每段之間有重疊得到的功率譜')；grid%Welchclearall;clc;%N=500fs=500;t=0:1/fs:1;F=0:1:fs;xn=sin(2*pi*60*t)+2*sin(2*pi*110*t)+randn(1,length(t));%截取時間段上的離散信號樣點序列Pxx1=abs(1/fs*fft(xn))「2;w=hamming(167)';w=w*sqrt(167/sum(w.*w));Pxx2=(abs(3/fs*fft(w.*xn(1:167)))."2+abs(3/fs*fft(w.*xn(83:249)))."2+abs(3/fs*fft(w.*xn(168:334)))."2+abs(3/fs*fft(w.*xn(250:416)))「2+abs(3/fs*fft(w.*xn(335:501)))「2)/5;subplot(211);plot(F,10*log10(Pxx1));xlabel('頻率(Hz)');ylabel('功率譜(dB)');title('周期圖法得到的500點功率譜')；subplot(212);plot(0:3:fs,10*log(Pxx2))xlabel('頻率(Hz)');ylabel('功率譜(dB)');title('采用welch平均修正周期圖法得出的功率譜估計')；grid附錄V%平滑周期圖法，即在自相關(guān)函數(shù)上加窗的方法。clearall;clc;Fs=1;N=1000;n=0:N-l;xn=sqrt(20)*sin(2*pi*0.2*n)+sqrt(2)*sin(2*pi*0.213*n)+randn(size(n);r二xcorr(xn,'unbiased');nfft=1024;window1=boxcar(length(r));window2=hanning(length(r));window3=hamming(length(r));a1二r.*window1';a2=r.*window2';a3=r.*window3';%求功率譜pxx1=10*log10(fft(a1,nfft));pxx2=10*log10(fft(a2,nfft));pxx3=10*log10(fft(a3,nfft));f1=(0:length(pxx1)-1)*Fs/length(pxx1);f2=(0:length(pxx2)-1)*Fs/length(pxx2);f3=(0:length(pxx3)T)*Fs/length(pxx3);subplot(311)plotf1,pxx1)xlabel('頻率(Hz)');ylabel('功率譜(dB)');title('矩形窗’)subplot(312)plot(f2,pxx2)xlabel('頻率(Hz)');ylabel('功率譜(dB)');title('hanning窗')subplot(313)plot(f3,pxx3)xlabel('頻率(Hz)');ylabel('功率譜(dB)');title('hamming窗')grid%clearall;clc;lam=1;t= [0 : 0.01 : 2*pi];d=0