非參數(shù)譜估計的算法與實現(xiàn)內(nèi)容摘要功率譜估計是隨機(jī)信號分析中的一個重要內(nèi)容，主要研究信號在頻域中的各種特征，目的是根據(jù)有限觀測數(shù)據(jù)在頻域內(nèi)提取被噪聲污染的有用信號。本設(shè)計針對非參數(shù)譜估計的算法與實現(xiàn)展開研究。首先討論了一種最常見的譜估計法周期圖法。由于周期圖法是有偏估計，其方差過大，分辨率較低，因此，需要研究基于周期圖的各種改進(jìn)方法。主要的改進(jìn)途徑包括改善窗口形狀、對數(shù)據(jù)進(jìn)行平均和平滑等。為此，討論了基于周期圖的其他幾種非參數(shù)譜估計算法（即BT法、平均周期圖法，包括BARTLETT法和WELCH法）及其通過MATLAB的仿真實現(xiàn)。此外，信號來波方向估計是功率譜估計算法在空間領(lǐng)域的推廣應(yīng)用，具體就是以多元天線陣結(jié)合現(xiàn)代數(shù)字信號處理為基礎(chǔ)的新型測向技術(shù)。為此，本設(shè)計還討論比較了采用CAPON法、MUSIC算法實現(xiàn)波達(dá)方向估計的實現(xiàn)。關(guān)鍵詞周期圖法平均周期圖法BT譜估計CAPON法MUSIC算法波達(dá)方向ALGORITHMSANDIMPLEMENTATIONSOFNONPARAMETERSPECTRALESTIMATIONSABSTRACTPOWERSPECTRUMESTIMATEISANIMPORTANTPARTOFTHERANDOMSIGNALANALYSIS,WHICHPRIMARILYSTUDIESTHERANDOMSIGNALCHARACTERISTICSINTHEFREQUENCYDOMAINANDAIMSATRECOVERINGTHEUSEFULSIGNALINTHEFREQUENCYDOMAINCONTAMINATEDWITHNOISEBASEDONTHELIMITEDOBSERVEDDATATHEDESIGNFOCUSESONTHEALGORITHMANDIMPLEMENTATIONOFTHENONPARAMETERSPECTRALESTIMATIONSFIRST,THEMOSTGENERALSPECTRALESTIMATIONPERIODOGRAMMETHODISDISCUSSEDASTHEPERIODOGRAMMETHODISABIASEDESTIMATEWITHTHELARGEVARIANCEANDLOWRESOLUTION,THEIMPROVEDMETHODSBASEDONTHEPERIODGRAMMETHODNEEDTOBESTUDIEDTHEMAINIMPROVEMENTSINCLUDECHANGINGTHESHAPEOFWINDOWS,AVERAGINGANDSMOOTHINGTHEDATA,ETCASARESULT,SEVERALOTHERNONPARAMETERSPECTRALESTIMATIONALGORITHMSBASEDONTHEPERIODGRAMMETHODSUCHASTHEBTMETHOD,THEAVERAGINGPERIODGRAMMETHOD,INCLUDINGTHEBARTLETTMETHODANDWELCHMETHODANDTHEIRREALIZATIONSBYTHESOFTWAREOFMATLABAREGIVENMOREOVER,THEDIRECTIONOFARRIVALDOAESTIMATIONISTOAPPLYTHEPOWERSPECTRUMESTIMATIONALGORITHMSINTOTHESPACEDOMAIN,IE,DOAISTHENEWDIRECTIONFINDINGTECHNIQUEBASEDONTHEMULTIPLESENSORARRAYCOMBINEDWITHTHEMODERNDIGITALSIGNALPROCESSINGACCORDINGLY,THEREALIZATIONSOFDOAESTIMATIONALGORITHMSBASEDONTHEMETHODOFCAPONANDMUSICAREDISCUSSEDANDCOMPAREDINTHISDESIGNKEYWORDSPERIODGRAMAVERAGINGPERIODGRAMBTSPECTRUMCAPONMETHODMUSICALGORITHMDIRECTIONOFARRIVALDOA目錄1緒論111引言1111非參數(shù)譜估計1112空間譜估計112本文的主要工作22幾種常用的非參數(shù)譜估計算法與實現(xiàn)221引言222周期圖法3221周期圖譜估計3222周期圖法的性質(zhì)3223周期圖法在MATLAB中的仿真實現(xiàn)423BT法5231BT法譜估計5232BT譜估計的性質(zhì)6233BT譜估計在MATLAB中的仿真實現(xiàn)724平均周期圖法7241BARTLETT方法及在MATLAB中的仿真實現(xiàn)8242WELCH方法及在MATLAB中的仿真實現(xiàn)1025小結(jié)133空間譜估計的算法與實現(xiàn)1431引言1432全向天線的波束下傾14321全向天線的波束下傾的基本原理14322全向天線的波束下傾的仿真實現(xiàn)1533天線的波達(dá)方向估計15331天線的波達(dá)方向估計的基本原理16332求波達(dá)方向的方法1734天線陣的波束形成2035小結(jié)234結(jié)論23致謝23參考文獻(xiàn)24附錄24非參數(shù)譜估計的算法與實現(xiàn)1緒論11引言111非參數(shù)譜估計信號的頻譜分析是研究信號特性的重要手段之一。對于確定性信號，可以用傅里葉變換來考察其頻譜性質(zhì)。在通信系統(tǒng)中，最常見的信號往往不是確定信號，而是具有某種統(tǒng)計特性的隨機(jī)信號。對于隨機(jī)信號，由于它一般不是周期的，持續(xù)時間無限長，具有無限長能量的功率信號，不滿足傅里葉變換條件，隨機(jī)信號也不存在解析表達(dá)式，因此，對隨機(jī)信號，我們不能像確定信號那樣進(jìn)行頻譜分析。然而，雖然隨機(jī)信號的頻譜不存在，但如果隨機(jī)信號是平穩(wěn)的，其相關(guān)函數(shù)就是確定的，那么相關(guān)函數(shù)的傅里葉變換就是它的功率譜密度函數(shù)，簡稱功率譜。通常可以通過求隨機(jī)信號的功率譜來對這類信號進(jìn)行頻譜分析。功率譜是頻率的函數(shù)，其物理單位是W/HZ,反映了隨機(jī)信號各頻率成份功率能量的分布情況，可以揭示信號中隱含的周期性及相鄰譜峰等有用信息。功率譜估計的應(yīng)用極其廣泛，例如，在語音信號識別、雷達(dá)雜波分析、地震勘探信號處理、水聲信號處理、系統(tǒng)辨識中非線性系統(tǒng)識別、物理光學(xué)中透鏡干涉、流體力學(xué)的內(nèi)波分析、太陽黑子活動周期研究等許多領(lǐng)域，發(fā)揮了重要作用。然而，實際應(yīng)用中的平穩(wěn)隨機(jī)信號通常是有限長的，只能根據(jù)有限長信號估計原信號的真實功率譜，這就是功率譜估計問題。功率譜估計分為非參數(shù)模型譜估計和參數(shù)模型譜估計，本文主要研究非參數(shù)模型譜估計非參數(shù)模型譜估計又稱經(jīng)典譜估計，主要方法有周期圖法，平均周期圖法，BT法及改進(jìn)的周期圖法，如BARTLETT法和WELCH法。112空間譜估計電磁波到達(dá)方向DOA估計在無線電通信、雷達(dá)、聲納和導(dǎo)航等領(lǐng)域有著廣泛的應(yīng)用。由于空間信號存在多源可能，加上電磁波傳播的多徑效應(yīng)，將會給波束法測向帶來較大測量誤差，并受干擾信號的影響嚴(yán)重。利用陣列天線進(jìn)行測向，其空間分辨率因陣列孔徑尺寸的選擇受經(jīng)典瑞利限的限制，在瑞利限以內(nèi)的空間目標(biāo)則不能分辨。為了突破瑞利限約束，提高角度分辨力，新的測向技術(shù)不斷涌現(xiàn)。在波束形成技術(shù)、零點技術(shù)和時域譜估計技術(shù)的基礎(chǔ)上發(fā)展起來的空間譜估計技術(shù)是基于陣列信號處理技術(shù)的一種超分辨測向技術(shù)，用于提高在處理帶寬內(nèi)空問信號角度的估計精度和分辨力。空間譜分析采用了類似時域譜估計中的非線性處理，從而產(chǎn)生了一些特殊的算法，如多重信號分類法MUSIC、旋轉(zhuǎn)不變技術(shù)的參數(shù)估計法ESPRIT、最小內(nèi)積法JVW、投影矩陣法和矩陣分解法等。這些方法都巧妙地利用接收信號的相關(guān)矩陣的特征結(jié)構(gòu)。其中，MUSIC算法將陣列輸出數(shù)據(jù)的自相關(guān)矩陣進(jìn)行特征分解，然后利用這兩個子空間的正交性來估計信號的參數(shù)。這一算法的提出開創(chuàng)了空間譜估計算法研究的新時代，成為空間譜估計理論體系中的標(biāo)志性算法。這種算法也是本文會研究到的?？臻g譜估計通過陣列信號處理，估計出陣列信號的空間譜，即信號源在空間的分布情況。采用該技術(shù)不僅可以提高測向精度和分辨力，而且加快了信號處理的運算速度。12本文的主要工作第一章，扼要介紹了譜估計，非參數(shù)譜估計以及空間譜估計等概念。第二章，首先研究周期圖法的算法與實現(xiàn)，用MATLAB對其仿真，分析仿真結(jié)果，總結(jié)周期圖法的優(yōu)缺點，再引出其他幾種對周期圖法的改進(jìn)方法，如，BT法，平均周期圖法BARTLETT法以及WELCH法。從各種估計法的算法著手，結(jié)合MATLAB仿真，分析仿真結(jié)果，與周期圖法作比較，總結(jié)各種算法對周期圖法作出了哪些改進(jìn)。第三章，首先引入什么是空間譜估計，簡介其主要應(yīng)用。采用CAPON算法和MUSIC算法仿真實現(xiàn)了陣列側(cè)向以及波束形成。第四章，總結(jié)了本課題研究過程中的主要工作、理論結(jié)論。2幾種常用的非參數(shù)譜估計算法與實現(xiàn)21引言譜估計的非參數(shù)方法完全依賴于PSD的定義式KTIER121LIMNITINYE2來實現(xiàn)功率譜的估計，這些方法構(gòu)成了PSD估計的“經(jīng)典手段”。本章主要討論幾種非參數(shù)譜估計方法、這些方法的性質(zhì)以及這些方法在MATLAB中的仿真。這里首先介紹一種最常用的譜估計算法，即由2式直接導(dǎo)出的周期圖法。對于足夠長的數(shù)據(jù)長度，周期圖法能夠獲得滿意的分辨率，但是不是可靠的譜估計算法，這是因為它的方差很大，而且并不隨著數(shù)據(jù)長度的增加而降低。由于周期圖法的估計方差很大，促使人們不斷地研究新的改進(jìn)方法，試圖以降低分辨率為代價來獲得較小的估計方差。目前已經(jīng)提出了一些改進(jìn)的方法，我們僅討論其中最常見的幾種?？梢宰C明，在大量數(shù)據(jù)樣本的條件下，所有這些方法在性質(zhì)和性能方面或多或少是等價的。22周期圖法221周期圖譜估計一種信號功率譜密度估計方法。它的特點是為得到功率譜估值，先取信號序列的離散傅里葉變換，然后取其幅頻特性的平方并除以序列長度，由于序列的離散傅里葉變換具有周期性，因而這種功率譜也具有周期性，常稱為周期圖。SCHUSTER于1899年首先提出周期圖法，也稱直接法，取平穩(wěn)隨機(jī)信號YT的有限個觀察值YO，Y1，YT1對功率譜進(jìn)行估計，周期圖方法依賴于PSD的定義2式。略去只有信號樣本信息NITY1的條件下不能實現(xiàn)的期望和極限運算，我們得到21NITIPEY（周期圖）3用來確定時間序列中可能存在“隱周期”，是周期圖譜估計子3式最早的應(yīng)用之一，這也許可以認(rèn)為這種方法的由來。222周期圖法的性質(zhì)分析P的統(tǒng)計特性的重要性在于，證明了周期圖作為PSD估計的不可靠性，此外，還提供了一些如何改進(jìn)周期圖以便獲得更好的譜估計的內(nèi)部機(jī)理。我們用兩個部分來分析周期圖方法的性質(zhì)偏差分析和方差分析。通常用偏差和方差這兩種基本度量來刻劃估計子的性能，其主要原因是一個估計的總平方誤差就是它的偏差的平方與方差之和。為了說明這個問題，設(shè)A為任意一個帶估計量，A是啊的一個估計。估計的均方差（MSE）為MSE22AEEAEAARE222VRBIS4通過分別考慮MSE的偏差和方差分量，我們可以獲得誤差來源的深入了解，并得到減少誤差的方法。（一）估計的均值WE215其中W是窗函數(shù)N的傅里葉變換。當(dāng)N時，E,是無偏估計，當(dāng)N為有限值時，是有偏估計，偏差為DWBIAS216（二）估計的方差EDDN20021VAR7其中，0D是矩形窗0ND，1的傅里葉變換，可見，不是的一致估計。隨著N增大，譜估計的起伏增大，N時，2VAR。（三）估計的分辨率除了通過偏差和方差來分析周期圖法的性質(zhì)，我們還可以通過分辨率來分析。數(shù)據(jù)窗為長度為N的矩形窗時，NS2890RE，增大時，分辨率提高，但會使Y的起伏加劇，可見方差與分辨率是一對矛盾。周期圖法應(yīng)用比較廣泛，主要是由于它與序列的頻譜有對應(yīng)關(guān)系，可以采用FFT快速算法來計算。但是，這種方法需要對無限長的平穩(wěn)序列進(jìn)行截斷，相當(dāng)于對其加矩形窗，使之成為有限長數(shù)據(jù)。同時，這也意味著對自相關(guān)函數(shù)加三角窗，使功率譜與窗函數(shù)卷積，從而產(chǎn)生頻譜泄漏，容易使弱信號的主瓣被強(qiáng)信號的旁瓣所淹沒，造成頻譜的模糊和失真，使得譜分辨率較低。223周期圖法在MATLAB中的仿真實現(xiàn)仿真實例1用周期圖法對信號2130SIN20SINNWX進(jìn)行功率譜估計，其中NW為白噪聲。對信號X的采樣率為F1HZ，滿足取樣定理。程序代碼見附錄。N0,N1，當(dāng)N256和N1024時，MATLAB仿真得到的功率譜估計圖分別為00102030405060708091402002040三三三HZ三三三三三DB三256三三三三00102030405060708091402002040三三三HZ三三三三三DB三1024三三三三圖1周期圖法得出的功率譜估計隨著采樣點數(shù)增加，該估計是漸進(jìn)無偏的。從圖1中可以看出，采用周期圖法估計得到的功率譜很不平滑，也就是估計的協(xié)方差比較大，而且采用增加采樣點的方法也不能使周期圖變得更加平滑，這是周期圖法的缺點。周期圖法的優(yōu)點是能應(yīng)用離散傅里葉變換的快速算法來進(jìn)行估值。這種方法適用于長信號序列的情況，在有足夠的序列長度時，應(yīng)用改進(jìn)的周期圖法，可以得到較好的功率譜估值，因而應(yīng)用很廣。23BT法BT法即BLACKMANTURKEY法，這種方法是BLACKMAN與TURKEY兩人于1958年提出的，故以他們的名字命名。這一節(jié)，我們研究BT法，并與周期圖法作比較。231BT法譜估計正如我們已經(jīng)看到的，即使有很長的樣本長度，方差仍然很大時周期圖是這種譜估計子存在的主要問題。周期圖譜估計較差的統(tǒng)計質(zhì)量可以直觀地解釋為是由這樣兩個方面的原因同時引起的，即C中末端滯后NK的KR估計精度低，以及大量（即便是少量）的協(xié)方差估計誤差在C中累積相加。C的定義公式C1NKKIER（相關(guān)圖）8是基于相關(guān)的PSD定義1直接導(dǎo)出的相關(guān)圖譜估計算法。通過在C的定義式8式中截短求和的方法可以減少這兩方面的影響。這個思想就導(dǎo)出了BT估計子，1MKKIBTERW9其中N,K是一個偶函數(shù)（KW），且10，MKW0，同時W隨K平滑衰減到零。由于在式9中對樣本協(xié)方差序列的滯后進(jìn)行加權(quán)，故被稱為滯后窗。因為由這種方法求出的功率譜是通過自相關(guān)函數(shù)間接得到的，所以稱為間接法，又稱為自相關(guān)法或BT法。當(dāng)M較小時，式9的計算量不是很大，因此，該方法是在FFT問世之前常用的譜估計方法。當(dāng)NM時，BT法與周期圖法估計的功率譜是一樣的；當(dāng)NM時，BT法的偏差大于周期圖法，在窗函數(shù)滿足一定條件下是漸進(jìn)無偏估計；方差小于周期圖的方差；分辨率比周期圖法低，與窗函數(shù)選擇有關(guān)。令W為KW的DTFT,1MKKIKIEE10利用DTFT的性質(zhì)，我們可以寫成KIBTER,01,0MW和,NRR的積的DTFTKDTFKTF11記DTF0,1,0NRRP，可得DWWPBT212由于極大多數(shù)常用的窗函數(shù)的在0處都有一個相對較窄的主峰。由12式可得BT譜估計子9式就是一個局部加權(quán)的平均周期圖。232BT譜估計的性質(zhì)由于0，自然也要求0BT。選擇適當(dāng)?shù)拇昂瘮?shù)就能夠使估計普擁有這樣的性質(zhì)。如果滯后窗KW式正半定的（即0W），則加窗的協(xié)方差序列RK也是正半定的，這意味著對所有的，0BT13即BT譜估計具有非負(fù)性。證明既然KW對于是實對稱的，那么其DTFW是一個實的偶函數(shù)。進(jìn)一步是一個正半定序列，則對于所有，0。由于定義0P，由12式可知，0W即意味著BT。應(yīng)該注意的是，有些滯后窗，如三角窗，并不滿足13式中的假設(shè)，因此，這些窗函數(shù)的使用可能導(dǎo)致估計譜出現(xiàn)負(fù)值。233BT譜估計在MATLAB中的仿真實現(xiàn)仿真實例2用BT法對信號2130SIN20SINNWX進(jìn)行功率譜估計，其中NW為白噪聲。對信號X的采樣率為SF1HZ，滿足取樣定理，N0,N1。程序代碼見附錄。當(dāng)N512時，MATLAB仿真得到的功率譜估計圖為00102030405060708091505101520253035三三三HZ三三三三三DB三512三三三三三三圖2BT法得出的功率譜估計圖通過實驗仿真可以直觀地看出功率譜估計中BT法和周期圖法所得到的結(jié)果是一致的，其特點是離散性大，曲線粗糙，方差較大，但是分辨率較高。24平均周期圖法周期圖是信號功率譜的一個有偏估值；而且，當(dāng)信號序列的長度增大到無窮時，估值的方差不趨于零。因此，隨著所取的信號序列長度的不同，所得到的周期圖也不同，這種現(xiàn)象稱為隨機(jī)起伏。由于隨機(jī)起伏大，使用周期圖不能得到比較穩(wěn)定的估值。一些學(xué)者對此作了改進(jìn)。為了減小隨機(jī)起伏，MSBARTLETT特提出平均周期圖法，即先把信號序列分為若干段，對每段分別計算其周期圖，然后取各個周期圖的平均作為功率譜的估值。平均周期圖可以減小隨機(jī)起伏，但是，如果信號序列不是足夠長，由于每段序列長度變短，功率譜估值對不同頻率成分的分辨能力也隨之下降。另一種改進(jìn)方法是將周期圖與一個適當(dāng)?shù)念l域窗函數(shù)相褶積，從而對周期圖產(chǎn)生平滑作用，以減小隨機(jī)起伏。加窗處理的結(jié)果雖然可以使隨機(jī)起伏減小，但也會使周期圖的分辨能力下降。POWELCH提出一種把加窗處理與平均處理結(jié)合起來的方法。先把分段的數(shù)據(jù)乘以窗函數(shù)進(jìn)行加窗處理，分別計算其周期圖，然后進(jìn)行平均。WELCH方法是較常用的一種計算方法。為了得到較好的功率譜估值，加窗和平均處理均應(yīng)兼顧減小隨機(jī)起伏和保證有足夠的譜分辨率兩個方面。241BARTLETT方法及在MATLAB中的仿真實現(xiàn)（一）BARTLETT方法BARTLETT方法的基本思想很簡單為了減少周期圖的起伏，把長度為N的觀測樣本分割為長度為M的NL/個子樣本，然后對從這些子樣本得到的周期圖對于每一個取平均。將BARTLETT方法表示成如下的數(shù)學(xué)形式。令,1TJYTJLJM,114表示第J個子樣本的觀測值，并令21MITIJJEY15表示對應(yīng)的周期圖。BARTLETT譜估計由下式給出LJB116由于BARTLETT方法是在長度為M的數(shù)據(jù)段上計算的，故其分辨率應(yīng)該約為M/1。因此，與與原來周期圖方法比較，BARTLETT方法的分辨率降低了L倍。作為對分辨率降低的回報，可以預(yù)料，BARTLETT方法的方差減小了。事實上，可以發(fā)現(xiàn)BARTLETT方法把周期圖的方差減小了L倍。在選擇M（或L）時，分辨率和方差之間的折中是很明顯的。（二）BARTLETT方法的性質(zhì)正如我們所知道的，16式中的J可以寫成1MMKKIJJER17其中JR是第J個子樣本的協(xié)方差序列。把17式代入16式得KIMKLJBER1118我們發(fā)現(xiàn)B在形式上與使用矩形窗的BT估計很相似。子樣本協(xié)方差KRJ在J上的平均是ACSKR的一個估計。但是，在18式中ACS的估計并沒有有效地利用可以得到的數(shù)據(jù)滯后的乘積KTY，特別是當(dāng)K接近1M時。實際上，對于1MK，只有大約可以滯后乘積的/1被用來形成18式中ACS的估計?？梢灶A(yù)料，這些滯后的方差要大于BT估計中對應(yīng)的KR滯后，類似地，B方差要比BT的方差大。此外，由于BARTLETT方法使用了一個固定的矩形滯后窗，所以其分辨率泄漏折中的靈活性也不如BT法。由于上述原因，我們總結(jié)為在形式上，1416式定義的BARTLETT估計與使用長度為M的固定矩形滯后窗的BT估計類似，但其方差略大于后者。由此和以前討論的BT譜估計的性質(zhì)可知，與基本的周期圖方法比較，BARTLETT估計降低了分辨率，減少了方差（其倍數(shù)均為MNL/）。矩形窗的主瓣要比極大多數(shù)其他滯后窗的主瓣窄。所以，在BT類譜估計中，BARTLETT估計可望得到最少的譜峰兼并（即有最大的分辨率），但有最明顯的泄漏。（三）BARTLETT估計的MATLAB仿真仿真實例3用BARTLETT估計對信號TNTTX102SI602SIN進(jìn)行功率譜估計，其中TN為高斯白噪聲。對信號的采樣頻率為SF500HZ，滿足抽樣定理。N500，源程序代碼見附錄。得到的功率譜估計圖如下BARTLETT估計對周期圖法德改進(jìn)的思想是將信號分段進(jìn)行估計，并將這些估計結(jié)果進(jìn)行平均而減少估計的方差，使得估計功率譜圖變得平滑。如上例中將501點的信號分為三段，分別作周期圖法估計，然后加以平均。功率譜估計結(jié)果如圖3所示，將信號分段得到的估計值顯然平滑了許多。0501001502002503003504004505006040200三三三HZ三三三三三DB三三三三三三三三500三三三三0501001502002503003504004505001020304050三三三HZ三三三三三DB三三500三三三三三三三三三三三三圖3BARTLETT法頻率譜估計圖242WELCH方法及在MATLAB中的仿真實現(xiàn)一WELCH譜估計WELCH方法是BARTLETT方法經(jīng)過兩個方面的改進(jìn)得到的。首先，WELCH方法中的數(shù)據(jù)段允許重疊。其次，在計算在周期圖前，對每一個數(shù)據(jù)段加窗。為了用數(shù)學(xué)方式描述WELCH方法，令TKJYTJ1，SJM,119表示第J個數(shù)據(jù)段。在19式中，K是第J個序列起始點。若MK，則序列間沒有重疊（但是都是鄰接的），我們得到用于BARTLETT方法中的樣本分割（得到MNLS/個數(shù)據(jù)子樣本）。但是，在WELCH方法中的推薦值為2/，在這種情況下，得到NS/2個數(shù)據(jù)段（前后數(shù)據(jù)段之間重疊50）。計算與TYJ對應(yīng)的加窗周期圖21MITIJJEYTP20其中P表示時間窗T的功率21MIT21通過對20式中的加窗周期圖取平均得到PSD的WELCH估計SJW122要解釋上述對BARTLETT方法的改進(jìn)，從而導(dǎo)出WELCH方法的理由很簡單。通過允許數(shù)據(jù)段之間的重疊，就有更多的周期圖參與22式中的平均，故可以降低PSD估計的方差。在周期圖計算中引入窗函數(shù)，就能夠更有效地控制PSD估計中的偏差/分辨率性質(zhì)。除此之外，時間窗還可用來對位于沒個子樣本末端的數(shù)據(jù)給予較小的加權(quán)，這樣，盡管子樣本間是重疊的，但可以使相鄰的子樣本間的相關(guān)性減弱。通過22式中平均。這種去相關(guān)的重要作用能夠更有效地降低方差。（二）WELCH譜估計的MATLAB仿真仿真實例4用WELCH譜估計對信號TNTTX102SI602SIN進(jìn)行功率譜估計，其中TN為高斯白噪聲。對信號的采樣頻率為SF500HZ，滿足抽樣定理。N500，源程序代碼見附錄。得到的功率譜估計圖如下0501001502002503003504004505006040200三三三HZ三三三三三DB三三三三三三三三500三三三三050100150200250300350400450500806040200三三三HZ三三三三三DB三三500三三三三三三三三三三三三三三三三三三三三圖4采用樣本混疊增加分段數(shù)得到的功率譜050100150200250300350400450500806040200三三三HZ三三三三三DB三三三三三三三三500三三三三050100150200250300350400450500806040200三三三HZ三三三三三DB三三三WELCH三三三三三三三三三三三三三三三三圖5WELCH平均修正周期圖法得到的功率譜估計圖增加分段數(shù)可以進(jìn)一步降低估計的方差，然而若每段中的數(shù)據(jù)點數(shù)太少，就會使得估計的頻率分辨率下降很多。在樣本信號總點數(shù)一定的條件下，我們可以采用使分段相互重疊的方法來增加分段數(shù)，而保持每段中信號點數(shù)不變，這樣就在保證頻率分辨率的前提下進(jìn)一步降低估計的協(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方法對信號2130SIN20SINNWX進(jìn)行功率譜估計，其中NW為白噪聲。對信號X的采樣率為F1HZ，滿足取樣定理，N0,N1。程序代碼見附錄。當(dāng)N1000時，窗函數(shù)采用矩形窗、HANNING窗、HAMMING窗，MATLAB仿真得到的功率譜估計圖分別為0010203040506070809150050三三三HZ三三三三三DB三三三三0010203040506070809102040三三三HZ三三三三三DB三HANNING三0010203040506070809102040三三三HZ三三三三三DB三HAMMING三圖6WELCH法加不同的窗函數(shù)得出的功率譜估計圖仿真結(jié)果分析由圖6可以看出，使用不同的窗函數(shù)譜估計的質(zhì)量是不一樣的，矩形窗的主瓣較窄，分辨率較好，但方差較大，噪聲水平較高；而HANING窗和HAMMING窗的主瓣較寬，分辨率較低，但方差較小，噪聲水平較低。因此，在進(jìn)行譜分析時選擇何種窗函數(shù)，要視具體情況而定。如果強(qiáng)調(diào)高分辨率，能精確讀出主瓣頻率，而不關(guān)心幅度的精度，例如測量震動物體的自震頻率，可以選用主瓣寬度比較窄的矩形窗；對受到強(qiáng)干擾的窄帶信號，若干擾靠近信號，則可選用旁瓣幅度較小的窗函數(shù)，若離開通帶較遠(yuǎn)，則可選用漸近線衰減速度比較快的窗函數(shù)?？傊?，要針對不同的信號和不同的處理目的來選擇合適的窗函數(shù)，這樣才能得到良好的效果。25小結(jié)本章先研究周期圖譜估計，從均值，方差，分辨率討論周期圖法的性能。估計的均值得出當(dāng)N時，是無偏估計，當(dāng)N為有限值時，是有偏估計。估計的方差得出不是的一致估計。隨著N增大，譜估計的起伏增大。N增大時，分辨率提高，但會使Y的起伏加劇，可見方差與分辨率是一對矛盾。為此，提出了改進(jìn)周期圖法的幾種途徑，即改善窗口形狀、平均和平滑。BT法中涉及到兩次加窗一次是用于截取數(shù)據(jù)，叫數(shù)據(jù)窗；一次是用于截取相關(guān)函數(shù)，叫延遲窗，延遲窗常用三角窗。相關(guān)法的兩次加窗將造成泄漏效應(yīng)，影響到譜峰分辨率和方差。平均就是將截取的數(shù)據(jù)段再分成L個小段，分別計算功率譜后取功率譜的平均，這種方法使估計的方差減少，但偏差加大，分辨率下降。平滑是用一個適當(dāng)?shù)拇昂瘮?shù)與算出的功率譜進(jìn)行卷積，使譜線平滑。這種方法得出的譜估計是無偏的，方差也小，但分辨率下降?，F(xiàn)在比較常用的改進(jìn)方法是WELCH法，又叫加權(quán)交疊平均法，簡記為WOSA法，這種方法以加窗（加權(quán)）求取平滑，以分段重疊求得平均，因此集平均與平滑的優(yōu)點于一體，同時也不可避免帶有兩者的缺點，因此歸根到底是一種折中。WELCH法得出的是無偏譜估計，當(dāng)段數(shù)L增大時方差減小，趨于一致估計，分段平均減小了由數(shù)據(jù)樣本本身的隨機(jī)性帶來的方差，段數(shù)越多方差越小，但分辨率下降。非參數(shù)譜估計法的基礎(chǔ)是FFT，它們出現(xiàn)較早，運算量較小，物理概念明確，便于工程實現(xiàn)，仍是目前常用的譜估計方法，有的已被固化到FFT儀和譜分析儀中。但有一些難以克服的缺點譜的分辨率較低，它正比于（N是數(shù)據(jù)長度）。加窗的壞影響不可避免。較寬的主瓣降低分辨力，較大的旁瓣有可能掩蓋真實譜中較弱的成分，或是產(chǎn)生假的峰值。而沒有一個窗函數(shù)能使譜估計在方差、偏差和分辨率各方面同時改善的，使用窗函數(shù)只是一種折中的技巧，不是解決問題的根本辦法。不是真實譜的一致估計，且N增大時譜曲線起伏加劇?；谏鲜鲈?，非參數(shù)譜估計還有繼續(xù)研究的必要性。3空間譜估計的算法與實現(xiàn)31引言本章將討論使用一個M元無源陣元組成的陣列對N個輻射源進(jìn)行定位的問題。這里，從輻射源來的能量可能是聲音、電磁波等等，接收傳感器將該能量轉(zhuǎn)換為電信號。接收陣元可能為電磁天線、水聲傳感器或者地震檢波器。這類問題在雷達(dá)或聲納系統(tǒng)、通信、天體物理學(xué)、生物醫(yī)學(xué)、地震水下監(jiān)視以及其他許多領(lǐng)域具有廣泛的應(yīng)用。這一問題的基本內(nèi)容包括確定能量在空間（空氣、水下或地下）的分布，這里能量高度集中的點就是輻射源。這類問題稱為空間譜估計。之所以稱為譜估計還因為輻射源定位問題和前面討論的譜估計有著非常緊密的聯(lián)系。事實上我們將看到，前面所述的方法幾乎都可以用作推導(dǎo)輻射源定位問題的解。32全向天線的波束下傾321全向天線的波束下傾的基本原理廣播用的發(fā)射機(jī)天線建在天線塔上，希望覆蓋廣大的地域。在覆蓋區(qū)的邊沿，由于距離遠(yuǎn)，地球表面彎曲，電磁波場強(qiáng)衰減很快，因此為了有效地覆蓋既定的區(qū)域，保證區(qū)域內(nèi)場強(qiáng)不低于特定值，通常采用波束下傾的方法，將軸向排列的半波振子天線，通過調(diào)節(jié)天線的軸向間距、饋電的相位，使得軸向天線陣列的方向圖實現(xiàn)波束下傾。軸向排列的半振子天線結(jié)構(gòu)，陣元之間的距離是D，垂直軸線與電磁波輻射方向的夾角是，相鄰陣元饋入信號的相位差是，相鄰陣元發(fā)出電磁波到達(dá)同一地點的程差是COSD，相鄰陣元發(fā)出的電磁波到達(dá)同一地點的相位差參看公式（23）2COS1120DEEENJJJ（23）322全向天線的波束下傾的仿真實現(xiàn)根據(jù)以上討論，編寫出繪制軸向排列天線陣列的方向圖程序。源程序代碼見附錄。仿真結(jié)果如下圖24683021060240902701203001503301800圖7軸向陣列波束下傾仿真圖33天線的波達(dá)方向估計天線陣列信號處理作為信號處理的一個重要分支，在通信、雷達(dá)、地震勘測、射電天文等領(lǐng)域獲得了廣泛應(yīng)用和迅速發(fā)展。它包括信源定位確定陣列到信源的仰角和方位角，甚至距離（若信源位于近場）信源分離確定各個信源發(fā)射的信號波形。各個信源從不同方向到達(dá)陣列，即使它們在時域和頻域上是疊加的。信道估計確定信源與陣列之間的傳輸信道的參數(shù)（多徑參數(shù)）。陣列信號處理時改善蜂窩和個人通信服務(wù)系統(tǒng)質(zhì)量、覆蓋范圍和容量的一個強(qiáng)有力的工具。建立在波達(dá)方向估計、波束形成基礎(chǔ)上的智能天線的應(yīng)用，抑制了干擾信號，改善了欲接收信號的信噪比，降低了數(shù)字通信系統(tǒng)的誤碼率。將接收天線陣列用于公眾網(wǎng)的反向連接（客戶到基站）時，多個接收天線能夠收集更多的信號能量，若天線在空間足夠分離或極化各異，則多個天線能夠提供很好的分集接收，并抑制多徑傳輸引起的衰落。這些好處可以擴(kuò)大基站的覆蓋范圍，改善通信質(zhì)量。331天線的波達(dá)方向估計的基本原理全向天線不僅利用率不高，而且對各種信號不加區(qū)別地接收，降低了通信質(zhì)量。定點無線通信采用定向天線，大幅度地改善了通信質(zhì)量。面對眾多移動用戶的公眾通信網(wǎng)基站和專用移動通信網(wǎng)，采用天線指向即波束可變的天線（智能天線），可以使移動通信的質(zhì)量得到很大的改善。為使天線的波束指向可控，甚至形狀可控，采用陣列天線是合適的。在距離通信源足夠遠(yuǎn)的空間里，可以將到達(dá)的電磁波視為平面波。對于等距離直線陣天線，由于調(diào)制在載波上的基帶信號碼元寬度與波速的乘積遠(yuǎn)大于天線陣的尺寸，因此多個天線陣元上的信號的幅度可視為不變，而它們的載波的相位差則取決于其相互位置、尺寸、波長和到達(dá)方向。無線接收的無線信號中有許多成分，其中我們關(guān)心的是S信號。天線陣列各個陣元接收的電磁波信號因為陣元排列位置的不同帶來相位差。經(jīng)過特定參數(shù)的加權(quán)控制器W處理后，進(jìn)一步改變了各個陣元輸出信號的相位和幅度。處理的目標(biāo)是使得陣元輸出信號的和天線輸出Y中的S成分具有最大輸出。用S信號作為基準(zhǔn)信號，反饋控制單元的功能就是將輸出信號Y與基準(zhǔn)信號S的差值（即誤差信號），作為調(diào)節(jié)控制加權(quán)控制器W參數(shù)的依據(jù)。反饋控制是使減小，Y中的S成分加大，也就是說，天線陣列接收方向圖指向了信號的方向。以下通過兩個簡單的例子介紹智能天線中波達(dá)方向估計和波束形式的原理。等距離直線陣智能天線，我們把天線陣元沿X軸排列，從0到1M。若有一平面波以角（入射線與Z軸的夾角）和角（入射線與軸的夾角）入射到陣列上，第K號陣元上產(chǎn)生的信號為KX，它與0號陣元的相位差為COSIN2DA（24）式中，和D分別是入射波的波長和陣元的間距，KA亦稱陣因子。計入陣因子的影響，第K號陣元的輸出時KXA，即U。為了使天線陣的輸出滿足需要，在每個陣元上，用加權(quán)因子W進(jìn)行控制。這樣第號陣元上輸出的信號為KXAW，即KUW。若到達(dá)天線陣的信號是N個，則天線陣的輸出是N個信號在M個陣元上的輸出的疊加。將問題簡化為XY平面的二維問題（1SIN），并用解析式表達(dá)如下NXNXXN,21（25）111COS2EXPSDMJJA（26）其中，1為第一個信號的入射角。NA,21（27）MW（28）UAXTYHH（29）其中，HMUU110,。求多個信號到達(dá)的方向（波達(dá)方向）的方法有多種。下面討論其中兩種方法（CAPON法和MUSIC算法）及其仿真驗證結(jié)果。332求波達(dá)方向的方法1CAPON法CAPON法亦稱最小方差無畸變響應(yīng)MVDR（MINIMUMVARIANCEDISTORTIONLESSRESPONSE）。天線陣列中的陣元數(shù)決定了陣列方向圖設(shè)計中的自由度數(shù)。CAPON法將陣列中可控的自由度用來形成期望的波束形狀，達(dá)到對有用信號進(jìn)行提升和對無用信號進(jìn)行抑制的目的，并將其優(yōu)化問題表達(dá)為WRNYEUHWWMII2（30）其約束條件為10AWH可以證明上式的解為ARUHCAP1代入式（29），可以得到相應(yīng)到的功率為ARPUHCAP1（31）2MUSIC法MUSIC法亦稱多重信號分類（MULTIPLESIGNALCLASSIFICATION）。當(dāng)入射信號數(shù)N小于陣元數(shù)M時，的UR矩陣有與進(jìn)入陣列的信號數(shù)目相等的非零特征值及個為零的特征值。1NNVV2MV（NV是U相應(yīng)的噪聲特征矢量）因為A與的正交性，分母很小，峰值很大，這樣可以得出MUSIC法的空間譜為APHMUSIC（32）其中，HNV稱為噪聲子空間的正交投影估計。3仿真以下是對七單元線性天線陣在四信號輸入情況下了仿真仿真實例6設(shè)信號1從4/方向入射，信號2從3/方向入射，信號3從6/方向入射，信號4從/3方向入射。下面是用兩種方法，分別用兩種坐標(biāo)（直角坐標(biāo)和極坐標(biāo)）求波達(dá)方向估計。程序代碼見附錄。仿真結(jié)果如下49999999999999998000000000000009999999999999999600000000000000149999999999999990000000000000003021060240902701203001503301800圖8CAPON法作出的波達(dá)方向估計（極坐標(biāo)）020406080100120140160180100105101010151020102510301035圖9MUSIC法作出的波達(dá)方向估計（直角坐標(biāo)）500000100000015000003021060240902701203001503301800圖10MUSIC法作出的波達(dá)方向估計（極坐標(biāo)）020406080100120140160180101510101051001051010圖11CAPON法作出的波達(dá)方向估計（直角坐標(biāo)）從對以上四圖的分析可以得到如下結(jié)論兩種方法估計得都比較準(zhǔn)確。在4/3,/,6處有尖銳的方向圖線。MUSIC法方向圖像幅度更大。34天線陣的波束形成1原理我們以等距離圓陣為例來討論天線陣的波束形成。我們?nèi)A陣的圓心為O，在陣列圓上任取一點B，把天線陣元順序定為0到1M。若有一平面波以角入射到陣列上，第K號陣元上產(chǎn)生的信號為KX，它與到達(dá)陣元中心的波前的相位差是KKRJACOS2EXP；與R分別是入射波的波長與陣列圓的半徑，亦稱陣因子。為了使天線陣的輸出滿足需要，在每個陣元上加上加權(quán)因子K控制。這樣第號陣元上輸出的信號為KXA。若到達(dá)天線陣的信號是N個，天線陣輸出的是N個信號的M個陣元上的輸出的疊加。用解析表達(dá)式如下COS2EXP1KKRJA（33）HMKA,NNXNXX110（34）式中，K是陣元以O(shè)B為基準(zhǔn)順時針畫出的角度。NA,21（35）HMWW10（36）UWAXTYH（37）為了求得多個信號到達(dá)的方向（波達(dá)方向），可以采用上述的CAPON、MUSIC兩種方法）。波束形成可以采用下面的方法當(dāng)有多（N）個信號輸入時，其中1個信號時我們關(guān)心的，1N個信號是需要抑制的。方程組（38）描述了上述需求的約束條件（四個信號輸入，第一個信號是我們關(guān)心的，其余信號是需要抑制的）。THUW0,1（38）根據(jù)信號波達(dá)方向的知識AX及約束條件求解方程組（38），可以得到110,MW（39）代入式（37）可得到陣列輸出的方向特性。2仿真下面是四個輸入信號八單元圓陣列天線陣的波達(dá)方向估計和波束形成。仿真實例7設(shè)信號1從3/2方向入射，信號2從3/方向入射，信號3從2/3方向入射，信號4從685方向入射。下面是用兩種算法求波達(dá)方向估計及波束形成。源程序代碼見附錄。仿真結(jié)果如下050100150200250300350400100105101010151020102510301035圖12CAPON算法波達(dá)方向估計050100150200250300350400101510101051001051010圖13MUSIC算法波達(dá)方向估計051153021060240902701203001503301800圖14波束形成35小結(jié)本章在對非參數(shù)譜估計的理解基礎(chǔ)上，討論了其在空間領(lǐng)域進(jìn)行推廣，具體就是了采用CAPON算法和MUSIC算法求天線陣的波達(dá)方向，討論了天線陣的波束形成，及其MATLAB的仿真實現(xiàn)。4結(jié)論本設(shè)計的主要任務(wù)是研究非參數(shù)譜估計的算法與實現(xiàn)，另外對功率譜估計算法在空間領(lǐng)域的推廣應(yīng)用尤其是陣列側(cè)向問題進(jìn)行了討論。在經(jīng)典譜估計中，無論是周期圖法還是其改進(jìn)的方法，都存在著頻率分辨率低、方差性能不好的問題，原因是譜估計時需要對數(shù)據(jù)加窗截斷，用有限個數(shù)據(jù)或其自相關(guān)函數(shù)來估計無限個數(shù)據(jù)的功率譜，這其實是假定了窗以外的數(shù)據(jù)或自相關(guān)函數(shù)全為零，這種假定是不符合實際的，正是由于這些不符合實際的假設(shè)造成了經(jīng)典譜估計分辨率較差。另外，經(jīng)典譜估計的功率譜定義中既無求均值運算又無求極限運算，因而使得譜估計的方差性能較差，當(dāng)數(shù)據(jù)很短時，這個問題更為突出。如何選取最佳窗函數(shù)、提高頻譜分辨率，如何在短數(shù)據(jù)情況下提高信號譜估計質(zhì)量，還需要進(jìn)一步研究。CAPON法和MUSIC算法可以用于空間譜估計測向，充分利用了天線各個陣元從空間電磁場接收到的全部信息，而傳統(tǒng)測向方法僅僅利用了其中一部分信息，因而具有多方面的優(yōu)越性1可以實現(xiàn)對同一一信道中存在的多個信號同時測向；2超分辨率測向，即對處于天線陣固有波束寬度以內(nèi)的兩個以上方向的來波同時測向；3陣元位置可以隨意放置，各個陣元的方向特性也可以是任意的；4具有優(yōu)良的測向靈敏度和準(zhǔn)確度；5可以較好地克服測向模糊問題。參考文獻(xiàn)1PETRESTOICA,RANDOLPHLMOSESINTRODUCTIONTOSPECTRALANALYSIS,PRENTICEHALL,INC19972魏平，唐斌，陸明泉，唐海燕，談覓舟，肖先賜譜估計原理J電子科技大學(xué)研究生教材，19993胡廣書數(shù)字信號處理一理論、算法與實現(xiàn)M北京清華大學(xué)出版社，20034王曉峰，王炳和周期圖及其改進(jìn)方法中譜分析率的MATLAB分析J武警工程學(xué)院學(xué)報，2003675775姚武川，姚天任經(jīng)典譜估計方法的MATLAB分析J華中理工大學(xué)學(xué)報，2000，28445476劉剛，呂新華，攸陽陣列信號處理中基于MUSIC算法的空間譜估計J微計算機(jī)信息，2006，2243302303附錄附錄直接法CLEARALLCLCN256FS1N256N0N1XNSQRT20SIN2PI02NSQRT2SIN2PI0213NRANDNSIZENNFFT256PXX10LOG10ABSFFTXN,NFFT2/NF0LENGTHPXX1FS/LENGTHPXXSUBPLOT211PLOTF,PXXXLABEL頻率（HZ）YLABEL功率譜（DB）TITLE256點功率譜GRIDN1024FS1N1024N0N1XNSQRT20SIN2PI02NSQRT2SIN2PI0213NRANDNSIZENNFFT1024PXX10LOG10ABSFFTXN,NFFT2/NF0LENGTHPXX1FS/LENGTHPXXSUBPLOT212PLOTF,PXXXLABEL頻率（HZ）YLABEL功率譜（DB）TITLE1024點功率譜GRID附錄間接法，512點功率譜估計CLEARALLCLCFS1N512N0N1XNSQRT20SIN2PI02NSQRT2SIN2PI0213NRANDNSIZENNFFT512PXX10LOG10FFTXCORRXN,UNBIASED,NFFTF0LENGTHPXX1FS/LENGTHPXXPLOTF,PXXXLABEL頻率（HZ）YLABEL功率譜（DB）TITLE512點功率譜估計GRID附錄BARTLETTCLEARALLCLCN500FS500T01/FS1F01FSXNSIN2PI60T2SIN2PI110TRANDN1,LENGTHTPXX1ABS1/FSFFTXN2PXX2ABSFFTXN11672ABSFFTXN1683342ABSFFTXN3355012/3SUBPLOT211PLOTF,10LOG10PXX1XLABEL頻率（HZ）YLABEL功率譜（DB）TITLE周期圖法得到的500點功率譜SUBPLOT212PLOT03FS,10LOG10PXX2XLABEL頻率（HZ）YLABEL功率譜（DB）TITLE將500點信號分三段得到的功率譜附錄WELCHCLEARALLCLCN500FS500T01/FS1F01FSXNSIN2PI60T2SIN2PI110TRANDN1,LENGTHTPXX1ABS1/FSFFTXN2PXX3ABS3/FSFFTXN11672ABS3/FSFFTXN832492ABS3/FSFFTXN1683342ABS3/FSFFTXN2504162ABS3/FSFFTXN3355012/5SUBPLOT211PLOTF,10LOG10PXX1XLABEL頻