Hybrid Beamforming for Millimeter Wave

论文

SYSTEM MODEL

根据在毫米波 MIMO 无线传输中广泛使用的 Saleh-Valenzuela 信道模型，在 n = 1, 2, ., Nc 时，DFRC BS 与第 n 个通信用户之间的信道可表示为一个向量

$$h_n=\sqrt{\frac{N_{BS}}{L_n}}\sum _{l=1}^{L_n}{g}_l^{(n)}{\alpha }^H(N_{BS},{\theta }_l^{(n)})\text{\hspace{1em}(2)}$$

其中 Ln 是可解析信道路径的数量。$g_l^{(n)}$l 和${\theta }_l^{(n)}$l 分别表示第 l 个信道路径的信道增益和离开角（AoD），其中 l = 1，2，. . .，Ln。$\alpha (N_{BS},{\theta }_l^{(n)})$表示信道导向矢量，表示为

$$\alpha (N_{BS},{\theta }_l^{(n)})=\frac{1}{\sqrt{N_{BS}}}[1,e^{j\pi {\theta }_l^{(n)}},...,e^{j(N_{BS}-1)\pi {\theta }_l^{(n)}}{]}^T\text{\hspace{1em}(3)}$$

对于毫米波 MIMO 传输，DFRC BS 通常采用包括模拟波束成形和数字波束成形的混合波束成形架构。我们将模拟波束成形器和数字波束成形器分别记为 ${\mathbf{F}}_{RF}$ 和 ${\mathbf{F}}_{BB}$，其中 ${\mathbf{F}}_{RF}\in {\mathbb{C}}^{N_{BS}\times N_{RF}}$ 和 ${\mathbf{F}}_{BB}\in {\mathbb{C}}^{N_{RF}\times N_{RF}}$ 。在 DFRC BS 进行混合波束成形之前的传输信号记为 $x\in {\mathbb{C}}^{N_{RF}}$，满足 E{x} = 0 和 $Ex{x}^H=I_{N_{R}F}$。x 的前端部分，包括 x1, x2, ., xNc ，是要传输给 Nc 个通信用户的通信信号。x 的后部包括 $x_{Nc+1},.,x_{N_{RF}}$ , 是某一时刻的雷达波形。需要注意的是，在本文中，我们只考虑一个时间实例，重点是空间域的波束成形设计，而跳过时域中多个时间实例的波形设计。那么，Nc 个通信用户的接收信号（记为 $\text{y}\triangleq \text{ [}{\text{y}}_{\text{1}}\text{, }{\text{y}}_{\text{2}}\text{, }...\text{ , }{\text{y}}_{\text{Nc}}{\text{]}}^{\text{T}}\in {\mathbb{C}}^{N_c}$）可表示为

$$y=\mathbf{H}{\mathbf{F}}_{RF}{\mathbf{F}}_{BB}x+\eta \text{\hspace{1em}(4)}$$

其中$\mathbf{H}\triangleq [h_1^T,h_2^T,...,h_{N_c}^T{]}^T\in {\mathbb{C}}^{N_c\times N_{BS}}$是 DFRC BS 和 Nc 个通信用户之间的信道矩阵，并且 η 是满足$\eta \sim CN(0,{\sigma }^2{\mathbf{I}}_{N_C})$的加性白色高斯噪声（AWGN）向量。为了简化符号，我们定义$\mathbf{F}\triangleq [f_1,f_2,...,f_{N_{RF}}]\triangleq {\mathbf{F}}_{RF}{\mathbf{F}}_{BB}\in {\mathbb{C}}^{N_{BS}\times N_{RF}}$，其中 F 的第 n 列表示为 fn，对于 n = 1，2，. . .，Nc.然后第 n 个通信用户（n = 1，2，. . .，Nc）的 SINR 可以表示为

$${\gamma }_n=\frac{{\left|h_n{f}_n\right|}^2}{\sum _{i=1,i\ne n}^{N_{RF}}{\left|h_n{f}_i\right|}^2+{\sigma }^2}\text{\hspace{1em}(7)}$$

注意，不同通信用户可能具有不同的无线服务质量要求，如，视频服务需要比文本服务更高的 SINR。通过将第 n 个用户的 SINR 要求的阈值定义为 Γn，我们可以将 SINR 约束写为 γn ≥ Γn。 另一方面，为了检测多个目标，DFRC BS 需要产生各种雷达波束。作为一个简单的示例，DFRC BS 可以使用 DFT 码字扫描不同的空间角度。假设感兴趣的角度空间由 M 个点，φ1，φ2，. . .，φM，其中 −1 ≤ φ1 < φ2 <。. . < φM ≤ 1。较大的 M 导致角度空间的更精细采样和光束的更精确描述。然后，基于 M 个采样点的空间采样矩阵可以表示为 $\Phi =[\alpha (N_{BS},{\varphi }_1),\alpha (N_{BS},{\varphi }_2),...,\alpha (N_{BS},{\varphi }_M){]}^T$ 实际上，DFRC BS 的发射波束是 Σ NRF i=1 f i。则发射波束的波束图案为：|Φ ∑NRF i=1 f i|其实质是将发射波束投影到与采样点对应的信道导向矢量上，然后获得投影的绝对值。注意，向量的绝对值意味着获得向量的每个条目的绝对值以形成相同维度的向量。 给定一个目标雷达波束方向图$b\in {\mathbb{R}}^M$，其中 b 的每一项都是非负的，在满足通信用户信干噪比和 DFRC 基站总发射功率的约束下，以接近 b 为目标设计 DFRC 基站的发射波束。然后，发射波束设计问题可以公式化为

$$\underset{f_1,f_2,...,f_{N_{RF}}}{\min }{\Vert D(\left|\Phi \sum _{i=1}^{N_{RF}}{f}_i\right|-b)\Vert }_2\text{\hspace{1em}(9a)}$$ $$s.t.\sum _{i=1}^{N_{RF}}{\Vert f_i\Vert }_2^2⩽P_T\text{\hspace{1em}(9b)}$$ $${\gamma }_n⩾{\Gamma }_n,n=1,2,...,N_c\text{\hspace{1em}(9c)}$$

其中，PT 表示 DFRC BS 的总传输功率，(9c)是通信用户的 SINR 约束。D 是预定义的正对角矩阵，其中对角项是角度空间的对应采样点处的权重。D 中的权重越大，表示在相应的采样点处接近 b 的要求越高。

HYBRID BEAMFORMING DESIGN

为了给毫米波 MIMO ISAC 系统中的发射波束设计提供额外的自由度，我们将相位向量$p\in {\mathbb{R}}^M$引入到目标波束方向图 b，其中|Pi| = 1，其中 i = 1，2，· · ·，M。我们进一步定义了一个非负对角矩阵$A\in {\mathbb{R}}^{M\times M}$，其中 A 的对角元素来自于 b 的对应元素，即 A = diag{b}。则发射波束设计问题公式化可以重写为

$$\begin{array}{c}\underset{f_1,f_2,...,f_{N_{RF}},p}{\min }{\Vert D(\Phi \sum _{i=1}^{N_{RF}}{f}_i-\mathbf{A}p)\Vert }_2\\ s.t.\sum _{i=1}^{N_{RF}}{\Vert f_i\Vert }_2^2⩽P_T,\\ {\gamma }_n⩾{\Gamma }_n,n=1,2,...,N_c,\\ \left|p_i\right|=1,i=1,2,...,M\end{array}\text{\hspace{1em}(11)}$$

为了求解（11），我们提出了一种交替最小化方法。

1. 给定 p。f1,f2,...,fNRF 的优化可以表示为

   $$\begin{array}{c}\underset{f_1,f_2,...,f_{N_{RF}}}{\min }{\Vert D(\Phi \sum _{i=1}^{N_{RF}}{f}_i-\mathbf{A}p)\Vert }_2\\ s.t.\sum _{i=1}^{N_{RF}}{\Vert f_i\Vert }_2^2⩽P_T,\\ {\gamma }_n⩾{\Gamma }_n,n=1,2,...,N_c\end{array}\text{\hspace{1em}(12)}$$

   我们定义了一个辅助矩阵 $S\triangleq [\underset{N_{RF}}{\underbrace{{\mathbf{I}}_{N_{BS}},{\mathbf{I}}_{N_{BS}},...,{\mathbf{I}}_{N_{BS}}}}]\in {\mathbb{R}}^{N_{BS}\times (N_{RF}{N}_{BS})}$ ,它实质上是将 NRF 的标识矩阵 ${\mathbf{I}}_{N_{BS}}$ 并排组合在一起。 进一步定义 $f\triangleq [f_1^T,f_1^T,...,f_{N_{RF}}^T]\in {\mathbb{R}}^{N_{RF}{N}_{BS}}$，它将不同的波束成形向量串联在一起。那么，DFRC BS 的发射波束可重写为 $\sum _{i=1}^{N_{BS}}{f}_i=Sf$, 我们进一步定义 $S\triangleq [\underset{i-1}{\underbrace{0,...,0}},{\mathbf{I}}_{N_{BS}},\underset{N_{RF}-i}{\underbrace{0,...,0}}]\in {\mathbb{R}}^{N_{BS}\times (N_{RF}{N}_{BS})}$,注意，Si 是由 NRF - 1 零矩阵和一个标识矩阵 ${\mathbf{I}}_{N_{BS}}$ 组成的。那么通信用户的 SINR 约束可以改写为 ${\Vert t\Vert }_2⩽\sqrt{1+\frac{1}{{\Gamma }_n}}{h}_n{S}_{n}f,n=1,2,...,N_c$，其中，$t\triangleq [h_n{S}_{1}f,h_n{S}_{2}f,...,h_n{S}_{n}f,\sigma {]}^T\in {\mathbb{R}}^{N_{RF}+1}$ 则优化可以重写为

   $$\begin{array}{c}\underset{f}{\min }{\Vert D(\Phi Sf-\mathbf{A}p)\Vert }_2\\ s.t.{\Vert f\Vert }_2^2⩽P_T,{\Vert t\Vert }_2⩽\sqrt{1+\frac{1}{{\Gamma }_n}}{h}_n{S}_{n}f,n=1,2,...,N_c\end{array}$$

   这是一个二阶锥规划（SOCP）问题，可以使用现有的优化工具箱解决。

2. 给定 f1，f2，. . .，f NRF。p 的优化可以表示为

   $$\begin{array}{c}\underset{p}{\min }{\Vert D(\Phi \sum _{i=1}^{N_{RF}}{f}_i-\mathbf{A}p)\Vert }_2\\ s.t.\left|p_i\right|=1,i=1,2,...,M\end{array}\text{\hspace{1em}(20)}$$

   这是约束 LS 估计问题。在没有约束的情况下，无约束 LS 解是 $\stackrel{\sim }{p}=W\sum _{i=1}^{N_{RF}}{f}_i$， 其中，$W\triangleq (A^H{D}^{H}DA{)}^{-1}{A}^H{D}^{H}D\Phi$,是与优化过程无关的常数矩阵，并且可以在优化之前基于$\Phi$和 A 计算。

将 (20) 的可行解表示为 $\stackrel{\wedge }{p_i}$，其第 i 个条目是 $\stackrel{\wedge }{p_i}=\frac{\stackrel{\sim }{p_i}}{|\stackrel{\sim }{p_i}|},i=1,2,...,M$。我们交替执行上述 1 和 2 两个步骤，直到触发停止条件。停止条件可以简单地设定为达到预定的最大迭代次数。也可以设置为 (11) 的目标函数小于预定的阈值。假设$\underset{1}{\overset{\wedge }{f}},\underset{2}{\overset{\wedge }{f}},...,\stackrel{\wedge }{f_{N_{RF}}}$ 即为完成上述步骤后的结果。与 (6) 类似，我们定义 $\stackrel{\wedge }{\mathbf{F}}\triangleq [\underset{1}{\overset{\wedge }{f}},\underset{2}{\overset{\wedge }{f}},...,\stackrel{\wedge }{f_{N_{RF}}}]\in {\mathbb{C}}^{N_{BS}\times N_{RF}}$ 。基于 DFRC BS 的发射波束设计，现在我们考虑 FRF 和 FBB 的混合波束成形器设计。请注意，对于毫米波 MIMO 无线系统，模拟波束成形器通常由移相器网络实现，如图 1 所示。因此，我们对 FRF 的每个条目都有恒定的包络约束。那么，给定 (24) 中$\stackrel{\wedge }{\mathbf{F}}$，确定 FRF 和 FBB 的混合波束成形设计问题可表示为

$$\underset{F_{RF},F_{BB}}{\min }{\Vert \stackrel{\wedge }{\mathbf{F}}-F_{RF}{F}_{BB}\Vert }_F\text{\hspace{1em}(25c)}$$ $$s.t.{\Vert F_{RF}{F}_{BB}\Vert }_F^2⩽P_T\text{\hspace{1em}(25b)}$$ $$\left|F_{RF}(m,n)\right|=1,m=1,2,...,N_{BS},n=1,2,...,N_{RF}\text{\hspace{1em}(25c)}$$

其中 (25c) 是移相器导致的恒定包络约束，(25b) 是总传输功率约束。事实上，我们可以暂时忽略 (25b) 来求解 (25)，然后将得到的 F BB 归一化以满足 (25b)。我们仍然采用交替最小化法，迭代执行以下两个步骤。

1. 给定 FRF。(25) 中 FBB 的优化可表示为

   $$\underset{F_{BB}}{\min }{\Vert \stackrel{\wedge }{\mathbf{F}}-F_{RF}{F}_{BB}\Vert }_F\text{\hspace{1em}(26)}$$

   这是一个 LS 问题，其解为

   $$\overline{F_{BB}}=(F_{RF}^H{F}_{RF}{)}^{-1}{F}_{RF}^H\stackrel{\wedge }{\mathbf{F}}\text{\hspace{1em}(27)}$$

2. 给定 FBB。（25）中 FRF 的优化可表示为

   $$\underset{F_{RF}}{\min }{\Vert \stackrel{\wedge }{\mathbf{F}}-F_{RF}{F}_{BB}\Vert }_F\text{\hspace{1em}(28a)}$$ $$s.t.\left|F_{RF}(m,n)\right|=1,m=1,2,...,N_{BS},n=1,2,...,N_{RF}\text{\hspace{1em}(28b)}$$

   这是一个典型的黎曼流形优化问题，可用现有工具箱求解。 上述 1) 和 2) 两个步骤反复进行，直到触发停止条件。假设$\stackrel{\wedge }{F_{RF}},\stackrel{\wedge }{F_{BB}}$ 已得到，我们最终对$\stackrel{\wedge }{F_{BB}}$ 进行归一化， $\stackrel{\wedge }{F_{BB}}\leftarrow \frac{\sqrt{P_{\Gamma }}}{\Vert \stackrel{\wedge }{F_{RF}}\stackrel{\wedge }{F_{BB}}\Vert }\stackrel{\wedge }{F_{BB}}$来满足 (25b)。

代码：

%% main.m
clc;
clear;
close all;

%% initialize H
L_n = 3;

N_BS = 128;     % 128 antennas

% 信道增益 cn(0,1)   cn(0,0.01)
g_1 = [GD(0,0.5),GD(0,0.5),GD(0,0.5)];
g_2 = [GD(0,0.05),GD(0,0.05),GD(0,0.05)];
g_3 = [GD(0,0.05),GD(0,0.05),GD(0,0.05)];


% set aoa
theta_1 = SAGL(-70);
theta_2 = SAGL(-40);
theta_3 = SAGL(-10);

A_1 = zeros(N_BS,1);
A_2 = zeros(N_BS,1);
A_3 = zeros(N_BS,1);


m = 0:(N_BS-1);
A_1 = 1 / sqrt(N_BS) * exp(1j * m' * theta_1);% 没写pi
A_2 = 1 / sqrt(N_BS) * exp(1j * m' * theta_2);
A_3 = 1 / sqrt(N_BS) * exp(1j * m' * theta_3);


% BS和第n个用户之间的信道
h_1 = sqrt(N_BS / L_n) * (g_1(1) * A_1' + g_2(1) * A_1' + g_3(1) * A_1');
h_2 = sqrt(N_BS / L_n) * (g_1(2) * A_2' + g_2(2) * A_2' + g_3(2) * A_2');
h_3 = sqrt(N_BS / L_n) * (g_1(3) * A_3' + g_2(3) * A_3' + g_3(3) * A_3');

% 信道矩阵
H = [h_1.',h_2.',h_3.'].'; % 3x128

%% D

M = 400;    % sample points
D = eye(M);


%% Phi

x = 1:1:M;
omega = -1 + (2 .* x - 1)/M;   % sin(x) 1xM

n = 0:(N_BS-1);
Phi_trans = 1 / sqrt(N_BS) * exp(1j * n' .* pi .* omega); % N_BSxM

Phi = Phi_trans.'; % MxN_BS

%% b

P_T = 20;
N_RF = 3;       % 3 RF chains
result = sqrt(2*N_RF*P_T/(SAGL(30)-SAGL(10)+SAGL(60)-SAGL(40)));


t_ang = [SAGL(10);SAGL(30);SAGL(40);SAGL(60)];  % target angle

b = zeros(M,1);

% omega_1,omega_2,……,omega_M



for y = 1 : M
    if (omega(y)>= t_ang(1) && omega(y) <= t_ang(2)) || (omega(y)>= t_ang(3) && omega(y) <= t_ang(4))
        b(y) = result;
    else
        b(y) = 0;
    end

end

%% epsilon
epsilon = 0 ;

%% gamma1,2,3(N_c)


f_1 = rand(N_BS,1) + 1j*zeros(N_BS,1);
f_2 = rand(N_BS,1) + 1j*zeros(N_BS,1);
f_3 = rand(N_BS,1) + 1j*zeros(N_BS,1);

% SINR
gamma_1 = power(abs(h_1 * f_1),2) / (power(abs(h_2 * f_2),2) + power(abs(h_3 * f_3),2)) + epsilon^2;
gamma_2 = power(abs(h_2 * f_2),2) / (power(abs(h_1 * f_1),2) + power(abs(h_3 * f_3),2)) + epsilon^2;
gamma_3 = power(abs(h_3 * f_3),2) / (power(abs(h_1 * f_1),2) + power(abs(h_2 * f_2),2)) + epsilon^2;


%% p

p = exp(1j * linspace(0, 2*pi, M)).'; % Mx1

%% S and S_i

I_nbs = eye(N_BS);
S =  repmat(I_nbs,1,N_RF);


S_s_1 = zeros(N_BS, N_BS*N_RF); % S_1
S_s_2 = zeros(N_BS, N_BS*N_RF); % S_2
S_s_3 = zeros(N_BS, N_BS*N_RF); % S_3

for x = 1:3
    col_index = ((x-1)*N_BS + 1) : (x*N_BS);
    if x == 1
        S_s_1(:,col_index) = eye(N_BS);
    elseif x == 2
        S_s_2(:,col_index) = eye(N_BS);
    elseif x == 3
        S_s_3(:,col_index) = eye(N_BS);
    end
end

%%
F = [f_1,f_2,f_3]; % 128x3
f = [f_1.',f_2.',f_3.'].'; % 3x128

%%

% A = zeros(M,M);
A = diag(b);

%%

Tou = [30; 30; 30];


%%
% 通过while循环迭代，直到满足某个条件
epochs = 10;
epoch = 1;
while epoch <= epochs

    t = [h_1*S_s_1*f,h_2*S_s_2*f,h_3*S_s_3*f,epsilon].';

    cvx_begin

    num = N_RF * N_BS;

    variable f(num) complex
    minimize(norm(D * (Phi * S * f - A * exp(1j * p)) , 2))

    subject to

    norm(f, 2) <= sqrt(P_T);
    norm(t,2) <= sqrt(1+1/Tou(1))*t(1);
    norm(t,2) <= sqrt(1+1/Tou(2))*t(2);
    norm(t,2) <= sqrt(1+1/Tou(3))*t(3);

    cvx_end
    epoch = epoch + 1;


    f1_opt = f(1:128);
    f2_opt = f(129:256);
    f3_opt = f(257:384);


    W = pinv(A' * D' * D * A) * A' * D' * D * Phi;

    p_opt = W * (f1_opt + f2_opt + f3_opt);

    p_fiopt = zeros(M, 1);

    for pm=1:M
        if abs(p_opt(pm)) ~= 0
            p_fiopt(pm) = p_opt(pm) / abs(p_opt(pm));
        else
            p_fiopt(pm) = 0;
        end
    end

    p = p_fiopt;

end

%%

R = Phi * S * f;

abc = linspace(-90,90,400);
plot(abc,abs(R))

%% SAGL.m
function sin_angle = SAGL(angle)
    angle_in_degrees = angle;
    angle_in_radians = deg2rad(angle_in_degrees);
    sin_angle = sin(angle_in_radians);
end

%% GD.m
function complex_samples = GD(mu,variance)
    real_part = normrnd(mu, sqrt(variance));
    imag_part = normrnd(mu, sqrt(variance));
    complex_samples = complex(real_part, imag_part);
end