% Example computations related to one-dimensional deconvolution.
% Demonstration of generalized Tikhonov regularization.
%
% The following routines must be computed before this one:
% DC3_nocrimedata_comp.m
%
% Samuli Siltanen February 2011

% Regularization parameter
delta = 10;

% Plot parameters
lwidth = .5;
thickline = 2;
fsize = 20;
msize = 4;

% Load data computed using a very fine grid, avoiding inverse crime
load data/DC_nocrimedata xx Af A y n m mn1 mn2 sigma1 sigma2

% Choose data by commenting two unwanted lines among the three lines below. 
% There is no inverse crime in any of these choices.
%data = m;   % No added noise
data = mn1; % Small added noise
%data = mn2; % Large added noise

% Construct regularization matrix of size (n)x(n). This implements difference
% between consecutive values assuming periodic boundary conditions.
L = eye(n);
L = L-[L(:,end),L(:,1:end-1)];

% Construct original signal for comparison
xx = linspace(0,1,1024);
signal = DC_target(xx);

% Reconstruct from noisy data using Tikhonov regularization.
syst  = [A; sqrt(delta)*L];
m_aug = [data(:);zeros(n,1)];
rec   = syst\m_aug;

% Calculate relative error
truth  = DC_target(y);
relerr = round(norm(rec(:)-truth(:))/norm(truth(:))*100);

% Take a look at the reconstruction
figure(6)
clf
plot(xx,signal,'k','linewidth',lwidth)
hold on
plot(y,rec,'k','linewidth',thickline)
text(0.9,1.5,[num2str(relerr),'%'],'fontsize',fsize)

% Axis settings
axis([0 1 -.5 2])
set(gca,'xtick',[0 1/2 1],'fontsize',fsize)
set(gca,'ytick',[0  1 2],'fontsize',fsize)
set(gca,'PlotBoxAspectRatio' , [2 1 1])
box off
