% Section 6.5.4
% Boyd & Vandenberghe "Convex Optimization"
% Original by Lieven Vandenberghe
% Adapted for CVX by Argyris Zymnis - 11/27/2005
%
% Here we find a sparse basis for a signal y out of
% a set of Gabor functions. We do this by solving
%       minimize  ||A*x-y||_2 + ||x||_1
%
% where the columns of A are sampled Gabor functions.
% We then fix the sparsity pattern obtained and solve
%       minimize  ||A*x-y||_2
%
% NOTE: The file takes a while to run

clear

% Problem parameters
sigma = 0.05;  % Size of Gaussian function
Tinv  = 500;   % Inverse of sample time
Thr   = 0.001; % Basis signal threshold
kmax  = 30;    % Number of signals are 2*kmax+1
w0    = 5;     % Base frequency (w0 * kmax should be 150 for good results)

% Build sine/cosine basis
fprintf(1,'Building dictionary matrix...');
% Gaussian kernels
TK = (Tinv+1)*(2*kmax+1);
t  = (0:Tinv)'/Tinv;
A  = exp(-t.^2/(sigma^2));
ns = nnz(A>=Thr)-1;
A  = A([ns+1:-1:1,2:ns+1],:);
ii = (0:2*ns)';
jj = ones(2*ns+1,1)*(1:Tinv+1);
oT = ones(1,Tinv+1);
A  = sparse(ii(:,oT)+jj,jj,A(:,oT));
A  = A(ns+1:ns+Tinv+1,:);
% Sine/Cosine basis
k  = [ 0, reshape( [ 1 : kmax ; 1 : kmax ], 1, 2 * kmax ) ];
p  = zeros(1,2*kmax+1); p(3:2:end) = -pi/2;
SC = cos(w0*t*k+ones(Tinv+1,1)*p);
% Multiply
ii = 1:numel(SC);
jj = rem(ii-1,Tinv+1)+1;
A  = sparse(ii,jj,SC(:)) * A;
A  = reshape(A,Tinv+1,(Tinv+1)*(2*kmax+1));
fprintf(1,'done.\n');

% Construct example signal
a = 0.5*sin(t*11)+1;
theta = sin(5*t)*30;
b = a.*sin(theta);

% Solve the Basis Pursuit problem
disp('Solving Basis Pursuit problem...');
tic
cvx_begin
    variable x(30561)
    minimize(sum_square(A*x-b)+norm(x,1))
cvx_end
disp('done');
toc

% Reoptimize problem over nonzero coefficients
p = find(abs(x) > 1e-5);
A2 = A(:,p);
x2 = A2 \ b;

% Constants
M = 61; % Number of different Basis signals
sk = 250; % Index of s = 0.5

% Plot example basis functions;
%if (0) % to do this, re-run basispursuit.m to create A
figure(1); clf;
subplot(3,1,1); plot(t,A(:,M*sk+1)); axis([0 1 -1 1]);
title('Basis function 1');
subplot(3,1,2); plot(t,A(:,M*sk+31)); axis([0 1 -1 1]);
title('Basis function 2');
subplot(3,1,3); plot(t,A(:,M*sk+61)); axis([0 1 -1 1]);
title('Basis function 3');
%print -deps bp-dict_helv.eps

% Plot reconstructed signal
figure(2); clf;
subplot(2,1,1);
plot(t,A2*x2,'--',t,b,'-'); axis([0 1 -1.5 1.5]);
xlabel('t'); ylabel('y_{hat} and y');
title('Original and Reconstructed signals')
subplot(2,1,2);
plot(t,A2*x2-b); axis([0 1 -0.06 0.06]);
title('Reconstruction error')
xlabel('t'); ylabel('y - y_{hat}');
%print -deps bp-approx_helv.eps

% Plot frequency plot
figure(3); clf;

subplot(2,1,1);
plot(t,b); xlabel('t'); ylabel('y'); axis([0 1 -1.5 1.5]);
title('Original Signal')
subplot(2,1,2);
plot(t,150*abs(cos(w0*t)),'--');
hold on;
for k = 1:length(t);
  if(abs(x((k-1)*M+1)) > 1e-5), plot(t(k),0,'o'); end;
  for j = 2:2:kmax*2
    if((abs(x((k-1)*M+j)) > 1e-5) | (abs(x((k-1)*M+j+1)) > 1e-5)),
      plot(t(k),w0*j/2,'o');
    end;
  end;
end;
xlabel('t'); ylabel('w');
title('Instantaneous frequency')
hold off;
Building dictionary matrix...done.
Solving Basis Pursuit problem...
 
Calling sedumi: 61625 variables, 502 equality constraints
------------------------------------------------------------
SeDuMi 1.21 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
eqs m = 502, order n = 61125, dim = 61626, blocks = 2
nnz(A) = 7484105 + 0, nnz(ADA) = 252004, nnz(L) = 126253
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            7.63E+04 0.000
  1 :  -1.00E-01 3.44E+02 0.000 0.0045 0.9990 0.9990   3.00  1  1  5.1E-01
  2 :   5.06E+00 1.08E+02 0.000 0.3136 0.9000 0.9000   1.00  1  1  4.8E-01
  3 :   1.03E+01 4.82E+01 0.000 0.4468 0.9000 0.9000   1.00  1  1  4.1E-01
  4 :   1.19E+01 2.73E+01 0.000 0.5657 0.9000 0.9000   1.00  1  1  3.5E-01
  5 :   1.25E+01 1.52E+01 0.000 0.5567 0.9000 0.9000   1.00  1  1  2.7E-01
  6 :   1.27E+01 9.41E+00 0.000 0.6201 0.9000 0.9034   1.00  1  1  2.1E-01
  7 :   1.27E+01 5.70E+00 0.000 0.6057 0.9000 0.9068   1.00  1  1  1.5E-01
  8 :   1.28E+01 2.83E+00 0.000 0.4968 0.9272 0.9000   1.00  1  1  9.1E-02
  9 :   1.28E+01 1.33E+00 0.000 0.4707 0.9496 0.9000   1.00  1  1  4.7E-02
 10 :   1.28E+01 9.07E-01 0.000 0.6808 0.9000 0.9081   1.00  1  1  3.3E-02
 11 :   1.28E+01 5.50E-01 0.000 0.6068 0.9325 0.9000   1.00  1  1  2.1E-02
 12 :   1.28E+01 3.38E-01 0.000 0.6146 0.9170 0.9000   1.00  1  1  1.3E-02
 13 :   1.28E+01 1.99E-01 0.000 0.5880 0.9353 0.9000   1.00  1  1  7.6E-03
 14 :   1.28E+01 1.19E-01 0.000 0.6001 0.9052 0.9000   1.00  1  1  4.6E-03
 15 :   1.28E+01 7.74E-02 0.000 0.6485 0.9037 0.9000   1.00  1  1  3.0E-03
 16 :   1.28E+01 3.41E-02 0.000 0.4402 0.9000 0.9042   1.00  1  1  1.3E-03
 17 :   1.28E+01 1.26E-02 0.000 0.3690 0.9184 0.9000   1.00  2  2  4.9E-04
 18 :   1.28E+01 1.83E-03 0.000 0.1454 0.9190 0.9000   1.00  2  2  7.1E-05
 19 :   1.28E+01 4.35E-04 0.000 0.2377 0.9000 0.9070   1.00  2  2  1.7E-05
 20 :   1.28E+01 8.31E-05 0.000 0.1912 0.9011 0.9000   1.00  2  2  3.2E-06
 21 :   1.28E+01 2.17E-05 0.000 0.2607 0.9000 0.9018   1.00  2  2  8.5E-07
 22 :   1.28E+01 7.19E-06 0.000 0.3320 0.9000 0.9031   1.00  2  2  2.8E-07
 23 :   1.28E+01 2.79E-06 0.000 0.3884 0.9006 0.9000   1.00  3  3  1.1E-07
 24 :   1.28E+01 1.07E-06 0.000 0.3825 0.9079 0.9000   1.00  3  3  4.2E-08
 25 :   1.28E+01 3.84E-07 0.000 0.3589 0.9112 0.9000   1.00  3  3  1.5E-08
 26 :   1.28E+01 1.38E-07 0.000 0.3596 0.9076 0.9000   1.00  4  4  5.4E-09

iter seconds digits       c*x               b*y
 26     94.0   8.0  1.2845155824e+01  1.2845155686e+01
|Ax-b| =   1.6e-11, [Ay-c]_+ =   1.1E-11, |x|=  2.6e+00, |y|=  9.9e-01

Detailed timing (sec)
   Pre          IPM          Post
2.820E+00    9.395E+01    3.000E-01    
Max-norms: ||b||=1.497811e+00, ||c|| = 1,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 3306.09.
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +12.8452
done
Elapsed time is 86.389939 seconds.