close all; clear all; clc
% 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 = 60000;  % length of signal
t = (0:n)';

figure(1)
subplot(311)
temp = ones(ceil((n+1)/4),1);
exact= [temp; -temp; temp; -temp];
exact = exact(1:n+1) + 0.5*sin((2*pi/n)*t); %original was 2
plot(t,exact,'-');
axis([0 n+10 -2 2]);
ylabel('voltage','Interpreter','latex');
xlabel('sample index $$n$$','Interpreter','latex');
%title('signal'); 
title('signal, $x[n]$','Interpreter','latex')

exact_variation = sum(abs(exact(2:(n+1)) - exact(1:n)))

subplot(312)
noise = 0.2*randn(size(t));
plot(t,noise,'-');
axis([0 n+10 -2 2]);
ylabel('voltage','Interpreter','latex');
xlabel('sample index $$n$$','Interpreter','latex');
%title('noise');
title('noise, $w[n]$','Interpreter','latex')

subplot(313)
corrupt = exact+noise;
plot(t,corrupt,'-');
axis([0 n+10 -2 2]);
noisy_variation = sum(abs(corrupt(2:(n+1)) - corrupt(1:n)))
ylabel('voltage','Interpreter','latex');
xlabel('sample index $$n$$','Interpreter','latex');
%title('corrupted signal');
title('corrupted signal, $x_{\texttt{cor}}[n]=x[n]+w[n]$','Interpreter','latex')

figure_name=['L1_example_1_1.eps'];
   print('-depsc',figure_name)
%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
%      fprintf('tradeoff point %d\n',i);
%      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; %original was 10
   cvx_end
   plot(t,xrec','-');
   axis([0 n -2 2]);
   ylabel('voltage','Interpreter','latex');
   xlabel('sample index $$n$$','Interpreter','latex');

   title('estimate, ${\hat x}_{{\texttt{a}}}[n]$','Interpreter','latex')
   

   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('voltage','Interpreter','latex');
   xlabel('sample index $$n$$','Interpreter','latex');
   title('estimate, ${\hat x}_{\texttt{b}}[n]$','Interpreter','latex')

   subplot(313)
   cvx_begin quiet
    variable xrec(n+1)
    minimize(norm(xrec-corrupt))
    subject to
        norm(xrec(2:(n+1))-xrec(1:n),1) <= 7;
   cvx_end
   plot(t,xrec','-');
   axis([0 n -2 2]);
   xlabel('x');
  ylabel('voltage','Interpreter','latex');
   xlabel('sample index $$n$$','Interpreter','latex');
   title('estimate, ${\hat x}_{\texttt{c}}[n]$','Interpreter','latex')
   
   figure_name=['L1_example_1_2.eps'];
   print('-depsc',figure_name)

   done = 1
   %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);

% tradeoff curve with quadratic smoothing
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');
%print -deps tv_smooth_tradeoff.eps

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');
title('quadratic smoothing with ||xhat-xcor||_2=10');
subplot(312), plot(xrecon(:,2));
axis([0 n -2 2])
ylabel('yb');
title('quadratic smoothing with ||xhat-xcor||_2=7');
subplot(313), plot(xrecon(:,3));
axis([0 n -2 2])
xlabel('x');
ylabel('yc');
title('quadratic smoothing with ||xhat-xcor||_2=4');
%print -deps tv_smooth_tradeoff_examples.eps
% figure 6.12, page 316