mus_svd.m

适用于水下信号分析的矢量奇异值MUSIC算法源代码

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   空间谱估计中 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
close all
clc
%%%%%%%%%%%%%%%%%%%%%%%  generate signal  %%%%%%%%%%%%%%%%%
f0=5000;                      %入射信号频率
fs=1000;                      %采样频率
ts=1/fs;                      %采样间隔
M=8;                           %阵元数
L=100;                        %快拍数
SNR=20;                        %信噪比  %%%%%%                        
the0=[-10 20];
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;               %阵列采样数据协方差矩阵
[U S V]=svd(R);
e=U(:,1);       %最大特征值对应的特征矢量
x=X*X(1,:)'/L;  %数据矢量
%%%%%%%%%%%%%%%%%%%  构造矩阵Y  %%%%%%%%%%%%%%%%%%%
mm=4;p=5;
Y=zeros(mm,p);
for k=1:p
    Y1(:,k)=e(k:k+mm-1);
    Y2(:,k)=x(k:k+mm-1);
end
%%%%%%%%%%%% 解相干后的噪声子空间 %%%%%%%%%%%%%%%%
[U1 S1 V1]=svd(Y1);
Un1=U1(:,Q+1:mm);

[U2 S2 V2]=svd(Y2);
Un2=U2(:,Q+1:mm);
%%%%%%%%%%%  基于矢量奇异值算法的MUSIC %%%%%%%%%%%%%%%%%%%%
u=-90:0.5:90;                           %搜索范围
for ii=1:length(u)
    Au=exp(-j*2*pi*m'*d/lamda*sin(u(ii)*pi/180));   %搜索范围的导向矢量阵
    P_esvd(ii)=1/abs(Au(1:mm)'*Un1*Un1'*Au(1:mm));
    P_dsvd(ii)=1/abs(Au(1:mm)'*Un2*Un2'*Au(1:mm));    
end
plot(u,10*log10(P_esvd),'s-r',u,10*log10(P_dsvd),'x-g')%,u,10*log10(Pfb),'--b',u,10*log10(P),'-k');
legend('ESVD','DSVD,'),xlabel('入射角度'),ylabel('空间方位谱(dB)'),
title('矢量奇异值算法'),grid on;

%mm=4,p=5
%Q=3时， (1)[-5 0 5]无法完全分辨 ， 直至(2)[-10 0 10]才能有误差的分辨，  至(3)[-13 0 13]能无误差地分辨
%Q=2时，(1)[0 5]刚好能分辨 ，(2)[-10 20]