% Figures 6.8-6.10: Quadratic smoothing % 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