% Simple illustration of deconvolution in 1D, the method we use is
% Truncated Singular Value Decomposition
%
% Samuli Siltanen February 2017

%% Load real-world data
load data/deco03 mn mn2 m0

%% Construct computational model
% Build a Point Spread Function (PSF)
M = 9;
% Parameter determining the width of the PSF
a   = 10;
% Evaluation points
x = linspace(-a,a,2*M+1);
psf = PSF(x,a);
psf = psf/sum(psf(:));

% Construct "unknown signal" f
N = 256;
x = linspace(0,1,N);
fstep = zeros(N,1);
fstep(1:(end/2)) = 1;
fstep = 1-fstep;
fstep = fstep(:); % Force vector f to be vertical

% Construct the convolution matrix
A = convmtx(psf,N);
A = A(:,(M+1):(end-M));

% Simulate the "measurement"
m_simul = A*fstep;

%% Compute truncated SVD solution

% Determine the SVD of matrix A
[U,D,V] = svd(A);
svals   = diag(D);

% Compute reconstruction
r_alpha = 50;
Dp_alpha = zeros(size(A));
for iii = 1:r_alpha
    Dp_alpha(iii,iii) = 1/svals(iii);
end
f0 = V*(Dp_alpha.')*(U.')*m0(:);
fn = V*(Dp_alpha.')*(U.')*mn(:);


%% Save results to file

save data/TSVD fstep m0 m_simul A svals U V mn fn r_alpha f0

