一維DOA估計(jì)算法對(duì)比實(shí)驗(yàn)

Projectreportaboutone-dimensionalDOAestimation題目：考慮一個(gè)10陣元數(shù)的均勻線陣，現(xiàn)有三個(gè)信源入射，它們的波達(dá)方向(DOA)分別是10°,30°和50°,請(qǐng)用MUSIC算法，ESPRIT算法，Capon算法以及傳播算子(PM)來(lái)估計(jì)這些信源的波達(dá)方向。1.傳統(tǒng)DOA估計(jì)算法1.1信號(hào)接收模型考慮N個(gè)不同DOA0,n=1,2,…,N的窄帶遠(yuǎn)場(chǎng)信號(hào)s(t),在離散時(shí)間k入射n n有M個(gè)傳感器的均勻線陣(ULA)時(shí)。在k時(shí)刻收集L次快拍，則陣列接收到的信號(hào)可以表示為：其中N=[n,k1Z其中N=[n,k1Z=[z,z,k12A（0）=「a（0）,MxN 1n2,Z=A(0)S+N (kkk,n]tgCmxl表示與信號(hào)不相干的加性高斯白噪聲,M]gCmxL為離散時(shí)間k時(shí)1）zMa（0）,???, a（0）2 N」的陣列觀測(cè),矩陣s]tgCnxL表示信源矩陣。其中N2）a（0）=2）n九為波長(zhǎng)，九為波長(zhǎng)，d是陣元間距。傳播算子方法1.2MUSIC算法聲源定位方法常見的有Capon波束形成算法,MUSIC算法傳播算子方法以及ESPRIT方法。這些方法中,MUSIC算法是通過譜峰搜索來(lái)確定聲源波達(dá)方向,假設(shè)噪聲是獨(dú)立同分布的零均值的高斯過程,且信號(hào)和噪聲不相干。R為k陣列協(xié)方差矩陣，可以表示為：R=同5Zh}=A(0)RAh(0)+Q=U刀Uh+U￡Uh(3)k kk s kSSSNNN式中Q=也「NNH1=Q2IgCmxM表示噪聲協(xié)方差矩陣，E是期望操作，上標(biāo)Hk Lkk」 M表示共軛轉(zhuǎn)置，信號(hào)協(xié)方差矩陣R=1e{sSt}gCnxn。U表示(3)式特征分解后s kk SN個(gè)大特征值對(duì)應(yīng)的特征向量張成的空間，稱作信號(hào)子空間；U表示(3)式特征N分解后(M-N)個(gè)小特征值對(duì)應(yīng)的特征向量張成的空間，稱作噪聲子空間，且與

信號(hào)子空間正交。在搜集L次快拍下，R可以近似為：k1R沁一ZZh (4)kLkk在空間譜估計(jì)算法中，MUSIC算法通過譜估計(jì)式的峰值搜索，來(lái)確定入射的方向信息，這與粒子濾波中的似然函數(shù)具有相似的效果，MUSIC算法空間譜函數(shù)為P @)=MUSICP @)=MUSIC1a(0)hUUha(0)NN5)式中U是噪聲子空間，經(jīng)過式(4)特征分解可以得到。則DOA角度估計(jì)可以通N過譜峰搜索使得式(5)輸出值最大化的角度信息得到，即：二arg6)max P(0)二arg6)0g[-k/2,兀/2]k|表示復(fù)數(shù)的幅度值?，F(xiàn)考慮一個(gè)10陣元數(shù)的均勻線陣，有三個(gè)信源入射,它們的波達(dá)方向(DOA)分別是10°,30°和50°它們的波達(dá)方向(DOA)分別是10°,30°和50°，快拍數(shù)L=100。在信噪比相對(duì)較高的噪聲環(huán)境中，MUSIC譜函數(shù)能夠以圖1(a)所示的尖峰表示聲源DOA。然而，當(dāng)信噪比較低時(shí)，峰值可能被扭曲，估計(jì)的DOA可能偏離真實(shí)值，如圖1(b)所示。圖1(b)信噪比SNR=0圖1(a)信噪比SNR=10■■0024LO-1"T"T-1rmpjupnpUEftUJa-p-llPFMEmEUJD21.3Capon算法在空間譜估計(jì)算法中，Capon算法也是通過譜估計(jì)式的峰值搜索，來(lái)確定入射的方向信息，Capon算法的空間譜函數(shù)為P (oP (o)=Capon1a(0)hR-ia(0)k7)式中R是為陣列協(xié)方差矩陣，經(jīng)過式(4)計(jì)算得到。則DOA角度估計(jì)可以通過k譜峰搜索使得式(7)輸出值最大化的角度信息得到。在與1.2相同的條件下，Capon譜函數(shù)也能夠以圖2所示的尖峰表示聲源DOA。

m—-Bp〕AJPP.-E匚gE-34.Ess4Jpn_cceE10 303D40 90ao7U80m—-Bp〕AJPP.-E匚gE-34.Ess4Jpn_cceE10 303D40 90ao7U80曲figlt(dfigftt)圖2(a)信噪比SNR=1010 203D^10 50 60TO80抽§的fdfrgres)圖2(b)信噪比SNR=01.4ESPRIT算法ESPRIT算法基于空間旋轉(zhuǎn)不變性。ESPRIT算法作為一種子空間角估計(jì)方法，需要對(duì)接收數(shù)據(jù)的交叉譜矩陣進(jìn)行特征值分解(EVD)或奇異值分解(SVD)。該算法要求接收信號(hào)具有一定的不變性結(jié)構(gòu)，而不變性結(jié)構(gòu)通常依賴于已知的兩束相同子陣列之間的平移。該算法不需要對(duì)譜峰進(jìn)行搜索，可以獲得自動(dòng)配對(duì)的參數(shù)估計(jì)。該算法為DOA估計(jì)提供了封閉解。由方程(3)-(4)，我們可以得到接收信號(hào)的協(xié)方差矩陣，并可以表示為8)R=U刀Uh+U刀Uh8)k SSS NNN通過R特征分解得到信號(hào)子空間U，并可以表示為kU=ATgCmxnS其中A為方向矩陣,T為滿秩矩陣。且U和A都可以分成兩部分，即S一U■firstrow其中A為方向矩陣,T為滿秩矩陣。且U和A都可以分成兩部分，即S一U■firstrowS1—lastrowUS2US10)其中a=ae，21則結(jié)合(10)-(12)得A1lastrowfirstrow

A2e一j2ndsin0J九e一j2ndsin0N九11)12)U=AT,U=AT=A6TS11 S2 21U=UT-iQTS2 S1定義屮=T-iQT=U+U，則屮的特征值是?的對(duì)角元素，U+U=UT-iQTS2 S1定義屮=T-iQT=U+U，則屮的特征值是?的對(duì)角元素，U+是U的廣義逆矩SiS2 SiSi陣。接下來(lái)只要計(jì)算V的特征值九,九，…,九，再利用相位角函數(shù)angle和反三角1 2 N函數(shù)arcsin就可以得出角度的估計(jì)值，即O=arcsin\angle\k)九;'2兀d)n n115)現(xiàn)考慮一個(gè)10陣元數(shù)的均勻線陣，有三個(gè)信源入射，它們的波達(dá)方向(D0A)分別是10°,30°和50°,ESPRIT算法的估計(jì)效果如圖3所示。片hMjfl嚴(yán)吹薩吧SNR=5rSrupshnlFlW55SO453025201510a20100圖3信噪比SNR=5dB40 60TinWs1.5PM(譜峰搜索)方程(4)中陣列接收信號(hào)協(xié)方差矩陣R可以寫成兩部分kRk=[gh]16)其中G是由R的第1列到第N列組成，H是R剩下的(M-N)列組成。定義傳kk播矩陣P=(GhGLGhHeCnx(m-n)，Q=「Ph(-I)「L M-N」是單位矩陣，則空間譜函數(shù)為eC(M-N)xM,其中IM-NPpM(°)=a(O)hQhQa(O)17)DOA角度估計(jì)可以通過譜峰搜索使得式(15)輸出值最大化的角度信息得到。在與1.2相同的條件下,PM譜峰搜索算法的空間譜函數(shù)也能夠以圖4所示的尖峰表示聲源DOA。從圖4(b)可以看出，當(dāng)信噪比較低時(shí)，譜峰搜索的效果變低。圖4（a圖4（a）信噪比SNR=101.6PM（旋轉(zhuǎn)不變子空間）AaAb利用旋轉(zhuǎn)不變性的特性，基于旋轉(zhuǎn)不變性的PM實(shí)現(xiàn)了不需要搜索譜峰的角度估計(jì)。考慮到方向矩陣AAaAb18)其中AecNxN是滿秩矩陣，AwC（M-N川。定義線性算子P，使得滿足a b cA=PA，基于P，我們建立矩陣P，滿足bca c19)_I19)P=N

Pc因此，我們有PAaIPAaINPc20)那么，我們將矩陣P分成兩部分，P和P。其中P 是由矩陣P的第1121行到第（M-1）行構(gòu)成，PeC（MTkN是由矩陣P的第2行到第M行組成。則（20）2可寫成PAa一PPAa一P1A=_A-1ilastrowalastrow21)PAafirstrowA=firstrow=firstrowP2aA2A01其中AwC（PAafirstrowA=firstrow=firstrowP2aA2A01其中AwC（m-、n是方向矩陣a的第1行到第（M-1）行,iAWC（M-1人N是方向矩2陣A的第2行到第M行，0在式（12）中已經(jīng)定義。由（21）-（22）可知PA=A=A0=PA0，貝V2a211aP+P=P+PA0A-i=A0A-iwCnxn (23)1 2 1 1aa aa其中P+是P的廣義逆矩陣。定義申=P+P=A0A-1wCwn，與矩陣0有同樣的1 1 1 2 aa特征值，通過特征分解，我們可以得到N個(gè)不同大小的特征值九，…,九。最后,1 N角度估計(jì)可以由下式得到：=arcsin\-angle\k)九；2兀d丿n24)同式(16),R=[GH],通常情況下，P=[G+H]H。kc現(xiàn)考慮一個(gè)10陣元數(shù)的均勻線陣，有三個(gè)信源入射，它們的波達(dá)方向（DOA）分別是10°,30°和50°,PM的估計(jì)效果如圖5所示?！?1t)l3rOrj=3Q°10.=M<51SNR=5l5nap9tiot^100￡ 百20-15a 2Q 60 sd iaoTima/s圖5信噪比SNR=10dB2.算法性能比較實(shí)驗(yàn)考慮一個(gè)10陣元數(shù)的均勻線陣，有三個(gè)信源入射，它們的波達(dá)方向（DOA）分別是10°,30°和50°，初始快拍數(shù)L=50，信噪比SNR=10dB。在以下仿真實(shí)驗(yàn)中，估計(jì)性能的評(píng)價(jià)方式為均方根誤差（RMSE），定義為RMSE二丄另丄更(Q—0)NMC k,jk,j廠是角度真實(shí)值。k,jn=1 j=廠是角度真實(shí)值。k,j其中0是第j次蒙特卡洛試驗(yàn)的角度估計(jì)值，k,j實(shí)驗(yàn)一：信噪比SNR=5dB，快拍數(shù)從50到500變化，變化步長(zhǎng)50。運(yùn)行500次蒙特卡洛（MC）試驗(yàn)，上述幾種算法的對(duì)比仿真效果如圖6所示。101010Q 150 200105025D3TOMD4DD4W101010Q 150 200105025D3TOMD4DD4W500SnapsFiot/s-3 musice&lim^lion■ caponestimaliDh―8—espirrtesBirnation0 PM-spectrumestirnatiDnPM-rotaLienalegtimstioni圖6信噪比SNR=5dB實(shí)驗(yàn)二：快拍數(shù)L=50不變，信噪比SNR從5dB變化到20dB,變化步長(zhǎng)為5dB,運(yùn)行500次蒙特卡洛（MC）試驗(yàn)，上述幾種算法的對(duì)比仿真效果如圖7所示。10101015SNR/dB—10101015SNR/dB—?—musicestimation―?—caponestimation—BspifilEstimation0PM-spectrurriestimatioinPM-not^lionBI"estimation圖7快拍數(shù)L=50在上一節(jié)中，我們介紹了均勻線陣的一維角