Sparse Deconvolution
This numerical tour explores the use of sparse regularization to solve 1D deconvolution problems. We show applications to sparse spikes deconvolution (band pass filter applied to a train of Diracs) and piece-wise smooth signal deblurring (low pass filter).
Contents
Installing toolboxes and setting up the path.
You need to download the general purpose toolbox and the signal toolbox.
You need to unzip these toolboxes in your working directory, so that you have toolbox_general/ and toolbox_signal/ in your directory.
For Scilab user: you must replace the Matlab comment '%' by its Scilab counterpart '//'.
Recommandation: You should create a text file named for instance numericaltour.sce (in Scilabe) or numericaltour.m to write all the Scilab/Matlab command you want to execute. Then, simply run exec('numericaltour.sce'); (in Scilab) or numericaltour; (in Matlab) to run the commands.
Execute this line only if you are using Matlab.
getd = @(p)path(path,p); % scilab users must *not* execute this
Then you can add these toolboxes to the path.
% Add some directories to the path getd('toolbox_signal/'); getd('toolbox_general/');
Signal Deconvolution Problem on Synthetic Sparse Data
This section introduces the 1D deconvolution problem.
We generate a synthetic sparse signal, with only a small number of non zero coefficients.
n = 1024; % number of diracts p = round(n*.03); % location of the diracs sel = randperm(n); sel = sel(1:p); % signal f = zeros(n,1); f(sel) = 1; f0 = f .* sign(randn(n,1)) .* (1-.3*rand(n,1));
Display the signal as spike train.
clf; plot_sparse_diracs(f0);
We load a seismic filter, which is a second derivative of a Gaussian.
% dimension of the signal n = 1024; % width of the filter s = 5; % 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); % recenter the filter for fft use h1 = fftshift(h);
To ease implementation, we define a filtering function.
filter = @(u)real(ifft(fft(h1).*fft(u)));
Display the filter and its Fourier transform.
% Fourier transform (normalized) hf = real(fftshift(fft(h1))) / sqrt(n); % display q = 200; clf; subplot(2,1,1); plot(-n/2+1:n/2, h); axis([-q q min(h) max(h)]); title('Filter, Spacial (zoom)'); subplot(2,1,2); plot(-n/2+1:n/2, hf); axis([-q q 0 max(hf)]); title('Filter, Fourier (zoom)');
A blur is a linear operator corresponding to a convolution with filter h. It can be computed using a Fourier transform. Since the filter is symmetric, the operator is symmetric.
We compute blurry noisy measurements y=h*f0+w where w is the noise,.
% noise level sigma = .06; % noise w = sigma*randn(n,1); % blurring + noise y = filter(f0) + w;
Display signals and measurements.
clf; subplot(2,1,1); plot_sparse_diracs(f0); axis('tight'); title('Signal f0'); subplot(2,1,2); plot(y); axis('tight'); title('Measurements y=h*f0+w');
Pseudo-inverse Calculation
A naive way to retrieve a high definition signal is by inverting the blurring operator. It corresponds to the computation of the pseudo-inverse (which here is equal to the inverse because the blurring is invertible).
Since the blurring is diagonalized by the Fourier transform, it can be computed by inverting the Fourier transform of the filter.
q = 200; clf; subplot(2,1,1); plot(-n/2+1:n/2, hf); title('Fourier'); axis([-q q 0 max(hf)]); subplot(2,1,2); plot(-n/2+1:n/2, 1./hf); title('1/Fourier'); axis([-q q 0 100]);
We apply the inversion. The inverse is a high frequency boosting filter that is highly unstable.
ih1 = 1 ./ fft(h1); % avoid numerical explosion ih1(fft(h1)<1e-1) = 0; % do the inversion fPi = ifft( fft(y) .* ih1 );
One can sees that the noise totally blows up after inversion.
clf; subplot(2,1,1); plot(f0); axis('tight'); title('Signal'); subplot(2,1,2); plot(fPi); axis('tight'); title('Pseudo-inverse');
Sparse Deconvolution
This section introduces a sparse L1 regularization and uses iterative thresholdings to compute the solution.
To recover a good approximation of y one should performed a regularized inversion to avoid the noise to explode. This leads to consider the minimization of the L1 norm
min_f 1/2*|h \star f-y^2 + lambda*norm(f,1)|
Where lambda controls the amplitude of the regularization and \star is the convolution operator.
The energy is minimized by iterative thresholding. The iteration is
fSp = S_T( fSp - tau * h \star (h \star f - y) )
where S_T is a soft thresholding and the threshod is T=lambda*tau.
Set up parameters for the regularization. The gradient step size should be smaller than 2/norm(H)^2, where H is the filtering operator by h. This norm is the maximum of the Fourier transform of the filter.
% gradient descent step size tau = 1.95 / max(abs(fft(h)))^2; % regularization strenght lambda = .1; % initialization fSp = y;
For such an ill-posed inverse problem, the number of iterations needs to be very large.
niter = 2000;
The first step, like for pseudo inverse computations is a gradient descent step.
fSp = fSp - tau*filter( filter(fSp) - y );
The second step is a soft thresholding in a wavelet basis. The threshold must be set to lambda*tau.
fSp = perform_thresholding(fSp, lambda*tau, 'soft');
Exercice 1: (the solution is exo1.m) Compute the solution of L1 deconvolution. Keep track of the energy. Also keep track of the decay of the error err(i)=norm(fSp-fSpC)/norm(fSpC) where fSpC is the solution at convergence.
exo1;
Display the decay of the convergence error, in log plot. This shows the linear convergence speed.
clf; hh = plot(log10(err)); axis('tight'); set(hh, 'LineWidth', 2); set_label('Iteration #', 'log_{10}(error)');
Display the result.
clf; subplot(2,1,1); plot(f0); axis('tight'); title(['Signal f0']); subplot(2,1,2); plot(fSp); axis('tight'); title(['Sparse regularization, SNR=' num2str(snr(f0,fSp)) 'dB']);
Exercice 2: (the solution is exo2.m) Try the minimzation for several values of lambda, to understand the influence of this parameter on the sparsity of the solution, and on the values of spikes.
exo2;
Debiasing of Minimum L1 Solution
An issue with L1 minimization is that the value of the spikes underestimate the true value one would like to recover. This is due to the bias introduced by soft thresholding. As a result, the solution obtained by L1 minimization often has a quite poor SNR.
A natural way to remove this bias is to use the solution fSp of the L1 optimization as an indicator of the location of the spikes, and then optimize the value of the spikes at this location using another algorithm.
The location of the spikes is obtained by simple thresholding. We use a small threshold.
T = 1e-2;
The debiased value fDeb is obtained from fSp by minimizing norm(H(:,I)*fDeb(I) - y) where H is the matrix corresponding to the filtering by h, and I denote the indexes of the non-zero entries of fSp .
This quadratic minimization can be computed by an iterative algorithm, that use a projected gradient descent.
We initialize the debiased signal.
fDeb = fSp;
We perform a gradient descent.
fDeb = fDeb - tau*filter( filter(fDeb) - y );
We project on the support of the sparse L1 signal.
fDeb(abs(fSp)<T) = 0;
Exercice 3: (the solution is exo3.m) Implement the debiasing algorithm. Display the error err(i) at iteration i and the energy norm(H(:,I)*fDeb(I) - y).
exo3;
Display the decay of the convergence error, in log plot. This shows the linear convergence speed.
clf; hh = plot(log10(err)); axis('tight'); set(hh, 'LineWidth', 2); set_label('Iteration #', 'log_{10}(error)');
Display the results.
clf; subplot(2,1,1); plot(fSp); axis([1 n -.8 .8]); title('L1'); subplot(2,1,1); plot(fDeb); axis([1 n -.8 .8]); title('Debiased');
Exercice 4: (the solution is exo4.m) Try for several values of lambda to find the optimal denoising result fDeb0 by maximize the SNR. Each time, you need to solve the L1 minimization, and then perform the debiasing. You might want to use values of lambda between .2 and 0. For each value of lambda, you can re-use the result obtained at the previous iteration to initialize the new set of iterations. Record the optimal value lambda of the regularization parameter.
exo4;
Display results.
clf; subplot(2,1,1); plot(f0); axis([1 n -1 1]); title('Signal f0'); subplot(2,1,2); plot(fDeb0); axis('tight'); title(['lambda=' num2str(lambda), ', SNR=' num2str(snr(f0,fDeb0)) 'dB']);
Signal Deconvolution Problem on Piecewise Smooth Signals
This section use a more realistic signal.
First we load a 1D signal.
n = 1024;
f0 = load_signal('piece-regular', n);
f0 = rescale(f0);
We load a low pass Gaussian filter.
% width of the filter s = 6; % second derivative of Gaussian t = (-n/2:n/2-1)'; h = exp( -(t.^2)/(2*s^2) ); h = h/sum(h); % recenter the filter for fft use h1 = fftshift(h); filter = @(u)real(ifft(fft(h1).*fft(u)));
Display the filter and its Fourier transform.
% Fourier transform (normalized) hf = real(fftshift(fft(h1))) / sqrt(n); % display q = 100; clf; subplot(2,1,1); plot(-n/2+1:n/2, h); axis([-q q 0 max(h)]); title('Filter, Spacial (zoom)'); subplot(2,1,2); plot(-n/2+1:n/2, hf); axis([-q q 0 max(hf)]); title('Filter, Fourier (zoom)');
Then we compute blurry noisy measurements.
% noise level sigma = .02; % measurements y = filter(f0) + sigma*randn(n,1);
Display signals and measurements.
clf; subplot(2,1,1); plot(f0); axis('tight'); title('Signal f0'); subplot(2,1,2); plot(y); axis('tight'); title('Measurements y=h*f0+w');
Sparsity in Wavelets for Deconvolution
This section introduces a sparse regularization over wavelets coefficients and uses iterative thresholdings to compute the solution.
Real signals f0 are not directly sparse over the spacial domain, one should use the sparsity over a transformed domain, typically wavelets.
To recover a good approximation of y one should performed a regularized inversion to avoid the noise to explode. This leads to consider the minimization of
min_f 1/2*|h \star f-y^2 + lambda*norm(W*f,1)|
Where W represents an orthogonal wavelet transform and lambda controls the amplitude of the regularization.
The energy is minimized by iterative thresholding. The iteration is
fSp = S_T( fSp - tau * h \star (h \star f - y) )
where S_T is a wavelet soft thresholding and the threshod is T=lambda*tau.
Parameters for the wavelet transform.
% minimum scale Jmin = 1; % the wavelet transform is biorthogonal options.filter = '7-9'; % non redundant options.ti = 0;
Set up parameters for the regularization. The gradient step size should be smaller than 2/norm(H)^2 where H is the convolution operator.
% gradient descent step size tau = .95 / max(abs(fft(h)))^2; % regularization strenght lambda = .03; % number of iterations niter = 500; % initialization fSp = y;
The first step, like for pseudo inverse computations is a gradient descent step.
fSp = fSp - tau*filter( filter(fSp) - y);
The second step is a soft thresholding in a wavelet basis. The threshold must be set to lambda*tau.
% forward transform fSpW = perform_wavelet_transf(fSp, Jmin, +1, options); % thresholding fSpW = perform_thresholding(fSpW, lambda*tau, 'soft'); % forward transform fSp = perform_wavelet_transf(fSpW, Jmin, -1, options);
Exercice 5: (the solution is exo5.m) Compute the solution of L1 regularization in a wavelet basis. Keep track of the energy. Also keep track of the decay of the error err(i)=norm(fSp-fSpC)/norm(fSpC) where fSpC is the solution at convergence.
exo5;
Warning: the energy might not exactly strictly decaying. This is because the wavelet transform is not exactly orthogonal (7/9 wavelet transform), so that the iterative thresholding is not exactly minimizing the L1 norm of the wavelet coefficients.
Display the decay of the convergence error, in log plot. This shows the linear convergence speed.
clf; plot(log10(err)); axis('tight'); set_label('Iteration #', 'log_{10}(error)');
Display the result.
clf; subplot(2,1,1); plot(y); axis('tight'); title('Signal'); subplot(2,1,2); plot(fSp); axis('tight'); title(['Sparse regularization, SNR=' num2str(snr(f0,fSp)) 'dB']);
Exercice 6: (the solution is exo6.m) Try for several values of lambda to find the optimal denoising result fSp0 by maximize the SNR. You might want to use values of lambda between .01 and 0.
exo6;
Display results.
clf; subplot(2,1,1); plot(f0); axis('tight'); title('Signal f0'); subplot(2,1,2); plot(fSp0); axis('tight'); title(['Sparse regularization, SNR=' num2str(snr(f0,fSp0)) 'dB']);
