music.m

经典MUSIC算法包括相干和非相干信号

clear
close all
clc
f0=5000;                      %入射信号频率
fs=1000;                      %采样频率
ts=1/fs;                      %采样间隔
M=8;                           %阵元数
Q=2;                           %信号源数
L=100;                        %快拍数
SNR=5;                        %信噪比  %%%%%%                        
the0=[0 30];
theta1=the0(1)*pi/180;
theta2=the0(2)*pi/180;                   %z1,z2为信号入射方向
c=1500;
lamda=c/f0;
d=lamda/2;                  %c为水下声速，d为阵元间距
t=(1:L)*ts;                 % 采样时刻序列
K=sqrt(2*10^(SNR/10));      %信号的幅度,Pn=1，the load impedance is 1 ohm.
u0=5;  % s0=exp(j*2*pi*0.5*u0*t)
u1=10;
fai1=unidrnd(360,1,L)*pi/180;    %fai1 and fai2 are random numbers from discrete uniform distribution 0-360 degree with 1 row and L columns
s1=K*exp(j*(2*pi*0.5*u0*(t-1)+fai1));
fai2=unidrnd(360,1,L)*pi/180;
s2=K*exp(j*(2*pi*0.5*u0*(t-1)+fai2));
s3=K*exp(j*(2*(t-1)*pi*0.5*u1+pi/2));
s4=K*exp(j*(2*(t-1)*pi*0.5*u1+pi/3));
ss1=[s1;s2];                         %独立信号源矩阵  基带信号
ss2=[s3;s4];                        %相干信号源矩阵
tao1=d*sin(theta1)/c;   %阵元间的时延
tao2=d*sin(theta2)/c;
m=0:M-1;
a_theta1=exp(-j*2*pi*f0*(m)'*tao1);%阵列方向向量
a_theta2=exp(-j*2*pi*f0*(m)'*tao2);
A=[a_theta1 a_theta2];  %阵列方向矩阵，阵列流形    
Nn = wgn(M,L,0,'complex');%  M columns vector of length 100 containing complex white Gaussian noise, each component of which has a noise power of 0 dBW
X1=A*ss1+Nn;                            %阵列接收相互独立的信号，情况1
X2=A*ss2+Nn;                            %阵列接收的相干信号， 情况2
R1=X1*X1'/L;                              %阵列采样数据协方差矩阵
R2=X2*X2'/L;
[p1,r1]=eig(R1);                           %特征分解
[p2,r2]=eig(R2);                        %特征分解
mr1=fliplr(r1);            %reshaping eigenvalues diagnal matrix r1 to mr1
mr1=flipud(mr1);           %将r1的特征值位置交换，默认的是从小到大排，使之按从大到小
mp1=fliplr(p1);            %reshaping eigenvectors corresponding eigenvalues p1 to mp1改变相应的特征向量位置
Us1(:,[1,2])=mp1(:,[1,2]);                 %对R特征分解并取出最大特征值及其对应的特征向量为信号子空间Us
Un1(:,[1,2,3,4,5,6])=mp1(:,[3,4,5,6,7,8]); %对R特征分解并取出小特征值及其对应的特征向量为噪声子空间Un
mr2=fliplr(r2);            %reshaping eigenvalues diagnal matrix r2 to mr2
mr2=flipud(mr2);
mp2=fliplr(p2);            %reshaping eigenvectors corresponding eigenvalues p2 to mp2
Us2(:,[1,2])=mp2(:,[1,2]);                 %对R特征分解并取出最大特征值及其对应的特征向量为信号子空间Us
Un2(:,[1,2,3,4,5,6])=mp2(:,[3,4,5,6,7,8]); %对R特征分解并取出小特征值及其对应的特征向量为噪声子空间Un
u0=[1 0 0 0 0 0 0 0]';
theta=-90:0.1:90;
for ii=1:length(theta)
    a_theta=exp(-j*2*pi*d/lamda*m'*sin(theta(ii)/180*pi));%阵列方向向量
   P_music1(ii)=abs((a_theta'*a_theta)/(a_theta'*Un1*Un1'*a_theta)); 
   P_music2(ii)=abs((a_theta'*a_theta)/(a_theta'*Un2*Un2'*a_theta));
end
plot(theta,10*log10(P_music1),'r',theta,10*log10(P_music2),'b');
xlabel('入射角度(deg)'),ylabel('空间方位谱(dB)'),title('经典MUSIC算法'),grid on
legend('独立信号源','相干信号源');