%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   空間譜估計(jì)中 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
close all
clc
%%%%%%%%%%%%%%%%%%%%%%%  generate signal  %%%%%%%%%%%%%%%%%
f0=5000;                      %入射信號頻率
fs=1000;                      %采樣頻率
ts=1/fs;                      %采樣間隔
M=8;                           %陣元數(shù)
L=100;                        %快拍數(shù)
SNR=5;                        %信噪比  %%%%%%                        
the0=[-40 -20 10 20 40];
Q=length(the0);                           %信號源數(shù)

c=1500;
lamda=c/f0;
d=lamda/2;                  %c為水下聲速，d為陣元間距

t=(1:L)*ts;                 % 采樣時(shí)刻序列
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)));  %獨(dú)立信源
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];                         %獨(dú)立信號源矩陣  基帶信號
ss2=[s3;s4;s5;s6;s7];                        %相干信號源矩陣
%%%%%%%%%%%%%%%%%%%%% 生成導(dǎo)向矢量矩陣  %%%%%%%%%%%%%%%%%%%
m=0:M-1;
A=exp(-j*2*pi*m'*d/lamda*sin(the0*pi/180));
%%%%%%%%%%%%%%%%%%%%% 產(chǎn)生零均值方差為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;                            %陣列接收相互獨(dú)立的信號，情況1
X2=A*ss2+Nn;                            %陣列接收的相干信號， 情況2

% R1=X1*X1'/L;                              %陣列采樣數(shù)據(jù)協(xié)方差矩陣
R2=X2*X2'/L;
R=R2;

    %%%%%%%%%%%%空間平滑Music算法%%%%%%%%%%%%%%%%%%%%
    p=3;mm=6;%將均勻線陣分成相互交錯(cuò)的5個(gè)子陣，每個(gè)子陣陣元數(shù)為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);
    %%%%%%%%%%%%經(jīng)典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));   %搜索范圍的導(dǎo)向矢量陣
    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','經(jīng)典MUSIC'),xlabel('入射角度'),ylabel('空間方位譜(dB)'),
title('解相干空間平滑與經(jīng)典MUSIC算法'),grid on;