1、StOMP重構算法流程：

2、分段正交匹配追蹤(StOMP)Matlab代碼(CS_StOMP.m)

function [ theta ] = CS_StOMP( y,A,S,ts ) 
%CS_StOMP Summary of this function goes here 
%Version: 1.0 written by jbb0523 @2015-04-29 
%   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<y_columns 
        y = y‘;%y should be a column vector 
    end 
    [M,N] = size(A);%傳感矩陣A為M*N矩陣 
    theta = zeros(N,1);%用來存儲恢復的theta(列向量) 
    Pos_theta = [];%用來叠代過程中存儲A被選擇的列序號 
    r_n = y;%初始化殘差(residual)為y 
    for ss=1:S%最多叠代S次 
        product = A‘*r_n;%傳感矩陣A各列與殘差的內積 
        sigma = norm(r_n)/sqrt(M);%參見參考文獻第3頁Remarks(3) 
        Js = find(abs(product)>ts*sigma);%選出大於閾值的列 
        Is = union(Pos_theta,Js);%Pos_theta與Js並集 
        if length(Pos_theta) == length(Is) 
            if ss==1 
                theta_ls = 0;%防止第1次就跳出導致theta_ls無定義 
            end 
            break;%如果沒有新的列被選中則跳出循環 
        end 
        %At的行數要大於列數，此為最小二乘的基礎(列線性無關) 
        if length(Is)<=M 
            Pos_theta = Is;%更新列序號集合 
            At = A(:,Pos_theta);%將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列空間上的正交投影 
        r_n = y - At*theta_ls;%更新殘差 
        if norm(r_n)<1e-6%Repeat the steps until r=0 
            break;%跳出for循環 
        end 
    end 
    theta(Pos_theta)=theta_ls;%恢復出的theta 
end
 

3、StOMP單次重構測試代碼

以下測試代碼基本與OMP單次重構測試代碼一樣，除了調用CS_StOMP之外，一定要註意這裏的測量矩陣Phi =randn(M,N)/sqrt(M)，一定一定！！！

%壓縮感知重構算法測試 
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*theta 
Phi = randn(M,N)/sqrt(M);%測量矩陣為高斯矩陣 
A = Phi * Psi;%傳感矩陣 
y = Phi * x;%得到觀測向量y 
%% 恢復重構信號x 
tic 
theta = CS_StOMP(y,A); 
x_r = Psi * theta;% x=Psi * theta 
toc 
%% 繪圖 
figure; 
plot(x_r,‘k.-‘);%繪出x的恢復信號 
hold on; 
plot(x,‘r‘);%繪出原信號x 
hold off; 
legend(‘Recovery‘,‘Original‘) 
fprintf(‘\n恢復殘差：‘); 
norm(x_r-x)%恢復殘差 
 

4、門限參數ts、測量數M與重構成功概率關系曲線繪制例程代碼

clear all;close all;clc; 
%% 參數配置初始化 
CNT = 1000;%對於每組(K,M,N)，重復叠代次數 
N = 256;%信號x的長度 
Psi = eye(N);%x本身是稀疏的，定義稀疏矩陣為單位陣x=Psi*theta 
ts_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) 
tic 
for 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; 
    end 
end 
toc 
save 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)‘]); 
end 
for 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 
 

5、結語

[轉]壓縮感知重構算法之分段正交匹配追蹤(StOMP)