1、基于MUSIC算法的測向性能仿真 2013 年 1 月 16 日 摘 要 隨著移動通信技術(shù)的飛速發(fā)展，智能天線技術(shù)研究的不斷深入，來波方向(DOA)估計技術(shù)逐漸成為研究的熱點之一，而MUSIC算法是智能天線技術(shù)的典型算法。本文在對MUSIC算法進行分析的基礎(chǔ)上，設(shè)計了MUSIC算法的仿真程序，對不同情況下該算法的性能進行了仿真分析。仿真結(jié)果表明該算法在不同陣列結(jié)構(gòu)、信號入射角度時具有不同的性能。關(guān)鍵詞：智能天線；DOA；MUSIC；陣元目 錄摘 要I引 言1一、MUSIC算法介紹11.1 MUSIC算法的提出11.2波達(dá)方向估計問題中的陣列信號數(shù)學(xué)模型21.3陣列協(xié)方差矩陣的特征分解41.4

2、MUSIC算法的原理及實現(xiàn)61.5 MUSIC算法的實現(xiàn)步驟：8二、MUSIC算法的DOA估計仿真82.1MUSIC算法的基本仿真82.2 MUSIC算法DOA估計與陣元數(shù)的關(guān)系92.3 MUSIC算法DOA估計與陣元間距的關(guān)系102.4 MUSIC算法DOA估計與快拍數(shù)的關(guān)系112.5 MUSIC算法DOA估計與信噪比的關(guān)系122.6 MUSIC算法DOA估計與信號入射角度差的關(guān)系13三、MUSIC算法性能分析小結(jié)15參考文獻15附 錄16附錄一：MUSIC算法的基本仿真源代碼16附錄二：MUSIC算法DOA估計與不同陣元數(shù)關(guān)系仿真源代碼17附錄三： MUSIC算法DOA估計與陣元間距的關(guān)系

3、仿真源代碼18附錄四：MUSIC算法DOA估計與快拍數(shù)的關(guān)系仿真源代碼21附錄五： MUSIC算法DOA估計與信噪比的關(guān)系仿真源代碼22附錄六：MUSIC算法DOA估計與信號入射角度差的關(guān)系仿真源代碼24圖目錄圖1-1 等距線陣與遠(yuǎn)場信號2圖2-1 MUSIC算法的DOA估計譜9圖2-2陣元數(shù)不同時MUSIC算法的DOA估計譜10圖2-3 陣元間距不同時MUSIC算法的DOA估計譜11圖2-4 快拍數(shù)不同時MUSIC算法的DOA估計譜12圖2-5 信噪比不同時MUSIC算法的DOA估計譜13圖2-6 角度間隔不同時MUSIC算法的DOA估計譜14引 言智能天線技術(shù)是當(dāng)前無線移動通信領(lǐng)域頗為關(guān)注

4、和研究的熱點領(lǐng)域之一，可將無線電的信號導(dǎo)向到具體的方向上，產(chǎn)生空間定向波束，使天線主波束對準(zhǔn)用戶信號到達(dá)方向，旁瓣或零陷對準(zhǔn)干擾信號的到達(dá)方向，起到充分高效利用移動用戶信號并刪除或抑制干擾信號的目的。而波束形成的關(guān)鍵是要準(zhǔn)確知道信號的到達(dá)方向，即波達(dá)方向，所以波達(dá)角估計(DOA)是波束形成的基礎(chǔ)。本文著重分析了用于DOA估計的典型算法-MUSIC(Multiple Signal Classification)算法，然后對不同的條件下MUSIC算法的性能進行了Matlab的仿真和分析。一、MUSIC算法介紹1.1 MUSIC算法的提出多重信號分類（MUSIC）算法是Schmidt等人在1979年

5、提出的。這一算法的提出開創(chuàng)了空間譜估計算法研究的新時代，促進了特征結(jié)構(gòu)類算法的興起和發(fā)展，該算法已成為空間譜估計理論體系中的標(biāo)志性算法。此算法提出之前的有關(guān)算法都是針對陣列接收數(shù)據(jù)協(xié)方差矩陣進行直接處理，而MUSIC算法的基本思想則是對任意陣列輸出數(shù)據(jù)的協(xié)方差矩陣進行特征分解，從而得到與信號分類相對應(yīng)的信號子空間和與信號分量相正交的噪聲子空間，然后利用這兩個子空間的正交性構(gòu)造空間譜函數(shù)，通過譜峰搜索，檢測信號的DOA。正是由于MUSIC算法在特定的條件下具有很高的分辨力、估計精度及穩(wěn)定性，從而吸引了大量的學(xué)者對其進行深入的研究和分析。總的來說，它用于陣列的波達(dá)方向估計有以下一些突出的優(yōu)點：(1

6、)多信號同時測向能力(2)高精度測向(3)對天線波束內(nèi)的信號的高分辨測向(4)可適用于短數(shù)據(jù)情況(5)采用高速處理技術(shù)后可實現(xiàn)實時處理1.2波達(dá)方向估計問題中的陣列信號數(shù)學(xué)模型 為了分析推導(dǎo)的方便，現(xiàn)將波達(dá)方向估計問題中的數(shù)學(xué)模型作理想狀態(tài)的假設(shè)如下：(1)各待測信號源具有相同的極化、且互不相關(guān)的。一般考慮信號源為窄帶的，且各信號源具有相同的中心頻率。待測信號源的個數(shù)為D。(2)天線陣列是由M(MD)個陣元組成的等間距直線陣，各陣元特性相同，各向同性，陣元間隔為d，并且陣元間隔不大于最高頻率信號半波長。(3)天線陣列處于各信號源的遠(yuǎn)場中，即天線陣列接收從各信號源傳來的信號為平面波。(4)各陣元

7、上有互不相關(guān)，與各待測信號也不相關(guān)，方差為的零均值高斯白噪聲。(5)各接收支路具有完全相同的特性。1d2 3 M圖1-1 等距線陣與遠(yuǎn)場信號設(shè)由第k（k=1，2，D）個信號源輻射到天線陣列的波前信號為，前面已假設(shè)為窄帶信號，則可以表示為以下形式： （1.1） 式中是的復(fù)包絡(luò)，是信號的角頻率。前面已經(jīng)假設(shè)D個信號具有相同的中心頻率，所以有： （1.2） 式中c是電磁波波速，是公用的信號波長。設(shè)電磁波通過天線陣列尺寸所需的時間為，則根據(jù)窄帶假設(shè)，有如下近似： （1.3） 故延遲后的波前信號為： （1.4）所以，若以第一個陣元為參考點，則t時刻等間距直線陣中的第m(m=1，2，M)個陣元對第k個信號

8、源的感應(yīng)信號為： （1.5） 其中，為第m個陣元對第k個信號源的影響，前面以假設(shè)各陣元無方向性，所以可取。為第k個信號源的方位角，表示由第m個陣元與第1個陣元間的波程差所引起的信號相位差。計及測量噪聲和所有信號源來波，第m個陣元的輸出信號為： （1.6） 其中是測量噪聲，所有標(biāo)號為m表示該量屬于第m個陣元，所有標(biāo)號為k表示該量屬于第k個信號源。設(shè) （1.7） 為第m個陣元對第k個信號源的響應(yīng)函數(shù)。則第m個陣元的輸出信號為： （1.8）其中是第k個信號源在陣元上的信號強度。運用矩陣的定義，可以得到更為簡潔的表達(dá)式： X=AS+N （1.9）式中 （1.10） （1.11） = （1.12） （1

9、.13） （1.14） 對進行N點采樣，要處理的問題就變成了通過輸出信號的采樣估計出信號源的波達(dá)方向角。由此，可以很自然的將陣列信號看作是噪聲干擾的若干空間諧波的疊加，從而將波達(dá)方向估計問題與譜估計聯(lián)系起來。1.3陣列協(xié)方差矩陣的特征分解對陣列輸出x作相關(guān)處理，得到其協(xié)方差矩陣： （1.15） 其中，H表示矩陣共軛轉(zhuǎn)置。前面已假設(shè)信號與噪聲互不相關(guān)、且噪聲為零均值白噪聲，因此將式(1.9)代入式(1.15)，可以得到： = = （1.16）式中 （1.17） 稱為信號的相關(guān)矩陣。 （1.18） 是噪聲的相關(guān)矩陣，是噪聲功率，I是M*M階的單位矩陣。實際應(yīng)用中，通常無法直接得到，能使用的只有樣本

10、的協(xié)方差矩陣： （1.19）是的最大似然估計，當(dāng)采樣數(shù)時，它們是一致的，但實際情況中將由于樣本數(shù)有限而造成誤差。根據(jù)矩陣特征分解的理論，可以對陣列協(xié)方差矩陣進行特征分解。首先考慮理想情況，即無噪聲的情況： （1.20） 對于均勻線陣，矩陣A是由式(1.12)所定義的范德蒙德矩陣，只要滿足： （1.21） 則，它的各列相互獨立，這樣，若為非奇異矩陣（，各信號源兩兩不相干），且MD，則有： （1.22） 由于，所以有： （1.23） 即是Hermite矩陣，它的特征值都是實數(shù)。又由于是正定的，因此矩陣是半正定的，它有D個正特征值和M-D個零特征值。再考慮有噪聲存在的情況 （1.24） 由于0，為滿

11、秩陣，所以有M個正實特征值，分別對應(yīng)于M個特征向量。又由于是Hermite矩陣，所以各特征向量是相互正交的，即： （1.25） 與信號有關(guān)的特征值只有D個，分別等于矩陣的各特征值與之和，其余的M-D個特征值為，也就是說，是R的最小特征值，它是M-D維的。對應(yīng)的特征向量，i=1,2，M中，也有D個是與信號有關(guān)的，另外M-D個是與噪聲有關(guān)的，在下一節(jié)里，將利用以上這些特征分解的性質(zhì)求出信號源的波達(dá)方向。1.4 MUSIC算法的原理及實現(xiàn)通過對陣列協(xié)方差矩陣的特征分解，可以得到如下結(jié)論：將矩陣的特征值進行從小到大的排序，即 （1.26） 其中D個較大的特征值對應(yīng)于信號，M-D個較小的特征值對應(yīng)于噪聲

12、。矩陣的屬于這些特征值的特征向量也分別對應(yīng)于信號和噪聲，因此，可以把的特征值(特征向量)劃分為信號特征值(特征向量)與噪聲特征值(特征向量)。設(shè)是矩陣的第i個特征值，是與個相對應(yīng)的特征向量，則有： （1.27） 再設(shè)是的最小特征值 i=D+1，D+2，M （1.28） 將 （1.29）代入上式，可得： （1.30） 將上式右邊展開與左邊比較，可得： （1.31） 因是D*D維的滿秩矩陣，存在；而同樣存在，則上式兩邊同乘以后變成： （1.32） 于是有 i=D+1，D+2，M （1.33） 上式表明：噪聲特征值所對應(yīng)的特征向量(稱噪聲特征向量)，與矩陣A的列向量正交，而A的各列是與信號源的方向相

13、對應(yīng)的。這就是利用噪聲特征向量求解信號源方向的出發(fā)點。用各噪聲特征向量為列，構(gòu)造一個噪聲矩陣： （1.34） 定義空間譜： = （1.35） 該式中分母是信號向量和噪聲矩陣的內(nèi)積，當(dāng)和的各列正交時，該分母為零，但由于噪聲的存在，它實際上為一最小值，因此有一尖峰。由該式，使變化，通過尋找波峰來估計到達(dá)角。1.5 MUSIC算法的實現(xiàn)步驟：(1)根據(jù)N個接收信號矢量得到下面協(xié)方差矩陣的估計值： （1.36） 對上面得到的協(xié)方差矩陣進行特征值分解 （1.37） (2)按特征值的大小順序，把與信號個數(shù)D相等的特征值和對應(yīng)的特征向量看作信號部分空間，把剩下的M-D個特征值和特征向量看作噪聲部分空間。得到

14、噪聲矩陣： i=D+1，D+2，M （1.38） （1.39） (3)使變化，按照式 （1.40）來計算譜函數(shù)，通過尋求峰值來得到波達(dá)方向的估計值。二、MUSIC算法的DOA估計仿真2.1MUSIC算法的基本仿真模擬2個獨立窄帶信號分別以20，70的方向入射到均勻線陣上，信號間互不相關(guān)，與噪聲相互獨立，噪聲為理想高斯白噪聲，陣元間距為入射信號波長的12，快拍數(shù)為200，信噪比為20dB，陣元數(shù)分別為10，50，100。其仿真結(jié)果如圖2-1所示：圖2-1 MUSIC算法的DOA估計譜由圖2-1可以看出在符合假設(shè)的前提下，采用MUSIC算法能構(gòu)造出針狀的譜峰，可以很好的估計出入射信號的個數(shù)和方向，

15、能有效的估計出獨立信號源的DOA,并且在模型準(zhǔn)確的前提下，對DOA的估計可以達(dá)到任意精度，克服了傳統(tǒng)測向定位方法精度低的缺點 ,可以有效解決密集信號環(huán)境中多個輻射源的高分辨率、高精度測向定位問題。可以看出超分辨率的 MUSIC算法具有測向準(zhǔn)確度、靈敏度高的特點且具有潛在分辨多信號的能力，具有較好的性能和較高的效率，能提供高分辨率及漸近無偏的到達(dá)角估計，這對實際中的應(yīng)用具有十分重要的意義。2.2 MUSIC算法DOA估計與陣元數(shù)的關(guān)系模擬2個獨立窄帶信號分別以20，60的方向入射到均勻線陣上，信號間互不相關(guān)，與噪聲相互獨立，噪聲為理想高斯白噪聲，陣元間距為入射信號波長的12，快拍數(shù)為200，信噪

16、比為20dB，陣元數(shù)分別為5，10，20。其仿真結(jié)果如圖2-2所示：圖2-2陣元數(shù)不同時MUSIC算法的DOA估計譜 由圖2-2可以看出，其他條件不變的情況下，隨著陣元數(shù)的增加，DOA估計譜的波束寬度變窄，陣列的指向性變好，也就是說陣列分辨空間信號的能力增強。由此可以看出，要得到更加精確的DOA估計譜，可以增加陣元數(shù)量，但陣元數(shù)量越多，需要處理的數(shù)據(jù)越多，運算量越大，運行速度越慢。由上圖可以看出陣元數(shù)大到一定數(shù)量時，波形變化不會很明顯。因此，在實際應(yīng)用中，可根據(jù)具體條件適當(dāng)選取陣元數(shù)量，在確保估計譜準(zhǔn)確的前提下，盡量減少資源浪費，加快運行的速度，提高工作效率。2.3 MUSIC算法DOA估計與

17、陣元間距的關(guān)系 模擬2個獨立窄帶信號分別以20，60的方向入射到均勻線陣上，信號間互不相關(guān)，與噪聲相互獨立，噪聲為理想高斯白噪聲，陣元數(shù)為10，快拍數(shù)為200，信噪比為20dB，陣元間距分別為/6、/2、。其仿真結(jié)果如圖2-3所示：圖2-3 陣元間距不同時MUSIC算法的DOA估計譜 由圖2-3可以看出，在其他條件不變的前提下，當(dāng)陣元間距不大于半波長時，隨著陣元間距的增加，DOA估計譜的波束寬度變窄，陣列的指向性變好，也就是說MUSIC算法的分辨力隨著陣元間距的加大相應(yīng)提高，但當(dāng)陣元間距大于半波長時，估計譜除了信號源方向外在其他方向出現(xiàn)了虛假譜峰，也就失去了估計的準(zhǔn)確性。可見，在實際應(yīng)用中，要

18、十分注意陣元間的距離，可以適當(dāng)增加陣元間距但絕不能超過半波長，這一點非常重要，最好是將陣元間距設(shè)為半波長。2.4 MUSIC算法DOA估計與快拍數(shù)的關(guān)系模擬2個獨立窄帶信號分別以20，60的方向入射到均勻線陣上信號間互不相關(guān)，與噪聲相互獨立，噪聲為理想高斯白噪聲，陣元數(shù)為10，陣元間距為入射信號波長的12，信噪比為20dB，快拍數(shù)分別為5，50，200。其仿真結(jié)果如圖2-4所示：圖2-4 快拍數(shù)不同時MUSIC算法的DOA估計譜 由圖2-4可以看出，在其他條件不變的情況下，隨著快拍數(shù)的增加，DOA估計譜的波束寬度變窄，陣列的指向性變好，陣列分辨空間信號的能力增強，MUSIC算法的估計精度增加。

19、由此可見，可通過增加采樣快拍數(shù)來增加DOA估計的精確度，但是采樣快拍數(shù)越多，需要處理的數(shù)據(jù)就越多，MUSIC算法的運算量就越大，速度就越慢，所以在實際應(yīng)用中要合理的選取采樣快拍數(shù)，在確定DOA估計譜準(zhǔn)確的前提下，盡量減少運算量，加快工作速度，節(jié)省人力物力，節(jié)約資源。2.5 MUSIC算法DOA估計與信噪比的關(guān)系模擬2個獨立窄帶信號分別以20，60的方向入射到均勻線陣上，信號間互不相關(guān)，與噪聲相互獨立，噪聲為理想高斯白噪聲，陣元數(shù)為10，陣元間距為入射信號波長的12，快拍數(shù)為200，信噪比分別為-10dB，0dB，20dB。其仿真結(jié)果如圖2-5所示：圖2-5 信噪比不同時MUSIC算法的DOA估

20、計譜由圖2-5可以看出，在其他條件不變的情況下，隨著信噪比的增加，DOA估計譜的波束寬度變窄，陣列的指向性變好，MUSIC算法的分辨力增加，信噪比的高低直接影響著超分辨方位估計算法的性能。在低信噪比時，MUSIC算法的性能會急劇下降，因而提高算法在低信噪比條件下的估計性能是超分辨DOA算法的研究重點。有學(xué)者提出了一種基于多級維納濾波器(MSWF：Multistage Weiner Filtering)的信號波達(dá)方向(DOA)估計算法，該算法在信號可能入射方向用MSWF 估計信號子空間，并在 MSWF 分解后互相關(guān)函數(shù)最小以及信號子空間估值與噪聲子空間正交時判定估計有效，進而構(gòu)造空間譜來實現(xiàn)信號

21、 DOA 估計。已證明在低信噪比條件下，該算法比子空間類算法有更好的分辨率和誤差性能。低信噪比條件下對DOA的精確估計還有很大的發(fā)展改進空間，有待進一步研究。2.6 MUSIC算法DOA估計與信號入射角度差的關(guān)系模擬2個獨立窄帶信號入射到均勻線陣上，信號間互不相關(guān)，與噪聲相互獨立，噪聲為理想高斯白噪聲，陣列的陣元數(shù)為10，快拍數(shù)為200。陣元間距為入射信號波長的12，信噪比為20dB，信號入射角度差分別為5，10，40。其仿真結(jié)果如圖2-6所示：圖2-6 角度間隔不同時MUSIC算法的DOA估計譜圖2-6說明在其他條件不變的情況下，隨著信號入射角度差的增加，DOA估計譜的波束寬度變窄，陣列的指

22、向性變好，MUSIC算法的分辨力增加。當(dāng)信號來波方向間隔角度很小時，不能準(zhǔn)確估計信號源數(shù)。通常的陣列信號源數(shù)估計方法, 都是在信號來波方向角度差較大情況下進行的估計，當(dāng)信號來波方向的角度差比較小時, 這些方法估計都要失效。已有學(xué)者提出了平方根修正的 Gerschgorin半徑估計方法，對信號來波方向角度差小時,也能很好的估計信號源數(shù)。目前提出的一些信號源數(shù)估計方法，大都存在一定的應(yīng)用條件，這在一定程度上限制了DOA算法的實際運用，因此，研究符合實際應(yīng)用環(huán)境的實時、穩(wěn)健的信號源數(shù)與DOA的聯(lián)合估計等仍具有十分重要的現(xiàn)實意義。三、MUSIC算法性能分析小結(jié)通過上述幾組仿真可以看出超分辨率的 MUS

23、IC算法具有較好的性能和較高的效率，能提供高分辨率及漸近無偏的到達(dá)角估計。而且陣元數(shù)越多，快拍數(shù)越多，信噪比越高，信號入射角度差越大 MUSIC算法的分辨率越高，當(dāng)陣元間距不大于載波半波長時， MUSIC算法的分辨力隨著陣元間距的加大相應(yīng)提高，但當(dāng)陣元間距大于/2時，空間譜除了信號源方向外在其他方向出現(xiàn)虛假譜峰。在小信噪比和小的角度間隔時，MUSIC算法的估計性能下降，已有學(xué)者提出了一些的改進算法，但這些問題仍是當(dāng)前研究的熱點。此外，有學(xué)者提出了很多對MUSIC算法的改進。可見，MUSIC算法還有很大的發(fā)展空間，值得我們進一步研究。參考文獻1Eli張賢達(dá)，保錚通信信號處理M北京：國防工業(yè)出版社

24、，20002張賢達(dá)現(xiàn)代信號處理M北京：清華大學(xué)出版社，20023鄭洪MUSIC算法與波達(dá)方向估計研究D成都：四川大學(xué)，20054張娟智能天線中DOA估計算法的研究D哈爾濱：哈爾濱工業(yè)大學(xué)，20065張玲華，鄭寶玉隨機信號處理M北京：清華大學(xué)出版社，20036Simon Haykin自適應(yīng)濾波器原理M北京：電子工業(yè)出版社，20037樓順天，李博菡. 基于MATLAB的系統(tǒng)分析與設(shè)計信號處理M. 西安: 西安電子科技大學(xué)出版社, 2000.8王正林，劉明.精通MATLAB7M.北京：電子工業(yè)出版社，2006.9王進，趙擁軍，王志剛.低信噪比條件下的高分辨DOA估計算法J.計算機工程，2005.Vo

25、l.35.No.3.附 錄附錄一：MUSIC算法的基本仿真源代碼clc clear all format long %將數(shù)據(jù)顯示為長整型科學(xué)計數(shù)N=200;%快拍數(shù)doa=20 60/180*pi; %信號到達(dá)角w=pi/4 pi/3;%信號頻率M=10;%陣元數(shù)P=length(w); %信號個數(shù)lambda=150;%波長d=lambda/2;%陣元間距 snr=20;%信噪比B=zeros(P,M); %創(chuàng)建一個P行M列的0矩陣for k=1:P B(k,:)=exp(-j*2*pi*d*sin(doa(k)/lambda*0:M-1); %矩陣賦值end B=B; xx=2*exp(j

26、*(w*1:N); %仿真信號x=B*xx; x=x+awgn(x,snr);%加入高斯白噪聲R=x*x; %數(shù)據(jù)協(xié)方差矩陣 U,V=eig(R); %求R的特征值和特征向量UU=U(:,1:M-P); %估計噪聲子空間theta=-90:0.5:90; %譜峰搜索for ii=1:length(theta) AA=zeros(1,length(M); for jj=0:M-1 AA(1+jj)=exp(-j*2*jj*pi*d*sin(theta(ii)/180*pi)/lambda);end WW=AA*UU*UU*AA; Pmusic(ii)=abs(1/ WW); end Pmusic

27、=10*log10(Pmusic/max(Pmusic); %空間譜函數(shù)plot(theta,Pmusic,-k) xlabel(角度 theta/degree) ylabel(譜函數(shù)P(theta) /dB) title(MUSIC算法的DOA估計譜)grid on 附錄二：MUSIC算法DOA估計與不同陣元數(shù)關(guān)系仿真源代碼clc clear all format long %將數(shù)據(jù)顯示為長整型科學(xué)計數(shù)N=200; %快拍數(shù)doa=20 60/180*pi; %信號到達(dá)角w=pi/4 pi/3; %信號頻率M1=5; %陣元數(shù)M2=10;M3=20;P=length(w); %信號個數(shù)lam

28、bda=150; %波長d=lambda/2; %陣元間距snr=20; %信噪比B1=zeros(P,M1); B2=zeros(P,M2);B3=zeros(P,M3);for k=1:P B1(k,:)=exp(-j*2*pi*d*sin(doa(k)/lambda*0:M1-1); %矩陣賦值B2(k,:)=exp(-j*2*pi*d*sin(doa(k)/lambda*0:M2-1); B3(k,:)=exp(-j*2*pi*d*sin(doa(k)/lambda*0:M3-1); end B1=B1; B2=B2; B3=B3; xx=2*exp(j*(w*1:N); %仿真信號x

29、1=B1*xx; x2=B2*xx; x3=B3*xx; x1=x1+awgn(x1,snr); %加入高斯白噪聲x2=x2+awgn(x2,snr);x3=x3+awgn(x3,snr);R1=x1*x1; %數(shù)據(jù)協(xié)方差矩陣R2=x2*x2;R3=x3*x3;U1,V1=eig(R1); %求R的特征值和特征向量U2,V2=eig(R2); U3,V3=eig(R3); UU1=U1(:,1:M1-P); ; %估計噪聲子空間UU2=U2(:,1:M2-P); UU3=U3(:,1:M3-P); theta=-90:0.5:90; %譜峰搜索for ii=1:length(theta) AA

30、1=zeros(1,length(M1); for jj=0:M1-1 AA1(1+jj)=exp(-j*2*jj*pi*d*sin(theta(ii)/180*pi)/lambda);end WW1=AA1*UU1*UU1*AA1; Pmusic1(ii)=abs(1/ WW1); end for ii=1:length(theta) AA2=zeros(1,length(M2); for jj=0:M2-1 AA2(1+jj)=exp(-j*2*jj*pi*d*sin(theta(ii)/180*pi)/lambda);end WW2=AA2*UU2*UU2*AA2; Pmusic2(ii

31、)=abs(1/ WW2); end for ii=1:length(theta) AA3=zeros(1,length(M3); for jj=0:M3-1 AA3(1+jj)=exp(-j*2*jj*pi*d*sin(theta(ii)/180*pi)/lambda);end WW3=AA3*UU3*UU3*AA3; Pmusic3(ii)=abs(1/ WW3); end Pmusic1=10*log10(Pmusic1/max(Pmusic1); %空間譜函數(shù)Pmusic2=10*log10(Pmusic2/max(Pmusic2);Pmusic3=10*log10(Pmusic3/m

32、ax(Pmusic3);plot(theta,Pmusic1,-k,LineWidth,2.0)hold onplot(theta,Pmusic2,r,LineWidth,2.0) hold onplot(theta,Pmusic3,-b,LineWidth,2.0) hold offxlabel(角度 theta/degree) ylabel(譜函數(shù)P(theta) /dB) title(MUSIC算法的DOA估計譜)grid on 附錄三： MUSIC算法DOA估計與陣元間距的關(guān)系仿真源代碼clc clear all format long %將數(shù)據(jù)顯示為長整型科學(xué)計數(shù)N=200;%快拍數(shù)

33、doa=20 60/180*pi; %信號到達(dá)角w=pi/4 pi/3;%信號頻率M=10;%陣元數(shù)P=length(w); %信號個數(shù)lambda=150;%波長d=lambda/6;%陣元間距 snr=20;%信噪比B=zeros(P,M); %創(chuàng)建一個P行M列的0矩陣for k=1:P B(k,:)=exp(-j*2*pi*d*sin(doa(k)/lambda*0:M-1); %矩陣賦值end B=B; xx=2*exp(j*(w*1:N); %仿真信號x=B*xx; x=x+awgn(x,snr);%加入高斯白噪聲R=x*x; %數(shù)據(jù)協(xié)方差矩陣 U,V=eig(R); %求R的特征值

34、和特征向量UU=U(:,1:M-P); %估計噪聲子空間theta=-90:0.5:90; %譜峰搜索for ii=1:length(theta) AA=zeros(1,length(M); for jj=0:M-1 AA(1+jj)=exp(-j*2*jj*pi*d*sin(theta(ii)/180*pi)/lambda);end WW=AA*UU*UU*AA; Pmusic(ii)=abs(1/ WW); end Pmusic=10*log10(Pmusic/max(Pmusic); %空間譜函數(shù)plot(theta,Pmusic ,-k,linewidth,2.0)hold onclc

35、 clear all format long %將數(shù)據(jù)顯示為長整型科學(xué)計數(shù)N=200;%快拍數(shù)doa=20 60/180*pi; %信號到達(dá)角w=pi/4 pi/3;%信號頻率M=10;%陣元數(shù)P=length(w); %信號個數(shù)lambda=150;%波長d=lambda/2;%陣元間距 snr=20;%信噪比B=zeros(P,M); %創(chuàng)建一個P行M列的0矩陣for k=1:P B(k,:)=exp(-j*2*pi*d*sin(doa(k)/lambda*0:M-1); %矩陣賦值end B=B; xx=2*exp(j*(w*1:N); %仿真信號x=B*xx; x=x+awgn(x,s

36、nr);%加入高斯白噪聲R=x*x; %數(shù)據(jù)協(xié)方差矩陣 U,V=eig(R); %求R的特征值和特征向量UU=U(:,1:M-P); %估計噪聲子空間theta=-90:0.5:90; %譜峰搜索for ii=1:length(theta) AA=zeros(1,length(M); for jj=0:M-1 AA(1+jj)=exp(-j*2*jj*pi*d*sin(theta(ii)/180*pi)/lambda);end WW=AA*UU*UU*AA; Pmusic(ii)=abs(1/ WW); end Pmusic=10*log10(Pmusic/max(Pmusic); %空間譜函

37、數(shù)plot(theta,Pmusic,r,linewidth,2.0) hold onclc clear all format long %將數(shù)據(jù)顯示為長整型科學(xué)計數(shù)N=200;%快拍數(shù)doa=20 60/180*pi; %信號到達(dá)角w=pi/4 pi/3;%信號頻率M=10;%陣元數(shù)P=length(w); %信號個數(shù)lambda=150;%波長d=lambda;%陣元間距 snr=20;%信噪比B=zeros(P,M); %創(chuàng)建一個P行M列的0矩陣for k=1:P B(k,:)=exp(-j*2*pi*d*sin(doa(k)/lambda*0:M-1); %矩陣賦值end B=B; x

38、x=2*exp(j*(w*1:N); %仿真信號x=B*xx; x=x+awgn(x,snr);%加入高斯白噪聲R=x*x; %數(shù)據(jù)協(xié)方差矩陣 U,V=eig(R); %求R的特征值和特征向量UU=U(:,1:M-P); %估計噪聲子空間theta=-90:0.5:90; %譜峰搜索for ii=1:length(theta) AA=zeros(1,length(M); for jj=0:M-1 AA(1+jj)=exp(-j*2*jj*pi*d*sin(theta(ii)/180*pi)/lambda);end WW=AA*UU*UU*AA; Pmusic(ii)=abs(1/ WW); e

39、nd Pmusic=10*log10(Pmusic/max(Pmusic); %空間譜函數(shù)plot(theta,Pmusic,-.b,linewidth,2.0) hold offxlabel(角度 theta/degree) ylabel(譜函數(shù)P(theta) /dB) title(MUSIC算法的DOA估計譜)grid on 附錄四：MUSIC算法DOA估計與快拍數(shù)的關(guān)系仿真源代碼clc clear all format long %將數(shù)據(jù)顯示為長整型科學(xué)計數(shù)N1=5;%快拍數(shù)N2=50;N3=200;doa=20 60/180*pi; %信號到達(dá)角w=pi/4 pi/3;%信號頻率M=

40、10;%陣元數(shù)P=length(w); %信號個數(shù)lambda=150;%波長d=lambda/2;%陣元間距 snr=20;%信噪比B=zeros(P,M); %創(chuàng)建一個P行M列的0矩陣for k=1:P B(k,:)=exp(-j*2*pi*d*sin(doa(k)/lambda*0:M-1); %矩陣賦值end B=B; xx1=2*exp(j*(w*1:N1); %仿真信號xx2=2*exp(j*(w*1:N2);xx3=2*exp(j*(w*1:N3);x1=B*xx1; x2=B*xx2;x3=B*xx3;x1=x1+awgn(x1,snr);%加入高斯白噪聲x2=x2+awgn(

41、x2,snr);x3=x3+awgn(x3,snr);R1=x1*x1; %數(shù)據(jù)協(xié)方差矩陣 R2=x2*x2;R3=x3*x3;U1,V1=eig(R1); U2,V2=eig(R2); U3,V3=eig(R3); UU1=U1(:,1:M-P); UU2=U2(:,1:M-P); UU3=U3(:,1:M-P); theta=-90:0.5:90; %譜峰搜索for ii=1:length(theta) AA=zeros(1,length(M);for jj=0:M-1 AA(1+jj)=exp(-j*2*jj*pi*d*sin(theta(ii)/180*pi)/lambda);end

42、WW1=AA*UU1*UU1*AA; WW2=AA*UU2*UU2*AA; WW3=AA*UU3*UU3*AA; Pmusic1(ii)=abs(1/ WW1);Pmusic2(ii)=abs(1/ WW2); Pmusic3(ii)=abs(1/ WW3); end Pmusic1=10*log10(Pmusic1/max(Pmusic1);Pmusic2=10*log10(Pmusic2/max(Pmusic2);Pmusic3=10*log10(Pmusic3/max(Pmusic3);plot(theta,Pmusic1,-k,linewidth,2.0)hold onplot(the

43、ta,Pmusic2,-r,linewidth,2.0)hold onplot(theta,Pmusic3,-.b,linewidth,2.0)hold offxlabel(角度 theta/degree) ylabel(譜函數(shù)P(theta) /dB) title(MUSIC算法的DOA估計譜)grid on 附錄五： MUSIC算法DOA估計與信噪比的關(guān)系仿真源代碼clc clear all format long %將數(shù)據(jù)顯示為長整型科學(xué)計數(shù)N=200;%快拍數(shù)doa=20 60/180*pi; %信號到達(dá)角w=pi/4 pi/3;%信號頻率M=10;%陣元數(shù)P=length(w); %

44、信號個數(shù)lambda=150;%波長d=lambda/2;%陣元間距 snr1=-10;%信噪比snr2=0;snr3=20;B=zeros(P,M); %創(chuàng)建一個P行M列的0矩陣for k=1:P B(k,:)=exp(-j*2*pi*d*sin(doa(k)/lambda*0:M-1); %矩陣賦值end B=B; xx=2*exp(j*(w*1:N); %仿真信號x=B*xx; x1=x+awgn(x,snr1);%加入高斯白噪聲x2=x+awgn(x,snr2);x3=x+awgn(x,snr3);R1=x1*x1; %數(shù)據(jù)協(xié)方差矩陣 R2=x2*x2;R3=x3*x3;U1,V1=e

45、ig(R1); U2,V2=eig(R2); U3,V3=eig(R3); UU1=U1(:,1:M-P); UU2=U2(:,1:M-P); UU3=U3(:,1:M-P); theta=-90:0.5:90; %譜峰搜索for ii=1:length(theta) AA=zeros(1,length(M); for jj=0:M-1 AA(1+jj)=exp(-j*2*jj*pi*d*sin(theta(ii)/180*pi)/lambda);end WW1=AA*UU1*UU1*AA; WW2=AA*UU2*UU2*AA; WW3=AA*UU3*UU3*AA; Pmusic1(ii)=abs(1/ WW1);Pmusic2(ii)=abs(1/ WW2); Pmusic3(ii)=abs(1/ WW3);end Pmusic1=10*log10(Pmusic1/max(Pmusic1);Pmusic2=10*log10(Pmusic2/max(Pmusic2);Pmusic3=10*log10(P