基于MUSIC的無線通信測(cè)向系統(tǒng)設(shè)計(jì)與仿真

1、東北大學(xué)秦皇島分校計(jì)算機(jī)與通信工程學(xué)院結(jié)課論文基于MUSIC的無線通信測(cè)向系統(tǒng)設(shè)計(jì)與仿真專業(yè)名稱班級(jí)學(xué)號(hào)學(xué)生姓名指導(dǎo)教師設(shè)計(jì)時(shí)間摘 要：本文主要是對(duì)DOA（波達(dá)方向）估計(jì)中傳統(tǒng)MUSIC 算法作了簡(jiǎn)要的介紹，然后通過仿真發(fā)現(xiàn)MUSIC算法不適用與相關(guān)信號(hào)。針對(duì)MUSIC算法的不足引出了空間平滑的MUSIC算法，很好的解決了相關(guān)信號(hào)的問題。關(guān)鍵詞：DOA 估計(jì)；MUSIC 算法；空間平滑一 引言波達(dá)方向(Directionof-Arrival)估計(jì)是陣列信號(hào)處理領(lǐng)域中的重要的研究方向，它是雷達(dá)、聲納、主動(dòng)防護(hù)系統(tǒng)、通訊系統(tǒng)以及智能天線等多個(gè)技術(shù)領(lǐng)域的共性問題?；陉嚵行盘?hào)處理的波達(dá)方向估計(jì)方法可

2、以同時(shí)對(duì)空間不同方向上的多個(gè)信號(hào)源實(shí)現(xiàn)高分辨率的方向估計(jì)。對(duì)波達(dá)方向的估計(jì)是空間譜估計(jì)研究的主要課題。最經(jīng)典的超分辨率空間譜估計(jì)方法是Schmidt在1979年提出的MUSIC(Mukiple Signal Classification)算法，在模型準(zhǔn)確的條件下，該算法能精確地估計(jì)空間上互不相關(guān)信號(hào)的波達(dá)方向。由于多徑傳播、電磁干擾等因素的影響，相干信源存在的電磁環(huán)境是經(jīng)常碰到的。當(dāng)空間存在相干源時(shí)，經(jīng)典的MUSIC 算法已經(jīng)失去了其高分辨性能優(yōu)勢(shì)，有時(shí)甚至不能正確地估計(jì)出信源的真實(shí)方位。因此，若將其用于相干源，首先對(duì)陣列輸出的協(xié)方差矩陣進(jìn)行各種去相干處理， 本文采用空間平滑算法，保持了在相干

3、信號(hào)下較高的分辨率。二陣列信號(hào)處理統(tǒng)計(jì)模型在無線通信中我們通過天線對(duì)電磁波進(jìn)行發(fā)射和接收。為了增加電磁波的利用率和電磁波的波束形狀可控，一般采用陣列天線。在一般情況下，將一組傳感器按一定的方式設(shè)置在空間不同的位置上組成傳感器陣列，此傳感器陣列能夠接收空間的傳播信號(hào)，然后對(duì)所接收到的信號(hào)經(jīng)過適當(dāng)?shù)奶幚聿⑻崛∷璧男盘?hào)源和信號(hào)屬性等信息，包括信號(hào)輻射源輻射信號(hào)的數(shù)目、方向、幅度等。一般來說，構(gòu)成陣列的陣元可以按照任意的方式進(jìn)行排列，但是通常是按照直線等距、圓周等距或平面等距排列的，并且取向相同。為了簡(jiǎn)化天線陣列的分析，通常作如下假設(shè)：1. 窄帶假設(shè)：這樣可以保證所有陣元幾乎同時(shí)接收到該信號(hào)，即陣元

4、接收之間的信號(hào)包絡(luò)沒有變化；2. 信號(hào)的統(tǒng)計(jì)特性：假設(shè)入射到陣列的信號(hào)為平穩(wěn)且各態(tài)歷經(jīng)，這樣可以用時(shí)間平均來代替統(tǒng)計(jì)平均。噪聲為互不相關(guān)的白噪聲，方差為。3. 忽略陣元之間的互耦；4. 信號(hào)的數(shù)目要小于陣元的數(shù)目，并且陣列接收到得所有信號(hào)的波達(dá)方向互不相同，信號(hào)之間互不相關(guān)；5. 平面波假設(shè)：假設(shè)信源到陣列的距離遠(yuǎn)大于陣列的口徑，從而所有入射到陣列的信號(hào)波前金額以近似為平面波。假設(shè)在天線陣的原唱存在個(gè)信號(hào)源，則所有到達(dá)陣列的波前可近似為平面波。若天線陣由個(gè)全向天線組成，將第一個(gè)陣元設(shè)為參考陣元，則到達(dá)參考陣元的第個(gè)信號(hào)為： （1）式中，為第個(gè)信號(hào)的復(fù)包絡(luò)，包含信號(hào)信息。為空間信號(hào)的載波。由于信

5、號(hào)滿足窄帶假設(shè)條件，則,那么經(jīng)過傳播延遲后的信號(hào)可以表示為： （2）則理想情況下第個(gè)陣元接收到的信號(hào)可以表示為: （3）式中，為第個(gè)陣元到達(dá)第個(gè)陣元時(shí)相對(duì)于參考陣元的時(shí)延，為第陣元上的加性噪聲。根據(jù)式（2）和（3）可得，整個(gè)天線陣接收到得信號(hào)為： （4）式中，為信號(hào)的方向向量，為陣列流形，為信號(hào)矩陣，為加性噪聲矩陣，表示矩陣轉(zhuǎn)置。M321圖1 陣列信號(hào)模型圖三MUSCI算法原理 Music算法是由ROSchmidt于1979年提出來，1986年重新發(fā)表的。它是最早的也是最經(jīng)典的超分辨DOA估計(jì)方法，它利用了信號(hào)子空間和噪聲子空間的正交性，構(gòu)造空間譜函數(shù)，通過譜峰搜索，檢測(cè)信號(hào)的DOA。接收信號(hào)

6、的協(xié)方差矩陣為： （5）由于假設(shè)信號(hào)與噪聲是不相關(guān)的，且噪聲為平穩(wěn)的加性高斯白噪聲，因此式（5）中的二，三項(xiàng)為零，且有。則式（5）簡(jiǎn)化為式（6）： （6） 式（6）中的是有用信號(hào)的協(xié)方差矩陣。由于假設(shè)信號(hào)源之間互不相關(guān)，因此為滿秩矩陣，其秩為。而為維的矩陣，其秩也是，并且是Hermite半正定矩陣，其秩也是。因此，令的特征值為，那么的個(gè)特征值為： 它們對(duì)應(yīng)的特征向量分別為，其中前個(gè)對(duì)應(yīng)大特征值，后個(gè)對(duì)應(yīng)小特征值。由此可以看出，協(xié)方差矩陣經(jīng)過特征值分解后可以產(chǎn)生個(gè)較大的特征值和個(gè)較小的特征值，并且這個(gè)小特征值非常接近。所以當(dāng)這些小特征值的重?cái)?shù)確定了，那么信號(hào)的個(gè)數(shù)就可以由式（7）估計(jì)出來： （7

7、） 對(duì)于與個(gè)最小特征值對(duì)應(yīng)的特征向量，有：，即： ，因?yàn)闈M秩，非奇異，因此：或這表明與個(gè)最小特征值對(duì)應(yīng)的特征向量，和個(gè)信號(hào)特征值對(duì)應(yīng)的方向向量正交，即信號(hào)子空間和噪聲子空間正交。因此，我們構(gòu)造維的噪聲子空間： 并定義Music空間譜為： （8） 或 （9） 由于信號(hào)子空間和噪聲子空間正交，所以當(dāng)?shù)扔谛盘?hào)的入射角時(shí)，Music空間譜將產(chǎn)生極大值。因此當(dāng)對(duì)Music空間譜搜索時(shí)，其個(gè)峰值將對(duì)應(yīng)個(gè)信號(hào)的入射方向，這就是Music算法。四 空間平滑算法1前向空間平滑算法將M個(gè)陣元的均勻線陣，分成相互交錯(cuò)的P個(gè)子陣，每個(gè)子陣包含的陣元數(shù)為m個(gè)，即滿足Mp+m-1。信號(hào)源數(shù)為N。圖2 前向空間算法原理圖如

8、圖 3所示，取第一個(gè)子陣(最左邊的子陣)為參考子陣，那么各個(gè)子陣的輸出矢量分別為: （10） 對(duì)于第k個(gè)子陣有： （11）其中： （12） 那么該子陣的數(shù)據(jù)協(xié)方差矩陣為： （13）其中，是一個(gè)m×p的參考子陣(通常取第一個(gè)子陣)的導(dǎo)向矢量矩陣， ，為信號(hào)的協(xié)方差矩陣，。 前向空間平滑技術(shù)是通過求各個(gè)子陣協(xié)方差矩陣的均值來實(shí)現(xiàn)的，即取前向平滑修正的協(xié)方差矩陣為: （14）可以證明，當(dāng)滿足m>N, p>N時(shí)，前向空間平滑數(shù)據(jù)協(xié)方差矩陣是滿秩的。即可以通過特征分解求得相應(yīng)的信號(hào)子空間和噪聲子空間。2 前后向空間平滑算法如果按照?qǐng)D4劃分陣列，即稱為后向平滑的方法劃分子陣，那么各個(gè)

9、子陣的輸出矢量為:圖3 后向空間平滑算法原理圖 (15)那么，第k個(gè)子陣的數(shù)據(jù)矢量為： （16）比較前向平滑和后向平滑的數(shù)據(jù)矢量，可以得到前向平滑中第k個(gè)子陣與后向平滑中第p-k+1個(gè)子陣之間存在如下關(guān)系： （17）其中J為m的交換矩陣。 ，所以后向平滑第p-k+1個(gè)子陣的數(shù)據(jù)協(xié)方差矩陣為: (18) 那么后向空間平滑修正的數(shù)據(jù)矩陣為: (19) 取前向平滑和后向平滑數(shù)據(jù)協(xié)方差矩陣的平均，即前后向空間平滑的數(shù)據(jù)矩陣，即 (20)同樣可以證明，當(dāng)滿足m>N, p>N時(shí)，后向空間平滑數(shù)據(jù)協(xié)方差矩陣是滿秩的。五測(cè)向系統(tǒng)組成空間譜估計(jì)測(cè)向系統(tǒng)是利用信號(hào)來波在各天線陣元上感應(yīng)產(chǎn)生的電壓的幅度

10、和相位與來波方向有關(guān)這一特性來實(shí)現(xiàn)對(duì)空間多個(gè)信號(hào)同時(shí)測(cè)向的。系統(tǒng)由天線陣、電子開關(guān)、放大器、接收機(jī)、信號(hào)采集與控制、終端計(jì)算機(jī)、校準(zhǔn)源和電源等單元組成?；緲?gòu)成如圖：圖4 測(cè)向系統(tǒng)組成六matlab仿真流程Music算法的步驟歸納如下：（1）收集信號(hào)樣本，其中為采樣點(diǎn)數(shù)，估計(jì)協(xié)方差函數(shù)：（2）對(duì)進(jìn)行特征值分解： 式中為特征值對(duì)角陣，且從大到小順序排列是對(duì)應(yīng)的特征向量。（3）利用最小特征值的重?cái)?shù)，估計(jì)信號(hào)數(shù)，并構(gòu)造噪聲子空間。（4）搜索Music空間譜，找出個(gè)峰值，得到DOA估計(jì)值。圖5 MUSIC 算法軟件流程七. 仿真分析（1）當(dāng)入射信號(hào)為非相干信號(hào)時(shí)，設(shè)輸入信號(hào)數(shù)為2，入射角度分別為-20

11、 40，陣元數(shù)為4，快拍數(shù)為512，陣元間距為 /2（ 表示波長(zhǎng)）。采用普通MUSIC算法可以準(zhǔn)確的得出結(jié)果，結(jié)果見圖(6)。圖（6） 非相干信號(hào)仿真結(jié)果(2) 當(dāng)入射信號(hào)為相干信號(hào)時(shí)(s1=2s2)，設(shè)輸入信號(hào)數(shù)為2，入射角度分別為-20 40，陣元數(shù)為4，快拍數(shù)為512，陣元間距為 /2（ 表示波長(zhǎng)）。采用普通MUSIC算法不能得出結(jié)果，見圖（7）。圖（7） 相干信號(hào)仿真結(jié)果(3) 當(dāng)入射信號(hào)為相干信號(hào)時(shí)(s1=2s2)，設(shè)輸入信號(hào)數(shù)為2，入射角度分別為-20 40，陣元數(shù)為4，快拍數(shù)為512，陣元間距為 /2（ 表示波長(zhǎng)）。采用空間平滑MUSIC算法可以準(zhǔn)確得出結(jié)果，見圖（8）。圖（8）

12、 經(jīng)典MUSIC與前后向平滑MUSIC對(duì)比八結(jié)語采用MUSIC算法能構(gòu)造出針狀的譜峰，可以很好的估計(jì)出入射信號(hào)的個(gè)數(shù)和方向，能有效的估計(jì)出獨(dú)立信號(hào)源的DOA,克服了傳統(tǒng)測(cè)向定位方法精度低的缺點(diǎn) ,可以有效解決密集信號(hào)環(huán)境中多個(gè)輻射源的高分辨率、高精度測(cè)向定位問題。但若存在相干信源時(shí)，陣列輸出信號(hào)協(xié)方差的秩 ，對(duì)信號(hào)協(xié)方差矩陣進(jìn)行特征值分解后，得到的較大的特征值個(gè)數(shù)小于P，而特征值為 的個(gè)數(shù)將大于M P。與此相對(duì)應(yīng)的信號(hào)子空間的向量也少于P，即特征向量展開的信號(hào)子空間的維數(shù)少于的列數(shù)。對(duì)某些相干源的方向矢量， 將不正交于噪聲子空間，不出現(xiàn)零點(diǎn)，所以，有些源在空間譜曲線中將不呈現(xiàn)峰值，造成譜估計(jì)的

13、漏報(bào)。因此，我們要采用空間平滑算法對(duì)陣列輸出的協(xié)方差矩陣首先進(jìn)行預(yù)處理，使其陣列協(xié)方差矩陣的秩恢復(fù)為信號(hào)元數(shù)，然后再采用MUSIC算法，能夠很好的測(cè)出信號(hào)方向。附錄:matlab仿真源碼（1）普通MUSIC算法（非相關(guān)信號(hào)）clear all;clc;p=2; %入射信號(hào)數(shù)目M=4; %陣元個(gè)數(shù)fc=1e9; %入射信號(hào)中頻為1GDOA=-20,40/180*pi; %信號(hào)入射DOAfs=3*fc; %采樣頻率N=512; %采樣個(gè)數(shù)snr=10; %信噪比T=1/fs; %采樣時(shí)間間隔c=3e8; %波速d=c/fc*0.5; %陣元間距t=0:T:(N-1)*T; %采樣時(shí)間區(qū)間s1=sq

14、rt(2)*cos(2*pi*fc*t); %信號(hào)數(shù)據(jù)s2=sqrt(2)*cos(2*pi*(fc+5e8)*t);ss=s1;s2;s=ss(1:p,:);%計(jì)算陣列流形矩陣AA=zeros(M,p);for k=1:p for kk=1:M A(kk,k)=exp(-j*2*pi*fc*(kk-1)*d*sin(DOA(k)/c); endend%陣列接收數(shù)據(jù)y=A*s;y=awgn(y,snr);R=y*y'/N;%計(jì)算噪聲子空間v,dd=eig(R);if(dd(1,1)>dd(2,2) Un=v(:,p+1:M);else Un=v(:,1:(M-p);enddo=-

15、90:90;pu=zeros(1,length(do);kg=1;for k=-90:90 a=zeros(M,1); for kk=1:M a(kk,1)=exp(-j*2*pi*fc*(kk-1)*d*sin(k/180*pi)/c); end pu(1,kg)=1/(a'*Un*Un'*a); kg=kg+1;endplot(do,10*log10(abs(pu),'-r','linewidth',2);grid on;title('MUSIC測(cè)向');xlabel('波達(dá)方向');ylabel('M

16、USIC譜');disp('MUSIC測(cè)向結(jié)果:');for k=1:p k1,k2=max(pu); DOA_guji(k)=(k2-1)-90; pu(k2)=0;endDOA_guji（2）前后向空間平滑MUSIC算法（相關(guān)信號(hào)）clear all;clc;p=2; %入射信號(hào)數(shù)目M=8; %陣元個(gè)數(shù)L=5; %將陣列劃分為相互重疊的L個(gè)子陣m=M-L+1; %每個(gè)子陣中的陣元個(gè)數(shù)fc=1e9; %入射信號(hào)中頻為1GDOA=-20,40/180*pi; %信號(hào)入射DOAfs=3*fc; %采樣頻率N=512; %采樣個(gè)數(shù)snr=10; %信噪比T=1/fs; %采

17、樣時(shí)間間隔c=3e8; %波速d=c/fc*0.5; %陣元間距t=0:T:(N-1)*T; %采樣時(shí)間區(qū)間s1=sqrt(2)*cos(2*pi*fc*t); %信號(hào)數(shù)據(jù)s2=2*s1;ss=s1;s2;s=ss(1:p,:);%計(jì)算陣列流形矩陣AA=zeros(M,p);for k=1:p for kk=1:M A(kk,k)=exp(-j*2*pi*fc*(kk-1)*d*sin(DOA(k)/c); endendy=A*s;y=awgn(y,snr);%陣列接收數(shù)據(jù)R=y*y'/N;v,dd=eig(R); %對(duì)協(xié)方差矩陣進(jìn)行特征值分解%計(jì)算噪聲子空間if(dd(1,1)>

18、;dd(2,2) Un=v(:,p+1:m);else Un=v(:,1:(L-p);enddo=-90:90;pu1=zeros(1,length(do);kg=1;for k=-90:90 a=zeros(M,1); for kk=1:M a(kk,1)=exp(-j*2*pi*fc*(kk-1)*d*sin(k/180*pi)/c); end pu1(1,kg)=abs(1/(a'*Un*Un'*a); kg=kg+1;endhold on;rf=zeros(L,L);rb=zeros(L,L);Z=zeros(L,N);X=zeros(L,N);Rf=zeros(L,L);Rb=zeros(L,L);%前向平滑數(shù)據(jù)協(xié)方差矩陣for i=1:m for k=1:N Z=y(i:i+L-1,:) rf=rf+Z(:,k)*Z(:,k)' end Rf=Rf+rf;endRf=Rf/(N*m);%后向平滑數(shù)據(jù)協(xié)方差矩陣for i=1:m for k=1:N X=conj(y(M-i+1:-1:M-L-i+2,:); %取復(fù)數(shù)共軛 rb