关于SAMP算法的实现,大多数采用的是莱斯大学的工具箱中的samp或者CSDN彬彬有礼中的CS_SAMP算法。
rice大学的工具箱中的samp算法实现如下:
function [xr iter_num] = SAMP(y, Phi, step_size, sigma);
% SAMP: Sparsity Adaptive Matching Pursuit algoritm for compressed sensing.
% For theoretical analysis, please refer to the paper :
% Thong. T. Do, Lu Gan and Trac D. Tran ,"Sparsity Adaptive Matching
% Purusit for practical compressed sensing" available at http://dsp.ece.rice.edu/cs
% Written by Thong Do(thongdo@jhu.edu)
% Updated on July, 26th 2008
% parameter usage:
% y: Mx1 observation vector 表测量值
% Phi: MxN measurement matrix 感知矩阵,H
% step_size: any positive integer value not larger than sparsity 初始步长(0,k)
% sigma: noise energy when sensing 重构误差
% xr: reconstructed sparse signal 重构的稀疏表示
% iter_num: number of iterations 迭代次数
% Initialization
iter_num = 0; %迭代次数
actset_size = step_size; %步长
active_set = []; %支撑集
res = y; %残差
stg_idx = 1; % stage index 阶段1 分阶段进行的
while (norm(res)>sigma)
% candidate list 候选集
[val idx] = sort(abs(Phi'*res), 'descend'); %计算内积并进行降序排序;并记录其对应的列号
candidate_set = [active_set; idx(1:actset_size)]; % 候选集=上一次的支撑集+此次内积的 前 step_size 个对应的列号 ,(即与y相关度最高的step_size个对应的H的列号)
% finalist
[val idx] = sort(abs(pinv(Phi(:,candidate_set))*y), 'descend'); %候选集中的列进行最小二乘(即求X的稀疏估计 X^),并进行降序排序;
new_active_set = candidate_set(idx(1:actset_size)); %选出X^的前actset_size(当前步长)个作为支撑集
new_res = y-Phi(:,new_active_set)*pinv(Phi(:,new_active_set))*y; % 选出的支撑集 最小二乘进行新的估计并更新残差
if (norm(new_res) >= norm(res))
% shift into a new stage
stg_idx = stg_idx + 1;
actset_size = stg_idx*step_siz