稀疏信号恢复算法SOMP（只给出算法，不给实验）

通过看了一篇LEO的信道估计文章，找到了这篇用于多测量的OMP算法，也叫SOMP。当然有些人实现它通过多加一层for循环，把多测量问题转换成多次单测量。缺点就是要知道稀疏度，当然在有噪声的情况下，也有该算法的变体。我们在这里就不再介绍，当然论文在代码里有注释。

function [Omega,X_hat] = MMV_SOMP(Y, PHI, s)
    % SOMP:同时正交匹配追踪 simultaneous orthogonal matching pursuit
    %         论文：J. Determe, J. Louveaux, L. Jacques and
    %         F. Horlin, \On the Noise Robustness of Simultaneous Orthogonal Matching Pursuit,"
    %         in IEEE Transactions on Signal Processing, vol. 65, no. 4, pp. 864-875, Feb. 15
    %         2017. doi: 10.1109/TSP.2016.2626244. http://ieeexplore.ieee.org/document/7738592/
    % write:zts
    % date：23,11,17
    % version:1.0
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % 参数：
    % Y: 观测矩阵，每列是一个观测向量
    % PHI: 字典矩阵
    % s: 稀疏度
    
    % 获取矩阵大小 
    [~, N] = size(PHI);
    [~, K] = size(Y);
	% 1 初始化残差和索引
    r = Y;
	% 初始化索引集
    Omega = [];
    % 初始化重构矩阵
    X_hat=zeros(N,K);
	% 2 3
    for t = 1:s
	
		% 初始化原子投影
		proj=[];
		% 4
		for j=1:N
			proj = [proj,norm(r'*PHI(:,j),1)];
		end
        % 选取最大投影对应的原子索引
        [~, index] = max(proj);
        Omega = [Omega, index];
        % 重构X的每行元素
        for k=1:K
             X_hat(Omega,k)=pinv(PHI(:,Omega))*Y(:,k);
        end
		%更新测量向量投影
		Y_tplus1 = PHI(:,Omega)*pinv(PHI(:,Omega))*Y;
        %更新残差
        r = Y - Y_tplus1;

    end
    Omega = sort(Omega);
end

下面是多次单测量以达到同样目的的OMP

 function[x_omp]=	MMV_OMP(y,H,L_seq,SP)

     [N,K]=size(H);                             %传感矩阵H为N*K矩阵
     L=L_seq;                                   %时隙长度
     

     x_omp=zeros(K,L);                          %用来存储恢复的x(K×L的矩阵)
     for i=1:L_seq                              %每次循环恢复x的一个时隙
         x_omp_col=zeros(K,1);                      %存放每个时隙恢复的x的列
         At=zeros(N,SP);                            %用来存储迭代过程中H被选择的列，最多SP列
         Pos_theta=zeros(1,SP);                     %用来迭代过程中存储H被选择的列序号
         y_col=y(:,i);                              %取对应的y的列
         r_n=y_col;                                 %初始化残差(residual)为y的列
         
        for ii=1:SP                                  %迭代S次，S为输入参数

             product=H'*r_n;                       %传感矩阵H各列与残差的内积，得到N行1列
             
            [val,pos]=max(abs(product));           %找到最大内积绝对值，即与残差最相关的列
                                                   %val得到值，pos得到列序号
                                                   
             At(:,ii)=H(:,pos);                     %存储这一列

             Pos_theta(ii)=pos;                     %存储这一列的序号

             %H(:,pos)=zeros(N,1);                   %清零H的这一列，其实此行可以不要，因为它与残差正交

             theta_ls=(At(:,1:ii)'*At(:,1:ii))^(-1)*At(:,1:ii)'*y_col;       %最小二乘解pinv(At(:,1:ii))*y_row

             r_n=y_col-At(:,1:ii)*theta_ls;             %更新残差

        end

         x_omp_col(Pos_theta)=theta_ls;                   %恢复出的theta

         x_omp(:,i)=x_omp_col;
     end
                                               

 end