index.html

信号处理系列导航

set_rand_seeds(123456,21);

%% 
% First we load a seismic filter, which is a second derivative of a
% Gaussian.

% dimension of the signal
n = 1024;
% width of the filter
s = 13;
% second derivative of Gaussian
t = -n/2:n/2-1;
h = (1-t.^2/s^2).*exp( -(t.^2)/(2*s^2) );
h = h-mean(h);
% normalize it
h = h/norm(h);
% recenter the filter for periodic boundary conditions
h1 = fftshift(h);

%% 
% We compute the filtering matrix. To stabilize the recovery, we sub-sample
% by a factor of 2 the filtering.

% sub-sampling (distance between wavelets)
sub = 2;
% number of atoms in the dictionary
p = n/sub;
% the dictionary, with periodic boundary conditions
[Y,X] = meshgrid(1:sub:n,1:n);
D = reshape( h1(mod(X-Y,n)+1), [n p]);

%%
% We generate a sparse signal to recover.

% spacing min and max between the spikes.
m = 5; M = 40;
k = floor( (p+M)*2/(M+m) )-1;
spc = linspace(M,m,k)';
% location of the spikes
sel = round( cumsum(spc) );
sel(sel>p) = [];
% randomization of the signs and values
x = zeros(p,1);
si = (-1).^(1:length(sel))'; si = si(randperm(length(si)));
% creating of the sparse spikes signal.
x(sel) = si;
x = x .* (1-rand(p,1)*.5);
% sparsity of the solution
M = sum(x~=0);

%%
% Now we perform the measurements.

% noise level
sigma = .06*max(h);
% noise
w = randn(n,1)*sigma;
% measures
y = D*x + w;

%%
% Display.

clf;
subplot(2,1,1);
plot_sparse_diracs(x);
title('Signal x');
subplot(2,1,2);
plot(y);
set_graphic_sizes([], 20);
axis('tight');
title('Noisy measurements y=D*x+w');

%% Homotopy Method for L1 Regularization
% Instead of using a greedy matching pursuit to recover an approximation of
% |x|, we now turn to more global minimization method. L1 minimization is a
% convex regularization that approximates the L0 pseudo-norm regularization
% (that would be ideal to find highly sparse Diracs train).
% The minimization reads

%%
% |x(lambda) = argmin_x 1/2*norm(y-D*x)^2 + lambda*norm(x,1)|

%%
% The Homotopy method starts with a large |lambda| so that |x(lambda)=0|
% and progressively decays the value of |lambda| so that the sparsity of
% |x(lambda)| only increases/decreases by 1 at each step.

%%
% The intial solution is 0.

xbp = zeros(p,1);

%%
% The solution |xbp| is update |xbp = xbp + gamma*d| where the direction of
% update is supported on the set of vectors that maximally correlate with
% the residual |y-D*xbp|

% compute the correlation with the residual
C = D'*(y-D*xbp);
lambda = max(abs(C));
% find the locations that maximally correlates
S = find( abs(abs(C/lambda)-1)<1e-9);
% compute the complementary set
I = ones(p,1); I(S)=0;
Sc = find(I);

%%
% The direction of descent |d(S)| on |S| is computed so that its image by |D(:,S)*d(S)| correlates
% as +1 or -1 with the atoms in |S| (same speed on all the coefficients).
% It is zero outside |S|.

d = zeros(p,1);
d(S) = (D(:,S)'*D(:,S)) \ sign( C(S) );

%%
% Now we compute the value |gamma| so that either

%%
% * 1) |xbp + gamma*d| correlates as much as |lambda| with one atom outside |S|.
% * 2) |xbp + gamma*d| correlates as much as |-lambda| with one atom outside |S|.
% * 3) one of the coordinates |xbp + gamma*d| in |S| becomes 0.

%%
% In situations 1) and 2), the number of non-zero coordinates of |xbp|
% (its sparsity) increases by 1. In situation 3), its sparsity stays
% constant because one coefficients appears and another deasapears. 

%% 
% Compute minimum |gamma| so that situation 1) is in force.

v = D(:,S)*d(S);
w = ( lambda-C(Sc) ) ./ ( 1 - D(:,Sc)'*v );
gamma1 = min(w(w>0));


%% 
% Compute minimum |gamma| so that situation 2) is in force.

w = ( lambda+C(Sc) ) ./ ( 1 + D(:,Sc)'*v );
gamma2 = min(w(w>0));

%%
% Compute minimum |gamma| so that situation 3) is in force.

w = -xbp(S)./d(S);
gamma3 = min(w(w>0));

%%
% Compute minimum |gamma| so that 1), 2) or 3) is in force, and update the
% solution.

gamma = min([gamma1 gamma2 gamma3]);

%%
% We can display the residual correlation before/after the update. You see
% that a new maximally correlating atoms has appeared, and that the overall correlation has decayed a little.

clf;
subplot(2,1,1);
plot( abs(D'*(y-D*xbp)), '.-' ); axis tight;
title('Correlation before update');
subplot(2,1,2);
plot( abs(D'*(y-D*(xbp+gamma*d))), '.-' ); axis tight;
title('Correlation after update');

%%
% Update the solution.
xbp = xbp + gamma*d;

%%
% _Exercice 1:_ (the solution is <../private/sparsity_seismic_bp/exo1.m exo1.m>)
% Implements the homotopy algorithm by applying several times (let's say up
% to |1.5*M| times) the previous update rule. Keep track of the evolution of |lambda| and of the
% evolution of the sparsity. 

exo1;

%%
% Display the solution computed by L1.
clf;
subplot(2,1,1);
plot_sparse_diracs(x);
title('Signal x');
subplot(2,1,2);
plot_sparse_diracs(xbp);
title('Recovered by L1 minimization');


%%
% The solution computed by L1 as a tendency to under-estimate the true
% values of the coefficients. This is bias inherent to L1 minimization. To
% remove this bias, one can perform a back projection that performs an L2
% best fit on the support selected by L1.

% find the support
sel = find(xbp~=0);
% perform the fit
xproj = zeros(p,1);
xproj(sel) = D(:,sel) \ y;
% display
clf;
subplot(2,1,1);
plot_sparse_diracs(xbp);
title('Recovered by L1 minimization');
subplot(2,1,2);
plot_sparse_diracs(xproj);
title('Back-projection');

%%
% _Exercice 2:_ (the solution is <../private/sparsity_seismic_bp/exo2.m exo2.m>)
% In the same code, add a tracking of the error |norm(x-xproj)|, and
% return the solution that mimizes this error. Save the value of the
% optimal |lambda| in |lambda_opt| and the value |xbp| of the solution for this |lambda| 
% for a future use.

exo2;

%% 
% Display the best solution

err_bp = norm(x-xproj)/norm(x);
clf;
subplot(2,1,1);
plot_sparse_diracs(x);
title('Signal x');
subplot(2,1,2);
plot_sparse_diracs(xproj);
title( strcat(['L1 recovery, error = ' num2str(err_bp, 3)]));

%% Iterative Soft Thresholding for L1 Minimization
% If the value of |lambda| is known, one can use an iterative algorithm to
% find the solution |xbp| of the L1 minimization.

%%
% We use the optimal value of |lambda| already computed.
% We also save the true solution computed by homotopy.

lambda = lambda_opt;
xbp_opt = xbp;

%%
% The gradient descent step size depends on the conditionning of the
% matrix.

tau = 2/norm(D).^2;

%%
% The iterative algorithm starts with the zero vector.

xbp = zeros(p,1);

%%
% The gradient step updates the value of the solution by decaying the value of |norm(y-D*xbp)^2|.

xbp = xbp + tau*D'*( y-D*xbp );

%%
% The thresholding step improves the sparsity of the solution.

xbp_prev = xbp; 
xbp = perform_thresholding( xbp, tau*lambda, 'soft' );
% display
clf;
subplot(2,1,1);
plot(xbp_prev); axis('tight');
set_graphic_sizes([], 20);
title('Estimate before soft thresholding');
subplot(2,1,2);
plot(xbp); axis('tight');
set_graphic_sizes([], 20);
title('Estimate after soft thresholding');

%%
% _Exercice 3:_ (the solution is <../private/sparsity_seismic_bp/exo3.m exo3.m>)
% Implement the iterative soft thresholding algorithm by applying many
% times these two steps. Keep track of the variational energy minimized by
% this algorithm. Keep also track of the distance to the solution.

exo3;

##### SOURCE END #####
-->
   </body>
</html>