% 1D deblurring example. A given input signal xtrue is convolved with a
% Gaussian kerne and white noise is added. The problem is to estimate the
% input from the noisy blurred sample.

% Discretization of the interval [0,1] by N points

N = 160;
t = linspace(0,1,N);

% Defining the Gaussian blurring kernel

width = 0.05;
a = 1/sqrt(2*pi*width^2)*exp(-0.5*(1/width^2)*t.^2);
A = 1/N*toeplitz(a);

% Plotting the kernel at two different locations

figure(1)
plot(t,A(10,:),'k-','LineWidth',2)
text(0.2,0.7*max(A(10,:)),'n=10')
hold on
plot(t,A(100,:),'r-','LineWidth',2)
text(0.8,0.7*max(A(100,:)),'n=100')
hold off
pause

% Defining different input signals to claculate the data

% Boxcar input

%xtrue = zeros(N,1);
%n = 40;
%xtrue(50:50+n) = ones(n+1,1);
%figure(2)
%plot(t,xtrue,'k-','LineWidth',2)
%axis([0,1,-0.1,1.1*max(xtrue)])

% Smooth input

xtrue = (t'-0.5).*exp(-0.5*1e3*(t'-0.5).^2);
figure(2)
plot(t,xtrue,'k-','LineWidth',2)
axis([0,1,-0.02,0.02])

pause

% Claculating the blurred noisy singnal. The standard deviation of the
% additive noise is defined percentagewise of the maximum of the noiseless
% signal.

noise = 4;       % Noise level in percentages of the max. of noiseless signal
y0 = A*xtrue;    % Noiseless signal
eta = max(y0)*noise/100;                % STD of the additive noise
y = y0 + max(y0)*noise/100*randn(N,1);

figure(3)
plot(t,y,'k-','LineWidth',2)

pause

% MAP solution with different Gaussian priors

% White noise prior

%STD = 0.1;    % Reasonable values: Boxcar: STD = 0.5, Smooth: 0.001 ... 0.01
%L = 1/STD*eye(N);

% Second order smoothness prior

STD = 10;     % Reasonable vaues: Boxcar: STD = 20 ... 1000, Smooth: 1 ... 10
L = N^2*(2*eye(N) - diag(ones(N-1,1),1) - diag(ones(N-1,1),-1));
L = 1/STD*L;

xMAP = [(1/eta)*A;L]\[(1/eta)*y;zeros(N,1)];

figure(4)
plot(t,xMAP,'k-','LineWidth',2)
hold on
plot(t,xtrue,'r-','LineWidth',2)
hold off

