%% TV denoising demo
%
% Ivan Selesnick
% selesi@nyu.edu
% Revised 2014, 2017

%% Load data

clear

s = load('blocks.txt');          % blocks signal
y = load('blocks_noisy.txt');    % noisy blocks signal

N = length(s);        % N : length of signal
n = 1:N;

ax = [0 N -3 6];

figure(1)
clf
plot(s)
axis(ax)
title('Clean signal');
print -dpdf figures/demo_true

figure(2)
clf
plot(y)
axis(ax)
title('Noisy signal');
print -dpdf figures/demo_noisy


%% TV Denoising
% using L1 norm

lam = 2;
Nit = 100;
x_L1 = tvd_mm(y, lam, Nit);
% x_L1 = tvd_ncvx_ver1(y, lam, 'L1', Nit);   % Equivalent

figure(1)
clf
plot(x_L1)
axis(ax)
title('TV denoising');
print -dpdf figures/demo_TVD_L1norm


%% TV denoising using atan penalty

x_atan = tvd_ncvx_ver1(y, lam, 'atan', Nit);

figure(1)
clf
plot(x_atan)
axis(ax)
title('TV denoising using atan penalty');
print -dpdf figures/demo_TVD_atan

%% TV denoising using MC penalty

x_mcp = tvd_ncvx_ver1(y, lam, 'mcp', Nit);

figure(1)
clf
plot(x_mcp)
axis(ax)
title('TV denoising using MCP penalty');
print -dpdf figures/demo_TVD_mcp

%% Denoising error

figure(1)
clf
plot(1:N, x_L1 - s, 1:N, x_mcp - s)
legend('L1 norm', 'MC penalty')
xlim([0 N])
title('Denoising error')
print -dpdf figures/demo_compare


