music_smooth.m

适用于水下信号分析的空间平滑MUSIC算法源代码

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   空间谱估计中 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
close all
clc
%%%%%%%%%%%%%%%%%%%%%%%  generate signal  %%%%%%%%%%%%%%%%%
f0=5000;                      %入射信号频率
fs=1000;                      %采样频率
ts=1/fs;                      %采样间隔
M=8;                           %阵元数
L=100;                        %快拍数
SNR=5;                        %信噪比  %%%%%%                        
the0=[-40 -20 10 20 40];
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)+2*pi*rand(1,L)));
% s2=K*exp(j*(2*pi*u1*(t-1)+2*pi*rand(1,L)));  %独立信源
s3=K*exp(j*(2*pi*u1*(t-1)+pi/6));
s4=K*exp(j*(2*pi*u1*(t-1)+pi/3));    %相干信源 
s5=K*exp(j*(2*pi*u1*(t-1)+pi/2));
s6=K*exp(j*(2*pi*u1*(t-1)+2*pi/3));
s7=K*exp(j*(2*pi*u1*(t-1)+5*pi/6));

% ss1=[s1;s2];                         %独立信号源矩阵  基带信号
ss2=[s3;s4;s5;s6;s7];                        %相干信号源矩阵
%%%%%%%%%%%%%%%%%%%%% 生成导向矢量矩阵  %%%%%%%%%%%%%%%%%%%
m=0:M-1;
A=exp(-j*2*pi*m'*d/lamda*sin(the0*pi/180));
%%%%%%%%%%%%%%%%%%%%% 产生零均值方差为1的 M*L高斯白噪声%%%%%%%%%%%       
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;
R=R2;

    %%%%%%%%%%%%空间平滑Music算法%%%%%%%%%%%%%%%%%%%%
    p=3;mm=6;%将均匀线阵分成相互交错的5个子阵，每个子阵阵元数为4
    Rf=zeros(mm);
    Rb=zeros(mm);
    J=fliplr(eye(mm));
    for k=1:p
        Rk=R(k:mm+k-1,k:mm+k-1);
        Rf=Rf+Rk;
        Rbk=J*conj(Rk)*J;
        Rb=Rb+Rbk;
    end
    Rf=Rf/p;
    Rb=Rb/p;
    Rw=(Rf+Rb)/2;
    [V1 D1]=eig(Rf);
    [V2 D2]=eig(Rb);
    [V3 D3]=eig(Rw);
    Un1=V1(:,1:mm-Q);
    Un2=V2(:,1:mm-Q);
    Un3=V3(:,1:mm-Q);
    %%%%%%%%%%%%经典Music算法%%%%%%%%%%%%%%%%%%%%
    [V D]=eig(R);
    Un=V(:,1:M-Q);
%……………………………谱峰搜索…………………………………………………………    
u=-90:0.5:90;                           %搜索范围
for ii=1:length(u)
    Au=exp(-j*2*pi*m'*d/lamda*sin(u(ii)*pi/180));   %搜索范围的导向矢量阵
    Pf(ii)=1/abs(Au(1:mm)'*Un1*Un1'*Au(1:mm));
    Pb(ii)=1/abs(Au(1:mm)'*Un2*Un2'*Au(1:mm));
    Pfb(ii)=1/abs(Au(1:mm)'*Un3*Un3'*Au(1:mm));
    P(ii)=1/abs(Au'*Un*Un'*Au);
end
plot(u,10*log10(Pf),'s-r',u,10*log10(Pb),'x-g',u,10*log10(Pfb),'--b',u,10*log10(P),'-k');
legend('前向空间平滑MUSIC','后向空间平滑MUSIC','双向空间平滑MUSIC','经典MUSIC'),xlabel('入射角度'),ylabel('空间方位谱(dB)'),
title('解相干空间平滑与经典MUSIC算法'),grid on;