www.gusucode.com > 基于Matlab编写压缩感知重建算法集 > 基于Matlab编写压缩感知重建算法集/CS_recovery_algorithms_OMP_SP_IHT/software/Demo_CS_IRLS.m
function Demo_CS_IRLS() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % the DCT basis is selected as the sparse representation dictionary % instead of seting the whole image as a vector, I process the image in the % fashion of column-by-column, so as to reduce the complexity. % Author: Chengfu Huo, roy@mail.ustc.edu.cn, http://home.ustc.edu.cn/~roy % Reference: R. Chartrand and W. Yin, “Iteratively Reweighted Algorithms % for Compressed Sensing,” 2008. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %------------ read in the image -------------- img=imread('lena.bmp'); % testing image img=double(img); [height,width]=size(img); %------------ form the measurement matrix and base matrix --------------- Phi=randn(floor(height/3),width); % only keep one third of the original data Phi = Phi./repmat(sqrt(sum(Phi.^2,1)),[floor(height/3),1]); % normalize each column mat_dct_1d=zeros(256,256); % building the DCT basis (corresponding to each column) for k=0:1:255 dct_1d=cos([0:1:255]'*k*pi/256); if k>0 dct_1d=dct_1d-mean(dct_1d); end; mat_dct_1d(:,k+1)=dct_1d/norm(dct_1d); end %--------- projection --------- img_cs_1d=Phi*img; % treat each column as a independent signal %-------- recover using omp ------------ sparse_rec_1d=zeros(height,width); Theta_1d=Phi*mat_dct_1d; for i=1:width column_rec=cs_irls(img_cs_1d(:,i),Theta_1d,height); sparse_rec_1d(:,i)=column_rec'; % sparse representation end img_rec_1d=mat_dct_1d*sparse_rec_1d; % inverse transform %------------ show the results -------------------- figure(1) subplot(2,2,1),imagesc(img),title('original image') subplot(2,2,2),imagesc(Phi),title('measurement mat') subplot(2,2,3),imagesc(mat_dct_1d),title('1d dct mat') psnr = 20*log10(255/sqrt(mean((img(:)-img_rec_1d(:)).^2))) subplot(2,2,4),imagesc(img_rec_1d),title(strcat('1d rec img ',num2str(psnr),'dB')) disp('over') %************************************************************************% function hat_x=cs_irls(y,T_Mat,m) % y=T_Mat*x, T_Mat is n-by-m % y - measurements % T_Mat - combination of random matrix and sparse representation basis % m - size of the original signal % the sparsity is length(y)/4 %hat_x_tp=(T_Mat'*T_Mat)^(-1)*T_Mat'*y; % initialization, leading to warning hat_x_tp=T_Mat'*y; epsilong=1; p=1; % solution for l-norm p times=1; while (epsilong>10e-9) && (times<length(y)/4) weight=(hat_x_tp.^2+epsilong).^(p/2-1); Q_Mat=diag(1./weight,0); hat_x=Q_Mat*T_Mat'*inv(T_Mat*Q_Mat*T_Mat')*y; if(norm(hat_x-hat_x_tp,2) < sqrt(epsilong)/100) epsilong=epsilong/10; end hat_x_tp=hat_x; times=times+1; end