分块压缩感知程序matlab

1、%本程序实现图像LENA的压缩传感%程序作者：沙威，香港大学电气电子工程学系，wshaeee.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_OMPclcclear%读文件X=i

2、mread(lena256.bmp);X=double(X);a,b=size(X);size_kuai=16*4;X2=zeros(size_kuai); %恢复矩阵X3=zeros(a,b); %恢复矩阵%小波变换矩阵生成ww=DWT(size_kuai);%随机矩阵生成M=12*4;R=randn(M,size_kuai);ticfor i_x=1:ceil(a/size_kuai)for i_y=1:ceil(b/size_kuai)XX=X(i_x-1)*size_kuai+1:i_x*size_kuai,(i_y-1)*size_kuai+1:i_y*size_kuai);%小波变

3、换让图像稀疏化(注意该步骤会耗费时间，但是会增大稀疏度)X1=ww*sparse(XX)*ww;X1=full(X1);%测量Y=R*X1;% OMP算法for i=1:size_kuai %列循环rec=omp_fenkuai(Y(:,i),R,size_kuai);X2(:,i)=rec;endX3(i_x-1)*size_kuai+1:i_x*size_kuai,(i_y-1)*size_kuai+1:i_y*size_kuai)=ww*sparse(X2)*ww;%小波反变换end endX3=full(X3);use_time=toc%原始图像figure(1);imshow(uin

4、t8(X);title(原始图像);%压缩传感恢复的图像figure(2);imshow(uint8(X3);title(分块恢复的图像);%误差(PSNR) errorx=sum(sum(abs(X3-X).A2); psnr=10*log10(255*255/(errorx/a/b)endhat_y(pos_array)=aug_y;% MSE误差% PSNR% OMP的函数% s-测量；T-观测矩阵；N-向量大小function hat_y=omp_fenkuai(s,T,N)Size=size(T);M=Size(1);hat_y=zeros(1,N);Aug_t=;r_n=s;for

5、 times=1:M/4for col=1:Nproduct(col)=abs(T(:,col)*r_n);数(内积值)endval,pos=max(product);Aug_t=Aug_t,T(:,pos);T(:,pos)=zeros(M,1);为了简单我把它置零)aug_y=(Aug_t*Aug_t)A(-1)*Aug_t*s;r_n=s-Aug_t*aug_y;pos_array(times)=pos;if (norm(r_n)40)break;end%观测矩阵大小%测量% 待重构的谱域（变换域）向量%增量矩阵（初始值为空矩阵）%残差值%迭代次数（稀疏度是测量的1/4）%恢复矩阵的所有

6、列向量%恢复矩阵的列向量和残差的投影系%最大投影系数对应的位置%矩阵扩充%选中的列置零（实质上应该去掉,%最小二乘，使残差最小%残差%纪录最大投影系数的位置%残差足够小%重构的向量%程序作者：沙威，香港大学电气电子工程学系，wshaeee.hku.hk%参考文献：小波分析理论与MA TLAB R2007实现，葛哲学，沙威，第20章 小波变换在矩阵方程求解中的应用（沙威、陈明生编写）%构造正交小波变换矩阵，图像大小N*N , N=2AP , P是整数。function ww=DWT(N)%矩阵构造for jj=rank_min:rank_maxnn=2Ajj；%构造向量p1_0=sparse(h

7、,zeros(1,nn-L);p2_0=sparse(g,zeros(1,nn-L);%向量圆周移位for ii=1:nn/2p1(ii,:)=circshift(p1_0,2*(ii-1);p2(ii,:)=circshift(p2_0,2*(ii-1);end%构造正交矩阵w1=p1;p2;mm=2Arank_max-length(w1);w=sparse(w1,zeros(length(w1),mm);zeros(mm,length(w1),eye(mm,mm);ww=ww*w;clear p1;clear p2;endh,g= wfilters(sym8,d);%分解低通和高通滤波器%