% Example computations related to one-dimensional deconvolution.
% Demonstration of truncated SVD solution.
%
% Samuli Siltanen September 2008

% Plot parameters
lwidth = .2;
redlwidth = .7;
fsize = 8;
msize = 7;
sigma = .02; % Noise level

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create data on finer grid to avoid inverse crime

% Construct signal
Nf  = 1499;
sf  = linspace(0,1,Nf);
Dsf = sf(2)-sf(1);
xf  = zeros(size(sf));
xf((sf>(1/4))&(sf<(1/2))) = 1;
xf((sf>(3/4))&(sf<(7/8))) = 2;

% Construct normalized point spread function
a    = .1;
nPSFf = ceil(a/Dsf); 
sPSFf = [-nPSFf:nPSFf]*Dsf;
PSFf = zeros(size(sPSFf));
PSFf((sPSFf>(-a))&(sPSFf<a)) = 1/(2*a);
PSFf = PSFf/(Dsf*sum(PSFf));

% Construct convolved signal
mf = Dsf*conv2(xf,PSFf,'same');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Construct objects needed for reconstruction
N  = 64;
s  = linspace(0,1,N);
Ds = s(2)-s(1);

% Construct normalized point spread function
nPSF = ceil(a/Ds); 
sPSF = [-nPSF:nPSF]*Ds;
PSF = zeros(size(sPSF));
PSF((sPSF>(-a))&(sPSF<a)) = 1/(2*a);
PSF = PSF/(Ds*sum(PSF));

% Construct measurement matrix A
A = convmtx(PSF,N);
A = Ds*A(:,nPSF+1:end-nPSF);

% Construct measurements with and without noise
m = interp1(sf,mf,s,'linear');
m = m(:);
mn = m + sigma*randn(size(m));

% Construct continuum quantities for visual comparison
% Construct signal
Nc = 5000;
sc = linspace(-2,2,2*Nc+1);
Xc = zeros(size(sc));
Xc((sc>(1/4))&(sc<(1/2))) = 1;
Xc((sc>(3/4))&(sc<(7/8))) = 2;
PSFc = zeros(size(sc));
PSFc((sc>(-a))&(sc<a)) = 1/(2*a);
Mc = (sc(2)-sc(1))*conv2(Xc,PSFc,'full');
Mc = Mc(Nc+[1:length(sc)]);

% Plot measurement
figure(1)
clf
plot(sc,Mc,'r','linewidth',redlwidth)
hold on
plot(s,mn,'k','linewidth',lwidth)
plot(s,mn,'k.','markersize',msize)
axis([0 1 -.2 2.2])
set(gca,'xtick',[0 1/4 1/2 3/4 7/8 1],'fontsize',fsize)
set(gca,'xticklabel',{'0' '1/4' '1/2' '3/4' '7/8' '1'},'fontsize',fsize)
set(gca,'ytick',[0  1 2],'fontsize',fsize)
axis square
box off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute SVD
[U,D,V] = svd(A);

% Take a look at matrix A and its SVD
figure(2)
clf
subplot(1,2,1)
spy(A)
[row,col] = size(A);
axis([1 col 1 row])
set(gca,'xtick',[1 round(col/2) col])
set(gca,'ytick',[1 round(row/2) row])
axis square
title('Nonzero elements of matrix A')
subplot(1,2,2)
semilogy(diag(D))
xlim([1 N])
set(gca,'xtick',[1 round(N/2) N])
axis square
title('Singular values of matrix A')

% Take a look at singular vectors (columns of matrix V)
figure(4)
clf
subplot(2,3,1)
plot(1:N,V(:,1),'k','linewidth',lwidth)
hold on
plot(1:N,V(:,1),'k.','markersize',msize)
title('First singular vector')
axis([1 N -.2 .4])
set(gca,'ytick',[0 1 2])
set(gca,'xtick',[1 round(N/2) N])
set(gca,'xticklabel',{})
set(gca,'yticklabel',{})
axis square
subplot(2,3,2)
plot(1:N,V(:,2),'k','linewidth',lwidth)
hold on
plot(1:N,V(:,2),'k.','markersize',msize)
title('Second singular vector')
axis([1 N -.2 .4])
set(gca,'ytick',[0 1 2])
set(gca,'xtick',[1 round(N/2) N])
set(gca,'xticklabel',{})
set(gca,'yticklabel',{})
axis square
subplot(2,3,3)
plot(1:N,V(:,10),'k','linewidth',lwidth)
hold on
plot(1:N,V(:,10),'k.','markersize',msize)
title('Tenth singular vector')
axis([1 N -.2 .4])
set(gca,'ytick',[0 1 2])
set(gca,'xtick',[1 round(N/2) N])
set(gca,'xticklabel',{})
set(gca,'yticklabel',{})
axis square
subplot(2,3,4)
plot(1:N,V(:,end-10),'k','linewidth',lwidth)
hold on
plot(1:N,V(:,end-10),'k.','markersize',msize)
title('Tenth-to-last singular vector')
axis([1 N -.2 .4])
set(gca,'ytick',[0 1 2])
set(gca,'xtick',[1 round(N/2) N])
set(gca,'xticklabel',{})
set(gca,'yticklabel',{})
axis square
subplot(2,3,5)
plot(1:N,V(:,end-1),'k','linewidth',lwidth)
hold on
plot(1:N,V(:,end-1),'k.','markersize',msize)
title('Second-to-last singular vector')
axis([1 N -.2 .4])
set(gca,'ytick',[0 1 2])
set(gca,'xtick',[1 round(N/2) N])
set(gca,'xticklabel',{})
set(gca,'yticklabel',{})
axis square
subplot(2,3,6)
plot(1:N,V(:,end),'k','linewidth',lwidth)
hold on
plot(1:N,V(:,end),'k.','markersize',msize)
title('Last singular vector')
axis([1 N -.2 .4])
set(gca,'ytick',[0 1 2])
set(gca,'xtick',[1 round(N/2) N])
set(gca,'xticklabel',{})
set(gca,'yticklabel',{})
axis square

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Demonstration of inversion by truncated SVD

% Naive reconstruction from noisy data
recnaive = inv(A)*mn;

% Reconstruct from noisy data using only large singular values.
% The singular values smaller than tolerance are discarded.
Ns               = 30;
[row,col]        = size(D.');
Dplus            = sparse(row,col);
svals            = diag(D);
Dplus(1:Ns,1:Ns) = diag(1./svals(1:Ns));
rec              = V*Dplus*U.'*mn(:);

% Take a look at the reconstruction
figure(3)
clf
subplot(1,2,1)
plot(sc,Xc,'r','linewidth',redlwidth)
hold on
plot(s,recnaive,'k','linewidth',lwidth)
plot(s,recnaive,'k.','markersize',msize)
title('Naive reconstruction')
axis([0 1 -.2 2.2])
set(gca,'xtick',[0 1/4 1/2 3/4 7/8 1],'fontsize',fsize)
set(gca,'xticklabel',{'0' '1/4' '1/2' '3/4' '7/8' '1'},'fontsize',fsize)
set(gca,'ytick',[0  1 2],'fontsize',fsize)
axis square
box off
subplot(1,2,2)
plot(sc,Xc,'r','linewidth',redlwidth)
hold on
plot(s,rec,'k','linewidth',lwidth)
plot(s,rec,'k.','markersize',msize)
title(['SVD truncated at ', num2str(Ns)])
axis([0 1 -.2 2.2])
set(gca,'xtick',[0 1/4 1/2 3/4 7/8 1],'fontsize',fsize)
set(gca,'xticklabel',{'0' '1/4' '1/2' '3/4' '7/8' '1'},'fontsize',fsize)
set(gca,'ytick',[0  1 2],'fontsize',fsize)
axis square
box off


