隨機(jī)信號(hào)的DOA估計(jì)方法

1、 隨機(jī)信號(hào)的DOA估計(jì)方法一、 實(shí)驗(yàn)?zāi)康?、 掌握利用周期圖法、Capon 方法、MUSIC方法實(shí)現(xiàn)隨機(jī)信號(hào)DOA估計(jì)的方法。2、 觀察陣元數(shù)目、陣元間距、信噪比、入射方向等參數(shù)對(duì)角度譜估計(jì)性能的影響。3、 理解特征結(jié)構(gòu)類方法進(jìn)行DOA估計(jì)的優(yōu)點(diǎn)。二、 實(shí)驗(yàn)原理信號(hào)的來波方向（DOA）估計(jì)石陣列信號(hào)處理領(lǐng)域中的一個(gè)重要內(nèi)容。陣列信號(hào)模型如圖1所示，設(shè)均勻陣列中有M個(gè)陣元，陣元間距為D，記信號(hào)波長(zhǎng)為，則陣列等效孔徑為，表示入射到陣列上信號(hào)的來波方向（DOA），以信號(hào)傳播方向與陣列法線方向的夾角來表示（順時(shí)針方向?yàn)檎D 1陣列信號(hào)模型假設(shè)信號(hào)源位于遠(yuǎn)場(chǎng)，即信號(hào)在到達(dá)各個(gè)陣元時(shí)的波前為平面波，以

2、原點(diǎn)處的陣元為參考點(diǎn)，則個(gè)陣元接收到的信號(hào)為 （1）其中，為信號(hào)中心頻率，為波長(zhǎng)。對(duì)于窄帶解析復(fù)信號(hào)，有（2）其中為角頻率，則第i個(gè)陣元上收到的信號(hào)可以表示為（3）如果有d個(gè)入射源信號(hào)，它們的入射角分別為，則有（4）M個(gè)陣元接收到的信號(hào)用矩陣表示為（5）其中將矩陣寫成如下形式，這里為導(dǎo)向矢量。信號(hào)的DOA估計(jì)大多采用搜索夫人方法，通過對(duì)譜估計(jì)函數(shù)進(jìn)行峰值搜索，估計(jì)信號(hào)波數(shù)到達(dá)的方向。本實(shí)驗(yàn)將對(duì)周期圖法、Capon法以及MUSIC方法予以討論。1 周期圖法已知接收信號(hào)觀測(cè)樣本序列為有限長(zhǎng)序列，記，其自相關(guān)矩陣為。本次實(shí)驗(yàn)中根據(jù)各態(tài)歷經(jīng)假設(shè)，對(duì)次快拍求平均估計(jì)自相關(guān)矩陣，從而有。使用周期圖方法進(jìn)行

3、角度譜估計(jì)的結(jié)果為（6）（7）因此可以通過譜峰搜索估計(jì)信號(hào)的波達(dá)方向。2 Capon方法Capon方法是一種利用空域?yàn)V波實(shí)現(xiàn)DOA估計(jì)的方法，通過在期望方向形成波束，并利用剩余自由度在干擾方向形成零陷從而一致干擾和噪聲。為了保證輸出信號(hào)中包含期望方向的信號(hào)，通常保持其外那個(gè)方向的增益為1，在此基礎(chǔ)上使輸出信號(hào)的功率最小，輸出功率反映了信號(hào)的空間譜特性。Capon方法可以描述成一個(gè)優(yōu)化問題：在約束條件或的前提下，使輸出功率最小，采用拉各朗日乘子法可得最優(yōu)解為，對(duì)應(yīng)的角度譜估計(jì)結(jié)果為（8）其中。3 MUSIC方法（5）式的接收信號(hào)形式中，為入射信號(hào)，為白噪聲，如果個(gè)新號(hào)批次獨(dú)立，且與噪聲不相關(guān)，則

4、有（9）其中。自相關(guān)矩陣的特征值和特征向量為：（10）因此可將表示為：，將特征值按從大到小排序，即：其特征結(jié)構(gòu)性質(zhì)為：（1）（11）（2）（12）個(gè)大特征值對(duì)應(yīng)的特征向量被稱為信號(hào)特征向量，而個(gè)小特征值對(duì)應(yīng)的特征向量被稱為噪聲特征向量，它們張成的空間是信號(hào)特征向量張成的空間的正交補(bǔ)。MUSIC方法將接收信號(hào)的自相關(guān)矩陣進(jìn)行特征分解，得到信號(hào)子空間和噪聲子空間，利用噪聲特征向量和信號(hào)特征向量的正交關(guān)系，通過譜峰搜索可以估計(jì)到達(dá)波方向，從而得到譜估計(jì)結(jié)果為（13）其中是由噪聲特征向量張成的噪聲子空間。搜索尋找譜峰來估計(jì)到達(dá)角度。三、 實(shí)驗(yàn)步驟1、仔細(xì)閱讀有關(guān)角度譜估計(jì)的內(nèi)容，按照?qǐng)D2給出的流程圖編

5、寫使用周期圖法進(jìn)行角度譜估計(jì)程序。選擇參數(shù)M=8,d=2=0.5,N=100,SNR=10dB,中心頻率f=300MHz,兩個(gè)信號(hào)源的入射角度分別為-20度和30度，觀察譜估計(jì)圖如圖1所示。圖1：周期圖角度譜估計(jì)如圖1所示，譜峰位于-20°和30°處，與信號(hào)源的入射角度相符。2、觀察并記錄參數(shù)變化對(duì)角度譜估計(jì)性能的影響。（1） 改變陣元間距d=4=0.25, d=2=2，其他輸入同步驟 1，觀察譜峰位置是否正確，以及陣元間距對(duì)譜分辨率的影響。d=4 d=2圖2：周期圖角度譜估計(jì)（陣元間距）對(duì)于空間采樣，與奈奎斯特采樣定理類似，必須滿足d2，所以d=4時(shí)，譜峰位置正確，而d=

6、2時(shí)，類似欠采樣，除了在正確的角度上出現(xiàn)譜峰外，在其他多處地方也出現(xiàn)了譜峰，所以譜峰位置不正確。而由d=4與d=2時(shí)的譜峰對(duì)比可知，d=2時(shí)分辨率更高，因此在滿足d2的條件下，d越大，譜分辨率越高。（2）當(dāng)只有一個(gè)入射源信號(hào)時(shí)，改變=0°,45°,90°，其他輸入同步驟1，觀察入射方向?qū)亲V估計(jì)性能的影響。=0° =45°=90°圖3：周期圖角度譜估計(jì)（入射角度）由圖3可知，在=0°,45°時(shí)，角度譜上在相應(yīng)位置出現(xiàn)正確的譜峰，而在=90°時(shí)，在90°和-90°處都出現(xiàn)一個(gè)譜峰，影響D

7、OA估計(jì)心能。不過，在實(shí)際中，90°和-90°物理朝向相差明顯，可以人工判斷到底是90°還是-90°。（3）改變M=4，16，其它條件同步驟1，觀察陣元數(shù)目對(duì)譜估計(jì)性能的影響。M=4 M=16圖4：周期圖角度譜估計(jì)（陣元數(shù)目）由圖4可知，M=4時(shí)，相對(duì)于M=8的譜峰更寬，分辨率更低；M=16時(shí)，相對(duì)于M=8的譜峰更窄，分辨率更高，但旁瓣數(shù)量也更多。（4）改變信噪比 SNR=0dB,20 dB ，其他輸入同步驟 1，觀察信噪比對(duì)譜估計(jì)性能的影響。SNR=0dB SNR=20 dB圖5：周期圖角度譜估計(jì)（信噪比）由圖5可知，信噪比越大，所得角度譜譜峰比旁瓣要

8、越高，那么當(dāng)信噪比小到一定程度，旁瓣會(huì)越來越高，容易與主瓣混淆，影響譜估計(jì)性能。3、編寫Capon 最大似然法進(jìn)行角度譜估計(jì)的程序。運(yùn)行程序，輸入?yún)?shù)同步驟1，觀察譜峰位置是否正確，并重復(fù)步驟2 的內(nèi)容。步驟1：如圖6所示，譜峰位于-20°和30°處，與信號(hào)源的入射角度相符。圖6：Capon角度譜估計(jì)步驟2：（1） d=4 d=2圖7：Capon角度譜估計(jì)（陣元間距）對(duì)于空間采樣，與奈奎斯特采樣定理類似，必須滿足d2，所以d=4時(shí)，譜峰位置正確，而d=2時(shí)，類似欠采樣，除了在正確的角度上出現(xiàn)譜峰外，在其他多處地方也出現(xiàn)了譜峰，所以譜峰位置不正確。而由d=4與d=2時(shí)的譜峰對(duì)

9、比可知，d=2時(shí)分辨率更高，因此在滿足d2的條件下，d越大，譜分辨率越高。（2）=0° =45°=90°圖8：Capon角度譜估計(jì)（入射角度）由圖8可知，在=0°,45°時(shí)，角度譜上在相應(yīng)位置出現(xiàn)正確的譜峰，而在=90°時(shí)，在90°和-90°處都出現(xiàn)一個(gè)譜峰，影響DOA估計(jì)心能。在實(shí)際中，90°和-90°物理朝向相差明顯，可以人工判斷到底是90°還是-90°。（3）M=4 M=16圖9：Capon角度譜估計(jì)（陣元數(shù)目）由圖9可知，M=4時(shí)，相對(duì)于M=8的譜峰更寬，分辨率更低；

10、M=16時(shí)，相對(duì)于M=8的譜峰更窄，分辨率更高，Capon旁瓣數(shù)量比周期圖旁瓣數(shù)量要少得多。（4）SNR=0dB SNR=20 dB圖10：Capon角度譜估計(jì)（信噪比）由圖10可知，信噪比越大，所得角度譜譜峰比旁瓣要越高，那么當(dāng)信噪比小到一定程度，旁瓣會(huì)越來越高，容易與主瓣混淆，影響譜估計(jì)性能。4、編寫MUSIC方法進(jìn)行角度譜估計(jì)的程序。運(yùn)行程序，輸入?yún)?shù)同步驟1，觀察譜峰位置是否正確，并重復(fù)步驟2 的內(nèi)容。步驟1：如圖11所示，譜峰位于-20°和30°處，與信號(hào)源的入射角度相符。圖11：MUSIC角度譜估計(jì)步驟2：（1） d=4 d=2圖12：MUSIC角度譜估計(jì)（陣元

11、間距）對(duì)于空間采樣，與奈奎斯特采樣定理類似，必須滿足d2，所以d=4時(shí)，譜峰位置正確，而d=2時(shí)，類似欠采樣，除了在正確的角度上出現(xiàn)譜峰外，在其他多處地方也出現(xiàn)了譜峰，所以譜峰位置不正確。而由d=4與d=2時(shí)的譜峰對(duì)比可知，d=2時(shí)分辨率更高，因此在滿足d2的條件下，d越大，譜分辨率越高。（2） =0° =45°=90°圖13：MUSIC角度譜估計(jì)（入射角度）由圖13可知，在=0°,45°時(shí)，角度譜上在相應(yīng)位置出現(xiàn)正確的譜峰，而在=90°時(shí)，在90°和-90°處都出現(xiàn)一個(gè)譜峰，影響DOA估計(jì)心能。（3） M=4 M

12、=16圖14：MUSIC角度譜估計(jì)（陣元數(shù)目）由圖14可知，M=4時(shí)，相對(duì)于M=8的譜峰更寬，分辨率更低；M=16時(shí)，相對(duì)于M=8的譜峰更窄，分辨率更高，MUSIC幾乎沒有旁瓣。（4）SNR=0dB SNR=20 dB圖15：MUSIC角度譜估計(jì)（信噪比）由圖15可知，信噪比越大，所得角度譜譜峰比旁瓣要越高，那么當(dāng)信噪比小到一定程度，旁瓣會(huì)越來越高，容易與主瓣混淆，影響譜估計(jì)性能。5、改變兩個(gè)源信號(hào)入射角度分別為-5°,5°，其他輸入?yún)?shù)同步驟1，將周期圖法、Capon 最大似然法以及MUSIC方法的角度譜估計(jì)結(jié)果畫在同一坐標(biāo)，比較三種方法在分辨率上的差別。如圖16所示，入

13、射角度為-5°,5°，使用周期圖法根本無法區(qū)分，分辨率太低。Capon最大似然法能夠區(qū)分兩個(gè)角度，但兩峰有明顯交界處，分辨率比周期圖要高。MUSIC方法能夠明顯地區(qū)分兩個(gè)角度，且交界處比較少，分辨率最高。所以三種方法分辨率比較為：MUSIC > Capon >周期圖法。圖16：周期圖、Capon、MUSIC角度譜估計(jì)比較圖四、 實(shí)驗(yàn)結(jié)論通過實(shí)驗(yàn)可得：使用周期圖法、Capon法、MUSIC法都能在一定條件下實(shí)現(xiàn)DOA估計(jì)方法。陣元間距d的大小能夠影響到角度譜分辨率，d越大，分辨率越高，但必須滿足d2，否則由于采樣率過低導(dǎo)致結(jié)果錯(cuò)誤。 一般情況下，入射角對(duì)譜估計(jì)性能

14、沒有影響，但當(dāng)入射角接近90度或-90度，譜峰也會(huì)出現(xiàn)在-90度或90度處，影響判斷。陣元數(shù)目越多，分辨率越高，對(duì)于周期圖法，旁瓣越明顯。信噪比越高，所得譜峰越明顯，旁瓣越少。當(dāng)信噪比低于一定程度，旁瓣會(huì)影響到譜估計(jì)結(jié)果。三種方法的分辨率不一樣，其分辨率大小比較為：MUSIC > Capon >周期圖法。五、 思考題1 推導(dǎo)周期圖法進(jìn)行角度譜估計(jì)的矩陣表達(dá)形式：PBT=IM=1MEHRE2 當(dāng)入射源信號(hào)不是彼此獨(dú)立的隨機(jī)信號(hào)，而是頻率相同的正弦信號(hào)（相干信源）時(shí)，上述三種角度譜估計(jì)方法的性能有什么變化？六、 附錄DOA.mclose all;clc;clear;M = 8;lumb

15、t = 1;d = lumbt/2;N = 100;SNR =10;f = 300000000;th1 = -20/180*pi;th2 = 30/180*pi;c=3*108;% generate signalTable = (1/sqrt(2)*-1-1i -1+1i 1-1i 1+1i;T1 = randi(0,1,1,N*2);T_IQ1 = reshape(T1,2,N).'symbol1 = bi2de(T_IQ1,'left-msb');s1 = Table(symbol1+1);T2 = randi(0,1,1,N*2);T_IQ2 = reshape(

16、T2,2,N)'symbol2 = bi2de(T_IQ2,'left-msb');s2 = Table(symbol2+1);s = s1;s2;% generate white complex gaussian noisen = 1/(sqrt(2*10(SNR/10)*(randn(M,N) + 1i*randn(M,N);% Receive transformation matrixw = 2*pi*f;t1 = d*sin(th1)/c;a1 = exp(1i*w*(0:M-1)*t1);t2 = d*sin(th2)/c;a2 = exp(1i*w*(0:M

17、-1)*t2);A = a1;a2.'X = A*s + n;% PeriodogramR = X*X'/N;for k=0:1:180 th = pi/180*k-pi/2; tth = d*sin(th)/c; E = (exp(1i*w*tth*(0:M-1).' P_BT(k+1) = E'*R*E/M;endfigure(1),clfplot(-90:1:90,10*log10(abs(P_BT);% legend('Periodogram');%hold on% CaponR_inv = inv(R);for k=0:1:180 th

18、 = pi/180*k-pi/2; tth = d*sin(th)/c; E = (exp(1i*w*tth*(0:M-1).' P_Capon(k+1) = 1/(E'*R_inv*E);end figure(2),clfplot(-90:1:90,10*log10(abs(P_Capon);% legend('Capon');% MUSICMUSIC_N = n*n'/N;MUSIC_delt = mean(diag(MUSIC_N);MUSIC_S = s*s'/N;MUSIC_R = A*MUSIC_S*A' + MUSIC_N;V,Lbt = eig(MUSIC_R);Lbt1 = diag(Lbt);p = find(abs(Lbt1-MUSIC_delt)<0.1);for k=0:1:180 th =