function?mcyOMP2
clc;clear
X=imread（‘G:\Program?Files\MATLAB\R2012b\graduate\rm.jpg‘）;?%?讀文件
if?ndims（X）>2
????X=rgb2gray（X）;
end
X=imresize（X[256?256]）;
X=double（X）;
[ab]=size（X）;
ww=DWT（a）;?%?小波變換矩陣生成
X1=ww*sparse（X）*ww‘;?%?小波變換讓圖像稀疏化（注意該步驟會耗費時間，但是會增大稀疏度
X1=full（X1）;
M=128;?%?隨機矩陣生成
R=randn（Ma）;
Y=R*X1;?%?測量
%?OMP算法
X2=zeros（ab）;?%?恢復矩陣
for?i=1:b?%?列循環???????
????rec=omp（Y（:i）Ra）;
????X2（:i）=rec;
end
figure（1）;?%?原始圖像
imshow（uint8（X））;
title（‘原始圖像‘）;
figure（2）;?%?變換圖像
imshow（uint8（X1））;
title（‘小波變換后的圖像‘）;
figure（3）;?%?壓縮傳感恢復的圖像
X3=ww‘*sparse（X2）*ww;?%?小波反變換
X3=full（X3）;
imshow（uint8（X3））;
title（‘恢復的圖像‘）;
%?誤差（PSNR）
errorx=sum（sum（abs（X3-X）.^2））;????????%?MSE誤差
PSNR=10*log10（255*255/（errorx/a/b））???%?PSNR
MSE=errorx/a/b
match=1-sum（abs（X3-X）.^2）/sum（abs（X3+X）.^2）
%?OMP的函數
%?s-測量；T-觀測矩陣；N-向量大小
function?hat_y=omp（sTN）
Size=size（T）;?????????????????????????????????????%?觀測矩陣大小
M=Size（1）;???