clc;
clear;
%% 特征分解
N=200;
snr=[0 10 10];                                     %信噪比
signal_direction=[20 40 -60];                      %实际信号方向
signal_frequency=[1.2e3 0.8e3 1.6e3];              %信号频率  
fs=1e4;                                            %采样频率                                 
array_num=16;                                      %阵元数目                          
signal_num=3;                                      %信号数目
dlamda=0.5;                                        % d/lamda
D_number=3;                                        %干扰数目
%% 阵列流行
A=zeros(array_num,signal_num);
for i=1:signal_num
    A(:,i)=exp(j*2*pi*dlamda*sin(signal_direction(i)*pi/180)*[0:(array_num-1)].');       %实际信号的的阵列流行
end
p=10.^(snr/10);
s0=p(1).^0.5*exp((j*2*pi*(signal_frequency(1)/fs))*[1:N]');
s1=p(2).^0.5*exp((j*2*pi*(signal_frequency(2)/fs))*[1:N]');
s2=p(3).^0.5*exp((j*2*pi*(signal_frequency(3)/fs))*[1:N]');
S=[s0,s1,s2];                                      %信号空间                                         
Noise=randn(array_num,N)+j*randn(array_num,N);     %定义噪声

X=A*S.'+Noise;                                     %实际信号定义

R=X*X'./N;                                         %实际相关矩阵(我们只能得到真实的相关矩阵,但是我们还以为他是我们假设的那个)

W1=inv(R)*A(:,1)./(A(:,1)'*inv(R)*A(:,1));         %MVDR波束形成器的权向量


%%
theta=-90:0.1:90;
A1=zeros(array_num,signal_num);
for i=1:length(theta)
    A1(:,i)=exp((j*2*pi*dlamda)*sin(theta(i)*pi/180)*[0:(array_num-1)]');
    beam1(i)=abs(W1'*A1(:,i));
end
beam1=beam1/max(beam1);    % 归一化

beam1=10*log(beam1); %空间匹配滤波器


figure(2);
plot(theta,beam1);
xlabel('入射角度');ylabel('dB');
grid on