压缩感知之分段正交匹配追踪法StOMP

% 分段匹配追踪（StOMP）算法。
%  读文件
% X=imread（‘C:\Users\zkdn\Desktop\MATLABdata\I30\I30.jpg‘）;
% X=double（X）;
% [ab]=size（X）;
load C:\Users\zkdn\Desktop\MATLABdata\I1\I1.mat
X=I1;
background=imopen（Xstrel（‘disk‘15））;%获取背景信息
X=imsubtract（Xbackground）;%利用函数去除背景

X=double（X）;
[ab]=size（X）;
%  小波变换矩阵生成
ww=dwt2（a‘haar‘）;
%  小波变换让图像稀疏化（注意该步骤会耗费时间，但是会增大稀疏度）
X1=ww*sparse（X）*ww‘;
% X1=X;
X1=full（X1）;

%  随机矩阵生成
M=200;
R=randn（Ma）;
% R=randn（a）;
% R=mapminmax（R0255）;
% R=round（R）;

%  测量值
Y=R*X1;

%  OMP算法
%  恢复矩阵
X2=zeros（ab）; 
%  按列循环
for i=1:b 
    %  通过OMP，返回每一列信号对应的恢复值（小波域）
    rec=CS_StOMP（Y（:i）Ra）;%[ theta ] = CS_StOMP（ yASts ）
    %  恢复值矩阵，用于反变换
    X2（:i）=rec;
end
%%%%
% figure（4）;
% X2=full（X2）;
% imshow（uint8（X2））;
% title（‘OMP计算后图像‘）;

%  原始图像
figure（1）;
imshow（uint8（X））;
title（‘原始图像‘）;

%  变换图像
% figure（2）;
% imshow（uint8（X1））;
% title（‘小波变换后的图像‘）;

%  压缩传感恢复的图像
figure（3）;
%  小波反变换
X3=ww‘*sparse（X2）*ww;  
% X3=X2;
X3=full（X3）;
imshow（uint8（X3））;
title（‘STOMP恢复的图像‘）;

%  误差（PSNR）
%  MSE误差
errorx=sum（sum（abs（X3-X）.^2））;        
%  PSNR
psnr=10*log10（255*255/（errorx/a/b））   



function [ theta ] = CS_StOMP（ yAS ）  
%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  
%   underdete