Signal Denoising with Wavelets
This numerical tour uses wavelets to perform signal denoising using thresholding estimators.
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 Scilab) or numericaltour.m (in Matlab) to write all the Scilab/Matlab commands 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/');
Loading a Signal and Making Noise
Here we consider a simple additive Gaussian noise.
First we load the 1D signal.
% size n = 2048; % clean signal name = 'piece-regular'; x0 = load_signal(name, n); x0 = rescale(x0,.05,.95);
Then we add some Gaussian noise.
sigma = 0.07; x = x0 + randn(size(x0))*sigma;
Display.
clf; subplot(2,1,1); plot(x0); axis([1 n 0 1]); title('Clean signal'); subplot(2,1,2); plot(x); axis([1 n 0 1]); title('Noisy signal');
Hard Thresholding vs. Soft Thresholding
A thresholding is a 1D non-linear function applied to each wavelet coefficients. The most important thresholding are the hard thresholding (related to L0 minimization) and the soft thresholding (related to L1 minimization).
% threshold value T = 1; v = -linspace(-3,3,2000); % hard thresholding of the t values v_hard = v.*(abs(v)>T); % soft thresholding of the t values v_soft = max(1-T./abs(v), 0).*v; % display clf; hold('on'); plot(v, v_hard); plot(v, v_soft, 'r--'); axis('equal'); axis('tight'); legend('Hard thresholding', 'Soft thresholding'); hold('off');
Wavelet Denoising with Hard Thesholding
It is possible to perform non linear denoising by thresholding the wavelet coefficients. This allows to better respect the sharp features of the signal.
First we compute the wavelet coefficients of the noisy signal.
options.ti = 0; Jmin = 4; xW = perform_wavelet_transf(x,Jmin,+1,options);
Then we hard threshold the coefficients below the noise level. In practice a threshold of 3*sigma is close to optimal for piecewise regular signals.
T = 3*sigma; xWT = perform_thresholding(xW,T,'hard'); clf; subplot(2,1,1); plot_wavelet(xW,Jmin); axis([1 n -1 1]); title('Noisy coefficients'); set_axis(0); subplot(2,1,2); plot_wavelet(xWT,Jmin); axis([1 n -1 1]); title('Thresholded coefficients'); set_axis(0);
One can then reconstruct from these noisy coefficients.
xhard = perform_wavelet_transf(xWT,Jmin,-1,options);
We can display and measure performance using SNR.
clf; subplot(2,1,1); plot(x); axis([1 n 0 1]); title('Noisy'); subplot(2,1,2); plot(xhard); axis([1 n 0 1]); title(strcat(['Hard denoising, SNR=' num2str(snr(x0,xhard))]));
Wavelet Denoising with Soft Thesholding
Is is possible to replace the hard decision performed by hard thresholding by a softer one, using soft thresholding. For 1D signal, this does not improve the recovery quality. For natural images, on contrary, this tends to improve the results.
- First one must use a lower threshold because soft thresholding also lower the value of non-thresholded coeffcients. In practice, a threshold of 3/2*sigma works well.
- The low frequency part of the coefficients must not be thresholded.
T = 1.5*sigma; xWT = perform_thresholding(xW,T,'soft'); % re-inject the low frequencies xWT(1:2^Jmin) = xW(1:2^Jmin); % re-construct xsoft = perform_wavelet_transf(xWT,Jmin,-1,options);
We can display the results.
clf; subplot(2,1,1); plot(xhard); axis([1 n 0 1]); title(strcat(['Hard denoising, SNR=' num2str(snr(x0,xhard))])); subplot(2,1,2); plot(xsoft); axis([1 n 0 1]); title(strcat(['Soft denoising, SNR=' num2str(snr(x0,xsoft))]));
Exercice 1: (the solution is exo1.m) Determine the best threshold T for both hard and soft thresholding. Test several T values in [.8*sigma, 4.5*sigma]. Do not forget to keep the low frequency wavelet coefficients. What can you conclude from these results ? Test with another image.
exo1;
Translation Invariant Wavelet Transform
Orthogonal wavelet transforms are not translation invariant. It means that the processing of an image and of a translated version of the image give different results. A translation invariant wavelet transform is implemented by ommitting the sub-sampling at each stage of the transform. This correspond to the decomposition of the image in a redundant familly of N*(J+1) atoms where N is the number of pixel and J is the number of scales of the transforms.
For Scilab, we need to extend a little the available memory.
extend_stack_size(4);
The invariant transform is obtained using the same function, by activating the switch options.ti=1.
options.ti = 1; xW = perform_wavelet_transf(x0,Jmin,+1,options);
xW(:,:,1) corresponds to the low scale residual. Each xW(:,1,j) corresponds to a scale of wavelet coefficient, and has the same size as the original signal.
nJ = size(xW,3); clf; subplot(5,1, 1); plot(x0); axis('tight'); title('Signal'); i = 0; for j=1:3 i = i+1; subplot(5,1,i+1); plot(xW(:,1,nJ-i+1)); axis('tight'); title(strcat(['Scale=' num2str(j)])); end subplot(5,1, 5); plot(xW(:,1,1)); axis('tight'); title('Low scale');
Translation Invariant Wavelet Denoising
Orthogonal wavelet denoising does not performs very well because of its lack of translation invariance. A much better result is obtained by not sub-sampling the wavelet transform, which leads to a redundant tight-frame.
First we compute the translation invariant wavelet transform
options.ti = 1; xW = perform_wavelet_transf(x,Jmin,+1,options);
Then we threshold the set of coefficients.
T = 3.5*sigma;
xWT = perform_thresholding(xW,T,'hard');
We can display some wavelets coefficients
J = size(xW,3)-2; clf; subplot(2,1,1); imageplot(xW(:,:,J)); title('Noisy coefficients'); subplot(2,1,2); imageplot(xWT(:,:,J)); title('Thresholded coefficients');
We can now reconstruct.
xti = perform_wavelet_transf(xWT,Jmin,-1,options);
Display the results.
clf; subplot(2,1,1); plot(xhard); axis([1 n 0 1]); title(strcat(['Hard orthogonal, SNR=' num2str(snr(x0,xhard))])); subplot(2,1,2); plot(xti); axis([1 n 0 1]); title(strcat(['Hard invariant, SNR=' num2str(snr(x0,xti))]));
Exercice 2: (the solution is exo2.m) Determine the best threshold T for both hard and soft thresholding, but now in the translation invariant case. Do not forger not to threshold the low scale wavelets coefficients. What can you conclude ?
exo2;
