上一篇《阵列信号处理及MATLAB仿真-----阵列信号绪论》里面说了阵列信号处理研究的四个主要问题：波束形成技术、空间谱估计、信号源定位、信源分离 。接下来我们就波束形成来做一个详细的学习。
一、波束形成的定义：
首先说一下它的物理意义，阵列天线的方向图是全方向的，但是阵列的输出经过加权求和后，可以被调整到阵列接收的方向，也就是说的增益聚集在一个方向，相当于形成一个“波束”。
那么波束形成技术的基本思想是：通过将各个阵元输出进行加权求和，在一段时间内将天线阵列波束“导向”到一个方向，对期望信号得到最大输出功率的导向位置，及给出了波达方向估计。
二、常用波束形成算法：
波束形成原理：
为什么要用波束形成：利用阵元直接相干叠加而获得的输出，它的缺点在于只有在垂直于阵列平面方向的入射波在阵列的输出端才能同相叠加，以致形成方向图中的主瓣的极大值。反过来说，要是阵列可以围绕它的中心轴旋转，那么当阵列输出最大时，空间波一定是由垂直与阵列平面的方向入射的。但是，有些天线阵列很大、无法转动。所以设计一种相控阵天线法（也就是说的常规波束形成法）。
波束形成方法中，阵列输出选取一个适当的加权向量来补偿各个阵元的传播延时，让在某一期望方向上阵列输出可以同相叠加，进而使阵列在该方向上产生一个主瓣波束，在其他方向上产生较小的响应，用这种方法对整个空间进行波束扫描就能够确定待检测信号的方位。
以一维M元等距线阵为例，波束形成算法结构图如图所示，假设空间信号为窄带信号，每一个通道用一个复加权系数来调整通道的幅度和相位。
阵列的输出可表示为：
若采用向量表示各个阵元输出及加权系数：
则阵列的输出用向量表示为：
为了在某一个方向θ上补偿各个阵元之间的时延从而形成一个主瓣，常规波束形成器在期望方向上的加权向量的构成为：
根据这个加权向量，若空间仅有一个来自方向θ的信号，它的方向向量a(θ)的形式和这个权向量一样，那么我们就得到了：
根据上述公式可以得到常规的波束形成器的输出功率可表示为：
其中，矩阵R为阵列输出x(t)的协方差矩阵，R=E｛x(t)xH(t)｝
接下来分析一下常规波束形成的角分辨率的问题。一般来说，当空间有两个同频率信号投射到阵列时，假如他们的空间方位角的间隔小于阵列主瓣波束宽度，那么不仅不能无法分辨。还会影响系统的正常工作，也就是对于阵列远场中的两个信号源，仅当他们之间的角度分离大于阵元间隔（阵列孔径）的倒数时，他们才能被分辨出来，这就是瑞利准则，瑞利准则说明了常规波束形成方法存在的缺点是角度分辨率低。
波束形成的最佳权向量：
上面第一节我们了解了波束形成的原理，从原理可以看出其中最重要的是里面提到的权向量。那么波束形成的最佳权向量又是什么呢？咱们一块来接着学习。
首先我们要有一个概念，波束形成的“导向”作用是通过调整加权系数完成的。令权向量为w=[w1,...,wM]T,那么输出可以为：
对于不同的权向量，公式5对不同方向的回波有不同的响应，从而导致不同方向的空间波束。通常采用移相器进行加权处理，也就是只调整信号相位，不改变信号幅值，因为信号在任一瞬间各个阵元上的幅度是相同的。
终上所述，假设空间只有一个来自方向θk的回波，它的方向向量为a(θk),那么当权向量w取a(θk)时，输出y(n)=a(θk)Ha(θk)=M最大，从而实现导向定位作用。这时，各路的加权信号为相干叠加，这一结果叫做空域匹配滤波。
但是匹配滤波只有在白噪声背景下才是最佳的，假如存在干扰信号需要另外考虑。接下来我们看一下复杂情况下的波束形成。
假设在远距离有一个感兴趣的信号d(t)(也叫期望信号，它的波达方向为θd)和J个不感兴趣的信号ij(t),j=1,2,...,J(干扰信号，波达方向为θij)。令每个阵元上的加性白噪声为nk(t)（方差为δ2）。那么在这些条件下，第k个阵元上的接收信号可表示为：
公式6右边分别表示信号、干扰和噪声。若采用矩阵形式表示，那么有
简单化简后：
那么，N个快拍的波束形成器输出y(t)=wHx(t)(t=1,2,3,...,N)的平均功率为：
我们忽略不同用户之间的相互作用（也就是感兴趣信号、不感兴趣信号和噪声信号三者之间的相互作用）那么输出的平均功率最后可化简为：
（8）
当N趋于无穷大时，上式可以写成：
上式中R=E{x(t)xH(t)}是阵列输出的协方差矩阵。另外当N趋近于无穷大时，公式8可表示为
公式10中，假设了各加性噪声具有相同的方差σn2
为了保证来自方向θd期望信号的正确接收，并同时抑制其他J个干扰，很容易根据公式10得到权向量的约束条件。
上述约束条件也叫做“置零条件”，原因是强迫接收阵列波束方向图的“零点”指向所有J个干扰信号，那么公式10可以化简为：
但是从提高信噪比的方面来说，干扰置零不是最佳的方法，因为虽然选定的权值使得干扰输出为0.但有可能使得噪声变大。所以抑制干扰和噪声应一起考虑进来，因此波束形成器最佳权向量的确定可以叙述为：在权向量约束条件下，求满足下式的权向量w。
这个问题的解可以通过拉格朗日乘子法去求，具体求解过程就不一一写出来了。我们通过最后的求解得到接收的来自方向θd期望信号的波束形成器的最佳权向量为：
其中μ为一个比例常数，θd期望信号的波达方向。这样就能够决定J+1个发射信号的波束形成的最佳权向量。此时，波束形成器将只接收来自方向θd的期望信号，并抑制所有其他波达方向的非期望信号。
我们注意到约束条件wHa(θd)=1,可以等价于aH(θd)w=1，那么我们把公式（11）两边分别乘以aH(θd)。那么最后我们可以得到常数μ。
根据上述阵列处理的基本问题可以看出，空域处理和时域处理截然不同。传统的时域处理主要提取的是信号的包络信息，作为载体的载波在完成传输任务后不再有用；而传统的空域处理则为了区别波达方向，主要是利用载波在不同阵元间的相位差，包络反而不起作用，并且利用窄带信号的复包络在各个阵元之间的延迟可以忽略的特点简化计算。
从公式11和12中可以看出，波束形成器的最佳权向量w取决于阵列方向向量，但对于移动用户来说它的方向向量一般是未知的，需要估计出来。因此在计算波束形成的最佳权向量之前必须在已经知道的阵列几何结构的前提下先估计期望信号的波达方向。
三、波束形成的准则
波束形成算法是在一定准则下综合各个输入信息来计算最优权值的数学方法。下面说一些最常见、最重要的准则。
最大信噪比准则（MSNR）：使得期望信号分量功率和噪声分量功率之比最大，但是必须知道噪声的统计量和期望信号的波达方向。
最大信干噪比准则（MSINR）:使得期望信号功率与干扰功率及噪声分量功率的和的比值最大。
最小均方误差准则（MMSE）：在非雷达应用中，阵列协方差矩阵中通常都是含有期望信号的。基于这种情况提出了该准则。使得阵列输出与某种期望响应的均方误差最小，不需要知道期望信号的波达方向。
最大似然比准则（MLH）：在对有用信号完全先验未知的情况下，参考信号无法设置，因此，在干扰噪声背景下，首先要取得对有用信号的最大似然估计。
线性最小方差准则（LFMV）：对有用信号形式和来向完全已知，在某种约束条件下使阵列输出的方差最小。
可以证明，在理想情况下这几种准则得到的权是等价的，且可写成通式wopt=R-1Ha(θd),通常为维纳解。其中，a(θd)是期望信号的方向函数，亦称约束导向向量，RH是不含期望信号的阵列协方差。
比较了MMSE、MSNR和LCMV这3种统计最佳约束形成技术的准则、代价函数、最佳解及具有的优缺点。
四、MMSE、MSNR和LCMV程序MATLAB仿真（里面的一些程序我也不太确定，如果大佬发现有哪些不足和错误的地方，还望赐教）
1、MMSE
clc
clear all
close all
M=18; %天线数?
L=100; %蹇媿鏁?
thetas=20; %入射角度
thetai=[-30 30]; %干扰角度
n=[0:M-1]';
vs=exp(-j*pi*n*sin(thetas/180*pi)); %信号方向向量
vi=exp(-j*pi*n*sin(thetai/180*pi)); %干扰方向向量
f=16000; %载波频率
t=[0:1:L-1];
snr=10; %信噪比
inr=10; %信干比
A=[vs vi];
rvar=10^(snr/20);
de=rvar*exp(j*2*pi*f*t); %期望信号
v1 = sin(2*pi*2*f*t/(8*f)); %干扰信号1
v2 = sin(2*pi*4*f*t/(8*f)); %干扰信号2
st=[de;v1;v2];
xi=sqrt(10^(inr/10/2))*vi*[randn(length(thetai),L)+j*randn(length(thetai),L)]; %干扰信号
noise=[randn(M,L)+j*randn(M,L)]/sqrt(2); %噪声
Xt=A*st+noise; %对比书上的公式
X=xi+noise; %噪声和信号
R=X*X'/(L-1); %求协方差
r=X*conj(de)'/(L-1); %X和de的协方差
wop1=inv(R)*r; %计算权值
sita=48*[-1:0.001:1]; %探测范围
v=exp(-j*pi*n*sin(sita/180*pi)); %不同角度的方向矢量
B=abs(wop1'*v);
figure(1)
plot(sita,20*log10(B/max(B)),'k');
title('MMSE波束图')
xlabel("角度/degree")
ylabel("功率/dB")
grid on
axis([-48 48 -50 0]);
hold off
2、LCMV
clc
clear all
close all
M=18; %天线数
L=100; %快扫个数
thetas=20; %入射角度
thetai=[-30 30]; %干扰角度
n=[0:M-1]';
vs=exp(-j*pi*n*sin(thetas/180*pi)); %信号方向矢量
vi=exp(-j*pi*n*sin(thetai/180*pi)); %干扰方向矢量
f=1600; %载波频率
t=[0:1:L-1]/200;
snr=10; %信噪比
inr=10; %干扰比
A=[vs vi];
rvar=10^(snr/20);
%de=rvar*exp(j*2*pi*f*t/(8*f)); %期望信号
de=sin(2*pi*f*t/(8*f)); %期望信号
%xi1=[randn(length(thetai),L)+j*randn(length(thetai),L)];
xi=sqrt(10^(inr/10/2))*vi*[randn(length(thetai),L)+j*randn(length(thetai),L)]; %干扰信号
noise=[randn(M,L)+j*randn(M,L)]/sqrt(2); %噪声
X=xi+noise; %信号加噪声
R=X*X'/(L-1); %求协方差矩阵
wop1=inv(R)*vs/(vs'*inv(R)*vs);
sita=50*[-1:0.001:1]; %探测范围
v=exp(-j*pi*n*sin(sita/180*pi)); %不同探测角度的方向矢量
B=abs(wop1'*v);
B1=abs(wop2'*v);
figure(1)
plot(sita,20*log10(B/max(B)),'k');
title('LCMV波束图')
xlabel("角度/degree")
ylabel("功率/dB")
grid on
axis([-50 50 -50 0]);
hold off
3、MSNR
clc
clear all
close all
M=18; %天线�
L=100; %快扫个数
thetas=20; %信号方向
thetai=[-30 -15 0 15 30]; %干扰方向
n=[0:M-1]';
vs=exp(-j*pi*n*sin(thetas/180*pi)); %信号方向矢量
vi=exp(-j*pi*n*sin(thetai/180*pi)); %干扰信号
f=16000; %频率
t=[0:1:L-1]/200;
snr=10; %信噪比
inr=10; %信干比�
rvar=10^(snr/20);
de=rvar*exp(j*2*pi*f*t); %期望信号
xs=sqrt(10^(snr/10))*vs*exp(j*2*pi*f*t);
xi1=[randn(length(thetai),L)+j*randn(length(thetai),L)];
xi=sqrt(10^(inr/10/2))*vi*[randn(length(thetai),L)+j*randn(length(thetai),L)]; %干扰信号的方位矢量
noise=[randn(M,L)+j*randn(M,L)]/sqrt(2); %噪声
X=xi+noise; %信号加噪声
A=[vs vi];
St=[de;xi1];
Xt=A*St+noise;
R=X*X'/(L-1); %求信号的协方差矩阵
Rn=noise*noise'/(L-1); %求噪声的协方差矩阵�
MAX_lamda=max(eig(inv(R)*Rn)); %最大特征值��
[m1,n1]=eig(inv(R)*Rn);
[x,y]=find(n1==max(max(n1)));
wop2=m1(:,y); %计算权值
sita=48*[-1:0.001:1]; %探测范围
v=exp(-j*pi*n*sin(sita/180*pi));
B1=abs(wop2'*v);
figure(1)
plot(sita,20*log10(B1/max(B1)),'k');
title('MSNR波束图')
xlabel("角度/degree")
ylabel("功率/dB")
grid on
axis([-48 48 -50 0]);
hold off
4、网上一位博主的LCMV程序，他的这个程序中求协方差矩阵的方法，是按照书上的公式得到Xt，然后求的协方差，虽然它的程序图效果很好，但是多次运行之后会出现效果很差的波束图。具体原因有待进一步研究。大家可以去看一下波束形成算法：LCMV实现（MATLAB）
clc;
clear;
M = 18; % 天线个数�
lambda = 10;
d = lambda / 2;
L = 100; %蹇媿鏁�
thetas = [10]; % 信号方向
thetai = [-30 30]; % 干扰方向
n = [0:M-1]';
vs = exp(-1j * 2 * pi * n * d * sind(thetas) / lambda); % 信号
vn = exp(-1j * 2 * pi * n * d * sind(thetai) / lambda); % 干扰信号
f = 1600; % 频率
t = [0:L-1];
di = sin(2*pi*f*t/(8*f)); % 期望信号
vn1 = sin(2*pi*2 * f*t/(8*f)); % 干扰信号1
vn2 = sin(2*pi*4 * f*t/(8*f)); % 干扰信号2
A = [vs vn];
St = [di;vn1;vn2];
Xt = A*St + randn(M,L); % 对标书上的公式
R_x = 1/L * (Xt * Xt');
R_x_inv = inv(R_x);
W_opt = R_x_inv * vs / (vs' * R_x_inv * vs);
% 探测范围
sita = 50 * [-1:0.001:1];
% 探测信号�
v = exp(-1i*2*pi*n* d*sind(sita)/lambda);
B = abs(W_opt' * v);
plot(sita,20*log10(B/max(B)),'k')
title('波束图�')
xlabel('角度/degree')
ylabel('功率/dB')
grid on
axis([-50 50 -50 0]);
主要参考：张小飞.阵列信号处理及MATLAB实现[M].电子工业出版社