Denoising Signals and Images
This example discusses the problem of signal recovery from noisy data. The general denoising procedure involves three steps. The basic version of the procedure follows the steps described below:
Decompose: Choose a wavelet, choose a level N. Compute the wavelet decomposition of the signal at level N.
Threshold detail coefficients: For each level from 1 to N, select a threshold and apply soft thresholding to the detail coefficients.
Reconstruct: Compute wavelet reconstruction using the original approximation coefficients of level N and the modified detail coefficients of levels from 1 to N.
Two points must be addressed in particular:
how to choose the threshold,
and how to perform the thresholding.
Soft or Hard Thresholding?
Thresholding can be done using the function
wthresh which returns soft or hard thresholding of the input signal. Hard thresholding is the simplest method but soft thresholding has nice mathematical properties. Let
thr denote the threshold.
thr = 0.4;
Hard thresholding can be described as the usual process of setting to zero the elements whose absolute values are lower than the threshold. The hard threshold signal is x if x >
thr, and is 0 if x
y = linspace(-1,1,100); ythard = wthresh(y,'h',thr);
Soft thresholding is an extension of hard thresholding, first setting to zero the elements whose absolute values are lower than the threshold, and then shrinking the nonzero coefficients towards 0. The soft threshold signal is sign(x)(x-
thr) if x >
thr and is 0 if x
ytsoft = wthresh(y,'s',thr); subplot(1,3,1) plot(y) title('Original') subplot(1,3,2) plot(ythard) title('Hard Thresholding') subplot(1,3,3) plot(ytsoft) title('Soft Thresholding')
As can be seen in the figure above, the hard procedure creates discontinuities at x = t, while the soft procedure does not.
Threshold Selection Rules
Recalling step 2 of the denoise procedure, the function
thselect performs a threshold selection, and then each level is thresholded. This second step can be done using
wthcoeff, directly handling the wavelet decomposition structure of the original signal. Four threshold selection rules are implemented in the function
thselect. Typically it is interesting to show them in action when the input signal is a Gaussian white noise.
rng default y = randn(1,1000);
Rule 1: Selection using principle of Stein's Unbiased Risk Estimate (SURE)
thr = thselect(y,'rigrsure')
thr = 2.0518
Rule 2: Fixed form threshold equal to sqrt(2*log(length(y)))
thr = thselect(y,'sqtwolog')
thr = 3.7169
Rule 3: Selection using a mixture of the first two options
thr = thselect(y,'heursure')
thr = 3.7169
Rule 4: Selection using minimax principle
thr = thselect(y,'minimaxi')
thr = 2.2163
Minimax and SURE threshold selection rules are more conservative and would be more convenient when small details of the signal lie near the noise range. The two other rules remove the noise more efficiently.
Let us use the "blocks" test data credited to Donoho and Johnstone as the first example. Generate original signal xref and a noisy version x adding a standard Gaussian white noise.
sqrt_snr = 4; % Set signal to noise ratio init = 2055615866; % Set rand seed [xref,x] = wnoise(1,11,sqrt_snr,init);
First denoise the signal using
wdenoise with default settings. Compare the result with the original and noisy signals.
xd = wdenoise(x); subplot(3,1,1) plot(xref) axis tight title('Original Signal') subplot(3,1,2) plot(x) axis tight title('Noisy Signal') subplot(3,1,3) plot(xd) axis tight title('Denoised Signal - Signal to Noise Ratio = 4')
Denoise the noisy signal a second time, this time using soft heuristic SURE thresholding on detail coefficients obtained from the decomposition of x, at level 3 by
sym8 wavelet. Compare with the previous denoised signal.
xd2 = wdenoise(x,3,'Wavelet','sym8',... 'DenoisingMethod','SURE',... 'ThresholdRule','Soft'); figure subplot(3,1,1) plot(xref) axis tight title('Original Signal') subplot(3,1,2) plot(xd) axis tight title('First Denoised Signal: Default Settings') subplot(3,1,3) plot(xd2) axis tight title('Second Denoised Signal')
Since only a small number of large coefficients characterize the original signal, both denoised signals compare well with the original signal. You can use the Wavelet Signal Denoiser to explore the effects other denoising parameters have on the noisy signal.
Dealing with Non-White Noise
When you suspect a non-white noise, thresholds must be rescaled by a level-dependent estimation of the level noise. As a second example, let us try the method on the highly perturbed part of an electrical signal. Let us use
db3 wavelet and decompose at level 3. To deal with the composite noise nature, let us try a level-dependent noise size estimation.
load leleccum indx = 2000:3450; x = leleccum(indx); % Load electrical signal and select part of it.
Denoise the signal using soft fixed form thresholding and level-dependent noise size estimation.
xd = wdenoise(x,3,'Wavelet','db3',... 'DenoisingMethod','UniversalThreshold',... 'ThresholdRule','Soft',... 'NoiseEstimate','LevelDependent'); Nx = length(x); figure subplot(2,1,1) plot(indx,x) axis tight title('Original Signal') subplot(2,1,2) plot(indx,xd) axis tight title('Denoised Signal')
The result is quite good in spite of the time heterogeneity of the nature of the noise after and before the beginning of the sensor failure around time 2410.
The denoising method described for the one-dimensional case applies also to images and applies well to geometrical images. The two-dimensional denoising procedure has the same three steps and uses two-dimensional wavelet tools instead of one-dimensional ones. For the threshold selection, prod(size(y)) is used instead of length(y) if the fixed form threshold is used.
Generate a noisy image.
load woman init = 2055615866; rng(init); x = X + 15*randn(size(X));
In this case fixed form threshold is used with estimation of level noise, thresholding mode is soft and the approximation coefficients are kept.
[thr,sorh,keepapp] = ddencmp('den','wv',x); thr
thr = 107.9838
Denoise image using global thresholding option.
xd = wdencmp('gbl',x,'sym4',2,thr,sorh,keepapp); figure('Color','white') colormap(pink(255)) sm = size(map,1); image(wcodemat(X,sm)) title('Original Image')
figure('Color','white') colormap(pink(255)) image(wcodemat(x,sm)) title('Noisy Image')
image(wcodemat(xd,sm)) title('Denoised Image')