讀書期間，在上《無線通信基礎》這門課時，因為教材 [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