www.pudn.com > Wavelet_OMP.zip > Wavelet_OMP.m, change:2009-10-14,size:2499b


%  本程序实现图像LENA的压缩传感 
%  程序作者:沙威,香港大学电气电子工程学系,wsha@eee.hku.hk 
%  算法采用正交匹配法,参考文献 Joel A. Tropp and Anna C. Gilbert  
%  Signal Recovery From Random Measurements Via Orthogonal Matching 
%  Pursuit,IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 53, NO. 12, 
%  DECEMBER 2007. 
%  该程序没有经过任何优化 
 
function Wavelet_OMP 
 
clc;clear 
 
%  读文件 
X=imread('lena256.bmp'); 
X=double(X); 
[a,b]=size(X); 
 
%  小波变换矩阵生成 
ww=DWT(a); 
 
%  小波变换让图像稀疏化(注意该步骤会耗费时间,但是会增大稀疏度) 
X1=ww*sparse(X)*ww'; 
X1=full(X1); 
 
%  随机矩阵生成 
M=190; 
R=randn(M,a); 
 
%  测量 
Y=R*X1; 
 
%  OMP算法 
X2=zeros(a,b);  %  恢复矩阵 
for i=1:b  %  列循环        
    rec=omp(Y(:,i),R,a); 
    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 
 
 
%  OMP的函数 
%  s-测量;T-观测矩阵;N-向量大小 
function hat_y=omp(s,T,N) 
 
Size=size(T);                                     %  观测矩阵大小 
M=Size(1);                                        %  测量 
hat_y=zeros(1,N);                                 %  待重构的谱域(变换域)向量                      
Aug_t=[];                                         %  增量矩阵(初始值为空矩阵) 
r_n=s;                                            %  残差值 
 
for times=1:M/4;                                  %  迭代次数(稀疏度是测量的1/4) 
    for col=1:N;                                  %  恢复矩阵的所有列向量 
        product(col)=abs(T(:,col)'*r_n);          %  恢复矩阵的列向量和残差的投影系数(内积值)  
    end 
    [val,pos]=max(product);                       %  最大投影系数对应的位置 
    Aug_t=[Aug_t,T(:,pos)];                       %  矩阵扩充 
    T(:,pos)=zeros(M,1);                          %  选中的列置零(实质上应该去掉,为了简单我把它置零) 
    aug_y=(Aug_t'*Aug_t)^(-1)*Aug_t'*s;           %  最小二乘,使残差最小 
    r_n=s-Aug_t*aug_y;                            %  残差 
    pos_array(times)=pos;                         %  纪录最大投影系数的位置 
     
    if (norm(r_n)<9)                              %  残差足够小 
        break; 
    end 
end 
hat_y(pos_array)=aug_y;                           %  重构的向量