【毫米波信道协方差估计算法】COMP/DSOMP/DCOMP原文分析

Spatial Channel Covariance Estimation for mmWave Hybrid MIMO Architecture[1]

说明：本篇为基于上文所做的笔记，可能存在理解错误，若存在疑虑问题请参照原文，敬请指正！

---

场景描述

• 利用接收机基带的低维信道投影进行高维信道估计和模拟预编码的设计。 首先，因为在数模混合波束赋形系统中，RF链路数量必少于天线数量，即基带接收到的导频信号经过模拟预编码变成更低维的空间信息，想要从这样的低维空间信息恢复出原来的信道情况十分困难。

• 减轻通信链路构建和维护所需的信道估计花销。 在某些场景中与其采用实时性较强的瞬时CSI，不如采取统计CSI；也就是说，在一个时隙中勉强能够完成或者极其影响该时隙内通信质量的工作，不如联合多个时隙完成，保证每个时隙中的通信质量，同时惠及接下来的几个时隙。

• 充分利用毫米波空间信道协方差矩阵的各种性质，来降低矩阵估计的复杂度。 毫米波空间信道协方差矩阵有着稀疏的性质，更具体来说，将该矩阵假设为厄米特稀疏矩阵更为一般。毫米波空间信道协方差矩阵中每个元素均对应着一个确切的AoA，而通过利用AoA的相对时间不变性能够显著提升估计的精确度。

• 化繁为简，一步能够解决就一步到位。 目前很多毫米波空间信道协方差估计办法还是需要在第一步估计信道矩阵，不管其算法内估计有多么粗略，本来也就不需要显式地估计信道。

---

信道及系统建模

TDD系统中，假设基站配有$N_{\mathrm{A}\mathrm{N}\mathrm{T}}$副天线，$N_{\mathrm{R}\mathrm{F}}$条RF链路，它通过估计UL信道协方差矩阵用于DL混合预编码设计；移动站只配有单天线；信道存在$L$条路径，且AoA在$T$个时隙内保持不变。

$t$时隙的UL信道建模如下
$${\mathbf{h}}_t={\mathbf{A}\mathbf{g}}_t,t=1,2,\cdots ,T\text{\hspace{1em}(1)}$$

${\mathbf{g}}_t\in {\mathbb{C}}^{L\times 1}$ — 信道路径增益向量，设为i.i.d.高斯随机向量，即$\mathcal{C}\mathcal{N}(0,{\mathbf{R}}_g)$.

$\mathbf{A}\in {\mathbb{C}}^{N_{\mathrm{A}\mathrm{N}\mathrm{T}}\times L}$ — $L$个阵列响应向量组成的矩阵，可表示为$\mathbf{A}=[\mathbf{a}({\theta }_1)\cdots \mathbf{a}({\theta }_L)]$，其中${\theta }_l$为信道路径$l$对应的AoA.

假设训练信号为$s_t=1$, 那么$t$时隙的接收信号有
$${\mathbf{x}}_t={\mathbf{A}\mathbf{g}}_t+{\mathbf{n}}_t,t=1,2,\cdots ,T.\text{\hspace{1em}(2)}$$

${\mathbf{n}}_t\sim \mathcal{C}\mathcal{N}(0,{\sigma }^2\mathbf{I}).$

将$(2)$的所有时隙整合为矩阵形式有
$$\mathbf{X}=\mathbf{A}\mathbf{G}+\mathbf{N}\text{\hspace{1em}(3)}$$
假设接收时采取的模拟组合编码矩阵为$\mathbf{W}\in {\mathbb{C}}^{N_{\mathrm{R}\mathrm{F}}\times N_{\mathrm{A}\mathrm{N}\mathrm{T}}}$, 则基带接收信号可表示为
$$\mathbf{Y}=\mathbf{W}\mathbf{X}=\mathbf{W}\mathbf{A}\mathbf{G}+\mathbf{W}\mathbf{N}\text{\hspace{1em}(4)}$$
基带接收信号是我们实际可获得的，用以计算空间信道协方差矩阵，即${\mathbf{R}}_h=\mathbb{E}[{\mathbf{h}}_t{\mathbf{h}}_t^{\ast }]={\mathbf{A}\mathbf{R}}_g{\mathbf{A}}^{\ast }$.

---

SVM基线方法——OMP

从式$(4)$中拆分出$t$时隙的基带接收信号可得
$${\mathbf{y}}_t=\mathbf{W}\tilde{\mathbf{A}}{\tilde{\mathbf{g}}}_t+{\mathbf{W}\mathbf{n}}_t,t=1,2,\cdots ,T.\text{\hspace{1em}(5)}$$

$\tilde{\mathbf{A}}\in {\mathbb{C}}^{N_{\mathrm{A}\mathrm{N}\mathrm{T}}\times D}$ — 码典矩阵，列向量为预设AoA集合对应的阵列响应向量；

${\tilde{\mathbf{g}}}_t\in {\mathbb{C}}^{D\times 1}$ — 信道路径增益向量，其中仅有$L$个非零元素.

此时，设$\mathbf{\Phi }=\mathbf{W}\tilde{\mathbf{A}}$作为感知矩阵，最优化问题可设置为
$$\underset{{\tilde{\mathbf{g}}}_t}{\min }{\Vert {\tilde{\mathbf{g}}}_t\Vert }_0\mathrm{s}.\mathrm{t}.{\Vert {\mathbf{y}}_t-\mathbf{\Phi }{\tilde{\mathbf{g}}}_t\Vert }_2\le ϵ\text{\hspace{1em}(6)}$$

$||\cdot |{|}_0$，向量零范数，指代向量内非零元素个数。

上述问题的意义在于保证恢复${\tilde{\mathbf{g}}}_t$尽量稀疏的条件下，使得恢复信号重构后与接收信号的残差小于某个阈值。提到残差和稀疏信号恢复，当然此问题可通过OMP算法求解。对于每个时隙$t$均需估计${\tilde{\mathbf{g}}}_t$, 而后通过计算下式获取信道样本协方差矩阵
$${\hat{\mathbf{R}}}_{\mathbf{h}}=\tilde{\mathbf{A}}\left(\frac{1}{T}\sum _{t=1}^T{\tilde{\mathbf{g}}}_t{\tilde{\mathbf{g}}}_t^{\ast }\right){\tilde{\mathbf{A}}}^{\ast }\text{\hspace{1em}(7)}$$
这种单时隙SVM方法缺点也正在逐个求解${\tilde{\mathbf{g}}}_t$，没有联合多个时隙进行估计；解决这个问题的一个可行方法就是使用向量化技术一次就估计出所有信道向量。再将式$(5)$中的信号模型整合可得
$$\mathbf{Y}=\mathbf{\Phi }\tilde{\mathbf{G}}+\mathbf{W}\mathbf{N},\text{\hspace{1em}(8)}$$
对式$(8)$进行向量化可得
$$\mathrm{v}\mathrm{e}\mathrm{c}(\mathbf{Y})=({\mathbf{I}}_T\otimes \mathbf{\Phi })\mathrm{v}\mathrm{e}\mathrm{c}(\tilde{\mathbf{G}})+\mathrm{v}\mathrm{e}\mathrm{c}(\mathbf{W}\mathbf{N}),\text{\hspace{1em}(9)}$$
那么最优化问题可如下描述
$$\underset{\mathrm{v}\mathrm{e}\mathrm{c}(\tilde{\mathbf{G}})}{\min }{\Vert \mathrm{v}\mathrm{e}\mathrm{c}(\tilde{\mathbf{G}})\Vert }_0\mathrm{s}.\mathrm{t}.{\Vert \mathrm{v}\mathrm{e}\mathrm{c}(\mathbf{Y})-({\mathbf{I}}_T\otimes \mathbf{\Phi })\mathrm{v}\mathrm{e}\mathrm{c}(\tilde{\mathbf{G}})\Vert }_2\le ϵ.\text{\hspace{1em}(10)}$$
以上问题的解与问题$(6)$相比有着更长的长度，但同为向量，也可使用OMP算法进行求解，亦可用同样方法获取信道样本协方差矩阵。然而，尽管进行了联合求解，与单时隙SVM方法相比并没有什么特别的增益，这主要是因为在根本上没有利用$\tilde{\mathbf{G}}$的特征。我们之前假设在$T$个时隙中AoA并不发生变化，也就是说$\tilde{\mathbf{G}}$具有行稀疏的特性。

---

MMV基线方法——SOMP

使用MMV方法以利用$\tilde{\mathbf{G}}$的行稀疏性，可将问题建模为
$$\underset{\tilde{\mathbf{G}}}{\min }{\Vert \tilde{\mathbf{G}}\Vert }_{\mathrm{r}\mathrm{o}\mathrm{w},0}\mathrm{s}.\mathrm{t}.{\Vert \mathbf{Y}-\mathbf{\Phi }\tilde{\mathbf{G}}\Vert }_F\le ϵ.\text{\hspace{1em}(11)}$$
其中，$\Vert \tilde{\mathbf{G}}{\Vert }_{\mathrm{r}\mathrm{o}\mathrm{w},0}$意为行稀疏度。这种解为矩阵的最优化问题可以采用SOMP算法。在求解后，依然可以使用同样的方法求出信道样本协方差矩阵。可以看出这里利用了矩阵行稀疏的性质，如图1(b)所示。

以上SVM方法和MMV方法均作为本文中的基准方法与之后提出的算法进行对比。以上算法均采用两步走的方法，先估计信道向量，之后计算协方差矩阵，显得多此一举；另外，均没有利用空间信道协方差矩阵的厄米特性质。

---

COMP——利用协方差矩阵表现出来的稀疏厄米特矩阵性质

经过对比，我们需要一种空间信道协方差矩阵的估计方法，尽量简化步骤，一步求解；能够利用如图1(d)所示的稀疏厄米特矩阵的性质。

首先，我们发现以上方法均通过${\mathbf{y}}_t$作为测试样本去恢复${\tilde{\mathbf{g}}}_t$，而后计算获取${\mathbf{R}}_{\tilde{\mathbf{g}}}$和${\mathbf{R}}_{\mathbf{h}}$，这也是使用两步的原因，可不可以通过某种方法避免这个中间过程，直接获取${\mathbf{R}}_{\tilde{\mathbf{g}}}$呢？

先探究${\mathbf{R}}_{\mathbf{y}}$与${\mathbf{R}}_{\tilde{\mathbf{g}}}$的关系，如下
$${\mathbf{R}}_{\mathbf{y}}={\mathbf{\Phi }\mathbf{R}}_{\tilde{\mathbf{g}}}{\mathbf{\Phi }}^{\ast }+{\sigma }^2{\mathbf{W}\mathbf{W}}^{\ast },\text{\hspace{1em}(12)}$$
在假设${\mathbf{R}}_{\tilde{\mathbf{g}}}$具有稀疏厄米特矩阵性质之前，先看看图1(c)这种情况，即${\mathbf{R}}_{\tilde{\mathbf{g}}}$是一个只有几个非零对角元素的稀疏对角矩阵。这种情况对应着信道路径增益之间相互独立，是一种较为理想的情况。类似式$(9)$采用向量化的方法，有
$$\mathrm{v}\mathrm{e}\mathrm{c}({\mathbf{R}}_{\mathbf{y}})=({\mathbf{\Phi }}^C\circ \mathbf{\Phi })\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mathbf{R}}_{\tilde{\mathbf{g}}})+{\sigma }^2\mathrm{v}\mathrm{e}\mathrm{c}({\mathbf{W}\mathbf{W}}^{\ast }),\text{\hspace{1em}(13)}$$
其中，$\circ$是Khatri-Rao积，相当于基于列的克罗内克积，$\mathbf{\Phi }\in {\mathbb{C}}^{N_{\mathrm{R}\mathrm{F}}\times D}，{\mathbf{\Phi }}^C\circ \mathbf{\Phi }\in {\mathbb{C}}^{N_{\mathrm{R}\mathrm{F}}{N}_{\mathrm{R}\mathrm{F}}\times D}$；$\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mathbf{R}}_{\tilde{\mathbf{g}}})$为${\mathbf{R}}_{\tilde{\mathbf{g}}}$对角线元素组成的列向量。既然可以如上向量化表示，也就说明这种情况下也可使用OMP算法对$\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mathbf{R}}_{\tilde{\mathbf{g}}})$进行求解。

当然，对于这种对角稀疏的假设，显然信道路径增益之间存在相关性则以上方法不适用；在样本并不丰富，$T$较小的情况下，也不能保证最终计算出的样本协方差矩阵对角线以外元素全部趋近于0.

假设${\mathbf{R}}_{\tilde{\mathbf{g}}}$为更一般的稀疏厄米特矩阵，因为没有较为合适的向量表示形式，最优化问题可表示为矩阵形式，通过MMV方法求解
$$\underset{{\mathbf{R}}_{\tilde{\mathbf{g}}}}{\min }{\Vert {\mathbf{R}}_{\tilde{\mathbf{g}}}\Vert }_{\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{i}\mathrm{c}\mathrm{e},0}\mathrm{s}.\mathrm{t}.{\Vert {\mathbf{R}}_{\mathbf{Y}}-\mathbf{\Phi }{\mathbf{R}}_{\tilde{\mathbf{g}}}{\mathbf{\Phi }}^{\ast }\Vert }_F\le ϵ,\text{\hspace{1em}(14)}$$
其中，${\Vert {\mathbf{R}}_{\tilde{\mathbf{g}}}\Vert }_{\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{i}\mathrm{c}\mathrm{e},0}$表示${\mathbf{R}}_{\tilde{\mathbf{g}}}$中元素分布在多少个行与列（行列标号相同算作一个），则该最优化问题旨在残差满足收敛条件的情况下，让元素尽量分布在几行几列的交叉位置上，行列数尽可能少且相同（协方差矩阵的对称性）。

求解该问题，我们采用COMP方法，算法如下

Algorithm 1 Covariance OMP

---

​ Input: $\mathbf{\Phi },\mathbf{Y},L$

​ Initialization: ${\mathbf{R}}_{\mathbf{y}}=\frac{1}{T}{\mathbf{Y}\mathbf{Y}}^{\ast },\mathbf{V}={\mathbf{R}}_{\mathbf{y}},\mathcal{S}=\varnothing$

​ for $n=1:L$ do

​ $j=\arg {\max }_i|[\mathbf{\Phi }{]}_{:,i}^{\ast }\mathbf{V}[\mathbf{\Phi }{]}_{:,i}|$

​ S = S ∪ { j } \mathcal{S}=\mathcal{S}\cup\{j\} S=S∪{j}

​ ${\mathbf{R}}_{\tilde{\mathbf{g}}}=[\mathbf{\Phi }{]}_{:,\mathcal{S}}^{†}{\mathbf{R}}_{\mathbf{y}}([\mathbf{\Phi }{]}_{:,\mathcal{S}}^{†}{)}^{\ast }$

​ $\mathbf{V}={\mathbf{R}}_{\mathbf{y}}-[\mathbf{\Phi }{]}_{:,\mathcal{S}}{\mathbf{R}}_{\tilde{\mathbf{g}}}[\mathbf{\Phi }{]}_{:,\mathcal{S}}^{\ast }$

​ end for

​ Output: $\mathcal{S},{\mathbf{R}}_{\tilde{\mathbf{g}}}$

如图2(a)所示，信道及接收端模拟组合编码将信道路径增益组合成基带接收信号，通过对接收信号求解信号样本协方差矩阵，经过COMP算法获得${\mathbf{R}}_{\tilde{\mathbf{g}}}$，最终求得${\mathbf{R}}_{\mathbf{h}}$.

---

DSOMP&DCOMP——基于时变模拟组合编码进行的信道协方差估计

以上压缩感知算法有一个主要限制，即使噪声较弱，为完美恢复信号同样需要许多检测。SMV方法中，完美恢复信号所需检测量为$N_{\mathrm{R}\mathrm{F}}=\mathcal{O}(L\log \frac{D}{L})$，这需要足够的RF链路来支持算法。同时MMV方法能够通过多个时隙的检测补偿SMV的劣势，但是因为感知矩阵$\mathbf{\Phi }$不变，使得这种多时隙联合估计的方法不是特别的有效。

我们尝试让感知矩阵$\mathbf{\Phi }$时变起来，采取时变的模拟组合矩阵$\mathbf{W}$，但并不是每个时隙均发生变化。则如下
$$\left[\begin{array}{c}{\mathbf{y}}_1\\ ⋮\\ {\mathbf{y}}_T\end{array}\right]=\left[\begin{array}{c}{\mathbf{W}}_1\tilde{\mathbf{A}}\\ ⋮\\ {\mathbf{W}}_T\tilde{\mathbf{A}}\end{array}\right]\tilde{\mathbf{g}}=\left[\begin{array}{c}{\mathbf{\Phi }}_1\\ ⋮\\ {\mathbf{\Phi }}_T\end{array}\right]\tilde{\mathbf{g}}\text{\hspace{1em}(15)}$$
其中假设$\mathbf{g}$在估计过程中为常量，即信道为静态信道。

如果以片段为视角，也就是片段$s$聚合$T_s$个连续时隙，而对于每个片段改变模拟组合编码矩阵。那么对于片段$s$，$\mathbf{y}$的估计样本协方差为
$${\hat{\mathbf{R}}}_{\mathbf{y},s}=\frac{1}{T_s}\sum _{t=1}^{T_s}{\mathbf{y}}_{s,t}{\mathbf{y}}_{s,t}^{\ast }={\mathbf{W}}_s{\hat{\mathbf{R}}}_{\mathbf{h},s}{\mathbf{W}}_s^{\ast }+{\mathbf{W}}_s{\hat{\mathbf{R}}}_{\mathbf{n},s}{\mathbf{W}}_s^{\ast }\text{\hspace{1em}(16)}$$
当然，上式我们不能想当然地认为噪声的样本协方差就可以近似为${\sigma }^2\mathbf{I}$，同样也不能将片段$s$的信道样本协方差近似为该相干时间内的信道协方差矩阵，这都需要$T_s$足够大。那么如果$T_s$足够大，式$(16)$可以进行近似，近似经过向量化及整合后可得
$$\left[\begin{array}{c}\mathrm{v}\mathrm{e}\mathrm{c}({\hat{\mathbf{R}}}_{\mathbf{y},1})\\ ⋮\\ \mathrm{v}\mathrm{e}\mathrm{c}({\hat{\mathbf{R}}}_{\mathbf{y},S})\end{array}\right]\approx \left[\begin{array}{c}({\mathbf{W}}_1^C\otimes {\mathbf{W}}_1)\\ ⋮\\ ({\mathbf{W}}_S^C\otimes {\mathbf{W}}_S)\end{array}\right]\mathrm{v}\mathrm{e}\mathrm{c}({\mathbf{R}}_{\mathbf{h}})+{\sigma }^2{\mathbf{n}}_{\mathrm{e}\mathrm{f}\mathrm{f}}\text{\hspace{1em}(17)}$$
其中，${\mathbf{n}}_{\mathrm{e}\mathrm{f}\mathrm{f}}=[{\mathrm{v}\mathrm{e}\mathrm{c}}^T({\mathbf{W}}_1{\mathbf{W}}_1^{\ast })\cdots {\mathrm{v}\mathrm{e}\mathrm{c}}^T({\mathbf{W}}_S{\mathbf{W}}_S^{\ast }){]}^T$.

对于式$(17)$，若要通过最小二乘法得出$\mathrm{v}\mathrm{e}\mathrm{c}({\mathbf{R}}_{\mathbf{h}})$，其左乘矩阵需要保证行数大于列数，即$N_{\mathrm{R}\mathrm{F}}^{2}S\ge N_{\mathrm{A}\mathrm{N}\mathrm{T}}^2$，可见RF链路数量小于天线数量也是可以的。但是如果真如上述方法求解信道协方差矩阵，首先，没有利用毫米波信道的稀疏性，其中最小二乘法需要求解大型矩阵的逆，运算复杂度很高；其次，如果$T_s$没有足够大，无法近似。也就是说模拟组合编码在长时间需要固定。

再重新描述我们的目标，我们拟进行基于压缩感知的信道协方差估计，其中使用时变模拟组合编码，信道也是时变的。假设在每个时隙改变模拟组合编码以最大化检测的扩展效果，这可以应用在SOMP上，亦可应用在COMP上。

首先，重温SOMP。上文式$(11)$中描述问题可以等效为
$$\underset{{\mathbf{g}}_1,\cdots ,{\mathbf{g}}_T}{\min }{\Vert \sum _{t=1}^T{\left|{\tilde{\mathbf{g}}}_t\right|}_{\mathrm{e}\mathrm{l}\mathrm{e}}\Vert }_0\mathrm{s}.\mathrm{t}.\sum _{t=1}^T{\Vert {\mathbf{y}}_t-\mathbf{\Phi }{\tilde{\mathbf{g}}}_t\Vert }_2^2\le ϵ,\text{\hspace{1em}(18)}$$
其中，$|\mathbf{a}{|}_{\mathrm{e}\mathrm{l}\mathrm{e}}$定义为元素为向量$\mathbf{a}$内元素绝对值的向量。上式可理解为想要估计一组信道增益向量，使所有时隙的增益向量具有相同的稀疏结构（体现了对毫米波信道稀疏性的利用），即在相同位置上存在非零元素。进一步，考虑时变组合编码矩阵，可以将$\mathbf{\Phi }$替换为${\mathbf{\Phi }}_t$, 这也反映了感知矩阵是时变的，DSOMP可阐述为以下算法：

Algorithm 2 Dynamic SOMP(DSOMP)

---

Input: 时变感知矩阵 { Φ 1 , ⋯   , Φ T } \{\mathbf{\Phi}_1,\cdots,\mathbf{\Phi}_T\} {Φ1​,⋯,ΦT​}, $T$时隙基带接收信号矩阵$\mathbf{Y}$, 信道主要路径数量$L$.

Initialization: 残差矩阵$\mathbf{V}=\mathbf{Y}$, 选取组合编码向量编号集合$\mathcal{S}=\varnothing$

for $n=1:L$ do
$$j=\underset{i}{\mathrm{argmax}}\sum _{t=1}^T\left|{\left[{\mathbf{\Phi }}_t^{\ast }{\left[\mathbf{V}\right]}_{:,t}\right]}_i\right|$$

• 这里$i$实际指代选取每个感知矩阵的第$i$列，但是实际上除非感知矩阵的码本相当庞大，这种时变感知矩阵才能近似全局最优的。 将这个中间变量$i$存入$j$.

S = S ∪ { j } g t = [ Φ t ] : , S † [ Y ] : , t , ∀ t \begin{aligned} \mathcal{S}&=\mathcal{S}\cup\{j\}\\ \mathbf{g}_t&=\left[\mathbf{\Phi}_t\right]_{:,\mathcal{S}}^\dagger\left[\mathbf{Y}\right]_{:,t},\quad\forall t \end{aligned} Sgt​​=S∪{j}=[Φt​]:,S†​[Y]:,t​,∀t​

• 将刚才得到相关性最强的码字编号存入集合$\mathcal{S}$, 作者在这里没有明确解释，如果单从该式看的话，很容易理解为最终结果${\mathbf{g}}_t\in {\mathbb{C}}^{|\mathcal{S}|\times 1}$，但放到下一步显然不适合，最好理解为这里对于集合内存在的几列进行计算，最终在其余行对应增益向量的元素补零。

$$\mathbf{V}=\mathbf{Y}-\left[{\mathbf{\Phi }}_1{\mathbf{g}}_1\cdots {\mathbf{\Phi }}_T{\mathbf{g}}_T\right]$$

• 根据以上更新的各时隙的增益向量得出新残差，可以得知每进行一次循环至少能够找到一个主要路径对应的信道增益。

end for

Output: $\mathcal{S},\tilde{\mathbf{G}}=[{\mathbf{g}}_1\cdots {\mathbf{g}}_T]$.

同时，将时变组合编码应用于COMP，进行算法变型。之前我们将基带接收信号样本协方差矩阵是一个信道增益向量样本协方差矩阵的函数，如下表示
$${\hat{\mathbf{R}}}_{\mathbf{y}}=\frac{1}{T}\sum _{t=1}^T{\mathbf{y}}_t{\mathbf{y}}_t^{\ast }=\mathbf{\Phi }\left(\frac{1}{T}\sum _{t=1}^T{\tilde{\mathbf{g}}}_t{\tilde{\mathbf{g}}}_t^{\ast }\right){\mathbf{\Phi }}^{\ast }=\mathbf{\Phi }{\hat{\mathbf{R}}}_{\tilde{\mathbf{g}}}{\mathbf{\Phi }}^{\ast },\text{\hspace{1em}(19)}$$
但是如果感知矩阵随时间变化，我们不能够想当然的如上表示，因为
$${\hat{\mathbf{R}}}_{\mathbf{y}}=\frac{1}{T}\sum _{t=1}^T{\mathbf{y}}_t{\mathbf{y}}_t^{\ast }=\frac{1}{T}\sum _{t=1}^T{\mathbf{\Phi }}_t{\tilde{\mathbf{g}}}_t{\tilde{\mathbf{g}}}_t^{\ast }{\mathbf{\Phi }}_t^{\ast }.\text{\hspace{1em}(20)}$$

那么，我们就无法单纯用$\mathbf{y}$的样本协方差矩阵直接获取${\hat{\mathbf{R}}}_{\tilde{\mathbf{g}}}$, 取而代之，我们用单时隙${\mathbf{y}}_t{\mathbf{y}}_t^{\ast }$, 极端下可以视作单时隙的样本协方差。注意，即使在信道路径不相干的情况下，单个时隙中${\mathbf{y}}_t{\mathbf{y}}_t^{\ast }$不可视作对角矩阵。 而利用稀疏厄米特矩阵的性质，所有时隙的${\mathbf{y}}_t{\mathbf{y}}_t^{\ast }$都是稀疏厄米特矩阵，有着相同非零元素位置。那么有如下DCOMP算法：

Algorithm 3 Dynamic COMP(DCOMP)

---

Input: 时变感知矩阵 { Φ 1 , ⋯   , Φ T } \{\mathbf{\Phi}_1,\cdots,\mathbf{\Phi}_T\} {Φ1​,⋯,ΦT​}, $T$时隙基带接收信号矩阵$\mathbf{Y}$, 信道主要路径数量$L$.

Initialization: 单时隙样本协方差矩阵${\mathbf{R}}_{\mathbf{y},t}=[\mathbf{Y}{]}_{:,t}[\mathbf{Y}{]}_{:,t}^{\ast }$, 该时隙残差${\mathbf{V}}_t={\mathbf{R}}_{\mathbf{y},t},\forall t$, 码字编号集合$\mathcal{S}=\varnothing$.

for $n=1:L$ do
j = arg ⁡ max ⁡ i ∑ t = 1 T ∣ [ Φ t ] : , i ∗ V t [ Φ t ] : , i ∣ , S = S ∪ { j } , R g , t = [ Φ t ] : , S † R y , t ( [ Φ t ] : , S † ) ∗ , ∀ t V t = R y , t − [ Φ t ] : , S R g , t [ Φ t ] : , S ∗ , ∀ t \begin{aligned} j&=\underset{i}{\arg\max}\sum_{t=1}^T\left|\left[\mathbf\Phi_t\right]_{:,i}^*\mathbf{V}_t\left[\mathbf\Phi_t\right]_{:,i}\right|,\\ \mathcal{S}&=\mathcal{S}\cup\{j\},\\ \mathbf{R}_{\mathbf{g},t}&=\left[\mathbf\Phi_t\right]_{:,\mathcal{S}}^\dagger\mathbf{R}_{\mathbf{y},t}\left(\left[\mathbf\Phi_t\right]_{:,\mathcal{S}}^\dagger\right)^*,\quad\forall t\\ \mathbf{V}_t&=\mathbf{R}_{\mathbf{y},t}-\left[\mathbf\Phi_t\right]_{:,\mathcal{S}}\mathbf{R}_{\mathbf{g},t}\left[\mathbf\Phi_t\right]_{:,\mathcal{S}}^*,\quad\forall t \end{aligned} jSRg,t​Vt​​=iargmax​t=1∑T​∣∣∣​[Φt​]:,i∗​Vt​[Φt​]:,i​∣∣∣​,=S∪{j},=[Φt​]:,S†​Ry,t​([Φt​]:,S†​)∗,∀t=Ry,t​−[Φt​]:,S​Rg,t​[Φt​]:,S∗​,∀t​
end for

Output: $\mathcal{S},\frac{1}{T}{\sum }_{t=1}^T{\mathbf{R}}_{\mathbf{g},t}$.

之前说过，使用时变感知矩阵能够高效地进行压缩感知，也就是说所需时隙总数$T$的数目可以随之降低，也就是说可以使用更短的时隙来估计协方差矩阵。

---

仿真结果

性能对比指标设定为最佳$p$维波束形成矩阵效率
$$\eta =\frac{\mathrm{T}\mathrm{r}\left({\mathbf{U}}_{{\hat{\mathbf{R}}}_{\mathbf{h}}}^{\ast }{\mathbf{R}}_{\mathbf{h}}{\mathbf{U}}_{{\hat{\mathbf{R}}}_{\mathbf{h}}}\right)}{\mathrm{T}\mathrm{r}\left({\mathbf{U}}_{{\mathbf{R}}_{\mathbf{h}}}^{\ast }{\mathbf{R}}_{\mathbf{h}}{\mathbf{U}}_{{\mathbf{R}}_{\mathbf{h}}}\right)}$$
其中，${\mathbf{U}}_{{\hat{\mathbf{R}}}_{\mathbf{h}}}$为估计协方差${\hat{\mathbf{R}}}_{\mathbf{h}}$特征向量以列填充的矩阵，${\mathbf{U}}_{{\mathbf{R}}_{\mathbf{h}}}$为理想协方差${\mathbf{R}}_{\mathbf{h}}$特征向量以列填充的矩阵。该指标有$0\le \eta \le 1$, 且越大估计越精确。

从图3可以看出，在RF链路充足的情况下，DCOMP算法在较短时间内可以较好的估计协方差。OMP和SOMP作为基准算法，可以看出本文提出的COMP明显优于二者。当时隙长度足够的情况下，DCOMP和DSOMP算法性能相近，均优于感知矩阵固定的COMP。

从图4和图3的差别可见，在RF链路并不充分的情况下DCOMP算法明显优于DSOMP，SOMP竟劣于OMP。且该情况下时变感知矩阵带来的增益就更加明显。

---

References

[1] S. Park and R. W. Heath, “Spatial channel covariance estimation for mmWave hybrid MIMO architecture,” 2016 50th Asilomar Conference on Signals, Systems and Computers, 2016, pp. 1424-1428, doi: 10.1109/ACSSC.2016.7869611.