% In this routine we compute discrete convolution data using a convolution
% matrix. 
%
% Samuli Siltanen January 2011

% Choose level of discretization 
n  = 32;

% Choose relative amplitude of noise (0.01 equals 1 percent noise)
noiselevel = 0.01;

% Parameter for the point spread function. Should satisfy 0<a<1/2
a = 0.04;

% Construct sample points
x  = [0:(n-1)]/n;
Dx = x(2)-x(1);

% Construct normalized discrete point spread function
nPSF = ceil(a/Dx); 
xPSF = [-nPSF:nPSF]*Dx;
PSF  = zeros(size(xPSF));
ind   = abs(xPSF)<a;
PSF(ind) = DC_PSF(xPSF(ind),a);
Ca   = 1/(Dx*trapz(PSF));
PSF  = Ca*PSF;

% Construct f
f = DC_target(x);

% Construct convolution matrix
A = Dx*DC_convmtx(PSF,n);

% Reset random number generator to make it possible to repeat experiments
% with exactly same noise. If you want to experiment with different noise
% components, comment this line.
reset(RandStream.getDefaultStream);

% Compute ideal and noisy data 
m  = A*f(:);
mn = m + noiselevel*max(abs(f))*randn(n,1);

% Save results to disc
save data/DC_discr_data A x m mn PSF noiselevel

% Plot the results
DC2_discretedata_plot

