xianmusic.m

MUSIC（Multiple Signal Classification）算法线阵仿真程式

clc
clear
for i=1:10;
sm=6;                              %阵元数
snr=10;                            %信噪比dB
nd=120;                            %快拍数
snr=10.^(snr/20);                  %信噪比
f=4*10^3*10^6;                     %信号频率  
wl=3*10^8/f;                       %波长
d=0.5*wl;                          %阵元间距
dix=[0:sm-1]*d;                    %均匀线阵
doa=[55  50  -30 -50]*pi/180;      %波达方向
ll=length(doa);                    %信号数

for m=1:sm;
    for n=1:ll;
     a(m,n)=exp(-j*2*pi/wl*dix(m)*sin(doa(n)));
    end
end

for m=1:ll;
   for n=1:nd;
      p1=rand(1,1);
      p2=rand(1,1);      
      sr(m,n)=sqrt(-2*snr*snr*log(p1))*cos(2*pi*p2);
      si(m,n)=sqrt(-2*snr*snr*log(p1))*sin(2*pi*p2);
      s(m,n)=sr(m,n)+j*si(m,n);
   end 
end

for m=1:sm;
    for n=1:nd;
       p3=rand(1,1);
       p4=rand(1,1);
       wr(m,n)=sqrt(-2*log(p3))*cos(2*pi*p4);
       wi(m,n)=sqrt(-2*log(p3))*sin(2*pi*p4);
       w(m,n)=wr(m,n)+j*wi(m,n);
    end
end

x=(a*s+w);
r=x*x'/nd;

[z,d]=eig(r,'nobalance');
d1=diag(d) ;                      %d1为列向量，各个元素为d的主对角线元素值，diag为建立或提取对角阵
[lambtal,k]=sort(d1);             %将d1元素按照升序排列付给lambtal；k对应的是lambtal的元素对应在d1列向量中的位置
lambtal=(fliplr(lambtal'))' ;     %fliplr:左右方向反转数组；此语句的功能是将特征值按照降幂排列
z=z(:,k) ;                        %特征向量按照升幂特征值的位置对应
z=fliplr(z);                      %特征向量按照降幂特征值的位置对应
                                  %利用MDL准则判断信号源数目
                                  
for p=1:sm;
    c1=1;
    c2=0;
    for i=sm-p+1:sm;
        c1=c1*lambtal(i);
        c2=c2+lambtal(i);
    end
    c1=c1^(1.0/p);
    c2=c2/p;
    lr(p)=(c1/c2)^(nd*p/2);
    mdl(p)=(sm-p)*(sm+p+1)*log10(nd)/4-log10(lr(p));
end
[lambta_min ,nx]=min(mdl)
d=sm-nx
                                   %利用MDL准则判断信号源数目
en=z(:,d+1:sm);                    
en=en/norm(en);                    %信号子空间

threat=[-90:0.1:90]*pi/180;
num=length(threat);

for k=1:sm;
   for m=1:num;
    b(k,m)=exp(-j*2*pi/wl*dix(k)*sin(threat(m))); 
   end
end

for i=1:num;
    pmu(i)=(b(:,i))'*en*en'*b(:,i);
end

h=[-90:0.1:90];
pmu=1./pmu;
pmu=abs(pmu);
pmu=10*log(pmu);
figure(1)
plot(h,pmu);
xlabel('入射角/（度）','fontsize',10)
ylabel('Pmusic','fontsize',10)




n=0;
for a=2:(num-1);
   if pmu(a)>pmu(a-1);
       if pmu(a)>pmu(a+1);
            n=n+1;
            pm(n)=a;
            k(n)=pmu(a);
        end
    end
end
jiaodu=(pm-1)/10-90;
disp(jiaodu)
end
grid on;