基于MUSIC的无线通信测向系统设计与仿真

1、东北大学秦皇岛分校计算机与通信工程学院结课论文基于MUSIC的无线通信测向系统设计与仿真专业名称班级学号学生姓名指导教师设计时间摘 要：本文主要是对DOA（波达方向）估计中传统MUSIC 算法作了简要的介绍，然后通过仿真发现MUSIC算法不适用与相关信号。针对MUSIC算法的不足引出了空间平滑的MUSIC算法，很好的解决了相关信号的问题。关键词：DOA 估计；MUSIC 算法；空间平滑一 引言波达方向(Directionof-Arrival)估计是阵列信号处理领域中的重要的研究方向，它是雷达、声纳、主动防护系统、通讯系统以及智能天线等多个技术领域的共性问题。基于阵列信号处理的波达方向估计方法可

2、以同时对空间不同方向上的多个信号源实现高分辨率的方向估计。对波达方向的估计是空间谱估计研究的主要课题。最经典的超分辨率空间谱估计方法是Schmidt在1979年提出的MUSIC(Mukiple Signal Classification)算法，在模型准确的条件下，该算法能精确地估计空间上互不相关信号的波达方向。由于多径传播、电磁干扰等因素的影响，相干信源存在的电磁环境是经常碰到的。当空间存在相干源时，经典的MUSIC 算法已经失去了其高分辨性能优势，有时甚至不能正确地估计出信源的真实方位。因此，若将其用于相干源，首先对阵列输出的协方差矩阵进行各种去相干处理， 本文采用空间平滑算法，保持了在相干

3、信号下较高的分辨率。二阵列信号处理统计模型在无线通信中我们通过天线对电磁波进行发射和接收。为了增加电磁波的利用率和电磁波的波束形状可控，一般采用阵列天线。在一般情况下，将一组传感器按一定的方式设置在空间不同的位置上组成传感器阵列，此传感器阵列能够接收空间的传播信号，然后对所接收到的信号经过适当的处理并提取所需的信号源和信号属性等信息，包括信号辐射源辐射信号的数目、方向、幅度等。一般来说，构成阵列的阵元可以按照任意的方式进行排列，但是通常是按照直线等距、圆周等距或平面等距排列的，并且取向相同。为了简化天线阵列的分析，通常作如下假设：1. 窄带假设：这样可以保证所有阵元几乎同时接收到该信号，即阵元

4、接收之间的信号包络没有变化；2. 信号的统计特性：假设入射到阵列的信号为平稳且各态历经，这样可以用时间平均来代替统计平均。噪声为互不相关的白噪声，方差为。3. 忽略阵元之间的互耦；4. 信号的数目要小于阵元的数目，并且阵列接收到得所有信号的波达方向互不相同，信号之间互不相关；5. 平面波假设：假设信源到阵列的距离远大于阵列的口径，从而所有入射到阵列的信号波前金额以近似为平面波。假设在天线阵的原唱存在个信号源，则所有到达阵列的波前可近似为平面波。若天线阵由个全向天线组成，将第一个阵元设为参考阵元，则到达参考阵元的第个信号为： （1）式中，为第个信号的复包络，包含信号信息。为空间信号的载波。由于信

5、号满足窄带假设条件，则,那么经过传播延迟后的信号可以表示为： （2）则理想情况下第个阵元接收到的信号可以表示为: （3）式中，为第个阵元到达第个阵元时相对于参考阵元的时延，为第阵元上的加性噪声。根据式（2）和（3）可得，整个天线阵接收到得信号为： （4）式中，为信号的方向向量，为阵列流形，为信号矩阵，为加性噪声矩阵，表示矩阵转置。M321图1 阵列信号模型图三MUSCI算法原理 Music算法是由ROSchmidt于1979年提出来，1986年重新发表的。它是最早的也是最经典的超分辨DOA估计方法，它利用了信号子空间和噪声子空间的正交性，构造空间谱函数，通过谱峰搜索，检测信号的DOA。接收信号

6、的协方差矩阵为： （5）由于假设信号与噪声是不相关的，且噪声为平稳的加性高斯白噪声，因此式（5）中的二，三项为零，且有。则式（5）简化为式（6）： （6） 式（6）中的是有用信号的协方差矩阵。由于假设信号源之间互不相关，因此为满秩矩阵，其秩为。而为维的矩阵，其秩也是，并且是Hermite半正定矩阵，其秩也是。因此，令的特征值为，那么的个特征值为： 它们对应的特征向量分别为，其中前个对应大特征值，后个对应小特征值。由此可以看出，协方差矩阵经过特征值分解后可以产生个较大的特征值和个较小的特征值，并且这个小特征值非常接近。所以当这些小特征值的重数确定了，那么信号的个数就可以由式（7）估计出来： （7

7、） 对于与个最小特征值对应的特征向量，有：，即： ，因为满秩，非奇异，因此：或这表明与个最小特征值对应的特征向量，和个信号特征值对应的方向向量正交，即信号子空间和噪声子空间正交。因此，我们构造维的噪声子空间： 并定义Music空间谱为： （8） 或 （9） 由于信号子空间和噪声子空间正交，所以当等于信号的入射角时，Music空间谱将产生极大值。因此当对Music空间谱搜索时，其个峰值将对应个信号的入射方向，这就是Music算法。四 空间平滑算法1前向空间平滑算法将M个阵元的均匀线阵，分成相互交错的P个子阵，每个子阵包含的阵元数为m个，即满足Mp+m-1。信号源数为N。图2 前向空间算法原理图如

8、图 3所示，取第一个子阵(最左边的子阵)为参考子阵，那么各个子阵的输出矢量分别为: （10） 对于第k个子阵有： （11）其中： （12） 那么该子阵的数据协方差矩阵为： （13）其中，是一个m×p的参考子阵(通常取第一个子阵)的导向矢量矩阵， ，为信号的协方差矩阵，。 前向空间平滑技术是通过求各个子阵协方差矩阵的均值来实现的，即取前向平滑修正的协方差矩阵为: （14）可以证明，当满足m>N, p>N时，前向空间平滑数据协方差矩阵是满秩的。即可以通过特征分解求得相应的信号子空间和噪声子空间。2 前后向空间平滑算法如果按照图4划分阵列，即称为后向平滑的方法划分子阵，那么各个

9、子阵的输出矢量为:图3 后向空间平滑算法原理图 (15)那么，第k个子阵的数据矢量为： （16）比较前向平滑和后向平滑的数据矢量，可以得到前向平滑中第k个子阵与后向平滑中第p-k+1个子阵之间存在如下关系： （17）其中J为m的交换矩阵。 ，所以后向平滑第p-k+1个子阵的数据协方差矩阵为: (18) 那么后向空间平滑修正的数据矩阵为: (19) 取前向平滑和后向平滑数据协方差矩阵的平均，即前后向空间平滑的数据矩阵，即 (20)同样可以证明，当满足m>N, p>N时，后向空间平滑数据协方差矩阵是满秩的。五测向系统组成空间谱估计测向系统是利用信号来波在各天线阵元上感应产生的电压的幅度

10、和相位与来波方向有关这一特性来实现对空间多个信号同时测向的。系统由天线阵、电子开关、放大器、接收机、信号采集与控制、终端计算机、校准源和电源等单元组成。基本构成如图：图4 测向系统组成六matlab仿真流程Music算法的步骤归纳如下：（1）收集信号样本，其中为采样点数，估计协方差函数：（2）对进行特征值分解： 式中为特征值对角阵，且从大到小顺序排列是对应的特征向量。（3）利用最小特征值的重数，估计信号数，并构造噪声子空间。（4）搜索Music空间谱，找出个峰值，得到DOA估计值。图5 MUSIC 算法软件流程七. 仿真分析（1）当入射信号为非相干信号时，设输入信号数为2，入射角度分别为-20

11、 40，阵元数为4，快拍数为512，阵元间距为 /2（ 表示波长）。采用普通MUSIC算法可以准确的得出结果，结果见图(6)。图（6） 非相干信号仿真结果(2) 当入射信号为相干信号时(s1=2s2)，设输入信号数为2，入射角度分别为-20 40，阵元数为4，快拍数为512，阵元间距为 /2（ 表示波长）。采用普通MUSIC算法不能得出结果，见图（7）。图（7） 相干信号仿真结果(3) 当入射信号为相干信号时(s1=2s2)，设输入信号数为2，入射角度分别为-20 40，阵元数为4，快拍数为512，阵元间距为 /2（ 表示波长）。采用空间平滑MUSIC算法可以准确得出结果，见图（8）。图（8）

12、 经典MUSIC与前后向平滑MUSIC对比八结语采用MUSIC算法能构造出针状的谱峰，可以很好的估计出入射信号的个数和方向，能有效的估计出独立信号源的DOA,克服了传统测向定位方法精度低的缺点 ,可以有效解决密集信号环境中多个辐射源的高分辨率、高精度测向定位问题。但若存在相干信源时，阵列输出信号协方差的秩 ，对信号协方差矩阵进行特征值分解后，得到的较大的特征值个数小于P，而特征值为 的个数将大于M P。与此相对应的信号子空间的向量也少于P，即特征向量展开的信号子空间的维数少于的列数。对某些相干源的方向矢量， 将不正交于噪声子空间，不出现零点，所以，有些源在空间谱曲线中将不呈现峰值，造成谱估计的

13、漏报。因此，我们要采用空间平滑算法对阵列输出的协方差矩阵首先进行预处理，使其阵列协方差矩阵的秩恢复为信号元数，然后再采用MUSIC算法，能够很好的测出信号方向。附录:matlab仿真源码（1）普通MUSIC算法（非相关信号）clear all;clc;p=2; %入射信号数目M=4; %阵元个数fc=1e9; %入射信号中频为1GDOA=-20,40/180*pi; %信号入射DOAfs=3*fc; %采样频率N=512; %采样个数snr=10; %信噪比T=1/fs; %采样时间间隔c=3e8; %波速d=c/fc*0.5; %阵元间距t=0:T:(N-1)*T; %采样时间区间s1=sq

14、rt(2)*cos(2*pi*fc*t); %信号数据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%阵列接收数据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=-

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测向');xlabel('波达方向');ylabel('M

16、USIC谱');disp('MUSIC测向结果:');for k=1:p k1,k2=max(pu); DOA_guji(k)=(k2-1)-90; pu(k2)=0;endDOA_guji（2）前后向空间平滑MUSIC算法（相关信号）clear all;clc;p=2; %入射信号数目M=8; %阵元个数L=5; %将阵列划分为相互重叠的L个子阵m=M-L+1; %每个子阵中的阵元个数fc=1e9; %入射信号中频为1GDOA=-20,40/180*pi; %信号入射DOAfs=3*fc; %采样频率N=512; %采样个数snr=10; %信噪比T=1/fs; %采

17、样时间间隔c=3e8; %波速d=c/fc*0.5; %阵元间距t=0:T:(N-1)*T; %采样时间区间s1=sqrt(2)*cos(2*pi*fc*t); %信号数据s2=2*s1;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); endendy=A*s;y=awgn(y,snr);%阵列接收数据R=y*y'/N;v,dd=eig(R); %对协方差矩阵进行特征值分解%计算噪声子空间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);%前向平滑数据协方差矩阵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);%后向平滑数据协方差矩阵for i=1:m for k=1:N X=conj(y(M-i+1:-1:M-L-i+2,:); %取复数共轭 rb