mus_md.m

适用于水下信号分析的矩阵分解源代码

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   空间谱估计中 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
close all
clc
%%%%%%%%%%%%%%%%%%%%%%%  generate signal  %%%%%%%%%%%%%%%%%
f0=5000;                      %入射信号频率
fs=1000;                      %采样频率
ts=1/fs;                      %采样间隔
M=8;                           %阵元数
L=100;                        %快拍数
SNR=5;                        %信噪比  %%%%%%                        
the0=[-10 0 10];
Q=length(the0);                           %信号源数

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.

u1=10;
s1=K*exp(j*(2*pi*u1*(t-1)+pi/6));
s2=K*exp(j*(2*pi*u1*(t-1)+pi/3));     
s3=K*exp(j*(2*pi*u1*(t-1)+pi/2));   %相干信源

ss=[s1;s2];%;s3];                  %相干信号源矩阵 基带信号
%%%%%%%%%%%%%%%%%%%%% 生成导向矢量矩阵  %%%%%%%%%%%%%%%%%%%
m=0:M-1;
A=exp(-j*2*pi*m'*d/lamda*sin(the0*pi/180));
%%%%%%%%%%%%%%%%%%%%% 数据协方差矩阵%%%%%%%%%%%       
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
X=A*ss+Nn;               %阵列接收的相干信号， 情况2                        
R=X*X'/L;               %阵列采样数据协方差矩阵
%%%%%%%%%%%%重构矩阵%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mm=4;l0=3;
Rm=zeros(mm,(l0+1)*M);
Rmt=Rm;
J1=fliplr(eye(mm));
J2=fliplr(eye(M));
for l=0:l0
    Rm(:,l*M+1:(l+1)*M)=R(l+1:l+mm,:);         % MD算法的矩阵
    Rmt(:,l*M+1:(l+1)*M)=J1*conj(R(l+1:l+mm,:))*J2;
end
Rm1=[Rm Rmt];   % MMD算法的矩阵

p=5;
Rf=zeros(mm);
Rb=zeros(mm); 
for k=1:p
    Rk=R(k:mm+k-1,k:mm+k-1);
    Rf=Rf+Rk;
    Rbk=J1*conj(Rk)*J1;
    Rb=Rb+Rbk;
end
Rf=Rf/p;
Rb=Rb/p;
Rm2=[Rf Rb];   % SMD算法的矩阵

%%%%%%%%%%%%%%基于矩阵分解的解相干%%%%%%%%%%%%%%%
[U S V]=svd(Rm);
[U1 S1 V1]=svd(Rm1);
[U2 S2 V2]=svd(Rm2);
Un=U(:,Q+1:mm);
Un1=U1(:,Q+1:mm);
Un2=U2(:,Q+1:mm);

u=-90:0.25:90;                           %搜索范围
for ii=1:length(u)
    Au=exp(-j*2*pi*m'*d/lamda*sin(u(ii)*pi/180));   %搜索范围的导向矢量阵
    P_md(ii)=1/abs(Au(1:mm)'*Un*Un'*Au(1:mm));
    P_mmd(ii)=1/abs(Au(1:mm)'*Un1*Un1'*Au(1:mm));    
    P_smd(ii)=1/abs(Au(1:mm)'*Un2*Un2'*Au(1:mm)); 
end
plot(u,10*log10(P_md),'s-r',u,10*log10(P_mmd),'x-g',u,10*log10(P_smd),'--b')%,u,10*log10(P),'-k');
legend('MD','MMD','SMD'),xlabel('入射角度'),ylabel('空间方位谱(dB)'),
title('矩阵分解算法'),grid on;

%Q=2时，(1)[0 5]  （2）[0 10]
%Q=3时，(1)[-5 0 5]  （2）[-10 0 10] (3)[-20 0 20]