%----------------------------------
% Demo: Total variation (TV) reconstruction
% by Jiaxian Pan and Wing-Kin Ma, Jan 2013
%
% In figure 1, we show the original signal and the corrupted signal
% 
% In figure 2, we show the results of the following reconstruction
% formulation
%
% min \|x- x_corrupted\|_2^2 + lambda* \|D * x\|_2^2
%
% which has a closed-form solution (LS).
% We try three values of lambda : 0.1, 1, 10
%
% In figure 3, we show the results of the following reconstruction
% formulation, which is popularly known as TV reconstruction
%
% min \|x- x_corrupted\|_2^2 + lambda* \|D * x\|_1
%
% We try three values of lambda : 0.1, 1, 10


close all;
clear;
n_temp=400;

%temp=sawtooth(linspace(0,2*pi,100))';
temp=linspace(0,1,n_temp)';
x_true=[0.2+3*temp;-temp;-5*temp;3*temp]; %the groundtruth

n=length(x_true);

x_corrupted=x_true+0.2*randn(n,1); %corrupt the signal by noise

figure(1);
subplot(211);
plot(x_true);
axis([0 n -1.2*max(abs(x_true)) 1.2*max(abs(x_true))]);
title('Original Signal');
subplot(212);
plot(x_corrupted);
axis([0 n -1.2*max(abs(x_true)) 1.2*max(abs(x_true))]);
title('Corrupted signal');


%--------------------
% Generate the difference matrix
D=diag(-ones(n-2,1),0)+diag(2*ones(n-3,1),1)+diag(-ones(n-4,1),2);
D=[D zeros(n-2,2)];
D(end,end-1)=2; D(end,end)=-1;
D(end-1,end-1)=-1;

%---------------------
% Solve the following problem
% min \|x- x_corrupted\|_2^2 + lambda* \|D * x\|_2^2
%------------------

%lambda=0.1
figure(2)
subplot(311)

x=(eye(n)+0.1*D'*D)\x_corrupted; %the solution is LS

plot(x','-');
axis([0 n -1.2*max(abs(x_true)) 1.2*max(abs(x_true))]);
title('xhat with lambda=0.1');

%lambda=1
subplot(312)

x=(eye(n)+1*D'*D)\x_corrupted; %the solution is LS

plot(x','-');
axis([0 n -1.2*max(abs(x_true)) 1.2*max(abs(x_true))]);
title('xhat with lambda=1');

%lambda=10
subplot(313)

x=(eye(n)+10*D'*D)\x_corrupted; %the solution is LS
plot(x','-');
axis([0 n -1.2*max(abs(x_true)) 1.2*max(abs(x_true))]);
title('xhat with lambda=10');


%---------solve the following problem for different value of TV
%      min \|xrec - x_corrupted\|_2^2 + lambda* \|D x\|_1


figure(3)

%lambda=0.1
subplot(311)
cvx_begin
    variable x(n)
    minimize(square_pos(norm(x-x_corrupted)) + 0.1*norm(D*x,1))
cvx_end
plot(x','-');
axis([0 n -1.2*max(abs(x_true)) 1.2*max(abs(x_true))]);
title('xhat with lambda=0.1');

%lambda= 1
subplot(312)
cvx_begin
    variable x(n)
    minimize(square_pos(norm(x-x_corrupted)) + norm(D*x,1))
cvx_end
plot(x','-');
axis([0 n -1.2*max(abs(x_true)) 1.2*max(abs(x_true))]);
title('xhat with lambda=1');

%lambda= 10
subplot(313)
cvx_begin
    variable x(n)
    minimize(square_pos(norm(x-x_corrupted)) + 10*norm(D*x,1))
cvx_end
plot(x','-');
axis([0 n -1.2*max(abs(x_true)) 1.2*max(abs(x_true))]);
title('xhat with lambda=10');


   