musicyuanzhen1.m

矢量阵指向性程序

clear;
f0=1000;
fs=80*f0;
Ns=2048;
N1=8;
c=1500;
r1=1*c/f0;
r2=0.5*r1;
c=1500;
fai=0;
SNR=10;
source=1;
w=2*pi*f0;
ph1=-180;                   
ph2=180;                   
phn=1;                  
ax=(ph2-ph1)/phn+1;

for i=1:N1;
       a=2;
    if mod(i,a)~=0
       for k=1:Ns;
        p(i,k)=exp(j*w*(k-1)/fs-j*w*r1*cos(fai*pi/180-(i-1)*2*pi/N1)/c);
       end
    else
       for k=1:Ns;
        p(i,k)=exp(j*w*(k-1)/fs-j*w*r2*cos(fai*pi/180-(i-1)*2*pi/N1)/c);
       end
    end
end
for k=1:N1
    n(k,:)=randn(1,length(p(1,:)));          
end
for k=1:N1;
    A(k,:)=sqrt(10^(-SNR/10)*sum(abs(p(k,:)).^2)/sum(abs(n(k,:)).^2));
end
for k=1:N1;
    N(k,:)= A(k,:)*n(k,:);
end
  
for k=1:N1
    P2(k,:)=p(k,:)+N(k,:);    
end 
R=1/Ns*P2*P2'; 
R1=inv(R);
[em,zm]=eig(R);                            
[zm1,pos1]=max(zm);                         
for l=1:source;
    [zm2,pos2]=max(zm1);                    
   zm1(:,pos2)=[];                          
   em(:,pos2)=[];                           
end
for i=1:ax;     
    num=ph1+(i-1)*phn;
    fain=num*pi/180;
   for k=1:N1;
         a=2;
      if mod(k,a)~=0
       x(k,i)=exp(-j*w*r1*cos(fain-(k-1)*2*pi/N1)/c); 
      else
      x(k,i)=exp(-j*w*r2*cos(fain-(k-1)*2*pi/N1)/c);   
      end
   end
      v=x(:,i);
   pm=v'*em;
   pmusic(i)=1/(pm*pm');                    
end
p1=max(pmusic);
D1=pmusic/p1;
zz1=20*log10(D1);
figure;
i=1:ax;
plot(ph1+(i-1)*phn,zz1,'*-');
grid on 
title('信噪比为10dB，MUSIC')
xlabel('角度')
ylabel('dB')
% for i=1:ax;     
%     num=ph1+(i-1)*phn;
%     fain=num*pi/180;
%      Xp1=0;
%     Xv1=0;
%     for n=1:N1;
%         a=2;
%       if mod(n,a)~=0
%         Xp01(n,1,i)=p1(n,1)*exp(j*w*r1*cos(fain-(n-1)*2*pi/N1)/c); 
%         Xp1=Xp1+Xp01(n,1,i);
%         Xv01(n,1,i)=Xp01(n,1,i)*(cos(fai*pi/180)*cos(fain)+sin(fai*pi/180)*sin(fain));
%         Xv1=Xv1+Xv01(n,1,i);
%       else
%         Xp02(n,1,i)=p2(n,1)*exp(j*w*r2*cos(fain-(n-1)*2*pi/N1)/c); 
%         Xp1=Xp1+Xp02(n,1,i);
%         Xv02(n,1,i)=Xp02(n,1,i)*(cos(fai*pi/180)*cos(fain)+sin(fai*pi/180)*sin(fain));
%         Xv1=Xv1+Xv02(n,1,i);
%       end
%     end
%     Rp(i)=abs(Xp1*Xv1);
% %       Rp(i)=abs((Xp1+Xv1)*Xv1);
% %         Rp(i)=abs((Xp1+Xv1).^2);
% %           Rp(i)=abs((Xp1+2*Xv1)*Xv1);
% %             Rp(i)=abs((Xp1+2*Xv1).^2);
% %               Rp(i)=abs((Xp1+2*Xv1)*(Xp1+Xv1));
% end 
% i=1:ax;
% maxRp=max(Rp);
% Rpp=20*log10(Rp(i)/maxRp);
% figure
% plot(ph1+(i-1)*phn,Rpp);
% grid on
% hold on
% axis([-180 180 -100 0])
% xlabel('角度')
% ylabel('dB')