music_2.txt

这里有三个music算法用于来波方向估计的算法

clc
clear
close
signalnum=2;
SNR1=10;
SNR2=20;
theta1=10;
theta2=40;
SA1=sqrt(2*10^(SNR1/10));
SA2=sqrt(2*10^(SNR2/10));
theta1=theta1*pi/180;
theta2=theta2*pi/180;
jay=sqrt(-1);
c=1500;lamda=1;
M = 10;fa=c/lamda;
d=lamda/2;
len=2;
fs=3000;
f1=500;
f2=600;
t=[0:1/fs:len];
ss1=sin(2*pi*f1*t);
s1=SA1*ss1;
ss2=cos(2*pi*f2*t);
s2=SA2*ss2;
%阵列流型
tau1=([0:d:(M-1)*d]*sin(theta1))/c;
tau2=([0:d:(M-1)*d]*sin(theta2))/c;
a1=exp(-jay*2*pi*fa*tau1).';
a2=exp(-jay*2*pi*fa*tau2).';
SRr=[a1 a2]*[s1;s2]+randn(M,length(s1));
Rr=SRr*SRr'/length(t);
flag0=1;
if flag0==0
    [u v]=svd(Rr);
    for i=1:M-signalnum
        ee1=u(:,i+2);
        e1(:,i)=u(:,i+2)./min(real(ee1));
    end
    Enr=e1;
else
    [u v]=eig(Rr);
    [Order Index]=sort(diag(v));
    for i=1:M-signalnum
        e1(:,i)=u(:,Index(i));
    end
    Enr=e1;
end
theta=[-90:0.1:90];
for j=1:length(theta)
    thetaa=theta(j)*pi/180;
    for i=1:M
        ae(i)=exp(-jay*2*pi*fa*(i-1)*d*sin(thetaa)/c);
    end
    a=ae.';
    Pr(j)=1/(a'*Enr*Enr'*a);
end
figure
plot(theta,10*log10(abs(Pr)/max(Pr)),'-r')
grid on
title('music算法估计声源方位谱')