让波束照进现实：从相控阵到波束赋形

读书期间，在上《无线通信基础》这门课时，因为教材 [1] 里都是关于数学公式的推导，包含矩阵论、随机过程等知识，我发觉自己很难在课堂上集中注意力。学期结束时，对一些概念和理论的理解也不是很深刻。后来，在从事了多年 LTE 和 5G 等无线通信的算法和系统设计后，我开始明白理解基本的概念和理论是如此重要，因为它们是不同应用和标准的基础，是进一步创新的源泉。

MATLAB 中包含了大量的示例和参考应用，可以帮助大家理解复杂的无线通信理论和概念。本文将通过一个相关示例展示 MATLAB 将如何助力您学习和理解什么是相控阵，波束赋形，预编码等概念，以及如何利用相控阵处理来获得的分集和空间复用增益，以提供无线通信系统的性能。

什么是相控阵？

波束赋形是一种在 5G 甚至 6G 通信系统、雷达、声纳、医疗成像和音频阵列系统中广泛使用的技术。要理解波束赋形，首先需要了解什么是波束。

简而言之，波束是由相控阵形成。相控阵中包含多个测量或者发射波的单元。这些单元按某种结构排列，并共同作用以产生所需的波束。其中，每个单元信号的相位都是可以调节的，从而来控制波束的方向。想要了解相控阵的更多细节，您可以观看 Brian Douglas 的视频 "What Are Phased Arrays?" [3]，他使用了大量的动画来解释这个概念。

在无线通信中，相控阵的单元一般是天线单元。从下面 Douglas 视频的快照中，您可以看到由全向天线所组成的均匀线性阵列 (Uniform Linear Array，ULA) 所产生的波束。在某些方向上，来自不同天线单元的信号被叠加，因此信号强度被增强；而在其他方向上，天线信号相互抵消，因此信号强度变弱。

波束的方向取决于天线阵列的设置，包括天线单元的数量、单元之间的间距、阵列的几何形状以及单个天线的方向图。MATLAB 中的 Sensor Array Analyzer[4] 应用程序为您提供了一个交互式工具，用于配置天线阵列，获得其基本性能特性，并且在 2D 或 3D 空间中可视化其波束的方向图。

在 MATLAB 中安装了 Phased Array Toolbox 后，您可以在工具栏的 “APPs” 选项卡的 “信号处理与通信” 下，找到这个 App 的图标，也可以在 MATLAB 中使用以下命令打开 Sensor Array Analyzer 应用程序。

sensorArrayAnalyzer

下面的动画演示如何使用 Sensor Array Analyzer[4]，您可以注意到这个应用程序能帮助我们通过交互式的界面，设置和调整相控阵的参数，可视化相控阵的波束，以及将相控阵对准某个方向。基于不同的导向矢量，我们就可以在不移动天线阵列物理位置的情况下，移动波束的方向啦！

接下来，将向您展示如何对天线阵列进行波束赋形，以提高无线链路的信噪比 (Signal-to-Noise Ratio，SNR)。

怎样做波束赋形?

天线阵列已经成为 5G 标准配置的一部分。这种无线通信系统通常被称为多输入多输出 (Multiple-Input-and-Multiple-Output, MIMO) 系统。为了展示多天线的优势，我们假设一个无线链路部署在 60 GHz，这是 5G 中新增的毫米波波段。

rng(0)     % set random seed
c = 3e8;    % propagation speed
fc = 60e9;   % carrier frequency

不失一般性，您可以将发射机放置在原点，将接收机放置在大约 1.6 公里外。

txcenter = [0;0;0];   % center of the transmitter
rxcenter = [1500;500;0]; % center of the receiver

然后根据发射机和接收机的中心位置，计算出射角（Angle-of-Departure, AoD）和到达角（Angle-of-Arrival, AoA）。

[~,txang] = rangeangle(rxcenter,txcenter); % AoD
[~,rxang] = rangeangle(txcenter,rxcenter); % AoA

下面将展示在两种常见无线信道下，如何利用天线阵列提高无线链路的性能。并通过链路仿真比较 MIMO 波束赋形相比于 SISO 的优势。

视距信道

首先，考虑视距 (Line-of-Sight，LOS) 传播。这是一种最简单的无线信道，通常发生在广阔的农村地区。在 LOS 信道下，发射机信号能直接到达接收机。在这种情况下，如何使用天线阵列提高接收机的信噪比，从而提高链路的误码率 (Bit Error Rate，BER) 呢?

作为对比，您可以先如下所示仿真 LOS 信道下 SISO 链路的误码率性能，其中 scatteringchanmtx 函数用于为不同的发射和接收阵列配置生成基于散射体的信道矩阵。该函数通过模拟发射阵列和接收阵列之间的多个散射体，假设信号首先从发射阵列传输到所有散射体，然后从散射体反射到接收阵列。

在这种情况下，每个散射体定义了发射阵列和接收阵列之间的一个信号路径，因此所得到的信道矩阵其实描述了一个多径环境。您可以通过选中该函数，然后右键单击以打开该函数来阅读其详细信息和代码实现。

考虑一个 SISO 无线通信链路，发射机和接收机都只有一个位于节点中心的天线单元，并且环境中存在一个从发射机到接收机的直接路径。这样的 LOS 信道可以视为多路径环境的一种特殊情况。

txsipos = [0;0;0]; % antenna position to the transmitter    
rxsopos = [0;0;0]; % antenna position to the receiver


g = 1; % path gain
sisochan = scatteringchanmtx(txsipos,rxsopos,txang,rxang,g); % generate a SISO LOS channel


Nsamp = 1e6;       % sample rate
x = randi([0 1],Nsamp,1); % data source


ebn0_param = -10:2:10;  % SNRs, unit in dB
Nsnr = numel(ebn0_param); % total number of SNRs


ber_siso = helperMIMOBER(sisochan,x,ebn0_param)/Nsamp; % caculate the bit error rate (BER)

至此，您已经获得了 SISO LOS 链路的误码率性能。通过查看文末的辅助函数 helperMIMOBER，大家可以注意到这个示例中考虑的是 BPSK 调制。

现在您可以创建一个 MIMO 链路。假设发射机和接收机都是包含半波长间距的 4 个天线单元的 ULA。

Ntx = 4; % number of transmitting antennas
Nrx = 4; % number of receiving antennas


txarray = phased.ULA('NumElements',Ntx,'ElementSpacing',lambda/2); % antenna array of the transmitter
txmipos = getElementPosition(txarray)/lambda; 


rxarray = phased.ULA('NumElements',Nrx,'ElementSpacing',lambda/2); % antenna array of the receiver 
rxmopos = getElementPosition(rxarray)/lambda;

有了接收天线阵列，阵列单元的接收信号是相干的，因此可以将接收阵列的波束方向预对准发射机方向，以提高信噪比。此外，发射机也可以将其天线阵列的主波束预先对准接收机，以达到进一步的性能提升。

txarraystv = phased.SteeringVector('SensorArray',txarray,...
  'PropagationSpeed',c); % the steering vector of transmit array
rxarraystv = phased.SteeringVector('SensorArray',rxarray,...
  'PropagationSpeed',c); % the steering vector of receive array


wt = step(txarraystv, fc, txang)';
wr = conj(step(rxarraystv, fc, rxang));


mimochan = scatteringchanmtx(txmipos,rxmopos,txang,rxang,g);
ber_mimo = helperMIMOBER(mimochan,x,ebn0_param,wt,wr)/Nsamp;


helperPlotSpatialMIMOScene(txmipos,rxmopos,txcenter,rxcenter,NaN,wt,wr)

波束赋形意味着调制信号在发射端发送前需要乘上一个权重向量 wt，而接收信号在接收端解调前需要乘上一个合并向量 wr。使用helperPlotSpatialMIMOScene 函数对场景进行可视化之后，您可以在上图中看到，发射天线和接收天线的主波束是指向对方的。

让我们比较在 LOS 传播信道下， SISO 和 MIMO 链路的误码率。

helperBERPlot(ebn0_param,[ber_siso(:) ber_mimo(1,:).']);
legend('SISO LOS','MIMO LOS');

上图的误码率 vs. 信噪比曲线显示，发射阵列和接收阵列的波束赋形一共带来了约 12 dB 的阵列增益。这里需要注意的是，增益是在假设发射机需要知道接收机的方向或位置的情况下实现的，同时接收机也需要知道信号的传入方向，通常这个角度可以利用波达估计算法得到。

多径信道

在大多数情况下，视距 (Line-of-Sight，LOS) 传播是比较困难的，无线通信通常发生在多径衰落环境中。本示例的其余部分将探讨在这种情况下相控阵是如何提高无线通信系统性能的。

不失一般性，假设信道中有 10 个随机放置的散射体，那么从发射机到接收机将有 10 条路径，如下图所示。

Nscat = 10; % number of scatterers
[~,~,scatg,scatpos] = helperComputeRandomScatterer(txcenter,rxcenter,Nscat);
helperPlotSpatialMIMOScene(txmipos,rxmopos,txcenter,rxcenter,scatpos,NaN,NaN)

简单起见，假设沿着所有路径传播的信号都在同一个符号周期内到达，因此信道是频率平坦的，而不是频率选择性的。

在这样的多径衰落信道中，我们需要仿真信道在多个时隙上进行变化，从而来获得误码率曲线。假设仿真时间为 1000 帧，每帧有包含 10000 比特。SISO 链路在多径信道下的误码率曲线可以如下得到。

Nframe = 1e3;    % number of frames
Nbitperframe = 1e4; % number of bits per frame
Nsamp = Nframe*Nbitperframe; % total number of data samples


x = randi([0 1],Nbitperframe,1); % generated the source data


nerr = zeros(1,Nsnr);
for m = 1:Nframe
  sisompchan = scatteringchanmtx(txsipos,rxsopos,Nscat);
  wr = sisompchan'/norm(sisompchan);
  nerr = nerr + helperMIMOBER(sisompchan,x,ebn0_param,1,wr);
end
ber_sisomp = nerr/Nsamp;

在 SISO 情况下，如果您比较 LOS 下的信道矩阵和多径下的信道矩阵，可以发现它只包含一个元素，而且由于多径衰落，元素的值从实数变为复数。在这种情况下，可以在接收端使用 wr = sisompchan'/norm(sisompchan) 计算接收信号所需要右乘的权值来提高信噪比。

现在将展示在 MIMO 情况下，如何通过仿真获得误码率曲线。其中，diagbfweights 函数被用来基于信道矩阵计算发射机的预编码权值 wp 和接收机的合并权值 wc。您可以右键单击该函数并打开它，查看其详细信息。

x = randi([0 1],Nbitperframe,Ntx);


nerr = zeros(Nrx,Nsnr);
for m = 1:Nframe
  mimompchan = scatteringchanmtx(txmipos,rxmopos,Nscat);
  [wp,wc] = diagbfweights(mimompchan);
  nerr = nerr + helperMIMOBER(mimompchan,x,ebn0_param,wp,wc);
end
ber_mimomp = nerr/Nsamp;

要理解预编码和合并权重是如何计算的，首先需要知道什么是奇异值分解 (Singular Value Decomposition，SVD)。使用以下命令，在MATLAB 查看 SVD 的文档。

help svd

简单来说，SVD ([U,S,V] = SVD (X)) 产生了与 X 相同维数且非负对角元素递减的对角矩阵 S，以及酉阵 U 和 V，是得 X = U*S*V'。尝试 在MATLAB 中输入以下命令：

bf = wp*mimompchan*wc
[u,s,v] = svd(mimompchan)

比较 bf 和 s，您可以发现 diagbfweights 函数是基于 SVD 的。wp*mimompchan*wc 是一个对角矩阵。在包含多个散射体的环境下，MIMO 多径信道的信道矩阵通常是满秩的。这意味这样的 MIMO 多径信道将可以类似于原始信道内的多个正交子信道。

因此，您可以通过 MIMO 多径信道同时发送多个数据流，这个策略被称为空间复用。空间复用的基本思想是使发射阵列中不同单元发送的数据流可以从接收信号中独立恢复。它的目标不全是提高信噪比，而主要是为了提高信息吞吐量。

现在绘制 SISO 多径情况下，以及 MIMO 多径情况下前两个数据流的误码率曲线。为了对比，您可以加上在上一节得到的 SISO LOS 链路的误码率曲线。

helperBERPlot(ebn0_param,[ber_siso(:) ...
ber_sisomp(:) ber_mimomp(1,:).' ber_mimomp(2,:).']);
legend('SISO LOS','SISO Multipath','MIMO Multipath Stream 1','MIMO Multipath Stream 2');

从上图中可以看到，与 SISO LOS 链路的误码率曲线相比，SISO 多径链路的误码率曲线由于多径传播引起的衰落而下降得慢得多。同时，在 MIMO 多径情况下，第一个子信道对应着主要的发射和接收方向，因此曲线斜率上的增益，也就是分集增益没有损失。虽然第二个流使用了一个不占主导地位的子信道，因此不能提供像第一个流那样高的增益，但由于现在可以同时传输多个数据流，整个链路的信息吞吐量可以得到提升。

总结

本文解释了如何使用阵列信号处理提高 MIMO 无线通信系统的性能。根据信道特性，阵列既可以通过阵列增益，也可以通过分集增益提高信噪比，或者通过空间复用提高信息容量。

您已经了解了 MIMO 无线通信的一些基础知识，比如什么是相控阵、波束成形、预编码、分集增益和空间复用。这个示例当中所使用的辅助函数将附在文末以帮助您理解实现细节。

除了直观地在各发射天线单元之间均匀分配功率外，您还可以通过注水算法进一步提高信道容量。具体请查看和运行本文的参考示例 “Improve SNR and Capacity of Wireless Communication Using Antenna Arrays” [2]。

最后，如果您想继续了解如何在 5G 标准下实现波束赋形，请观看教程视频 "Beamforming for MU-MIMO in 5G New Radio" [5]。同时，MATLAB 的无线通信文档[6] 中也包含有大量的示例和参考应用，可以帮助您更容易和清楚地理解相关的概念和理论。各种 AI forWireless[7] 的示例也将给您提供更多如何在无线通信领域应用AI技术的参考和支持。

推荐您观看最新一期的小迈步之通信系统设计——从基础到 AI+：

点开展开阅读附录：辅助函数。

附录：辅助函数

第一个函数 helperPlotSpatialMIMOScene 用于在 2D 区域绘制 MIMO 场景。SISO 场景可以被认为是 MIMO 场景的特殊情况。

function helperPlotSpatialMIMOScene(txarraypos,rxarraypos,...
txcenter,rxcenter,scatpos,wt,wr)
  rmax = norm(txcenter-rxcenter);
  
  if size(txarraypos,2) == 1 && size(rxarraypos,2) == 1
    spacing_scale = 0;
  elseif size(txarraypos,2) == 1
    spacing_scale = rmax/20/mean(diff(rxarraypos(2,:)));
  else
    spacing_scale = rmax/20/mean(diff(txarraypos(2,:)));
  end
  
  txarraypos_plot = txarraypos*spacing_scale+txcenter;
  rxarraypos_plot = rxarraypos*spacing_scale+rxcenter;
  
  clf;
  figure
  hold on;
  
  plot(txarraypos_plot(1,:),txarraypos_plot(2,:),'kv','MarkerSize',10,...
    'MarkerFaceColor','k');
  text(txcenter(1)-85,txcenter(2)-15,'TX');
  plot(rxarraypos_plot(1,:),rxarraypos_plot(2,:),'kv','MarkerSize',10,...
    'MarkerFaceColor','k');
  text(rxcenter(1)+35,rxcenter(2)-15,'RX');
  if isnan(scatpos)
    set(gca,'DataAspectRatio',[1,1,1]);
    line([txcenter(1) rxcenter(1)],[txcenter(2) rxcenter(2)]);
  else
    hscat = plot(scatpos(1,:),scatpos(2,:),'ro');
    for m = 1:size(scatpos,2)
      plot([txcenter(1) scatpos(1,m)],[txcenter(2) scatpos(2,m)],'b');
      plot([rxcenter(1) scatpos(1,m)],[rxcenter(2) scatpos(2,m)],'b');
    end
    legend(hscat,{'Scatterers'},'Location','SouthEast');
  end
  
  nbeam = rmax/5;
  if ~isnan(wt)
    txbeam_ang = -90:90;
    txbeam = abs(wt*steervec(txarraypos,txbeam_ang)); % wt row
    txbeam = txbeam/max(txbeam)*nbeam;
    [txbeampos_x,txbeampos_y] = pol2cart(deg2rad(txbeam_ang),txbeam);
    plot(txbeampos_x+txcenter(1),txbeampos_y+txcenter(2),'k');
  end
  if ~isnan(wr)
    rxbeam_ang = [90:180 -179:-90];
    rxbeam = abs(wr.'*steervec(rxarraypos,rxbeam_ang)); % wr column
    rxbeam = rxbeam/max(rxbeam)*nbeam;
    [rxbeampos_x,rxbeampos_y] = pol2cart(deg2rad(rxbeam_ang),rxbeam);
    plot(rxbeampos_x+rxcenter(1),rxbeampos_y+rxcenter(2),'k');
  end
  axis off;
  hold off;
end

第二个函数 helperComputeRandomScatterer 用于计算散射点位置等信息。

function [txang,rxang,g,scatpos] = ...
  helperComputeRandomScatterer(txcenter,rxcenter,Nscat)
  ang = 90*rand(1,Nscat)+45;
  ang = (2*(rand(1,numel(ang))>0.5)-1).*ang;
  
  r = 1.5*norm(txcenter-rxcenter);
  scatpos = phased.internal.ellipsepts(txcenter(1:2),rxcenter(1:2),r,ang);
  scatpos = [scatpos;zeros(1,Nscat)];
  g = ones(1,Nscat);
  
  [~,txang] = rangeangle(scatpos,txcenter);
  [~,rxang] = rangeangle(scatpos,rxcenter);
end

第三个函数 helperMIMOBER 用于计算采用 BPSK MIMO 链路的误码率。SISO 链路可以被认为是 MIMO 链路的特殊情况。

function nber = helperMIMOBER(chan,x,snr_param,wt,wr) 
  Nsamp = size(x,1);
  Nrx = size(chan,2);
  Ntx = size(chan,1);
  
  if nargin < 4
        wt = ones(1,Ntx);
    end
    
    if nargin < 5
        wr = ones(Nrx,1);
    end    
    
    xt = 1/sqrt(Ntx)*(2*x-1)*wt; % map to bpsk
    nber = zeros(Nrx,numel(snr_param),'like',1); % real
    
    for m = 1:numel(snr_param)
        n = sqrt(db2pow(-snr_param(m))/2)*... 
(randn(Nsamp,Nrx)+1i*randn(Nsamp,Nrx));
        y = xt*chan*wr+n*wr;
        xe = real(y)>0;
    nber(:,m) = sum(x~=xe);
  end
end

第四个函数 helperBERPlot 用于绘制误码率曲线。

function helperBERPlot(ebn0,ber) 
  figure
  h = semilogy(ebn0,ber);
  set(gca,'YMinorGrid','on','XMinorGrid','on');
  markerparam = {'x','+','^','v','.'};
  for m = 1:numel(h)
    h(m).Marker = markerparam{m};
  end
  xlabel('E_b/N_0 (dB)');
  ylabel('BER');
  lim([9e-6 1]);
end

审核编辑：汤梓红