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

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

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

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

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點(diǎn)功率譜圖1周期圖法得出的功率譜估計(jì))Bd(譜率功隨著采樣點(diǎn)數(shù)增加，該估計(jì)是漸進(jìn)無(wú)偏的。從圖1中可以看出，采用周期圖法估計(jì)得到的功率譜很不平滑，也就是估計(jì)的協(xié)方差比較大，而且采用增加采樣點(diǎn)的方法也不能使周期圖變得更加平滑，這是周期圖法的缺點(diǎn)。周期圖法的優(yōu)點(diǎn)是能應(yīng)用離散傅里葉變換的快速算法來(lái)進(jìn)行估值。這種方法適用于長(zhǎng)信號(hào)序列的情況，在有足夠的序列長(zhǎng)度時(shí)，應(yīng)用改進(jìn)的周期圖法，可以得到較好的功率譜估值，因而應(yīng)用很廣。BT法BT法即Blackman-Turkey法，這種方法是Blackman與Turkey兩人于1958年提出的，故以他們的名字命名。這一節(jié)，我們研究BT法，并與周期圖法作比較。BT法譜估計(jì)正如我們已經(jīng)看到的，即使有很長(zhǎng)的樣本長(zhǎng)度，方差仍然很大時(shí)周期圖是這

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

由這樣兩個(gè)方面的原因同時(shí)引起的，即e6）中末端滯后?沁n）的血）估計(jì)精度

c低，以及大量（即便是少量）的協(xié)方差估計(jì)誤差在e6）中累積相加。ce6）e6）的定義公式ce6）=匸r（k）e-iOkck=-（N-1）（相關(guān)圖）是基于相關(guān)的psd定義（1）直接導(dǎo)出的相關(guān)圖譜估計(jì)算法。通過(guò)在e6）的定義式c?6)=￡w(k)r(k)e-iOkBT⑼>M),（8）式中截短求和的方法可以減少這兩方面的影響。這個(gè)思想就導(dǎo)出了BT估計(jì)子,其中M<N,t（?6)=￡w(k)r(k)e-iOkBT⑼>M),同時(shí)w（k）隨k平滑衰減到零。由于w（k）在式（9）中對(duì)樣本協(xié)方差序列的滯后進(jìn)行加權(quán)，故被稱為滯后窗。因?yàn)橛蛇@種方法求出的功率譜是通過(guò)自相關(guān)函數(shù)間接得到的，所以稱為間接法，又稱為自相關(guān)法或BT法。當(dāng)M較小時(shí)，式（9）的計(jì)算量不是很大，因此，該方法是在FFT問(wèn)世之前常用的譜估計(jì)方法。當(dāng)M=N時(shí)，BT法與周期圖法估計(jì)的功率譜是一樣的；當(dāng)M<N時(shí)，BT法的偏差大于周期圖法，在窗函數(shù)滿足一定條件下是漸進(jìn)無(wú)偏估計(jì)；方差小于周期圖的方差；分辨率比周期圖法低，與窗函數(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處都有一個(gè)相對(duì)較窄的主峰。由（12）式可得BT譜估計(jì)子（9）式就是一個(gè)局部加權(quán)的平均周期圖。BT譜估計(jì)的性質(zhì)由于elw]>0，自然也要求06）>0。選擇適當(dāng)?shù)拇昂瘮?shù)就能夠使估計(jì)普擁BT有這樣的性質(zhì)。

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

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

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

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

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

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

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