主要内容：

1. StOMP的算法流程
2. StOMP的MATLAB实现
3. 一维信号的实验与结果
4. 门限参数Ts、测量数M与重构成功概率关系的实验与结果

一、StOMP的算法流程

分段正交匹配追踪(Stagewise OMP)也是由OMP改进而来的一种贪心算法，与CoSaMP、SP算法类似，不同之处在于CoSaMP、SP算法在迭代过程中选择的是与信号内积最大的2K或K个原子，而StOMP是通过门限阈值来确定原子。此算法的输入参数中没有信号稀疏度K，因此相比于ROMP及CoSaMP有独到的优势（这句话存在疑问）。

StOMP的算法流程：

二、StOMP的MATLAB实现(CS_StOMP.m)

function [ theta ] = CS_StOMP( y,A,S,ts )%   CS_StOMP%   Detailed explanation goes here%   y = Phi * x%   x = Psi * theta%    y = Phi*Psi * theta%   令 A = Phi*Psi, 则y=A*theta%   S is the maximum number of StOMP iterations to perform%   ts is the threshold parameter%   现在已知y和A，求theta%   Reference:Donoho D L，Tsaig Y，Drori I，Starck J L．Sparse solution of%   underdetermined linear equations by stagewise orthogonal matching %   pursuit[J]．IEEE Transactions on Information Theory，2012，58(2)：1094—1121    if nargin < 4        ts = 2.5; %ts范围[2,3],默认值为2.5    end    if nargin < 3        S = 10; %S默认值为10    end    [y_rows,y_columns] = size(y);    if y_rows
     
      ts*sigma); %选出大于阈值的列        Is = union(pos_num,Js); %pos_num与Js并集        if length(pos_num) == length(Is)            if ss==1                theta_ls = 0; %防止第1次就跳出导致theta_ls无定义            end            break; %如果没有新的列被选中则跳出循环        end        %At的行数要大于列数，此为最小二乘的基础(列线性无关)        if length(Is)<=M            pos_num = Is; %更新列序号集合            At = A(:,pos_num); %将A的这几列组成矩阵At        else %At的列数大于行数，列必为线性相关的,At'*At将不可逆            if ss==1                theta_ls = 0; %防止第1次就跳出导致theta_ls无定义            end            break; %跳出for循环        end        %y=At*theta，以下求theta的最小二乘解(Least Square)        theta_ls = (At'*At)^(-1)*At'*y; %最小二乘解        %At*theta_ls是y在At列空间上的正交投影        res = y - At*theta_ls; %更新残差        if norm(res)<1e-6 %Repeat the steps until r=0            break; %跳出for循环        end    end    theta(pos_num)=theta_ls; %恢复出的thetaend
     

三、一维信号的实验与结果

%压缩感知重构算法测试clear all;close all;clc;M = 64; %观测值个数N = 256; %信号x的长度K = 12; %信号x的稀疏度Index_K = randperm(N);x = zeros(N,1);x(Index_K(1:K)) = 5*randn(K,1); %x为K稀疏的，且位置是随机的Psi = eye(N); %x本身是稀疏的，定义稀疏矩阵为单位阵x=Psi*thetaPhi = randn(M,N)/sqrt(M); %测量矩阵为高斯矩阵A = Phi * Psi; %传感矩阵y = Phi * x; %得到观测向量y%% 恢复重构信号xtictheta = CS_StOMP(y,A);x_r = Psi * theta; % x=Psi * thetatoc%% 绘图figure;plot(x_r,'k.-'); %绘出x的恢复信号hold on;plot(x,'r'); %绘出原信号xhold off;legend('Recovery','Original')fprintf('\n恢复残差：');norm(x_r-x) %恢复残差

四、门限参数ts、测量数M与重构成功概率关系的实验与结果

clear all;close all;clc;%% 参数配置初始化CNT = 1000;%对于每组(K,M,N)，重复迭代次数N = 256;%信号x的长度Psi = eye(N);%x本身是稀疏的，定义稀疏矩阵为单位阵x=Psi*thetats_set = 2:0.2:3;K_set = [4,12,20,28,36];%信号x的稀疏度集合Percentage = zeros(N,length(K_set),length(ts_set));%存储恢复成功概率%% 主循环，遍历每组(ts,K,M,N)ticfor tt = 1:length(ts_set)    ts = ts_set(tt);    for kk = 1:length(K_set)        K = K_set(kk);%本次稀疏度        %M没必要全部遍历，每隔5测试一个就可以了        M_set=2*K:5:N;        PercentageK = zeros(1,length(M_set));%存储此稀疏度K下不同M的恢复成功概率        for mm = 1:length(M_set)           M = M_set(mm);%本次观测值个数           fprintf('ts=%f,K=%d,M=%d\n',ts,K,M);           P = 0;           for cnt = 1:CNT %每个观测值个数均运行CNT次                Index_K = randperm(N);                x = zeros(N,1);                x(Index_K(1:K)) = 5*randn(K,1);%x为K稀疏的，且位置是随机的                                Phi = randn(M,N)/sqrt(M);%测量矩阵为高斯矩阵                A = Phi * Psi;%传感矩阵                y = Phi * x;%得到观测向量y                theta = CS_StOMP(y,A,10,ts);%恢复重构信号theta                x_r = Psi * theta;% x=Psi * theta                if norm(x_r-x)<1e-6%如果残差小于1e-6则认为恢复成功                    P = P + 1;                end           end           PercentageK(mm) = P/CNT*100;%计算恢复概率        end        Percentage(1:length(M_set),kk,tt) = PercentageK;    endendtocsave StOMPMtoPercentage1000 %运行一次不容易，把变量全部存储下来%% 绘图for tt = 1:length(ts_set)    S = ['-ks';'-ko';'-kd';'-kv';'-k*'];    figure;    for kk = 1:length(K_set)        K = K_set(kk);        M_set=2*K:5:N;        L_Mset = length(M_set);        plot(M_set,Percentage(1:L_Mset,kk,tt),S(kk,:));%绘出x的恢复信号        hold on;    end    hold off;    xlim([0 256]);    legend('K=4','K=12','K=20','K=28','K=36');    xlabel('Number of measurements(M)');    ylabel('Percentage recovered');    title(['Percentage of input signals recovered correctly(N=256,ts=',...        num2str(ts_set(tt)),')(Gaussian)']);endfor kk = 1:length(K_set)    K = K_set(kk);    M_set=2*K:5:N;    L_Mset = length(M_set);    S = ['-ks';'-ko';'-kd';'-kv';'-k*';'-k+'];    figure;    for tt = 1:length(ts_set)        plot(M_set,Percentage(1:L_Mset,kk,tt),S(tt,:));%绘出x的恢复信号        hold on;    end    hold off;    xlim([0 256]);    legend('ts=2.0','ts=2.2','ts=2.4','ts=2.6','ts=2.8','ts=3.0');    xlabel('Number of measurements(M)');    ylabel('Percentage recovered');    title(['Percentage of input signals recovered correctly(N=256,K=',...        num2str(K),')(Gaussian)']);    end

1、门限参数ts分别为2.0,2.2,2.4,2.6,2.8,3.0时，不同稀疏信号下，测量值M与重构成功概率的关系：

2、稀疏度为4,12,20,28,36时，不同门限参数ts下，测量值M与重构成功概率的关系：

结论：

通过对比可以看出，总体上讲ts=2.4或ts=2.6时效果较好，较大和较小重构效果都会降低，这里由于没有ts=2.5的情况，但我们推测ts=2.5应该是一个比较好的值，因此一般默认取为2.5即可。

六、参考文章