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

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

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

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

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

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

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

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

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

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

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

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

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

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

14、%信號數(shù)據(jù)s2=sqrt(2)*cos(2*pi*(fc+5e8)*t);ss=s1;s2;s=ss(1:p,:);%計算陣列流形矩陣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;%計算噪聲子空間v,dd=eig(R);if(dd(1,1)dd(2,2) Un=v(:,p+1:M);else Un=v(:,1:(M-p);enddo=-90:90;pu=zeros(1,length(do);kg=

15、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測向);xlabel(波達方向);ylabel(MUSIC譜);disp(MUSIC測向結(jié)果:);for k=1:p k1,k2=max(pu); DOA_guji(k)=(k2-1)-90; pu(k2)=0;endDO

16、A_guji（2）前后向空間平滑MUSIC算法（相關(guān)信號）clear all;clc;p=2; %入射信號數(shù)目M=8; %陣元個數(shù)L=5; %將陣列劃分為相互重疊的L個子陣m=M-L+1; %每個子陣中的陣元個數(shù)fc=1e9; %入射信號中頻為1GDOA=-20,40/180*pi; %信號入射DOAfs=3*fc; %采樣頻率N=512; %采樣個數(shù)snr=10; %信噪比T=1/fs; %采樣時間間隔c=3e8; %波速d=c/fc*0.5; %陣元間距t=0:T:(N-1)*T; %采樣時間區(qū)間s1=sqrt(2)*cos(2*pi*fc*t); %信號數(shù)據(jù)s2=2*s1;ss=s1;s

17、2;s=ss(1:p,:);%計算陣列流形矩陣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); %對協(xié)方差矩陣進行特征值分解%計算噪聲子空間if(dd(1,1)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);

18、 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é)方差矩陣f