yuanmusic.m

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

%%%6 sensor. 4-8GHz%%%
clc
clear
tic
ccc1=cputime;
sm=7;%阵元数
snr=20;
nd=100;%快拍数
md=181; %?
snr=10.^(snr/20);
f=4*10^3*10^6;  
wl=3*10^8/f;
r=wl/2
diz=[0                0                0                0                0                0                0               ];
dix=[r*cos(2*pi/sm*0) r*cos(2*pi/sm*1) r*cos(2*pi/sm*2) r*cos(2*pi/sm*3) r*cos(2*pi/sm*4) r*cos(2*pi/sm*5) r*cos(2*pi/sm*6)]; 
diy=[r*sin(2*pi/sm*0) r*sin(2*pi/sm*1) r*sin(2*pi/sm*2) r*sin(2*pi/sm*3) r*sin(2*pi/sm*4) r*sin(2*pi/sm*5) r*sin(2*pi/sm*6)];
%参数都为距离不是角度 单位是米

fw=[270 200];     %fangwei%
fy=[60  50];      %fuyang%


rky=fw*pi/180;   %fangwei%
rkx=fy*pi/180;   %fuyang%
ll=length(rkx);
   p31=rand(1,1);
  for m=1:sm;
    for n=1:ll;
    %yxl  a(m,n)=exp(-j*2*pi/wl*((dix(m)-dix(1))*sin(rkx(n))*cos(rky(n))+(diy(m)-diy(1))*sin(rkx(n))*sin(rky(n))+diz(m)*cos(rkx(n))));
       a(m,n)=exp(-j*2*pi/wl*((dix(m))*sin(rkx(n))*cos(rky(n))+(diy(m))*sin(rkx(n))*sin(rky(n))+diz(m)*cos(rkx(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
en=z(:,d+1:sm);
en=en/norm(en);
en=en*en';


x=linspace(0,360,md);
rka=x*pi/720;  %fuyang
y=linspace(0,360,md);
rkb=y*pi/180;  %fangwei


for k=1:sm;
   for m=1:md;
      for n=1:md;
%yxl       b(k,m,n)=exp(-j*2*pi/wl1*((dix(k)-dix(1))*sin(rka(m))*cos(rkb(n))+(diy(k)-diy(1))*sin(rka(m))*sin(rkb(n))+diz(k)*cos(rka(m))));
        b(k,m,n)=exp(-j*2*pi/wl*((dix(k))*sin(rka(m))*cos(rkb(n))+(diy(k))*sin(rka(m))*sin(rkb(n))+diz(k)*cos(rka(m))));    %b的维数为6*181*181;
   end
   end
end
for i=1:md;
    for j=1:md;
      pmu(i,j)=(b(:,i,j))'*en*b(:,i,j);
   end
 end
pmu=1./pmu;
pmu=abs(pmu);
pmu=10*log(pmu);
figure(1)
plot3(x,y/4,pmu);
% mesh(x,y/4,pmu);
grid on;
% hidden off;
 meshc(x,y/4,pmu);
 
xlabel('方位角度（0\circ\sim 360\circ）')
ylabel('俯仰角度（0\circ\sim 90\circ）')

%figure(2)%
%contour(x,y/4,pmu);
%title('6阵元下的MUSIC','fontsize',12)
%xlabel('方位角度（0\circ\sim 360\circ）','fontsize',12)
%ylabel('俯仰角度（0\circ\sim 90\circ）','fontsize',12)
%grid on;
d=length(fy);
qx=0;
qy=0;
 [ii,yy]=max(pmu);
 [iii,yyy]=max(ii); 
 cc1=yy(yyy);
if pmu(1,1)>=pmu(1,2);
    if pmu(1,1)>=pmu(2,1);
        qx=qx+1;
        qy=qy+1;
%         pm(qx,qy)=pmu(i,j);
        pm(qx,qy)=pmu(1,1);
        pda(qx)=1;
        pdb(qy)=1;
    end
end  %

for k=2:md-1;
    if (pmu(1,k-1)<=pmu(1,k));
        if (pmu(1,k)>=pmu(1,k+1));
            if (pmu(1,k)>=pmu(2,k));
            qx=qx+1;
            qy=qy+1;
            pm(qx,qy)=pmu(1,k);
            pda(qx)=1;
            pdb(qy)=k;
            end
        end
    end
 end
    
for i=2:md-1;
    if (pmu(i,1)>=pmu(i-1,1));
        if (pmu(i,1)>=pmu((i+1),1));
            if (pmu(i,1)>=pmu(i,2));
                qx=qx+1;
                qy=qy+1;
                pm(qx,qy)=pmu(i,1);
                pda(qx)=i;
                pdb(qy)=1;
            end
        end
    end
end
for i=2:md-1;
  for j=2:md-1;
     if(pmu(i-1,j)<=pmu(i,j));
        if (pmu(i,j)>=pmu(i+1,j));
           if(pmu(i,j-1)<=pmu(i,j));
              if (pmu(i,j)>=pmu(i,j+1));
                 qx=qx+1;
                 qy=qy+1;
                 pm(qx,qy)=pmu(i,j);
                 pda(qx)=i;
                 pdb(qy)=j;
              end
           end
       end 
     end
 end 
end
aa=qx;
bb=qy;
pm=diag(pm);
[qq]=find(pm>-500);
pda=pda(qq);
pdb=pdb(qq);
pm1=pm(qq);
[pm2,k]=sort(pm1);      %将pm1元素按照升序排列付给pm2；k对应的是pm1的元素对应在pm2列向量中的位置
pm3=(fliplr(pm2'))';    %fliplr:左右方向反转数组；此语句的功能是将特征值按照降幂排列 
pm4=pm3(1:d);
pda=pda(:,k) ;          %pda按照升幂值的位置对应
pda=fliplr(pda);
pda=pda(1:(d));
pdb=pdb( :,k) ;           %pda按照升幂值的位置对应
pdb=fliplr(pdb)  
pdb=pdb(1:(d));
pm4;
pda=(pda-1)/2;
pdb=(pdb-1)*2;
disp('俯仰角度')
disp(pda)
disp('方位角度')
disp(pdb)
[fy,b]=sort(pda);
fw=pdb(:,b);
ccc2=cputime;
cc=ccc2-ccc1
toc