自适应波束形成算法及MATLAB仿真算法（RLS和LMS）

一、概述

　　自适应研究的重点一直都是自适应算法，经典的自适应波束形成算法可分为闭环算法（反馈控制算法）和开环算法（也称直接求解方法）。

　　一般而言，闭环算法比开环算法要简单，实现方便，但其收敛速率受到系统稳定性要求的限制。闭环算法包括最小均方（LMS）算法、差分最陡下降（DSD）算法、加速梯度（AG）算法及3种算法的变形。

　　后来更多的集中在开环算法的研究上。开环算法是一种直接求解方法，不存在收敛问题，可以提供更快的暂态响应性能，但同时也受到处理精度和阵列协方差矩阵求逆运算量的控制。事实上，开环算法可认为是实现自适应处理的最佳途径，现在被广泛应用，但开环运算量大。鉴于该问题，人们想到了采用自适应处理技术，它减轻了自适应算法的实时计算负荷，并且能较快产生较快的自适应响应。

二、自适应波束形成的最佳权向量：

　　传统自适应波束形成的结构如图1所示，波束形成的权重通过自适应信号处理获得。假定阵元m的输出为连续基带（复包络）信号xm(t)，经过A/D转换后，变成了离散基带信号xm(k),m=0,1,…,M-1，并且以阵元0为参考点。另外，假定共有Q个信源存在，wq(k)表示在时刻k对第q个信号解调所加的权向量，其中q=1,…,Q。权向量用某种准则确定，来使解调出来的第q个信号的质量在某种意义上最优。

 图1 传统自适应波束形成的结构

　　在最佳波束形成中，权向量通过代价函数的最小化来确定。在典型情况下，这种代价函数越小，阵列输出信号的质量就越好，因此当代价函数最小时，自适应阵列输出信号的质量最好。

　　代价函数有两种常用的形式：最小均方误差（MMSE）和最小二乘方法(LS)

　　1、MMSE方法

　　MMSE方法是在波形估计、信号检测和系统参数辨识等信号处理中广泛应用的一种优化方法。所以MMSE方法就是使得估计误差y(k)-dq(k)的均方值最小化，即代价函数取：

式中，x(k)=[x0(k),x1(k),…xM-1(k)]T。代价函数为第q个信号的阵列输出与该信号在时刻k的期望形式之间的平方误差的数学期望值，将上式展开可得。

由上式可得：

上式中Rx是数据向量x(k)的自相关矩阵，即

 而rxd是数据向量x(k)与期望信号dq(k)的互相关向量，即：

 令

，那么得到：

这就是MMSE方法下的最佳阵列权向量，它是维纳滤波理论中最佳滤波器的标准。

　　2、LS方法

　　在MMSE方法中，代价函数定义为阵列输出与第q个用户期望响应之间误差平方的总体平均（均方差），实际数据向量总是有限长的，假如直接定义代价函数为其误差平方，则得到LS方法。

　　假定有N个快拍的数据向量x(k),k=1,…,N，定义为代价函数；

 则求出的梯度为：

令梯度等于0，可以得到：

 这就是最小二乘意义下针对第q个用户的波束形成器的最佳权向量，上式中X和dq分别表示数据向量和期望信号向量，其值为：

 　　上面介绍的MMSE和LS的核心问题是，在对第q个用户进行波束形成时，需要在接收端使用该用户的期望响应。为了提供这个期望响应，必须周期性的发送对发射机和接收机二者都知道的训练序列。但训练序列会占用频谱资源，这是MMSE和LS共有的问题。

　　一种可以代替训练序列的方法是采用决策指向更新对期望响应进行学习。在决策指向更新中，期望信号样本的估计根据阵列输出和信号解调器输出重构。由于期望信号是在接收端产生的，不需要发射数据的知识，因此不用训练序列。

三、权向量更新的自适应算法

　　上面提到的自适应阵列的最佳权向量的确定需要求解相关方程，一般来说，不是很希望直接求解方程。理由如下：

　　（1）由于移动目标环境是实时变化的，所以权向量的解必须实时更新。

　　（2）由于估计最佳解需要的数据是含有噪声的，所以希望使用一种更新的技术，能够利用已经求得的权向量求平滑最佳响应的估计，从而减小噪声的影响。

　　所以希望使用自适应算法周期更新权向量。

　　自适应算法即可采用迭代模式也可采用分块模式。

　　（1）所谓迭代模式，就是在每个迭代步骤，n时刻的权向量加上一个校正量后，就组成了n+1时刻的权向量，用这个新的权向量来逼近最佳权向量。

　　（2）在分块模式下，权向量并不是实时更新的，而是间隔一定时间周期才更新；由于在一定时间周期内对应于一个数据块而不是一个数据点。所以这种更新又称为分块更新。

　　为了使阵列系统能够自适应工作，就必须将第二节介绍的方法归到自适应算法中。以MMSE为例子，看看如何将它变为一种自适应算法。

　　考虑到随机梯度算法，其更新权向量的一般公式为：

　　其中，

，μ称为收敛因子，它控制自适应算法的收敛速度，则：

　　上式中的数学期望用各自的瞬时值代替，就得到了k时刻的梯度估计值。

　　公式中，

，代表阵列输出与第q个用户期望响应dq(k)之间的瞬时误差。容易证明，梯度估计是真实梯度的无偏估计。

　　将梯度估计值公式带入更新权向量的一般公式中可得熟悉的LMS自适应算法。

　　MMSE方法可以用LMS算法实现，而LS方法的自适应算法可以递推最小二乘（RLS）算法，下表给出了3种自适应波束形成算法的比较，分别是LMS算法、RLS算法和Bussgang算法。从表中可以看出LMS和RLS算法需要使用训练序列，但Bussgang算法不需要训练序列。

　　注意Bussgang算法中，g[y(k)]是一个非线性的估计子，，它对解调器的输出信号y(k)起作用，并用g[y(k)]代替期望信号d(k),然后产生误差函数e(k)。

四、matlab仿真：

　　本文仅对上表中的最小均方算法和递推最小二乘算法进行仿真，如果后面有用到Bussgang算法，再来学习这个算法。需要注意一点的是在RLS算法中，遗忘因子的值越接近1波束形成效果越好，越接近0，波束形成的效果会很差。

 1、最小均方（LMS）算法matlab仿真

%最小均方算法（LMS）采用迭代模式，在每个迭代步骤n时刻的权向量加上一个校正量后
%就得到n+1时刻的权向量，用它来逼近最佳权向量。
clc
clear all
close all
M=16;                      %天线数
K=2;                       %信源数
theta=[0 30];              %信号入射角度
d=0.3;                     %天线间距
N=500;                    %采样点数
Mean_noise=0;              %噪声均值
varn_noise=1;              %噪声方差

SNR=10;                               %信噪比
INR=10;                               %干噪比

pp=zeros(100,500);
pp1=zeros(100,500);

rvar1=sqrt(varn_noise)*10^(SNR/20);          %信号功率
rvar2=sqrt(varn_noise)*10^(INR/20);          %干扰功率

%for q=1:100
s=[rvar1*exp(j*2*pi*(50*0.001*[0:N-1]));rvar2*exp(j*2*pi*(100*0.001*[0:N-1]+rand))];      %生成源信号
A=exp(-j*2*pi*d*[0:M-1]'.*sind(theta));              %来自波达方向theta的发射源的方向向量
noise=sqrt(varn_noise/2)*(randn(M,N)+j*randn(M,N));
Y=A*s+noise;                                             %接收信号
%LMS算法
L=200;                 %快拍数
w1=[];
beam1=[];
for i=1:length(theta)
    de=s(i,:);             %期望信号
    mu=0.0005;             %步长参数
    w=zeros(M,1);
    for k=1:N
        y(k)=w'*Y(:,k);               %预测下一个采样和误差
        e(k)=de(k)-y(k);              %误差
        w=w+mu*Y(:,k)*conj(e(k));     %更新权向量
    end
    w1=[w1 w];
    for i=1:L
        a=exp(-j*2*pi*d*[0:M-1]'.*sin(-pi/2+pi*(i-1)/L));
        beam(i)=20*log10(abs(w'*a));
        
    end
    beam1=[beam1;beam];
    
end
sum1=sum(beam1);
figure
angle=-90:180/200:(90-180/200);
plot(angle,beam1(1,:));
grid on;
xlabel("角度/degree")
ylabel("波束图/dB")
title("入射角度0，LMS波束图")
figure
angle=-90:180/200:(90-180/200);
plot(angle,beam1(2,:));
grid on;
xlabel("角度/degree")
ylabel("波束图/dB")
title("入射角度30，LMS波束图")

 2、递推最小二乘（RLS）算法仿真

clc
clear all
close all
M=16;    %天线数
K=2;     %信号源个数
theta=[0 30];      %信号入射角度
d=0.5;             %天线间距
N=500;             %采样点数

Mean_noise=0;              %噪声均值
varn_noise=1;              %噪声方差

SNR=10;                               %信噪比
INR=10;                               %干噪比

rvar1=sqrt(varn_noise)*10^(SNR/20);          %信号功率
rvar2=sqrt(varn_noise)*10^(INR/20);          %干扰功率

s=[rvar1*exp(j*2*pi*(50*0.001*[0:N-1]));rvar2*exp(j*2*pi*(100*0.001*[0:N-1]+rand))];      %生成源信号
A=exp(-j*2*pi*d*[0:M-1]'.*sind(theta));              %来自波达方向theta的发射源的方向向量
noise=sqrt(varn_noise/2)*(randn(M,N)+j*randn(M,N));
Y=A*s+noise;
beam1=[];
lamda=0.99;             %遗忘因子
L=200;                 %快拍数
for i=1:length(theta)
    %RLS算法  
    de=s(i,:);             %期望信号
    w=zeros(M,1);          %初始化权值
    P0=0.5*eye(M);
    P=P0;
    y=zeros(M,1);
    for k=1:N
        v(:,k)=P*Y(:,k);
        u=inv(lamda)*v(:,k)./(1+inv(lamda)*Y(:,k)'*v(:,k));
        a(k)=de(k)-w'*Y(:,k);
        w=w+u*conj(a(k));
        P1=inv(lamda)*[eye(M)-u.*Y(:,k)']*P;
        P=P1;
    end
    beam=zeros(1,L);
    for i=1:L
        a1=exp(-j*2*pi*d*[0:M-1]'.*sin(-pi/2+pi*(i-1)/L));
        beam(i)=20*log10(abs(w'*a1));
    end
    beam1=[beam1;beam];
end
sum1=sum(beam1);
figure
angle=-90:180/200:(90-180/200);
plot(angle,beam1(1,:));
grid on;
xlabel("角度/degree")
ylabel("波束图/dB")
title("入射角度0，RLS波束图")
figure
angle=-90:180/200:(90-180/200);
plot(angle,beam1(2,:));
grid on;
xlabel("角度/degree")
ylabel("波束图/dB")
title("入射角度30，RLS波束图")

五、参考内容

　　《阵列信号处理及MATLAB仿真》