Figures 6.21-6.23: Basis pursuit using Gabor functions
clear
sigma = 0.05;
Tinv = 500;
Thr = 0.001;
kmax = 30;
w0 = 5;
fprintf(1,'Building dictionary matrix...');
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,:);
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);
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');
a = 0.5*sin(t*11)+1;
theta = sin(5*t)*30;
b = a.*sin(theta);
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
p = find(abs(x) > 1e-5);
A2 = A(:,p);
x2 = A2 \ b;
M = 61;
sk = 250;
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');
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}');
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 = 30563
nnz(A) = 3742304 + 0, nnz(ADA) = 252004, nnz(L) = 126253
it : b*y gap delta rate t/tP* t/tD* feas cg cg prec
0 : 3.82E+04 0.000
1 : 4.12E-02 4.00E+02 0.000 0.0105 0.9900 0.9900 4.66 1 1 3.7E-01
2 : 3.74E+00 1.60E+02 0.000 0.3989 0.9000 0.9000 1.01 1 1 3.5E-01
3 : 8.68E+00 7.43E+01 0.000 0.4650 0.9000 0.9000 1.01 1 1 3.2E-01
4 : 1.16E+01 3.59E+01 0.000 0.4826 0.9000 0.9000 1.00 1 1 2.7E-01
5 : 1.23E+01 2.01E+01 0.000 0.5612 0.9000 0.9000 1.00 1 1 2.2E-01
6 : 1.26E+01 1.02E+01 0.000 0.5071 0.9139 0.9000 1.00 1 1 1.6E-01
7 : 1.27E+01 4.82E+00 0.000 0.4727 0.9311 0.9000 1.00 1 1 9.7E-02
8 : 1.28E+01 3.12E+00 0.000 0.6477 0.9051 0.9000 1.00 1 1 6.9E-02
9 : 1.28E+01 2.05E+00 0.000 0.6571 0.9000 0.9260 1.00 1 1 4.9E-02
10 : 1.28E+01 1.15E+00 0.000 0.5588 0.9095 0.9000 1.00 1 1 2.9E-02
11 : 1.28E+01 6.97E-01 0.000 0.6074 0.9000 0.9098 1.00 1 1 1.8E-02
12 : 1.28E+01 4.60E-01 0.000 0.6604 0.9000 0.9216 1.00 1 1 1.2E-02
13 : 1.28E+01 2.34E-01 0.000 0.5082 0.9039 0.9000 1.00 1 1 6.3E-03
14 : 1.28E+01 1.44E-01 0.000 0.6167 0.9000 0.9219 1.00 1 1 3.9E-03
15 : 1.28E+01 6.69E-02 0.000 0.4641 0.9000 0.9000 1.00 1 1 1.8E-03
16 : 1.28E+01 2.54E-02 0.000 0.3791 0.9049 0.9000 1.00 1 1 7.0E-04
17 : 1.28E+01 4.12E-03 0.000 0.1622 0.9095 0.9000 1.00 2 2 1.1E-04
18 : 1.28E+01 3.67E-04 0.000 0.0892 0.9900 0.9900 1.00 2 2 1.0E-05
19 : 1.28E+01 8.20E-05 0.000 0.2234 0.9000 0.9000 1.00 2 2 2.3E-06
20 : 1.28E+01 2.69E-05 0.000 0.3280 0.9000 0.9042 1.00 2 2 7.4E-07
21 : 1.28E+01 9.48E-06 0.000 0.3522 0.9000 0.9012 1.00 2 2 2.6E-07
22 : 1.28E+01 3.55E-06 0.000 0.3746 0.9000 0.9000 1.00 3 3 9.8E-08
23 : 1.28E+01 1.15E-06 0.249 0.3253 0.9000 0.9000 1.00 3 3 3.2E-08
24 : 1.28E+01 1.16E-07 0.494 0.1007 0.9290 0.9000 1.00 4 4 3.2E-09
iter seconds digits c*x b*y
24 246.1 8.0 1.2845155802e+01 1.2845155686e+01
|Ax-b| = 1.9e-11, [Ay-c]_+ = 3.8E-13, |x|= 3.5e+00, |y|= 9.9e-01
Detailed timing (sec)
Pre IPM Post
1.081E+01 2.461E+02 7.000E-02
Max-norms: ||b||=1.497811e+00, ||c|| = 1,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 3304.27.
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +12.8452
done
Elapsed time is 201.717853 seconds.