本文是 对 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 : L
ax = 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';
end
hr = 0;
for l = 1 : Lprime
ax = 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;
end
H = diag(hr') * G;
n = 0;
for t = 1 : T
w = 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
