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

1、0 / 41非參數(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)展開(kāi)研究。 首先討論了一種最常見(jià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ì)是功率譜

2、估計(jì)算法在空間領(lǐng)域的推廣應(yīng)用，具體就是以多元天線(xiàn)陣結(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á)方向Algorithms and Implementations of Non-Parameter Spectral EstimationsAbstract：Power spectrum estimate is an important part of the random signal analysis, which primarily s

3、tudies the random signal characteristics in the frequency domain and aims at recovering the useful signal in the frequency domain contaminated with noise based on the limited observed data.The design focuses on the algorithm and implementation of the non-parameter spectral estimations. First, the mo

4、st general spectral estimation periodogram method is discussed. As the periodogram method is a biased estimate with the large variance and low resolution, the improved methods based on the periodgram method need to be studied. The main improvements include changing the shape of windows, averaging an

5、d 1 / 41smoothing the data, etc. As a result, several other non-parameter spectral estimation algorithms based on the periodgram method (such as the BT method, the averaging periodgram method, including the Bartlett method and Welch method) and their realizations by the software of MATLAB are given.

6、Moreover, the direction of arrival (DOA) estimation is to apply the power spectrum estimation algorithms into the space domain, i.e., DOA is the new direction finding technique based on the multiple sensor array combined with the modern digital signal processing. Accordingly, the realizations of DOA

7、 estimation algorithms based on the method of Capon and MUSIC are discussed and compared in this design.Key words: Periodgram Averaging Periodgram BT Spectrum Capon method Music Algorithm Direction of arrival (DOA)目 錄1 緒 論 .11.1 引言.11.1.1 非參數(shù)譜估計(jì).11.1.2 空間譜估計(jì).11.2 本文的主要工作.22 幾種常用的非參數(shù)譜估計(jì)算法與實(shí)現(xiàn).22.1 引言.

8、22.2 周期圖法.32.2.1 周期圖譜估計(jì).32.2.2 周期圖法的性質(zhì).32.2.3 周期圖法在MATLAB中的仿真實(shí)現(xiàn).42.3 BT 法.52.3.1 BT法譜估計(jì).52.3.2 BT譜估計(jì)的性質(zhì).62.3.3 BT譜估計(jì)在MATLAB中的仿真實(shí)現(xiàn).72.4 平均周期圖法.72.4.1 Bartlett方法及在MATLAB中的仿真實(shí)現(xiàn).82.4.2 Welch方法及在MATLAB中的仿真實(shí)現(xiàn).102.5 小結(jié).133 空間譜估計(jì)的算法與實(shí)現(xiàn).143.1 引言.143.2 全向天線(xiàn)的波束下傾.143.2.1 全向天線(xiàn)的波束下傾的基本原理.143.2.2 全向天線(xiàn)的波束下傾的仿真實(shí)現(xiàn).1

9、53.3 天線(xiàn)的波達(dá)方向估計(jì).153.3.1 天線(xiàn)的波達(dá)方向估計(jì)的基本原理.163.3.2 求波達(dá)方向的方法.173.4 天線(xiàn)陣的波束形成.203.5 小結(jié).234 結(jié) 論.23致謝.23參考文獻(xiàn)：.24附錄.24非參數(shù)譜估計(jì)的算法與實(shí)現(xiàn)1 緒 論1.1 引言1.1.1 非參數(shù)譜估計(jì)信號(hào)的頻譜分析是研究信號(hào)特性的重要手段之一。對(duì)于確定性信號(hào)，可以用傅里葉變換來(lái)考察其頻譜性質(zhì)。在通信系統(tǒng)中，最常見(jiàn)的信號(hào)往往不是確定信號(hào)，而是具有某種統(tǒng)計(jì)特性的隨機(jī)信號(hào)。對(duì)于隨機(jī)信號(hào)，由于它一般不是周期的，持續(xù)時(shí)間無(wú)限長(zhǎng)，具有無(wú)限長(zhǎng)能量的功率信號(hào)，不滿(mǎn)足傅里葉變換條件，隨機(jī)信號(hào)也不存在解析表達(dá)式，因此，對(duì)隨機(jī)信號(hào)，

10、我們不能像確定信號(hào)那樣進(jìn)行頻譜分析。然而，雖然隨機(jī)信號(hào)的頻譜不存在，但如果隨機(jī)信號(hào)是平穩(wěn)的，其相關(guān)函數(shù)就是確定的，那么相關(guān)函數(shù)的傅里葉變換就是它的功率譜密度函數(shù)，簡(jiǎn)稱(chēng)功率譜。通?？梢酝ㄟ^(guò)求隨機(jī)信號(hào)的功率譜來(lái)對(duì)這類(lèi)信號(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í)中非線(xiàn)性系統(tǒng)識(shí)別、物理光學(xué)中透鏡干涉、流體力學(xué)的內(nèi)波分析、太陽(yáng)黑子活動(dòng)周期研究等許多領(lǐng)域，發(fā)揮了重要作用。然而，實(shí)際應(yīng)用中的平穩(wěn)隨機(jī)信

11、號(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ì)又稱(chēng)經(jīng)典譜估計(jì)，主要方法有：周期圖法，平均周期圖法，BT法及改進(jìn)的周期圖法，如Bartlett法和Welch法。1.1.2 空間譜估計(jì)電磁波到達(dá)方向(DOA)估計(jì)在無(wú)線(xiàn)電通信、雷達(dá)、聲納和導(dǎo)航等領(lǐng)域有著廣泛的應(yīng)用。由于空間信號(hào)存在多源可能，加上電磁波傳播的多徑效應(yīng)，將會(huì)給波束法測(cè)向帶來(lái)較大測(cè)量誤差，并受干擾信號(hào)的影響嚴(yán)重。利用陣列天線(xiàn)進(jìn)行測(cè)向，其空間分辨率因陣列孔徑尺寸的選擇受經(jīng)典瑞利限的限制，在瑞利限以?xún)?nèi)的空間目標(biāo)則不能

12、分辨。為了突破瑞利限約束，提高角度分辨力，新的測(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ì)精度和分辨力。空間譜分析采用了類(lèi)似時(shí)域譜估計(jì)中的非線(xiàn)性處理，從而產(chǎn)生了一些特殊的算法，如多重信號(hào)分類(lèi)法(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ù)。這一算法的提出開(kāi)創(chuàng)

13、了空間譜估計(jì)算法研究的新時(shí)代，成為空間譜估計(jì)理論體系中的標(biāo)志性算法。這種算法也是本文會(huì)研究到的。空間譜估計(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法，平均周期圖法(Bartlett法以及Welch法)。從各種估計(jì)法的算法著手，結(jié)合MATLAB仿真，分析仿真結(jié)果，與周期圖法作比較

14、，總結(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ù)方法完全依賴(lài)于PSD的定義式 (1) ktiekr (2)211lim)(NitiNetyNE來(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)

15、度，周期圖法能夠獲得滿(mǎn)意的分辨率，但是不是可靠的譜估計(jì)算法，這是因?yàn)樗姆讲詈艽?，而且并不隨著數(shù)據(jù)長(zhǎng)度的增加而降低。 由于周期圖法的估計(jì)方差很大，促使人們不斷地研究新的改進(jìn)方法，試圖以降低分辨率為代價(jià)來(lái)獲得較小的估計(jì)方差。目前已經(jīng)提出了一些改進(jìn)的方法，我們僅討論其中最常見(jiàn)的幾種?？梢宰C明，在大量數(shù)據(jù)樣本的條件下，所有這些方法在性質(zhì)和性能方面或多或少是等價(jià)的。2.2 周期圖法2.2.1 周期圖譜估計(jì)一種信號(hào)功率譜密度估計(jì)方法。它的特點(diǎn)是：為得到功率譜估值，先取信號(hào)序列的離散傅里葉變換，然后取其幅頻特性的平方并除以序列長(zhǎng)度，由于序列的離散傅里葉變換具有周期性，因而這種功率譜也具有周期性，常稱(chēng)為周期

16、圖。Schuster于1899年首先提出周期圖法，也稱(chēng)直接法，取平穩(wěn)隨機(jī)信號(hào)Y(t)的有限個(gè)觀察值 y(O)， y(1)， ，y(T-1)對(duì)功率譜進(jìn)行估計(jì)，周期圖方 法依賴(lài)于PSD的定義(2)式。略去只有信號(hào)樣本信息 的條件下不能Nity1實(shí)現(xiàn)的期望和極限運(yùn)算，我們得到 （周期圖） ( 211NitipetyN3)用來(lái)確定時(shí)間序列中可能存在 “隱周期”，是周期圖譜估計(jì)子 (3)式最早的應(yīng)用之一，這也許可以認(rèn)為這種方法的由來(lái)。2.2.2 周期圖法的性質(zhì) 分析 的統(tǒng)計(jì)特性的重要性在于，證明了周期圖作為 PSD 估計(jì)的不 p可靠性，此外，還提供了一些如何改進(jìn)周期圖以便獲得更好的譜估計(jì)的內(nèi)部機(jī)理。我們

17、用兩個(gè)部分來(lái)分析周期圖方法的性質(zhì)：偏差分析和方差分析。通常用偏差和方差這兩種基本度量來(lái)刻劃估計(jì)子的性能，其主要原因是一個(gè)估計(jì)的總平方誤差就是它的偏差的平方與方差之和。為了說(shuō)明這個(gè)問(wèn)題，設(shè)a 為任意一個(gè)帶估計(jì)量， 是啊的一個(gè)估計(jì)。估計(jì)的均方差（MSE）為a MSE22aaEaEaEaaE =aaEaEaEaaEaEaERe22 = (4)2varabiasa通過(guò)分別考慮 MSE 的偏差和方差分量，我們可以獲得誤差來(lái)源的深入了解，并得到減少誤差的方法。（一） 估計(jì)的均值： (5) WE21其中是窗函數(shù)的傅里葉變換。當(dāng)時(shí)，, W nN E是無(wú)偏估計(jì)，當(dāng) N 為有限值時(shí)，是有偏估計(jì)，偏差為： (6)

18、dWbias21（二） 估計(jì)的方差： (7) EdDDN20021var其中， 是矩形窗，的傅里葉變換，可見(jiàn)，不 0D 10nd1 Nn 是的一致估計(jì)。隨著 N 增大，譜估計(jì)的起伏增大，時(shí)， N。 2var（三）估計(jì)的分辨率：除了通過(guò)偏差和方差來(lái)分析周期圖法的性質(zhì)，我們還可以通過(guò)分辨率來(lái)分析。數(shù)據(jù)窗為長(zhǎng)度為 的矩形窗時(shí)，增大時(shí)，分辨N Ns289. 0ReN率提高，但會(huì)使的起伏加劇，可見(jiàn)方差與分辨率是一對(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í)，這

19、也意味著對(duì)自相關(guān)函數(shù)加三角窗，使功率譜與窗函數(shù)卷積，從而產(chǎn)生頻譜泄漏，容易使弱信號(hào)的主瓣被強(qiáng)信號(hào)的旁瓣所淹沒(méi)，造成頻譜的模糊和失真，使得譜分辨率較低。2.2.3 周期圖法在 MATLAB 中的仿真實(shí)現(xiàn)仿真實(shí)例仿真實(shí)例 1 1 用周期圖法對(duì)信號(hào)進(jìn)行功率譜估計(jì)，其中為白)()213. 02sin(2)2 . 02sin(20)(nwnnnx)(nw噪聲。對(duì)信號(hào)的采樣率為=1Hz，滿(mǎn)足取樣定理。程序代碼見(jiàn)附錄。 nxsfn=(0,N-1)，當(dāng)N=256和N=1024時(shí)，MATLAB仿真得到的功率譜估計(jì)圖分別為：01-40-2002040三 三 三

20、 Hz三三 三 三 三 dB三256三 三 三 三01-40-2002040三 三 三 Hz三三 三 三 三 dB三1024三 三 三 三圖 1 周期圖法得出的功率譜估計(jì)隨著采樣點(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)用很廣。2.3 BT

21、 法BT法即Blackman-Turkey法，這種方法是Blackman與Turkey兩人于1958年提出的，故以他們的名字命名。這一節(jié)，我們研究BT法，并與周期圖法作比較。2.3.1 BT 法譜估計(jì) 正如我們已經(jīng)看到的，即使有很長(zhǎng)的樣本長(zhǎng)度，方差仍然很大時(shí)周期圖是這種譜估計(jì)子存在的主要問(wèn)題。周期圖譜估計(jì)較差的統(tǒng)計(jì)質(zhì)量可以直觀地解釋為是由這樣兩個(gè)方面的原因同時(shí)引起的，即中末端滯后的估 cNk kr計(jì)精度低，以及大量（即便是少量）的協(xié)方差估計(jì)誤差在中累積相加。 c的定義公式 = （相關(guān)圖） c c 11NNkkiekr(8)是基于相關(guān)的PSD定義(1)直接導(dǎo)出的相關(guān)圖譜估計(jì)算法。通過(guò)在的定義 c

22、式(8)式中截短求和的方法可以減少這兩方面的影響。這個(gè)思想就導(dǎo)出了BT估計(jì)子， 11MMkkiBTekrkw(9)其中,是一個(gè)偶函數(shù)（），且，NM kw kwkw 10 w，同時(shí)隨平滑衰減到零。由于在式(9)中對(duì)樣本協(xié) Mkkw 0 kwk kw方差序列的滯后進(jìn)行加權(quán)，故被稱(chēng)為滯后窗。 因?yàn)橛蛇@種方法求出的功率譜是通過(guò)自相關(guān)函數(shù)間接得到的，所以稱(chēng)為間接法，又稱(chēng)為自相關(guān)法或BT法。當(dāng)較小時(shí)，式(9)的計(jì)算量不是很大，因M此，該方法是在FFT問(wèn)世之前常用的譜估計(jì)方法。 當(dāng)時(shí)，BT法與周期圖法估計(jì)的功率譜是一樣的；當(dāng) 時(shí)，BTNM NM 法的偏差大于周期圖法，在窗函數(shù)滿(mǎn)足一定條件下是漸進(jìn)無(wú)偏估計(jì)；方

23、差小于周期圖的方差；分辨率比周期圖法低，與窗函數(shù)選擇有關(guān)。令為的DTFT, W kw (10) W 11MMKkikiekekw利用DTFT的性質(zhì)，我們可以寫(xiě)成 kiBTekrk = 和,.0 , 0 , ,1,.,1, 0 , 0.,MwMw 的積的DTFT.0 , 0 ,1,.,1, 0 , 0.,NrNr = (11) kwDTFTkrDTFT記，可得DTFT.0 , 0 ,1,.,1, 0 , 0.,NrNr P (12) dWWPPBT21由于極大多數(shù)常用的窗函數(shù)的在處都有一個(gè)相對(duì)較窄的主峰。由 W0(12)式可得BT譜估計(jì)子(9)式就是一個(gè)局部加權(quán)的平均周期圖。2.3.2 BT 譜

24、估計(jì)的性質(zhì) 由于，自然也要求。選擇適當(dāng)?shù)拇昂瘮?shù)就能夠使估計(jì)普 0 0BT擁有這樣的性質(zhì)。 如果滯后窗式正半定的（即），則加窗的協(xié)方差序列 kw 0W也是正半定的，這意味著對(duì)所有的， krkw (13) 0BT即BT譜估計(jì)具有非負(fù)性。 證明證明 既然對(duì)于是實(shí)對(duì)稱(chēng)的，那么其 是一個(gè)實(shí)的 kw0kDTFT W偶函數(shù)。進(jìn)一步是一個(gè)正半定序列，則對(duì)于所有，。由于定義 kw 0W，由(12)式可知，即意味著。 0P 0W 0BT 應(yīng)該注意的是，有些滯后窗，如三角窗，并不滿(mǎn)足(13)式中的假設(shè)，因此，這些窗函數(shù)的使用可能導(dǎo)致估計(jì)譜出現(xiàn)負(fù)值。2.3.3 BT 譜估計(jì)在 MATLAB 中的仿真實(shí)現(xiàn) 仿真實(shí)例仿真

25、實(shí)例 2 2 用BT法對(duì)信號(hào)進(jìn)行功率譜估計(jì)，其中為白)()213. 02sin(2)2 . 02sin(20)(nwnnnx)(nw噪聲。對(duì)信號(hào)的采樣率為=1Hz，滿(mǎn)足取樣定理，n=(0,N-1)。程序代碼見(jiàn) nxsf附錄。當(dāng)N=512時(shí)，MATLAB仿真得到的功率譜估計(jì)圖為：01-505101520253035三 三 三 Hz三三 三 三 三 dB三512三 三 三 三 三 三圖2 BT法得出的功率譜估計(jì)圖通過(guò)實(shí)驗(yàn)仿真可以直觀地看出功率譜估計(jì)中BT法和周期圖法所得到的結(jié)果是一致的，其特點(diǎn)是離散性大，曲線(xiàn)粗糙，方差較大，但是分辨率較高。2.

26、4 平均周期圖法 周期圖是信號(hào)功率譜的一個(gè)有偏估值 ；而且，當(dāng)信號(hào)序列的長(zhǎng)度增大到無(wú)窮時(shí)，估值的方差不趨于零。因此，隨著所取的信號(hào)序列長(zhǎng)度的不同，所得到的周期圖也不同，這種現(xiàn)象稱(chēng)為隨機(jī)起伏。由于隨機(jī)起伏大，使用周期圖不能得到比較穩(wěn)定的估值。一些學(xué)者對(duì)此作了改進(jìn)。 為了減小隨機(jī)起伏， M.S. Bartlett 特提出平均周期圖法，即先把信號(hào)序列分為若干段，對(duì)每段分別計(jì)算其周期圖，然后取各個(gè)周期圖的平均作為功率譜的估值。平均周期圖可以減小隨機(jī)起伏，但是，如果信號(hào)序列不是足夠長(zhǎng)，由于每段序列長(zhǎng)度變短，功率譜估值對(duì)不同頻率成分的分辨能力也隨之下降。另一種改進(jìn)方法是將周期圖與一個(gè)適當(dāng)?shù)念l域窗函數(shù)相褶積

27、，從而對(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.1 Bartlett 方法及在 MATLAB 中的仿真實(shí)現(xiàn) （一）Bartlett方法Bartlett方法的基本思想很簡(jiǎn)單：為了減少周期圖的起伏，把長(zhǎng)度為的觀N測(cè)樣本分割為長(zhǎng)度為的個(gè)子樣本，然后對(duì)從這些

28、子樣本得到的周MMNL/期圖對(duì)于每一個(gè)取平均。將Bartlett方法表示成如下的數(shù)學(xué)形式。令 (14),1tMjytyjLjMt, 1, 1表示第個(gè)子樣本的觀測(cè)值，并令j (15) 211MitijjetyM表示對(duì)應(yīng)的周期圖。Bartlett譜估計(jì)由下式給出 (16) LjjBL11由于Bartlett方法是在長(zhǎng)度為的數(shù)據(jù)段上計(jì)算的，故其分辨率應(yīng)該約為。MM/1因此，與與原來(lái)周期圖方法比較，Bartlett方法的分辨率降低了倍。作為對(duì)分L辨率降低的回報(bào)，可以預(yù)料，Bartlett方法的方差減小了。事實(shí)上，可以發(fā)現(xiàn)Bartlett方法把周期圖的方差減小了倍。在選擇（或）時(shí)，分辨率和方差LML之間

29、的折中是很明顯的。 （二）Bartlett方法的性質(zhì)正如我們所知道的，(16)式中的可以寫(xiě)成 j (17) 11Mmkkijjekr其中是第個(gè)子樣本的協(xié)方差序列。把(17)式代入(16)式得 jrj (18) kiMMkLjjBekrL1111我們發(fā)現(xiàn)在形式上與使用矩形窗的BT估計(jì)很相似。子樣本協(xié)方差在 B krj上的平均是 ACS 的一個(gè)估計(jì)。但是，在(18)式中 ACS 的估計(jì)并沒(méi)有有j kr效地利用可以得到的數(shù)據(jù)滯后的乘積，特別是當(dāng)接近時(shí)。 ktytyk1M實(shí)際上，對(duì)于，只有大約可以滯后乘積的被用來(lái)形成(18)式中 1 MkM/1ACS 的估計(jì)?？梢灶A(yù)料，這些滯后的方差要大于BT估計(jì)中對(duì)

30、應(yīng)的滯后， kr類(lèi)似地，方差要比的方差大。此外，由于Bartlett方法使用了一個(gè) B BT固定的矩形滯后窗，所以其分辨率-泄漏折中的靈活性也不如BT法。由于上述原因，我們總結(jié)為：在形式上，(14)(16)式定義的Bartlett估計(jì)與使用長(zhǎng)度為的固定矩形滯M后窗的BT估計(jì)類(lèi)似，但其方差略大于后者。由此和以前討論的BT譜估計(jì)的性質(zhì)可知，與基本的周期圖方法比較，Bartlett估計(jì)降低了分辨率，減少了方差（其倍數(shù)均為）。矩形窗的主瓣要比極大多數(shù)其他滯后窗的主瓣窄。MNL/所以，在BT類(lèi)譜估計(jì)中，Bartlett估計(jì)可望得到最少的譜峰兼并（即有最大的分辨率），但有最明顯的泄漏。（三）Bartlet

31、t估計(jì)的MATLAB仿真仿真實(shí)例仿真實(shí)例 3 3 用Bartlett估計(jì)對(duì)信號(hào)進(jìn)行tntttx1102sin2602sin功率譜估計(jì)，其中為高斯白噪聲。對(duì)信號(hào)的采樣頻率為=500Hz，滿(mǎn)tntxsf足抽樣定理。N=500，源程序代碼見(jiàn)附錄。得到的功率譜估計(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ì)值顯然平滑了許多。050100150200250300350400450500-60-4

32、0-200三 三 三 Hz三三 三 三 三 dB三三 三 三 三 三 三 三 500三 三 三 三0501001502002503003504004505001020304050三 三 三 Hz三三 三 三 三 dB三三 500三 三 三 三 三 三 三 三 三 三 三 三圖3 Bartlett法頻率譜估計(jì)圖2.4.2 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方法，令 ， tKjytyj1S

33、jMt, 1, 1(19)表示第個(gè)數(shù)據(jù)段。在(19)式中，是第個(gè)序列起始點(diǎn)。若，則jkj1jMK 序列間沒(méi)有重疊（但是都是鄰接的），我們得到用于Bartlett方法中的樣本分割（得到個(gè)數(shù)據(jù)子樣本）。但是，在Welch方法中的推薦值為MNLS/K，在這種情況下，得到個(gè)數(shù)據(jù)段（前后數(shù)據(jù)段之間重疊50%）。2/MK NMS/2計(jì)算與對(duì)應(yīng)的加窗周期圖tyj (20) 211MitijjetytMP其中表示時(shí)間窗的功率：Pt (21)211MitMP通過(guò)對(duì)(20)式中的加窗周期圖取平均得到PSD的Welch估計(jì)： (22) SjjWS11要解釋上述對(duì)Bartlett方法的改進(jìn)，從而導(dǎo)出Welch方法的理

34、由很簡(jiǎn)單。通過(guò)允許數(shù)據(jù)段之間的重疊，就有更多的周期圖參與(22)式中的平均，故可以降低PSD估計(jì)的方差。在周期圖計(jì)算中引入窗函數(shù)，就能夠更有效地控制PSD估計(jì)中的偏差/分辨率性質(zhì)。除此之外，時(shí)間窗還可用來(lái)對(duì)位于沒(méi)個(gè)子樣本末端的數(shù)據(jù)給予較小的加權(quán)，這樣，盡管子樣本間是重疊的，但可以使相鄰的子樣本間的相關(guān)性減弱。通過(guò)(22)式中平均。這種去相關(guān)的重要作用能夠更有效地降低方差。（二）Welch譜估計(jì)的MATLAB仿真仿真實(shí)例仿真實(shí)例 4 4 用Welch譜估計(jì)對(duì)信號(hào)進(jìn)tntttx1102sin2602sin行功率譜估計(jì)，其中為高斯白噪聲。對(duì)信號(hào)的采樣頻率為=500Hz，tntxsf滿(mǎn)足抽樣定理。N=

35、500，源程序代碼見(jiàn)附錄。得到的功率譜估計(jì)圖如下：050100150200250300350400450500-60-40-200三 三 三 Hz三三 三 三 三 dB三三 三 三 三 三 三 三 500三 三 三 三050100150200250300350400450500-80-60-40-200三 三 三 Hz三三 三 三 三 dB三三 500三 三 三 三 三 三 三 三 三 三 三 三 三 三 三 三 三 三 三 三圖4 采用樣本混疊增加分段數(shù)得到的功率譜050100150200250300350400450500-80-60-40-200三 三 三 Hz三三 三 三 三 dB三三

36、 三 三 三 三 三 三 500三 三 三 三050100150200250300350400450500-80-60-40-200三 三 三 Hz三三 三 三 三 dB三三 三 welch三 三 三 三 三 三 三 三 三 三 三 三 三 三 三 三圖5 Welch 平均修正周期圖法得到的功率譜估計(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段，而且每

37、段之間有重疊，對(duì)每段分別求FFT后求平均值，得出功率估計(jì)值。功率譜估計(jì)結(jié)果如圖4所示，與圖3相比，顯然估計(jì)值又平滑了許多。另外，采用加窗法也可以降低估計(jì)的協(xié)方差，這也是窗函數(shù)應(yīng)用的一個(gè)方面。即在計(jì)算周期圖法之前，對(duì)數(shù)據(jù)分段并加非矩形窗（如海明窗、漢寧窗或凱瑟窗等），然后再采用分段長(zhǎng)度一半的混疊率，就能夠大大降低估計(jì)方差。這種方法稱(chēng)為Welch平均修正周期圖法，簡(jiǎn)稱(chēng)Welch 法。仿真實(shí)例4是采用海明窗的Welch 法。譜估計(jì)結(jié)果如圖5所示，與圖3，圖4相比，顯然Welch 法得出的功率譜估計(jì)顯然平滑性較好。（三）不同窗函數(shù)的Welch譜估計(jì)的MATLAB仿真仿真實(shí)例仿真實(shí)例 5 5 用Welc

38、h方法對(duì)信號(hào)進(jìn)行功率譜估計(jì)，其中為白)()213. 02sin(2)2 . 02sin(20)(nwnnnx)(nw噪聲。對(duì)信號(hào)的采樣率為=1Hz，滿(mǎn)足取樣定理，n=(0,N-1)。程序代碼見(jiàn) nxsf附錄。當(dāng)N=1000時(shí)，窗函數(shù)采用矩形窗、Hanning窗、hamming窗，MATLAB仿真得到的功率譜估計(jì)圖分別為：01-50050三 三 三 Hz三三 三 三 三 dB三三 三 三0102040三 三 三 Hz三三 三 三 三 dB三hanning三00.40.50

39、.102040三 三 三 Hz三三 三 三 三 dB三hamming三圖6 Welch法加不同的窗函數(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ù)，若離開(kāi)通

40、帶較遠(yuǎn)，則可選用漸近線(xiàn)衰減速度比較快的窗函數(shù)?？傊?，要針對(duì)不同的信號(hào)和不同的處理目的來(lái)選擇合適的窗函數(shù)，這樣才能得到良好的效果。2.5 小結(jié)本章先研究周期圖譜估計(jì)，從均值，方差，分辨率討論周期圖法的性能。估計(jì)的均值得出當(dāng)時(shí)，是無(wú)偏估計(jì)，當(dāng) N 為有限值時(shí)，是有偏估計(jì)。N估計(jì)的方差得出不是的一致估計(jì)。隨著 N 增大，譜估計(jì)的起伏增 大。增大時(shí)，分辨率提高，但會(huì)使的起伏加劇，可見(jiàn)方差與分辨率是N y一對(duì)矛盾。為此，提出了改進(jìn)周期圖法的幾種途徑，即改善窗口形狀、平均和平滑。BT 法中涉及到兩次加窗：一次是用于截取數(shù)據(jù)，叫數(shù)據(jù)窗；一次是用于截取相關(guān)函數(shù)，叫延遲窗，延遲窗常用三角窗。相關(guān)法的兩次加窗將造

41、成泄漏效應(yīng)，影響到譜峰分辨率和方差。平均就是將截取的數(shù)據(jù)段再分成 L 個(gè)小段，分別計(jì)算功率譜后取功率譜的平均，這種方法使估計(jì)的方差減少，但偏差加大，分辨率下降。平滑是用一個(gè)適當(dāng)?shù)拇昂瘮?shù)與算出的功率譜進(jìn)行卷積，使譜線(xià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ù)樣本本身的

42、隨機(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)生假的峰值。而沒(méi)有一個(gè)窗函數(shù)能使譜估計(jì)在方差、偏差和分辨率各方面同時(shí)改善的，使用窗函數(shù)只是一種折中的技巧，不是解決問(wèn)題的根本辦法。不是真實(shí)譜的一致估計(jì)，且 N 增大時(shí)譜曲線(xiàn)起伏加劇。基于上述原因，非參數(shù)譜估計(jì)還有繼續(xù)研究的必要性。3 空間譜

43、估計(jì)的算法與實(shí)現(xiàn)3.1 引言本章將討論使用一個(gè)元無(wú)源陣元組成的陣列對(duì)個(gè)輻射源進(jìn)行定位的問(wèn)mn題。這里，從輻射源來(lái)的能量可能是聲音、電磁波等等，接收傳感器將該能量轉(zhuǎn)換為電信號(hào)。接收陣元可能為電磁天線(xiàn)、水聲傳感器或者地震檢波器。這類(lèi)問(wèn)題在雷達(dá)或聲納系統(tǒng)、通信、天體物理學(xué)、生物醫(yī)學(xué)、地震水下監(jiān)視以及其他許多領(lǐng)域具有廣泛的應(yīng)用。這一問(wèn)題的基本內(nèi)容包括確定能量在空間（空氣、水下或地下）的分布，這里能量高度集中的點(diǎn)就是輻射源。這類(lèi)問(wèn)題稱(chēng)為空間譜估計(jì)。之所以稱(chēng)為譜估計(jì)還因?yàn)檩椛湓炊ㄎ粏?wèn)題和前面討論的譜估計(jì)有著非常緊密的聯(lián)系。事實(shí)上我們將看到，前面所述的方法幾乎都可以用作推導(dǎo)輻射源定位問(wèn)題的解。3.2 全向天

44、線(xiàn)的波束下傾3.2.1 全向天線(xiàn)的波束下傾的基本原理廣播用的發(fā)射機(jī)天線(xiàn)建在天線(xiàn)塔上，希望覆蓋廣大的地域。在覆蓋區(qū)的邊沿，由于距離遠(yuǎn)，地球表面彎曲，電磁波場(chǎng)強(qiáng)衰減很快，因此為了有效地覆蓋既定的區(qū)域，保證區(qū)域內(nèi)場(chǎng)強(qiáng)不低于特定值，通常采用波束下傾的方法，將軸向排列的半波振子天線(xiàn)，通過(guò)調(diào)節(jié)天線(xiàn)的軸向間距、饋電的相位，使得軸向天線(xiàn)陣列的方向圖實(shí)現(xiàn)波束下傾。軸向排列的半振子天線(xiàn)結(jié)構(gòu)，陣元之間的距離是，垂直軸線(xiàn)與電磁波輻d射方向的夾角是，相鄰陣元饋入信號(hào)的相位差是，相鄰陣元發(fā)出電磁波到達(dá)同一地點(diǎn)的程差是，相鄰陣元發(fā)出的電磁波到達(dá)同一地點(diǎn)的相位差參看公式cosd（23） 2cos1120deeeEEnjjj（

45、23）3.2.2 全向天線(xiàn)的波束下傾的仿真實(shí)現(xiàn)根據(jù)以上討論，編寫(xiě)出繪制軸向排列天線(xiàn)陣列的方向圖程序。源程序代碼見(jiàn)附錄。仿真結(jié)果如下圖： 2 4 6 83021060240902701203001503301800圖7 軸向陣列波束下傾仿真圖3.3 天線(xiàn)的波達(dá)方向估計(jì)天線(xiàn)陣列信號(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ù)）。

46、陣列信號(hào)處理時(shí)改善蜂窩和個(gè)人通信服務(wù)系統(tǒng)質(zhì)量、覆蓋范圍和容量的一個(gè)強(qiáng)有力的工具。建立在波達(dá)方向估計(jì)、波束形成基礎(chǔ)上的智能天線(xiàn)的應(yīng)用，抑制了干擾信號(hào)，改善了欲接收信號(hào)的信噪比，降低了數(shù)字通信系統(tǒng)的誤碼率。將接收天線(xiàn)陣列用于公眾網(wǎng)的反向連接（客戶(hù)到基站）時(shí)，多個(gè)接收天線(xiàn)能夠收集更多的信號(hào)能量，若天線(xiàn)在空間足夠分離或極化各異，則多個(gè)天線(xiàn)能夠提供很好的分集接收，并抑制多徑傳輸引起的衰落。這些好處可以擴(kuò)大基站的覆蓋范圍，改善通信質(zhì)量。3.3.1 天線(xiàn)的波達(dá)方向估計(jì)的基本原理全向天線(xiàn)不僅利用率不高，而且對(duì)各種信號(hào)不加區(qū)別地接收，降低了通信質(zhì)量。定點(diǎn)無(wú)線(xiàn)通信采用定向天線(xiàn)，大幅度地改善了通信質(zhì)量。面對(duì)眾多移動(dòng)

47、用戶(hù)的公眾通信網(wǎng)基站和專(zhuān)用移動(dòng)通信網(wǎng)，采用天線(xiàn)指向即波束可變的天線(xiàn)（智能天線(xiàn)），可以使移動(dòng)通信的質(zhì)量得到很大的改善。為使天線(xiàn)的波束指向可控，甚至形狀可控，采用陣列天線(xiàn)是合適的。在距離通信源足夠遠(yuǎn)的空間里，可以將到達(dá)的電磁波視為平面波。對(duì)于等距離直線(xiàn)陣天線(xiàn)，由于調(diào)制在載波上的基帶信號(hào)碼元寬度與波速的乘積遠(yuǎn)大于天線(xiàn)陣的尺寸，因此多個(gè)天線(xiàn)陣元上的信號(hào)的幅度可視為不變，而它們的載波的相位差則取決于其相互位置、尺寸、波長(zhǎng)和到達(dá)方向。無(wú)線(xiàn)接收的無(wú)線(xiàn)信號(hào)中有許多成分，其中我們關(guān)心的是信號(hào)。天線(xiàn)陣列S各個(gè)陣元接收的電磁波信號(hào)因?yàn)殛囋帕形恢玫牟煌瑤?lái)相位差。經(jīng)過(guò)特定參數(shù)的加權(quán)控制器處理后，進(jìn)一步改變了各個(gè)陣元

48、輸出信號(hào)的相位和幅度。處w理的目標(biāo)是使得陣元輸出信號(hào)的和天線(xiàn)輸出中的成分具有最大輸出。用YS信號(hào)作為基準(zhǔn)信號(hào)，反饋控制單元的功能就是將輸出信號(hào)與基準(zhǔn)信號(hào)的SYS差值（即誤差信號(hào)），作為調(diào)節(jié)控制加權(quán)控制器參數(shù)的依據(jù)。反饋控制是w使減小，中的成分加大，也就是說(shuō)，天線(xiàn)陣列接收方向圖指向了信號(hào)YSS的方向。以下通過(guò)兩個(gè)簡(jiǎn)單的例子介紹智能天線(xiàn)中波達(dá)方向估計(jì)和波束形式的原理。等距離直線(xiàn)陣智能天線(xiàn)，我們把天線(xiàn)陣元沿軸排列，從到。若有一平x01M面波以角（入射線(xiàn)與軸的夾角）和角（入射線(xiàn)與軸的夾角）入射到陣列zx上，第號(hào)陣元上產(chǎn)生的信號(hào)為，它與號(hào)陣元的相位差為：KKx0 cossin2dKAK（24）式中，和分

49、別是入射波的波長(zhǎng)和陣元的間距，亦稱(chēng)陣因子。計(jì)入陣因子dKA的影響，第號(hào)陣元的輸出時(shí)，即。為了使天線(xiàn)陣的輸出滿(mǎn)足需要，KKKxAKu在每個(gè)陣元上，用加權(quán)因子進(jìn)行控制。這樣第號(hào)陣元上輸出的信號(hào)為KwK，即。若到達(dá)天線(xiàn)陣的信號(hào)是個(gè)，則天線(xiàn)陣的輸出是個(gè)信KKKxAwKKuwNN號(hào)在個(gè)陣元上的輸出的疊加。將問(wèn)題簡(jiǎn)化為平面的二維問(wèn)題（），Mxy1sin并用解析式表達(dá)如下： nxnxnxnXN,21（25） 1111cos12expcos2exp1dMjdjA（26）其中，為第一個(gè)信號(hào)的入射角。1 NAAAA,21（27） MWWWW,21（28） UWAXWtyHHH（29）其中，。HMuuuU110,求

50、多個(gè)信號(hào)到達(dá)的方向（波達(dá)方向）的方法有多種。下面討論其中兩種方法（Capon法和MUSIC算法）及其仿真驗(yàn)證結(jié)果。3.3.2 求波達(dá)方向的方法1. Capon 法 Capon法亦稱(chēng)最小方差無(wú)畸變響應(yīng) MVDR（Minimum Variance Distortion-Less Response）。天線(xiàn)陣列中的陣元數(shù)決定了陣列方向圖設(shè)計(jì)中的自由度數(shù)。Capon法將陣列中可控的自由度用來(lái)形成期望的波束形狀，達(dá)到對(duì)有用信號(hào)進(jìn)行提升和對(duì)無(wú)用信號(hào)進(jìn)行抑制的目的，并將其優(yōu)化問(wèn)題表達(dá)為 WRWnyEUUHwwminmin2（30）其約束條件為.可以證明上式的解為 10AWH ARAARWUUHUUcap11代

51、入式（29），可以得到相應(yīng)到的功率為 ARAPUUHcap11（31） 2. MUSIC法 MUSIC法亦稱(chēng)多重信號(hào)分類(lèi)（Multiple Signal Classification）。當(dāng)入射信號(hào)數(shù)小于陣元數(shù)時(shí)，的矩陣有與進(jìn)入陣列的信號(hào)數(shù)目相等NMMM UUR的非零特征值及個(gè)為零的特征值。NM （是相應(yīng)的噪聲特征矢量）1NnvV2NvMvnVUUR 因?yàn)榕c的正交性，分母很小，峰值很大，這樣可以得出MUSIC法的空AnV間譜為 AAAAPHHmusic（32）其中，稱(chēng)為噪聲子空間的正交投影估計(jì)。HnnVV3. 仿真 以下是對(duì)七單元線(xiàn)性天線(xiàn)陣在四信號(hào)輸入情況下了仿真仿真實(shí)例仿真實(shí)例 6 6 設(shè)信號(hào)

52、1從方向入射，信號(hào)2從方向入射，信號(hào)3從4/3/6/方向入射，信號(hào)4從方向入射。下面是用兩種方法，分別用兩種坐標(biāo)（直4/3角坐標(biāo)和極坐標(biāo)）求波達(dá)方向估計(jì)。程序代碼見(jiàn)附錄。仿真結(jié)果如下： 4999999999999999800000000000000 9999999999999999600000000000000 149999999999999990000000000000003021060240902701203001503301800圖8 Capon法作出的波達(dá)方向估計(jì)（極坐標(biāo)）020406080100120140160180100105101010151020102510301035圖9

53、MUSIC法作出的波達(dá)方向估計(jì)（直角坐標(biāo)） 500000 1000000 15000003021060240902701203001503301800圖10 MUSIC法作出的波達(dá)方向估計(jì)（極坐標(biāo)）02040608010012014016018010-1510-1010-51001051010圖11 Capon法作出的波達(dá)方向估計(jì)（直角坐標(biāo)）從對(duì)以上四圖的分析可以得到如下結(jié)論：兩種方法估計(jì)得都比較準(zhǔn)確。在處有尖銳的方向圖線(xiàn)。MUSIC法方向圖像幅度更大。4/3 , 3/, 4/, 6/3.4 天線(xiàn)陣的波束形成1 原理我們以等距離圓陣為例來(lái)討論天線(xiàn)陣的波束形成。我們?nèi)A陣的圓心為，在陣列圓上任取

54、一點(diǎn)，把天線(xiàn)陣元順序定為OB到 。若有一平面波以角入射到陣列上，第號(hào)陣元上產(chǎn)生的信號(hào)為01MK，它與到達(dá)陣元中心的波前的相位差是；Kx kKrjAcos2exp與 分別是入射波的波長(zhǎng)與陣列圓的半徑，亦稱(chēng)陣因子。為了使天線(xiàn)陣的rKA輸出滿(mǎn)足需要，在每個(gè)陣元上加上加權(quán)因子控制。這樣第號(hào)陣元上輸出KK的信號(hào)為。若到達(dá)天線(xiàn)陣的信號(hào)是個(gè)，天線(xiàn)陣輸出的是個(gè)信號(hào)的KKKxANN個(gè)陣元上的輸出的疊加。用解析表達(dá)式如下：M )cos(2exp(1kKrjA（33） HMKKKAAAA,21 HNnxnxnxnX110,（34）式中，是陣元以為基準(zhǔn)順時(shí)針畫(huà)出的角度。KKOB NAAAA,21（35） HMwwwW

55、110,（36） UWAXWtyHH（37）為了求得多個(gè)信號(hào)到達(dá)的方向（波達(dá)方向），可以采用上述的Capon、MUSIC兩種方法）。波束形成可以采用下面的方法：當(dāng)有多（）個(gè)信號(hào)輸入時(shí)，其中1個(gè)信號(hào)時(shí)我們關(guān)心的，個(gè)信號(hào)是N1N需要抑制的。方程組（38）描述了上述需求的約束條件（四個(gè)信號(hào)輸入，第一個(gè)信號(hào)是我們關(guān)心的，其余信號(hào)是需要抑制的）。 THUW0 , 0 , 0 , 1（38）根據(jù)信號(hào)波達(dá)方向的知識(shí)及約束條件求解方程組（38），可以)(AXUU得到 110,MwwwW（39）代入式（37）可得到陣列輸出的方向特性。 2. 仿真 下面是四個(gè)輸入信號(hào)八單元圓陣列天線(xiàn)陣的波達(dá)方向估計(jì)和波束形成。仿

56、真實(shí)例仿真實(shí)例 7 設(shè)信號(hào)1從方向入射，信號(hào)2從方向入射，信號(hào)3從3/23/方向入射，信號(hào)4從方向入射。下面是用兩種算法求波達(dá)方向估計(jì)2/36/8 . 5及波束形成。源程序代碼見(jiàn)附錄。仿真結(jié)果如下：050100150200250300350400100105101010151020102510301035圖12 Capon算法波達(dá)方向估計(jì)05010015020025030035040010-1510-1010-51001051010圖13 MUSIC算法波達(dá)方向估計(jì) 0.5 1 1.53021060240902701203001503301800圖14 波束形成3.5 小結(jié) 本章在對(duì)非參數(shù)譜估

57、計(jì)的理解基礎(chǔ)上，討論了其在空間領(lǐng)域進(jìn)行推廣，具體就是了采用Capon算法和MUSIC算法求天線(xiàn)陣的波達(dá)方向，討論了天線(xiàn)陣的波束形成，及其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)

58、典譜估計(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è)向，充分利用了天線(xiàn)各個(gè)陣元從空間電磁場(chǎng)接收到的全部信息，而傳統(tǒng)測(cè)向方法僅僅利用了其中一部分信息，因而具有多方面的優(yōu)越性：1)可以實(shí)現(xiàn)對(duì)同一一信道中存在的多個(gè)信號(hào)同時(shí)測(cè)向；2)超分辨率測(cè)向，即對(duì)處于天線(xiàn)陣固有波束寬度以?xún)?nèi)的兩個(gè)以上方向的來(lái)波同時(shí)測(cè)向；3)陣元位置可以隨意放置，各個(gè)陣元的方向特性也可以是任意的；4)具有優(yōu)良的測(cè)向靈敏度和準(zhǔn)確度；5)

59、可以較好地克服測(cè)向模糊問(wèn)題。參考文獻(xiàn)：參考文獻(xiàn)：1 Petre Stoica,Randolph L.Moses. Introduction to spectral analysis, Prentice Hall,Inc.19972 魏平，唐斌，陸明泉，唐海燕，談覓舟，肖先賜. 譜估計(jì)原理J. 電子科技大學(xué)研究生教材，19993 胡廣書(shū)數(shù)字信號(hào)處理一理論、算法與實(shí)現(xiàn)M北京：清華大學(xué)出版社，20034 王曉峰，王炳和周期圖及其改進(jìn)方法中譜分析率的Matlab分析J武警工程學(xué)院學(xué)報(bào)，2003(6)：75775 姚武川，姚天任經(jīng)典譜估計(jì)方法的Matlab分析J華中理工大學(xué)學(xué)報(bào)，2000，28(4)：4

60、5476 劉剛，呂新華，攸陽(yáng)陣列信號(hào)處理中基于MUSIC算法的空間譜估計(jì)J微計(jì)算機(jī)信息，2006，22(4-3)：302 303附錄附錄附錄附錄 %直接法clear all;clc;%N=256Fs=1;N=256;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=256;Pxx=10*log10(abs(fft(xn,nfft).2)/(N);f=(0:length(Pxx)-1)*Fs/length(Pxx);subplot(211);plot(f,Pxx)xlabel(頻率（H