阵列式信号处理-MUSIC算法仿真

阵列式信号处理MUSIC算法仿真摘要：波达方向(DOA)估计是阵列信号处理研究的重要内容之一多重信号分类（MUSIC）算法是一种经典有效的DOA方法。概述了阵列信号处理研究的主要内容，详细介绍了多重信号分类算法的原理，并用MATLAB进行仿真实现该算法，研究了MUSIC算法的性质。关键词：阵列信号处理多重信号分类算法仿真1、实验目的（1）掌握MUSIC算法进行DOA估计的基本原理（2）了解影响DOA估计性能的因素2、实验内容假设阵元个数M=16，信号方向分别为0度、3度和-35度，噪声为高斯噪声仿真：（1）快拍次数一定，分辨力与SNR的关系（2）快拍次数一定，分辨力与孔径的关系；（3）SNR一定，分辨力与快拍的关系。3、实验原理窄带远场信号的DOA数学模型为阵列数据的协方差矩阵为由于信号与噪声相互独立，数据协方差矩阵可以分解为与信号、噪声相关的两部分，其中是信号的协方差矩阵，是信号部分。对进行特征分解有式中，是由大特征值对应的特征矢量张成的子空间也即信号子空间，而是由小特征值对应的特征矢量张成的子空间也即噪声子空间。根据前面所述的性质2可知，在理想的条件下数据空间中的信号子空间与噪声子空间是相互正交的，即信号空间中的导向矢量也与噪声子空间正交经典的MUSIC算法正是基于上述这个性质提出的，但考虑到实际接收数据矩阵是有限快拍数的，即数据协方差矩阵的最大似然估计为对进行特征分解可以计算得到噪声子空间特征矢量矩阵。由于噪声的存在，与并不能完全的正交。因此，实际上求DOA是以最小优化搜索实现的，即所以，MUSIC算法的谱估计公式为下面给出MUSIC算法的计算步骤：由阵列的接收数据得到数据协方差矩阵；对进行特征分解；由的特征值进行信号源数的判断；确定信号子空间与噪声子空间；根据信号参数的范围进行谱峰搜索；找出极大值点对应的角度就是信号入射方向。对于上述的MUSIC算法，还应该要注意以下几点：对于非理想情况下得到协方差矩阵的特征值满足下式：所以判断信号源数需要用到有关信号源估计的方法来进行信号源数的确定。线阵的信号参数搜索范围为，而面阵的搜索范围为。另外，一个确定阵列的导向矢量由阵元的位置唯一的确定；MUSIC算法的一种归一化形式，即在实际应用中，对于一维导向矢量有下式成立：4、实验仿真及结果分析4.1MUSIC随信噪比变化情况 信号数共3个，分别位于2°，30°和31°，它们功率相同，阵元数共16个。这里快拍数取500，以保证快拍数足够，可以较好得得出统计特性，以去除采样不足带来的分辨力低下。信噪比分别取为SNR=0,5,10,15,20dB，得到图2。从图中可以看到，随信噪比增大，MUSIC谱图旁瓣变低，主瓣变得更加尖锐，分辨力变得更好。Ri=(akj1*exp(1j*w1*i*Ts)+as2*akj2*exp(1j*w2*i*Ts)+as3*akj3*exp(1j*w3*i*Ts)+ri)*(akj1*exp(1j*w1*i*Ts)+as2*akj2*exp(1j*w2*i*Ts)+as3*akj3*exp(1j*w3*i*Ts)+ri)';R=R+Ri;endR=R/N;[V,D]=eigs(R,NN,'sm');fai=(-90:0.005:90)/180*pi;ps=zeros(1,length(fai));disp('¾­µämusicËã·¨¼ÆËãÀ´²¨½Ç¶È/¶È');fori=1:length(fai)%ÕóÁÐÁ÷ÐÎɨÃèki=k*w3/w1*sin(fai(i));aki=exp(-1j*(ki*m*d)).';ps(i)=real(1/(aki'*(V*V')*aki));endfori=2:length(fai)-1%Ñ°ÕÒ¼«´óÖµµã£¬¼´À´²¨·½Ïòifps(i)>ps(i+1)&&ps(i)>ps(i-1)&&ps(i)>max(ps)/20disp(fai(i)/pi*180);%À´²¨½Ç¶Èendendfigure;plot(fai/pi*180,10*log10(ps/max(ps)));%%%Çó¸ùËã·¨%%%%%symsz;Pz1=sym(zeros(1,M));Pz2=sym(zeros(1,M)).';fori=1:MPz1(i)=z^(1-i);Pz2(i)=z^(i-1);endtsym=solve(Pz1*(V*V')*Pz2);tdouble=double(tsym);r=abs(tdouble);rnum=length(r);%numsolve=zeros(1,K);%fori=1:K%Ñ°ÕÒºÏÊʵĸù%temp1=20;%temp2=1;%forj=1:length(r)%ifabs(r(j)-1)<temp1%temp1=abs(r(j)-1);%temp2=j;%end%end%numsolve(i)=tdouble(temp2);%r(temp2)=10;%endtsolve=zeros(1,rnum);j=0;fori=1:rnumif(abs(r(i)-1)<0.1)j=j+1;tsolve(j)=tdouble(i);endendnumsolve=tsolve(1:j);signal_position=asin(-wavelength/(2*pi*d)*angle(numsolve))/pi*180;real_position=signal_position(1);fori=2:length(signal_position)%È¥µôÖظùflag=0;forj=1:i-1ifabs(signal_position(i)-signal_position(j))<0.5flag=1;endendiffla