mmusic_new.m

修正MUSIC算法实现相干信号源的DOA估计算法

clear;
clc;

%%%%%%%%%%%%input%%%%%%%%%%%%
arg=[1.11,-1.17];d=1/2;
N=20;K=100;c=3*10^8;
ps = 0.9607; c = 3*10^8;f = 6*10^9;
lamda = c/f;
vv =  1.2190e-004;
ff = 2*pi*vv*10^3/lamda;
snr = 5;
%%%%%%%%%%%%input%%%%%%%%%%%%
xx=eye(N);
Jm=rot90(xx);%转移向量
fai=sin(arg*pi/180);
II=eye(2);
R=zeros(N,N);
w0=rand(1,K);
tic;
for k=1:K
    for n=1:N
        A(n,1)=exp(-i*2*pi*d*(n-1)*fai(1));
        A(n,2)=exp(-i*2*pi*d*(n-1)*fai(2));
    end
    s0=[1;ps*exp(i*ff)];
    s=10^(snr/20)*s0;
    y=A*s+randn(N,1)+i*randn(N,1);
    R=R+y*y';
end
R=R/K;
R1=R+Jm*conj(R)*Jm;
[U,S,V]=svd(R1);
for p=1:2
    S1(p)=S(p,p);
    UA(:,p)=U(:,p);
    VA(p,:)=V(p,:);
end
for m=1:N-2
    S2(m)=S(m+2,m+2);
    UO(:,m)=U(:,m+2);
    VO(m,:)=V(m+2,:);
end
SS=diag(S1);
SN=diag(S2);
sigma=trace(SN)/(N-2);
B0=SS-sigma*II;
for p=1:2
    B1(p)=1/B0(p,p);
end
B=diag(B1);
RA=UA*B*UA';
PO=UO*UO';

N1=-5;N2=5;t=0.005;
for nn=0:round((N2-N1)/t)
    fei=2*pi*d*sin((N1+nn*t)*pi/180);
    n=1:N;
    e0=exp(-i*fei*(n-1));
    e=e0.';
    num0=e'*e;
    num=e'*RA*e;
    den=e'*PO*e;
    fMUSIC=1/den;
    fM(nn+1)=0;
    fM(nn+1)=10*log10(fMUSIC);
end
toc;
pp=N1:t:N2;
plot(pp,fM)
title('MUSIC Spectrum');
xlabel('\theta/deg'); 
ylabel('Spectrum/dB');
grid on
% print -dtiff MMUISIC