svdmusic.m

来自「矢量阵指向性程序」· M 代码 · 共 73 行

%*****************  ******************
 format long     
 clc;
clear;
close all
c=1500;
f0=5000;
d=0.5*(c/f0);
Fs=5*f0;
lth=256;               % length of FFT
sita1=-5;
sita2=5;
sita3=15;
SNR=20;
count=16;
m=4;          %子阵阵元数
p=count-m+1;    %子阵数
ph1=-90;                   
ph2=90; 
source=3;
phn=1;                  
ax=(ph2-ph1)/phn+1;     
t01=d/c*sin(sita1*pi/180);
t02=d/c*sin(sita2*pi/180);
t03=d/c*sin(sita3*pi/180);
% ********产生信号********
for b=1:10;
for k=1:count
    for i=1:lth
        S(k,i)=exp(-j*(2*pi*f0*(i-1)/Fs-2*(k-1)*pi*f0*t01))+exp(-j*(2*pi*f0*(i-1)/Fs-2*(k-1)*pi*f0*t02))+exp(-j*(2*pi*f0*(i-1)/Fs-2*(k-1)*pi*f0*t03));     
    end  
     n(k,:)=randn(1,length(S(k,:))); 
     A(k,:)=sqrt(10^(SNR/10)/sum(abs(S(k,:)).^2)*lth);
     P2(k,:)=A(k,:).*S(k,:)+n(k,:); 
end 
X1=(1/lth)*P2*P2(1,:)';
F=zeros(m,p);
for i=1:p; 
    f1=[zeros(m,i-1),eye(m),zeros(m,count-m-i+1)];
    f2=[zeros(1,i-1),eye(1),zeros(1,count-m-i+1)];
    F1= f1*X1*f2;
    F=F+F1;
end
R=F/p;
[em,zm,cm]=svd(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;
    fai=num*pi/180;
   for k=1:m;
       x(k,i)=exp(j*(k-1)*2*pi*d*f0*sin(fai)/c); 
   end
      v=x(:,i);
   pm=v'*em;
   pmusic(b,i)=1/(pm*pm');                    
end
end
pmusic=sum(pmusic)/20;
p1=max(pmusic);
D1=pmusic/p1;
zz1=10*log10(D1);
figure;
i=1:ax;
plot(ph1+(i-1)*phn,zz1);
legend('SVDMUSIC',4);
title('信噪比为20dB，相干信号源')
grid on