% Example computations related to two-dimensional deconvolution.
% Here we study reconstruction using truncated SVD
%
% Samuli Siltanen September 2008
% Construct example signal (two-dimensional)
N = 32;
x = zeros(N,N);
x(round(N/4):round(N/2),round(N/4):round(N/2)) = 1;
% Take a look at the signal
figure(1)
clf
imagesc(x)
colormap gray
title('Signal x')
axis square
% Construction of the measurement matrix A
A = oblur(N,5);
% Compute SVD
[U,D,V] = svd(full(A));
% Take a look at matrix A and its SVD
figure(2)
clf
subplot(1,2,1)
spy(A)
axis square
title('Nonzero elements of matrix A')
subplot(1,2,2)
plot(diag(D),'r','linewidth',3)
axis([1 length(diag(D)) 0 max(diag(D))])
set(gca,'xtick',[1 round(length(diag(D))/2) length(diag(D))])
axis square
title(['Singular values of matrix A'])
% Construct ideal and noisy measurements m and mn
m = A*x(:);
m = reshape(m,N,N);
noiselevel = 0.02;
mn = m + noiselevel*randn(size(m));
% Take a look at the noisy measurement
figure(1)
clf
subplot(1,2,1)
imagesc(x)
colormap gray
title('Signal x')
axis square
axis off
subplot(1,2,2)
imagesc(mn)
colormap gray
title('Measurement')
axis square
axis off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Demonstration reconstruction using truncated SVD
% Reconstruct using truncated SVD (Ns singular values included)
Ns = 80;
[row,col] = size(D.');
Dplus = sparse(row,col);
svals = diag(D);
Dplus(1:Ns,1:Ns) = diag(1./svals(1:Ns));
recn = V*Dplus*U.'*mn(:);
recn = reshape(recn,N,N);
% Take a look at the reconstruction
figure(3)
clf
subplot(1,2,1)
imagesc(x)
colormap gray
title('Signal x')
axis square
axis off
subplot(1,2,2)
imagesc(recn)
colormap gray
title(['SVD truncation with Ns = ', num2str(Ns)])
axis square
axis off
% Take a look at singular vectors (columns of matrix V)
figure(4)
subplot(2,2,1)
imagesc(reshape(V(:,1),N,N))
title('First singular vector')
colormap gray
axis square
axis off
subplot(2,2,2)
imagesc(reshape(V(:,10),N,N))
title('Tenth singular vector')
colormap gray
axis square
axis off
subplot(2,2,3)
imagesc(reshape(V(:,end-10),N,N))
title('Tenth-to-last singular vector')
colormap gray
axis square
axis off
subplot(2,2,4)
imagesc(reshape(V(:,end),N,N))
title('Last singular vector')
colormap gray
axis square
axis off
figure(3)