```% Section 6.3.3
% Boyd & Vandenberghe "Convex Optimization"
% Original by Lieven Vandenberghe
% Adapted for CVX Argyris Zymnis - 10/2005
%
% Suppose we have a signal x, which is mostly smooth, but has several
% rapid variations (or jumps). If we apply quadratic smoothing on
% this signal (see SMOOTHREC_CVX) then in order to remove the noise
% we will not be able to preserve the signal's sharp transitions.
%
% We can instead apply total variation reconstruction on the signal
% by solving
%        minimize ||x_hat - x_cor||_2 + lambda*TV(x_hat)
%
% where TV(x) = sum(abs(x_(i+1)-x_i)) , for i = 1 to n-1.
% The parameter lambda controls the ''smoothness'' of x_hat.
%
% Figure 1 shows the original and corrupted signals.
% Figure 2 shows the tradeoff curve obtained when varying lambda
% and figure 3 shows three reconstructed signals with different
% total variation.
%
% Figure 4 is a tradeoff curve for quadratic smoothing, while figure 5
% shows three reconstructed signals with quadratic smoothing.
% Note how TV reconstruction does a better job of preserving the
% sharp transitions in the signal while removing the noise.

n = 2000;  % length of signal
t = (0:n)';

figure(1)
subplot(211)
temp = ones(ceil((n+1)/4),1);
exact= [temp; -temp; temp; -temp];
exact = exact(1:n+1) + 0.5*sin((2*pi/n)*t);
plot(t,exact,'-');
axis([0 n+10 -2 2]);
ylabel('ya');
title('signal');
exact_variation = sum(abs(exact(2:(n+1)) - exact(1:n)))

subplot(212)
noise = 0.1*randn(size(t));
corrupt = exact+noise;
plot(t,corrupt,'-');
axis([0 n+10 -2 2]);
noisy_variation = sum(abs(corrupt(2:(n+1)) - corrupt(1:n)))
ylabel('yb');
xlabel('x');
title('corrupted signal');
%print -deps tv_exact_corrupt.eps % figure 6.11, page 315

% tradeoff curve, total variation vs ||x-xcorr||_2
% figure 6.13 page 316
fprintf('computing 100 points on tradeoff curve ... \n');
nopts = 100;
TVs = linspace(0.01,.9*noisy_variation,nopts);

obj1 = [];  obj2 = [];
for i=1:nopts
cvx_begin quiet
variable xrec(n+1)
minimize(norm(xrec-corrupt))
subject to
norm(xrec(2:(n+1))-xrec(1:n),1) <= TVs(i);
cvx_end
obj1 = [obj1, TVs(i)];
obj2 = [obj2, norm(full(xrec-corrupt))];
end;
obj1 = [0 obj1 noisy_variation];
obj2 = [norm(corrupt) obj2 0];

figure(2)
plot(obj2,obj1,'-'); hold on
plot(0,noisy_variation,'o');
plot(norm(corrupt),0,'o');  hold off
xlabel('x');
ylabel('y');
title('||Dxhat||_1 versus ||xhat-x||_2');
%print -deps tv_tradeoff.eps % figure 6.13, page 316

figure(3)
subplot(311)
% solve total variation problem
cvx_begin quiet
variable xrec(n+1)
minimize(norm(xrec-corrupt))
subject to
norm(xrec(2:(n+1))-xrec(1:n),1) <= 10;
cvx_end
plot(t,xrec','-');
axis([0 n -2 2]);
ylabel('ya');
title('xhat with TV=10');

subplot(312)
cvx_begin quiet
variable xrec(n+1)
minimize(norm(xrec-corrupt))
subject to
norm(xrec(2:(n+1))-xrec(1:n),1) <= 8;
cvx_end
plot(t,xrec','-');
axis([0 n -2 2]);
ylabel('yb');
title('xhat with TV=8');

subplot(313)
cvx_begin quiet
variable xrec(n+1)
minimize(norm(xrec-corrupt))
subject to
norm(xrec(2:(n+1))-xrec(1:n),1) <= 5;
cvx_end
plot(t,xrec','-');
axis([0 n -2 2]);
xlabel('x');
ylabel('yc');
title('xhat with TV=5');

%print -deps tv_rec_10_8_5.eps % figure 6.14, page 317

% quadratic smoothing, figure 6.12, page 316
% In this case it is not a good idea to use CVX
% as the sparsity in the closed form solution
% makes it very easy to solve directly
A = sparse(n,n+1);
A(:,1:n) = -speye(n,n);  A(:,2:n+1) = A(:,2:n+1)+speye(n,n);

nopts = 100;
lambdas = logspace(-10,10,nopts);
obj1 = [];  obj2 = [];
for i=1:nopts

lambda = lambdas(i);
x = (A'*A+lambda*speye(n+1,n+1)) \ (lambda*corrupt);
obj1 = [obj1, norm(full(A*x))];
obj2 = [obj2, norm(full(x-corrupt))];
end;

figure(4)
plot(obj2,obj1,'-'); hold on
plot(0,norm(A*corrupt),'o');
plot(norm(corrupt),0,'o'); hold off
xlabel('x');
ylabel('y');
title('||Dxhat||_2 vs ||xhat-xcor||_2');

nopts = 3;
alphas = [10 7 4];
xrecon = [];

for i=1:3
alpha = alphas(i);
u = 10;  l = -10;  normx = Inf;
while (abs(normx-alpha) > 1e-3)
lambda = 10^((u+l)/2);
x = (A'*A+lambda*speye(n+1,n+1)) \ (lambda*corrupt);
normx = norm(x-corrupt);
if (normx > alpha), l = (u+l)/2; else u = (u+l)/2;  end;
end;
xrecon = [xrecon, x];

end;

figure(5)
subplot(311), plot(xrecon(:,1));
axis([0 n -2 2])
ylabel('ya');
subplot(312), plot(xrecon(:,2));
axis([0 n -2 2])
ylabel('yb');
subplot(313), plot(xrecon(:,3));
axis([0 n -2 2])
xlabel('x');
ylabel('yc');
% figure 6.12, page 316
```
```exact_variation =

7.9968

noisy_variation =

231.1383

computing 100 points on tradeoff curve ...