MVDR波束形成MATLAB实现

clc
clear
%%设置初始参数
v=3e8;
f=2.4e9;
lamda=v/f;
d=0.5*lamda;
w=2*pi*f;
m=10;%阵元数
n=1000;%快拍数
snr=0;%0dB
inr=20;%20dB
thetas=90/180*pi;
thetai=40/180*pi;
%%生成导向向量a（theta）
theta=0:180;
theta=theta/180*pi;
tao=d*cos(theta)/v;
atheta = zeros(181,10);%阵列流形向量
for ii=1:181
    for jj=1:10
        atheta(ii,jj)=exp(1i*w*tao(ii)*(jj-1));
    end
end
atheta=atheta.';
%期望和入射信号导向向量产生
a_theta_s=zeros(1,10);
a_theta_i=zeros(1,10);
for ii=1:10
    a_theta_s(ii)=exp(1j*w*tao(91)*(ii-1));
    a_theta_i(ii)=exp(1j*w*tao(41)*(ii-1));
end
a_theta_i=a_theta_i.';
a_theta_s=a_theta_s.';
    
%%入射信号及干扰信号的产生、
S=zeros(2,1000);
T=randi([0,1],[1,2*n]);%2000的行向量
T_IQ=reshape(T,2,n).';%reshape成两个1000的行向量的转置
symbol=bi2de(T_IQ,'left-msb'); %把10变成2
Table=(1/sqrt(2))*[-1-j -1+j 1-j 1+j];
S(1,:)=Table(symbol+1);%取矩阵中的第k行元素  期望信号
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T=randi([0,1],[1,2*n]);%2000的行向量
T_IQ=reshape(T,2,n).';%reshape成两个1000的行向量的转置
symbol=bi2de(T_IQ,'left-msb'); %把10变成2
Table=(1/sqrt(2))*[-1-j -1+j 1-j 1+j];
S(2,:)=Table(symbol+1);%取矩阵中的第k行元素  期望信号
Noise=sqrt(1/2)*(randn(m,n)+j*randn(m,n));
X=a_theta_s*10^(snr/10)*S(1,:)+a_theta_i*10^(inr/10)*S(2,:)+Noise;
R0=X*X'/n;
R1=a_theta_s*a_theta_s'+a_theta_i*a_theta_i'*100+eye(m);



w0=inv(R0)*a_theta_s/(a_theta_s'*inv(R0)*a_theta_s);
w1=inv(R1)*a_theta_s/(a_theta_s'*inv(R1)*a_theta_s);

H0=w0'*atheta;
H1=w1'*atheta;

H0=10*log(abs(H0));
H1=10*log(abs(H1));

angles=0:1:180;
subplot(2,1,1)
plot(angles,H0);
set(gca,'YLim',[-100,0])
title('H0')
subplot(2,1,2)
plot(angles,H1);
set(gca,'YLim',[-100,0])
title('H1')

H0为快拍1000次后的估计出的增益方向图
H1为已知信噪比情况下的理想增益方向图

从生成图的过程来看，只已知信号的期望方向，但是信噪比，干噪比均未知，干扰信号的方向也未知，却能通过快拍采样产生的估计协方差矩阵将干扰信号方向增益进行陷波。这也是该方法与最小二乘法的本质区别，是否已知波达方向和幅度。