学习笔记：OMP算法（Spatially Sparse Precoding in Millimeter wave MIMO Systems）

        自从开始研究毫米波通信之后，读了几篇经典论文，这其中就包括2014年这篇《Spatially Sparse Precoding in Millimeter wave MIMO Systems》。在把文章读的差不多之后，导师建议我自己尝试写一下这篇文章中的OMP算法。于是我在知乎、博客、书籍等渠道搜集了一下OMP算法的讲解，其中知乎上有一篇OMP讲解的十分透彻，用具体的数值将MP和OMP按照计算步骤进行对比，深入浅出，很适合新手入门，网址如下：https://zhuanlan.zhihu.com/p/52276805

        虽说论文中已经给出了OMP算法的伪代码，但是到自己上手的时候就犯了难：开始时会疑惑字典集如何选取，Fbb和Frf还有Fopt的矩阵维度未知等问题。搜集资料时发现matlab自带的OMP函数正是根据2014这篇论文来写的。在看matlab自带的OMP函数时恍然大悟：封装成函数后这些只需外部输入参数就好，而字典集选取是根据不同场景进行。读懂matlab自带的OMP函数也费了一番功夫，故通过本文记录学习代码的问题。

OMP算法就是简单来说就是给定压缩信号y和压缩矩阵A求解重构矩阵x。

论文中给出的OMP算法：

论文的优化目标是下式（1）

    （1）

算法可以分为以下几步：

1. 将信号Fopt赋值给压缩信号Fres

2.在字典集At中寻找对Fres贡献值最大的原子（可以通过矩阵相乘得出Φ=At’*Fres）。

3.根据Φ选出贡献值最大的原子，并将它放入新的压缩矩阵Frf，通过最小二乘法式（3）计算新压缩矩阵Frf对Fres的贡献值。

4.最小二乘法得到的贡献值将赋给重构信号矩阵Fbb。

5.删除所选矢量的贡献值，更新压缩信号Fres。

6.重复步骤1234，直到原压缩矩阵的所有列都被选完。

（2）

(2)的解为式(3)

（3）

其中

注：上式(2)和(2)的解式都摘自知乎文章https://zhuanlan.zhihu.com/p/52276805（可以把Fopt看做是y，Frf看做新的压缩矩阵Anew，Fbb看做λ。）

下面是matlab2019b版本工具箱中的OMP函数，博主给每句代码都加了中文注释，如有异见请在评论中写出，可以共同讨论一下。

function [Fbb,Frf] = helperOMPTransmitWeights(H,NtRF,Ns,At)
% helperOMPTransmitWeights is in support of
% MassiveMIMOHybridBeamformingExample. It may change in a future release.

%   Copyright 2017 The MathWorks, Inc.

%   Reference:
%   [1] Oma El Ayach, et al.,"Spatially Sparse Precoding in Millimeter wave
%   MIMO Systems", IEEE Transactions on Wireless Communications, Vol. 13,
%   No. 3, March 2014.

% use Comms convention
[~,~,v] = svd(H.');%正常的svd分解得到的是v'，要得到v只需求H'的svd分解，故此处是对H'进行svd分解。
Fopt = v(:,1:Ns);%最初发射信号维度是Ns*1，由于是转置后的v赋值给Fopt，所以是Ns列
Nt = size(H,1);%返回信道H的行数

% Algorithm 1: Spatially sparse precoding via Orthogonal Matching Pursuit
Frf = complex(zeros(Nt,NtRF));%生成Nt*NtRF维复数矩阵
Fres = Fopt;                  %将理想信号赋值给恢复矩阵
for m = 1:NtRF              
    Psi = At'*Fres;          %所有的原子贡献度，得到的Psi是一个列向量
    [~,k] = max(diag(Psi*Psi'));          %Psi*Psi'得到的对角线值是各个原子贡献值的平方，max输出列向量最大值的行号
    Frf(:,m) = At(:,k);                                % 将最大贡献值对应原子赋值给新的压缩矩阵
    Fbb = (Frf(:,1:m)'*Frf(:,1:m))\(Frf(:,1:m)'*Fopt); % 计算新压缩矩阵的最小二乘解
    temp = Fopt-Frf(:,1:m)*Fbb;
    Fres = temp/norm(temp,'fro');           %更新压缩信号
end
Fbb = sqrt(Ns)*Fbb/norm(Frf*Fbb,'fro');  %功率约束

% match H.'*Frf*Fbb*X.' to X*Fbb*Frf*H
Fbb = Fbb.';    %矩阵转置
Frf = Frf.';

% [EOF]