% 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 does not vary too rapidly
% and that x is corrupted by some small, rapidly varying noise v,
% ie. x_cor = x + v. Then if we want to reconstruct x from x_cor
% we should solve (with x_hat as the parameter)
%        minimize ||x_hat - x_cor||_2 + lambda*phi_quad(x_hat)
%
% where phi_quad(x) = sum(x_(i+1)-x_i)^2 , for i = 1 to n-1.
% The parameter lambda controls the ''smoothness'' of x_hat.
%
% The first figure which is generated shows the original and
% the corrupted signals. The second figure shows the tradeoff curve
% obtained when varying lambda and the third figure shows three
% reconstructed signals.
%
% NOTE: This is not a good problem to use CVX on. By exploiting
% the sparsity in this case, we can solve this problem much more
% efficiently using least squares.


randn('state',0);

n = 4000;  t = (0:n-1)';
exact = 0.5*sin((2*pi/n)*t).*sin(0.01*t);
corrupt = exact + 0.05*randn(size(exact));

figure(1)
subplot(211)
plot(t,exact,'-');
axis([0 n -0.6 0.6])
title('original signal');
ylabel('ya');

subplot(212)
plot(t,corrupt,'-');
axis([0 n -0.6 0.6])
xlabel('x');
ylabel('yb');
title('corrupted signal');
%print -deps smoothrec_signals.eps % figure 6.8, page 313

A = sparse(n-1,n);
A(:,1:n-1) = -speye(n-1,n-1);  A(:,2:n) = A(:,2:n)+speye(n-1,n-1);

% tradeoff curve, figure 6.9, page 313
nopts = 100;
lambdas = logspace(-10,10,nopts);

obj1 = [];  obj2 = [];

fprintf('computing 100 points on tradeoff curve ... \n');

for i=1:nopts

  lambda = lambdas(i);
  cvx_begin quiet
    variable x(n)
    minimize(norm(x-corrupt)+lambda*norm(x(2:n)-x(1:n-1)))
  cvx_end
  obj1 = [obj1, norm(full(A*x))];
  obj2 = [obj2, norm(full(x-corrupt))];

  fprintf('tradeoff point %d\n',i);
end;

figure(2)
plot(obj2,obj1,'-');  hold on;
plot(0,norm(A*corrupt),'o');
plot(norm(corrupt),0,'o');  hold off;
xlabel('x');
ylabel('y');
title('||xhat-xcorr||_2 vs. ||D xhat||_2');
%print -deps smoothrec_tradeoff.eps % figure 6.9, page 313

%three smooth signals, figure 6.10, page 314
nopts = 3;
alphas = [8 3 1];
xrecon = [];

for i=1:3
   fprintf(1,'Reconstructed Signals: %d of 3 \n',i)
   alpha = alphas(i);
   cvx_begin quiet
    variable x(n)
    minimize(norm(x(2:n)-x(1:n-1)))
    subject to
        norm(x-corrupt) <= alpha;
   cvx_end
   xrecon = [xrecon, x];

end

figure(3)
subplot(311), plot(xrecon(:,1));
axis([0 n -0.6 0.6])
ylabel('ya');
title('||xhat-xcorr||_2=8');
subplot(312), plot(xrecon(:,2));
axis([0 n -0.6 0.6])
ylabel('yb');
title('||xhat-xcorr||_2=3');
subplot(313), plot(xrecon(:,3));
axis([0 n -0.6 0.6])
xlabel('x');
ylabel('yc');
title('||xhat-xcorr||_2=1');
%print -deps smoothrec_results.eps % figure 6.10, page 314

computing 100 points on tradeoff curve ... 
tradeoff point 1
tradeoff point 2
tradeoff point 3
tradeoff point 4
tradeoff point 5
tradeoff point 6
tradeoff point 7
tradeoff point 8
tradeoff point 9
tradeoff point 10
tradeoff point 11
tradeoff point 12
tradeoff point 13
tradeoff point 14
tradeoff point 15
tradeoff point 16
tradeoff point 17
tradeoff point 18
tradeoff point 19
tradeoff point 20
tradeoff point 21
tradeoff point 22
tradeoff point 23
tradeoff point 24
tradeoff point 25
tradeoff point 26
tradeoff point 27
tradeoff point 28
tradeoff point 29
tradeoff point 30
tradeoff point 31
tradeoff point 32
tradeoff point 33
tradeoff point 34
tradeoff point 35
tradeoff point 36
tradeoff point 37
tradeoff point 38
tradeoff point 39
tradeoff point 40
tradeoff point 41
tradeoff point 42
tradeoff point 43
tradeoff point 44
tradeoff point 45
tradeoff point 46
tradeoff point 47
tradeoff point 48
tradeoff point 49
tradeoff point 50
tradeoff point 51
tradeoff point 52
tradeoff point 53
tradeoff point 54
tradeoff point 55
tradeoff point 56
tradeoff point 57
tradeoff point 58
tradeoff point 59
tradeoff point 60
tradeoff point 61
tradeoff point 62
tradeoff point 63
tradeoff point 64
tradeoff point 65
tradeoff point 66
tradeoff point 67
tradeoff point 68
tradeoff point 69
tradeoff point 70
tradeoff point 71
tradeoff point 72
tradeoff point 73
tradeoff point 74
tradeoff point 75
tradeoff point 76
tradeoff point 77
tradeoff point 78
tradeoff point 79
tradeoff point 80
tradeoff point 81
tradeoff point 82
tradeoff point 83
tradeoff point 84
tradeoff point 85
tradeoff point 86
tradeoff point 87
tradeoff point 88
tradeoff point 89
tradeoff point 90
tradeoff point 91
tradeoff point 92
tradeoff point 93
tradeoff point 94
tradeoff point 95
tradeoff point 96
tradeoff point 97
tradeoff point 98
tradeoff point 99
tradeoff point 100
Reconstructed Signals: 1 of 3 
Reconstructed Signals: 2 of 3 
Reconstructed Signals: 3 of 3