稀疏信号恢复SAMP（只给出自己复现的代码）

很多人已经公布过这个算法代码，不过还是再发一遍，针对多测量模型，用多次单测量的结果来实现。师兄复现的代码都没有论文标注。仅仅复现，具体实现还要看自己的系统是否符合条件。也有可能是我们的代码出现错误，毕竟调代码就是个麻烦活，谁也不想干。

function [x_SAMP] = SAMP(y,H,step_size,L_seq,sigma2)
% parameter usage:
%   y: MxL observation vector      测量值
%   H: NxK measurement matrix    感知矩阵，H
%   step_size: any positive integer value not larger than sparsity    初始步长（1-K）
%   sigma2: noise energy when sensing   重构误差
%   xr: reconstructed sparse signal   重构的稀疏表示
%   iter_num: number of iterations    迭代次数
%   步长=1
% Initialization
[N,K]=size(H);                             %传感矩阵H为N*K矩阵
x_SAMP = zeros(K,L_seq);
for i=1:L_seq
    active_set = [];          %支撑集
    actset_size = step_size;  %步长=1
    iter_num = 0;             %迭代次数（没用上）
    stg_idx = 1;              % stage index   阶段1 分阶段进行的
    y_col=y(:,i);
    res=y_col; 
    x_SAMP_col=zeros(K,1);
    while (norm(res)^2/N>sigma2)
        % candidate list   候选集
        [~, idx1] = sort(abs(H'*res), 'descend');       %计算内积并进行降序排序；并记录其对应的列号           
        candidate_set = [active_set; idx1(1:actset_size)];   %候选集=上一次的支撑集+此次内积的 前 step_size 个对应的列号 ，（即与y相关度最高的step_size个对应的H的列号） 
        if length(candidate_set)>=N   %候选集合大于了观测值，跳出循环
            break;
        end
        % finalist
        [~, idx] = sort(abs(pinv(H(:,candidate_set))*y_col), 'descend');   %候选集中的列进行最小二乘（即求X的稀疏估计 X^），并进行降序排序；
        new_active_set = candidate_set(idx(1:actset_size));                %选出X^的前actset_size(当前步长)个作为支撑集
        new_res = y_col-H(:,new_active_set)*pinv(H(:,new_active_set))*y_col;  % 选出的支撑集 最小二乘进行新的估计并更新残差
        %new_res = y_col-pinv(H(:,new_active_set))*y_col;  % 选出的支撑集 最小二乘进行新的估计并更新残差
        if (norm(new_res) >= norm(res))   
            % shift into a new stage
            stg_idx = stg_idx + 1;
            actset_size = stg_idx*step_size;
        else          
            % update residual and active set    
         res = new_res;
            active_set= new_active_set;
        end 
        iter_num = iter_num +1;%阶段计数
        
    end % loop
    % reconstruction    
    xr_active_set = pinv(H(:,active_set))*y_col;              %选好最终的支撑集进行最小二乘，得到原始信号的估计
    x_SAMP_col(active_set) = xr_active_set;
    x_SAMP(:,i)=x_SAMP_col;
end
end