基于特征子空间的波束形成算法原理介绍及MATLAB实现

原理介绍

期望信号的导向矢量属于导向矢量张成的子空间，因此可以将有误差的名义导向矢量向该子空间投影来修正导向矢量，进而减小误差。
对接收数据的协方差矩阵进行特征分解可以得到
$$\hat{R}=\sum _{m=1}^M{\lambda }_m{u}_m{u}_m^H$$
其中，${\lambda }_m$为其特征值，$u_m$为${\lambda }_m$对应的特征向量。假设${\lambda }_1\ge {\lambda }_2\ge \cdots {\lambda }_M$，可以选择前$K$较大特征值对应的特征向量张成的子空间为信号子空间，因此可以得到信号子空间为
$$U_s=\left[u_1,u_2,\cdots ,u_K\right]$$
进而得到投影矩阵
$$P_s=U_s{U}_s^H$$
通过对名义导向矢量进行投影可以得到修正后的导向矢量
$${\hat{a}}_0=P_s{\bar{a}}_0$$
那么可以得到最终的加权系数为
$$w=\frac{R^{-1}{\hat{a}}_0}{{\hat{a}}_0^H{R}^{-1}{\hat{a}}_0}$$
权系数也可以写作
$$w=U_s{\Lambda }_s^{-1}{U}_s^H{\bar{a}}_0$$
该算法的关键点在于重构投影矩阵，也就是怎么估计出准确的信号子空间。在实际应用中，基于大特征值所对应的特征矢量作为信号子空间的方法受限于信源数信息，及协方差矩阵的特征值容易受到低信噪比条件的影响导致大小非常接近，很难通过比较大小准确估计信号子空间。
针对该问题，研究者提出了一种改进的算法。
首先计算名义导向矢量${\bar{a}}_0$在特征矢量$u_m$，$m=1,2,\cdots ,M$上投影的模值
$$p\left(m\right)={\left|u_m^H{\bar{a}}_0\right|}^2$$
对模值进行排序
$$p_1\ge p_2\ge \cdots \ge p_M$$
对于给定的门限$\rho$，找到满足下面式子的最小$K^{\prime }$值
$$\left(p_1+p_2+\cdots p_{K^{\prime }}\right)/\sum _{m=1}^{M}p\left(m\right)>\rho$$
那么将$p_m$，$m=1,2,\cdots ,K^{\prime }$对应的特征向量作为新的矩阵$U_s$，构造新的投影矩阵
$$P=U_s{U}_s^H$$
通过对名义导向矢量进行投影可以得到修正后的导向矢量
$${\hat{a}}_0=P_s{\bar{a}}_0$$
那么可以得到最终的加权系数为
$$w=\frac{R^{-1}{\hat{a}}_0}{{\hat{a}}_0^H{R}^{-1}{\hat{a}}_0}$$
该算法期望通过相关性估计一个期望信号子空间，然而本质上只有期望信号子空间和干扰信号子空间组成的信号子空间，而期望信号子空间和干扰信号子空间并不是直和关系，除非导向矢量之间是正交的，即使是正交的，他们也是很难通过特征矢量分割的，而且该算法严重依赖门限值。同时，该算法在低信噪比的条件下较之前的特征子空间算法有较大的改善。

仿真参数设置

参数名称	参数值
阵元数	10
期望信号角度	$-5^{\circ }$
干扰信号角度	$-30^{\circ }$、$30^{\circ }$
SNR	10dB
INR	20
快拍数	60

仿真结果如下

这里采用的是固定信号子空间进行投影后的结果，当然也可以选择利用门限值确定信号子空间的方法，两者各有优缺点，不过都能改善波束形成器的性能。在阵列模型存在误差的条件下，更能体现该方法的优越性，这里就不再叙述。
代码如下：

clear;
close all;
clc;
warning off
%% 初始化
M = 10;             %阵元数
fs = 5000;          % 采样频率
f = 1000;           % 信号频率
snap = 60;         % 快拍数
T = 0.5;           %采样时间
t = 1/fs:1/fs:T;
c = 340;
lamda = c/f;              %波长
d = 0.5*lamda;          %阵元间距
theta0 =-5;                %期望信号角度
theta1 =-30;                %干扰角度
theta2 = 30;                %干扰角度
snr=10;                     %信噪比
inr1 =20;                   %干噪比
inr2 = 20;                   %干噪比
snr_noise = 0;              %噪声功率1,为0dBW

%% 导向矢量
a0 = exp(-1j*2*pi*d*sind(theta0)*(0:M-1)'/lamda);
a1 = exp(-1j*2*pi*d*sind(theta1)*(0:M-1)'/lamda);
a2 = exp(-1j*2*pi*d*sind(theta2)*(0:M-1)'/lamda);

%% 信号、干扰和噪声
tar_sig = wgn(1,length(t), snr);
inf1 = wgn(1,length(t),inr1); 
inf2 = wgn(1,length(t),inr2); 
noise = wgn(M,length(t),snr_noise);

%% 阵列接收信号
rec_sig = a0*tar_sig + a1*inf1 + a2*inf2 + noise; 
interference = a1*inf1 + a2*inf2;
sig = a0 * tar_sig;

%% 协方差矩阵
Rx = rec_sig(:,1:snap)*rec_sig(:,1:snap)'/snap; 
Rs = sig(:,1:snap)*sig(:,1:snap)'/snap; 
Ri = interference(:,1:snap)*interference(:,1:snap)'/snap; 
Rn = noise(:,1:snap)*noise(:,1:snap)'/snap; 

%% EIG算法
[U, V] = eig(Rx);
[V_new, index] = sort(diag(V),'descend');
U_new = U(:,index);
K = 3;
uu=U_new(:,1:K);
P = uu * uu';
a0 = P * a0;
w = inv(Rx)*a0 /(a0'*inv(Rx)*a0);
optimal_weight = w;
theta = -90:0.1:90; % scan angle
p = exp(-1j*2*pi*d*(0:M-1)'*sind(theta)/lamda);
y = optimal_weight'*p;
yy = 20*log10(abs(y)/max(abs(y)));
%% 绘图
figure(1);
plot(theta,yy,'linewidth', 2);
xlabel('角度(\circ)');ylabel('归一化增益(dB)')
grid on;
xlim([-90 90])