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