% denoise a noisy step function v by minimizing the functional
% J(u) = \int (u-v)^2 + mu \int |u'|^2

% start fresh
clear all;
clc;

% figure
h1 = figure(1);
close(h1);
h1 = figure(1);
set(h1,'Position',[10 10 500 500]);

% grid
N  = 128;
x  = (0:N)/N;
h  = 1/N;

% regularization
mu = 1.0e-3;
% noise level
ns = 0.05;

% exact step function
ux = ones(N+1,1);
ux(1:round(N/2)) = 0;

% data
randn('state',0);
v  = ux + ns*randn(N+1,1);

% system to solve
E = ones(N,1);
D = diag([E;0]+[0;E]) - diag(E,1) - diag(E,-1);
A = eye(N+1) + mu/h^2*D;

% solve system
u = A \ v;

% show result
plot(x,ux,'b',x,v,'r',x,u,'g');
legend('exact','noisy','denoised');
title('denoising');
xlabel('x');
ylabel('u');
