% One-dimensional de-blurring example

% Number of discretization points for 1D interval including end points
N = 128;
gamma = 0.02;
C = 1 / (sqrt(2*pi)*gamma);

K = zeros(N,N);
h = 1/(N-1);
x = linspace(0,1,N)';

for l = 1:N
    for k = 1:N
    	K(l,k) = h * C * exp(-(l-k)^2 * h^2 / (2 * gamma^2));
    end
end

p = sin(4*pi*x) + 0.5 * cos(24*pi*x) + 0.0 * cos(50*pi*x);
x0 = 0.4;
k = C * exp(-(x-x0).^2/(2*gamma^2));


d = K * p;

% s = svd(K);
% plot(s,'o')

% no noise
p_alpha = K\d;
figure;
plot(x,p,x,p_alpha,'Linewidth', 2), axis([0,1,-1.5,1.5]);
legend('exact data','reconstruction no noise, no regularizaton');

% noisy data
dn = d + 1e-4 * randn(N,1);

p_alpha = K\dn;
figure;
plot(x,p,x,p_alpha,'Linewidth', 2), axis([0,1,-1.5,1.5]);
legend('exact data','reconstruction tiny noise, no regularizaton');


% regularization parameter
alpha = 1e-8;
% solve Tikhonov system (this could also be done using a QR decomposition)
p_alpha = (K'*K + alpha * eye(N))\(K'*dn);
figure;
plot(x,p,x,p_alpha,'Linewidth', 2), axis([0,1,-1.5,1.5]);
legend('exact data', 'Tikhonov reconstruction');

return;
