智能发射面| 信道估计与代码复现：基于压缩感知

前言

本文是对 Compressed Channel Estimation and Joint Beamforming for Intelligent Reflecting Surface-Assisted Millimeter Wave Systems 一文的分析和复现。文章已在 IEEE Signal Processing Letter上发表，这里给出 arxiv链接：传送门

信道模型与稀疏表示

文章先考虑了一个单用户单天线场景。根据几何建模，基站-反射面（BS-IRS)信道可以表示为：
G = N M ρ ∑ l = 1 L ϱ l a r ( ϑ l , γ l ) a t H ( ϕ l ) \boldsymbol{G}=\sqrt{\frac{N M}{\rho}} \sum_{l=1}^{L} \varrho_{l} \boldsymbol{a}_{r}\left(\vartheta_{l}, \gamma_{l}\right) \boldsymbol{a}_{t}^{H}\left(\phi_{l}\right) G=ρNM l=1∑Lϱlar(ϑl,γl)atH(ϕl)
注意，作者认为， 智能反射面为平面，因此按UPA建模，而基站为ULA线天线，因此
a r ( ϑ l , γ l ) = a x ( u ) ⊗ a y ( v ) \boldsymbol{a}_{r}\left(\vartheta_{l}, \gamma_{l}\right)=\boldsymbol{a}_{x}(u) \otimes \boldsymbol{a}_{y}(v) ar(ϑl,γl)=ax(u)⊗ay(v)
也就是说，反射面端的响应同时与 ϑ l \vartheta_{l} ϑl 和 γ l \gamma_{l} γl 即仰角和水平角相关。其中，
a x ( u ) ≜ 1 M x [ 1 e j u … e j ( M x − 1 ) u ] T a y ( v ) ≜ 1 M y [ 1 e j v … e j ( M y − 1 ) v ] T \begin{array}{l} \boldsymbol{a}_{x}(u) \triangleq \frac{1}{\sqrt{M_{x}}}\left[1 e^{j u} \ldots e^{j\left(M_{x}-1\right) u}\right]^{T} \\ \boldsymbol{a}_{y}(v) \triangleq \frac{1}{\sqrt{M_{y}}}\left[\begin{array}{lll} 1 & e^{j v} & \ldots & \left.e^{j\left(M_{y}-1\right) v}\right]^{T} \end{array}\right. \end{array} ax(u)≜Mx 1[1eju…ej(Mx−1)u]Tay(v)≜My 1[1ejv…ej(My−1)v]T
很容易地，我们有：
G = ( F x ⊗ F y ) Σ F L H ≜ F P Σ F L H \boldsymbol{G}=\left(\boldsymbol{F}_{x} \otimes \boldsymbol{F}_{y}\right) \boldsymbol{\Sigma} \boldsymbol{F}_{L}^{H} \triangleq \boldsymbol{F}_{P} \boldsymbol{\Sigma} \boldsymbol{F}_{L}^{H} G=(Fx⊗Fy)ΣFLH≜FPΣFLH
也就是经典地信道的压缩感知写法，不再赘述了，其中 F P = ( F x ⊗ F y ) \boldsymbol{F}_{P}= \left(\boldsymbol{F}_{x} \otimes \boldsymbol{F}_{y}\right) FP=(Fx⊗Fy), F L \boldsymbol{F}_{L} FL 可以分别理解为UPA响应和ULA响应的码本。因此， Σ \boldsymbol{\Sigma} Σ 就是一个只有 L L L个非零元素的稀疏矩阵。第 l l l个非零元素的位置（第 i i i行，第 j j j列）分别揭示了第 l l l径对应于哪两个码字（也就是 a r \boldsymbol{a}_{r} ar和 a t \boldsymbol{a}_{t} at) ，以及其增益（也就是 ϱ l \varrho_{l} ϱl).

同理， IRS-UE的信道可以压缩写为：
h r = F P α \boldsymbol{h}_{r}=\boldsymbol{F}_{P} \boldsymbol{\alpha} hr=FPα
其中 α \boldsymbol{\alpha} α 是一个只有 L ′ L^\prime L′个非零值的稀疏向量。

信道估计

我们都知道，在IRS的信道估计中，其实就是估计级联信道，即：

H = diag ⁡ ( h r H ) G \boldsymbol{H}=\operatorname{diag}\left(\boldsymbol{h}_{r}^{H}\right) \boldsymbol{G} H=diag(hrH)G
基于上述的压缩表示，我们进一步推导：
H = diag ⁡ ( h r H ) G = ( a ) h r ∗ ∙ G \boldsymbol{H}=\operatorname{diag}\left(\boldsymbol{h}_{r}^{H}\right) \boldsymbol{G} \stackrel{(a)}{=} \boldsymbol{h}_{r}^{*} \bullet \boldsymbol{G} H=diag(hrH)G=(a)hr∗∙G
这里作者引出了一个概念： transposed Khatri-Rao product，这里维基上讲的很清楚，传送门。一张图就能说明了：

一言以蔽之，就是行间克罗内克积。那么也很容验证上式成立了。继续推导：

H = diag ⁡ ( h r H ) G = ( a ) h r ∗ ∙ G = ( b ) ( F P ∗ α ∗ ) ∙ ( F P Σ F L H ) = ( c ) ( F P ∗ ∙ F P ) ( α ∗ ⊗ ( Σ F L H ) ) = ( d ) ( F P ∗ ∙ F P ) ( α ∗ ⊗ Σ ) ( 1 ⊗ F L H ) = ( e ) D ( α ∗ ⊗ Σ ) F L H \begin{aligned} \boldsymbol{H} &=\operatorname{diag}\left(\boldsymbol{h}_{r}^{H}\right) \boldsymbol{G} \stackrel{(a)}{=} \boldsymbol{h}_{r}^{*} \bullet \boldsymbol{G} \\ & \stackrel{(b)}{=}\left(\boldsymbol{F}_{P}^{*} \boldsymbol{\alpha}^{*}\right) \bullet\left(\boldsymbol{F}_{P} \boldsymbol{\Sigma} \boldsymbol{F}_{L}^{H}\right) \\ & \stackrel{(c)}{=}\left(\boldsymbol{F}_{P}^{*} \bullet \boldsymbol{F}_{\boldsymbol{P}}\right)\left(\boldsymbol{\alpha}^{*} \otimes\left(\boldsymbol{\Sigma} \boldsymbol{F}_{L}^{H}\right)\right) \\ & \stackrel{(d)}{=}\left(\boldsymbol{F}_{P}^{*} \bullet \boldsymbol{F}_{P}\right)\left(\boldsymbol{\alpha}^{*} \otimes \boldsymbol{\Sigma}\right)\left(1 \otimes \boldsymbol{F}_{L}^{H}\right) \\ & \stackrel{(e)}{=} \boldsymbol{D}\left(\boldsymbol{\alpha}^{*} \otimes \boldsymbol{\Sigma}\right) \boldsymbol{F}_{L}^{H} \end{aligned} H=diag(hrH)G=(a)hr∗∙G=(b)(FP∗α∗)∙(FPΣFLH)=(c)(FP∗∙FP)(α∗⊗(ΣFLH))=(d)(FP∗∙FP)(α∗⊗Σ)(1⊗FLH)=(e)D(α∗⊗Σ)FLH
其中(b)步就是把一开始信道建模时的压缩表示代入，（c）就是利用了 transposed Khatri-Rao product 的性质：

( A ∙ B ) ( C ⊗ D ) = ( A C ) ∙ ( B D ) (\mathbf{A} \bullet \mathbf{B})(\mathbf{C} \otimes \mathbf{D})=(\mathbf{A} \mathbf{C}) \bullet(\mathbf{B} \mathbf{D}) (A∙B)(C⊗D)=(AC)∙(BD)

大家可以自己验证。 (d)继续使用了克罗内克积的性质，
( A B ) ⊗ ( C D ) = ( A ⊗ C ) ( B ⊗ D ) (\boldsymbol{A B )} \otimes(\boldsymbol{C D})=(\boldsymbol{A} \otimes \boldsymbol{C})(\boldsymbol{B} \otimes \boldsymbol{D}) (AB)⊗(CD)=(A⊗C)(B⊗D)
(e)就是简单的定义了新变量： D = ( F P ∗ ∙ F P ) \boldsymbol{D}=\left(\boldsymbol{F}_{P}^{*} \bullet \boldsymbol{F}_{P}\right) D=(FP∗∙FP). 但需要注意，这里作者提出了一个中啊哟的简化思路：

原文中作者省略了证明，但其实很简单，我以下面这个例子作更简洁的证明：

记 C ∙ D = E \mathbf{C} \bullet \mathbf{D}=\mathbf{E} C∙D=E, 那么 E \mathbf{E} E的前三列其实就是：
E ( : , 1 : 3 ) = d i a g ( C ( : , 1 ) ) D \mathbf{E}(:, 1:3) = \mathrm{diag}(\mathbf{C}(:,1))\mathbf{D} E(:,1:3)=diag(C(:,1))D
看出来了嘛？显然有：
E ( : , 4 : 6 ) = d i a g ( C ( : , 2 ) ) D \mathbf{E}(:, 4:6) = \mathrm{diag}(\mathbf{C}(:,2))\mathbf{D} E(:,4:6)=diag(C(:,2))D
也就是说， E \mathbf{E} E的所有列其实就相当于前三列进行简单的线性变换（乘上一个对角阵）就能得到。这就是文章中Proposition的结论。

D u = D ( : , 1 : M G ) \boldsymbol{D}_{u}=\boldsymbol{D}\left(:, 1: M_{G}\right) Du=D(:,1:MG)

现在，我们可以将式子进一步简化：

H = D ( α ∗ ⊗ Σ ) F L H = D u Λ F L H \boldsymbol{H}=\boldsymbol{D}\left(\boldsymbol{\alpha}^{*} \otimes \boldsymbol{\Sigma}\right) \boldsymbol{F}_{L}^{H}=\boldsymbol{D}_{u} \boldsymbol{\Lambda} \boldsymbol{F}_{L}^{H} H=D(α∗⊗Σ)FLH=DuΛFLH
这里给大家详细讲一下吧。（原来博客有误，经读者指出，现纠正）
首先注意到， D u = F P ∗ \mathbf{D}_u = \mathbf{F}_P^* Du=FP∗, 而 F P \mathbf{F}_P FP 是什么呢？作者没有严明，但代码中显示，这是一个DFT矩阵。这个也是标准的字典生成做法。那么DFT的好处是什么？DFT矩阵任意两列的哈达玛积必是DFT矩阵的一列。大家可以自己验证。

OK，有了这个结论，而 D \mathbf{D} D中不属于 D u \mathbf{D}_u Du的其他列呢？ 都可以看做是 D u \mathbf{D}_u Du中两列的哈达玛积，那么，必然等于 D u \mathbf{D}_u Du中的一列。也就是说 D \mathbf{D} D中的所有列去重之后，其实就是 D u \mathbf{D}_u Du。那么我们就可以把 D ( α ∗ ⊗ Σ ) \mathbf{D}\left(\boldsymbol{\alpha}^{*} \otimes \boldsymbol{\Sigma}\right) D(α∗⊗Σ)简写为 D u Λ \boldsymbol{D}_{u} \boldsymbol{\Lambda} DuΛ了。而且有， Λ \boldsymbol{\Lambda} Λ 的非零元素必定不超过 L × L ′ L \times L^\prime L×L′个。这是因为， ( α ∗ ⊗ Σ ) \left(\boldsymbol{\alpha}^{*} \otimes \boldsymbol{\Sigma}\right) (α∗⊗Σ)中的非零元素必不超过 L × L ′ L \times L^\prime L×L′个。

这是作者原文的描述：

受限于篇幅，我觉得这里简短的讲述可能不太好理解。我举个例子，比如现在某个 ( α ∗ ⊗ Σ ) \left(\boldsymbol{\alpha}^{*} \otimes \boldsymbol{\Sigma}\right) (α∗⊗Σ)中的非零元素在第 ( M G + 10 ) (M_G+10) (MG+10)行，那么他对应的是与 D \mathbf{D} D的 ( M G + 10 ) (M_G+10) (MG+10)列相乘，而根据上面的结论，这一列可以用 D u \mathbf{D}_u Du中的一列表示，那只需要把非零元素乘在那一列对应的 Λ \boldsymbol{\Lambda} Λ的行上即可了，结果是等效的。

所以(11)这个式子的推导核心就在于上面的这个结论： D u \mathbf{D}_u Du的任意一列都可以看做是 D u \mathbf{D}_u Du中两列的哈达玛积。

有了这个重要的结论，我们很快地：

y ( t ) = v H ( t ) H w ( t ) s ( t ) + ϵ ( t ) = ( a ) ( w T ( t ) ⊗ v H ( t ) ) vec ⁡ ( H ) + ϵ ( t ) = ( b ) ( w T ( t ) ⊗ v H ( t ) ) ( F L ∗ ⊗ D u ) vec ⁡ ( Λ ) + ϵ ( t ) = ( c ) ( w T ( t ) ⊗ v H ( t ) ) F ~ x + ϵ ( t ) \begin{aligned} y(t) &=\boldsymbol{v}^{H}(t) \boldsymbol{H} \boldsymbol{w}(t) s(t)+\epsilon(t) \\ & \stackrel{(a)}{=}\left(\boldsymbol{w}^{T}(t) \otimes \boldsymbol{v}^{H}(t)\right) \operatorname{vec}(\boldsymbol{H})+\epsilon(t) \\ & \stackrel{(b)}{=}\left(\boldsymbol{w}^{T}(t) \otimes \boldsymbol{v}^{H}(t)\right)\left(\boldsymbol{F}_{L}^{*} \otimes \boldsymbol{D}_{u}\right) \operatorname{vec}(\boldsymbol{\Lambda})+\epsilon(t) \\ & \stackrel{(c)}{=}\left(\boldsymbol{w}^{T}(t) \otimes \boldsymbol{v}^{H}(t)\right) \tilde{\boldsymbol{F}} \boldsymbol{x}+\epsilon(t) \end{aligned} y(t)=vH(t)Hw(t)s(t)+ϵ(t)=(a)(wT(t)⊗vH(t))vec(H)+ϵ(t)=(b)(wT(t)⊗vH(t))(FL∗⊗Du)vec(Λ)+ϵ(t)=(c)(wT(t)⊗vH(t))F~x+ϵ(t)

这个式子实在是过于简单了，不再赘述，这个不会推的可以参考混合波束成形| 蜂窝系统的信道估计和混合预编码:Channel Estimation and Hybrid Precoding for Millimeter Wave Cellular Systems 和混合波束成形| 宽带系统基于码本的信道估计《Channel Estimation for Hybrid Architecture-Based Wideband Millimete.

其中，
F ~ ≜ F L ∗ ⊗ D u and x ≜ vec ⁡ ( Λ ) \tilde{\boldsymbol{F}} \triangleq \boldsymbol{F}_{L}^{*} \otimes \boldsymbol{D}_{u} \text { and } \boldsymbol{x} \triangleq \operatorname{vec}(\boldsymbol{\Lambda}) F~≜FL∗⊗Du and x≜vec(Λ)
最后
y = Φ x + ϵ \boldsymbol{y}=\boldsymbol{\Phi} \boldsymbol{x}+\boldsymbol{\epsilon} y=Φx+ϵ
其中， where Φ ≜ W v F ~ and W v ≜ [ w T ( 1 ) ⊗ v H ( 1 ) ⋮ w T ( T ) ⊗ v H ( T ) ] \begin{array}{l} \text { where } \boldsymbol{\Phi} \triangleq \boldsymbol{W}_{v} \tilde{\boldsymbol{F}} \text { and } \\ \qquad \boldsymbol{W}_{v} \triangleq\left[\begin{array}{c} \boldsymbol{w}^{T}(1) \otimes \boldsymbol{v}^{H}(1) \\ \vdots \\ \boldsymbol{w}^{T}(T) \otimes \boldsymbol{v}^{H}(T) \end{array}\right] \end{array} where Φ≜WvF~ and Wv≜⎣⎢⎡wT(1)⊗vH(1)⋮wT(T)⊗vH(T)⎦⎥⎤

写成这个形式，由于 x x x稀疏，就可以用压缩感知的经典算法，如OMP进行求解了。

请注意：这个算法的重大漏洞在于： Λ \Lambda Λ这个矩阵，应该是块稀疏的形式，而不是任意稀疏的！ 但OMP算法中是不会考虑的，也因此会有损失。什么是块稀疏呢？比如 BS-IRS有三条径， IRS-UE有三条径，那么cascaded 信道的秩是多少？还是3. 但是用压缩感知算法，会算出秩是9的矩阵。这就是因为， Λ \Lambda Λ矩阵其实可以证明，非零值必定集中在 3行（对应BS-IRS的径数）和 3列（对应IRS-UE的径数）中，而不是任意的9行9列！这个在 *Channel Estimation for Reconfigurable Intelligent Surface Aided Multi-User MIMO Systems * 一文中提到了。
文章传送门： http://arxiv.org/abs/1912.03619

仿真代码

clear;
N = 16; %transmit antennas
Mx = 8; % rows of IRS
My = 8;  % columns of IRS
M = Mx * My; %total elements of IRS
NG = 64;  % number of codewords of IRS
MG_x = 32;
MG_y = 32;
MG = MG_x * MG_y;
L = 3;  % paths of BS-IRS
Lprime = 3;  % paths of IRS-UE
T = 1500;FL = 1 / sqrt(N) * exp(1j * (0 : N-1)' * (-1 : 2 / NG : 1 - 2/NG));
Fx = 1 / sqrt(Mx) * exp(1j * (0 : Mx-1)' * (-1 : 2 / MG_x : 1 - 2/MG_x));
Fy = 1 / sqrt(My) * exp(1j * (0 : My-1)' * (-1 : 2 / MG_y : 1 - 2/MG_y));Fp = kron(Fx, Fy);
f = conj(Fp(:, 1));
Du = diag(f) * Fp;v = ones(M, 1);
F_hat = kron(conj(FL), Du);
W =[];G = 0;
for l = 1 : Lax = 1 / sqrt(Mx) * exp(1j * (0 : Mx-1)' * (2 * rand() - 1));ay = 1 / sqrt(My) * exp(1j * (0 : My-1)' * (2 * rand() - 1));ar = kron(ax, ay);at = 1 / sqrt(N) * exp(1j * (0 : N-1)' * (2 * rand() - 1));G = G +  sqrt(2)/2 * ar * at';
endhr = 0;
for l = 1 : Lprimeax = 1 / sqrt(Mx) * exp(1j * (0 : Mx-1)' * (2 * rand() - 1));ay = 1 / sqrt(My) * exp(1j * (0 : My-1)' * (2 * rand() - 1));ar = kron(ax, ay);hr = hr + sqrt(M / Lprime) *  sqrt(2)/2 * ar;
endH = diag(hr') * G;
n = 0;for t = 1 : Tw = 1 / sqrt(N) * exp(1j * (0 : N-1)' *  (2 * rand() - 1));v = 1 / sqrt(N) * exp(1j * (0 : M-1)' *  (2 * rand() - 1));W = [W; kron(w.', v')];y(t) = v' * H * w + n;
end
phi = W * F_hat;[x,res] = omp(y.' ,phi, L * Lprime);
Sigma = reshape(x, MG, NG);
H_est = Du * Sigma * FL';
NMSE = (norm(H_est - H, 'fro') / norm(H,'fro'))^2