Automatic 1-D de-noising
XD = wden(X,TPTR,SORH,SCAL,N,
XD = wden(C,L,TPTR,SORH,SCAL,N,
XD = wden(W,'modwtsqtwolog',SORH,'mln',N,WNAME)
[XD,CXD] = wden(...)
[XD,CXD,LXD] = wden(...)
wden is a one-dimensional
wden performs an automatic
de-noising process of a one-dimensional signal using wavelets.
XD = wden(X,TPTR,SORH,SCAL,N, returns
a de-noised version
XD of input signal
by thresholding the wavelet coefficients.
TPTR character vector contains the threshold
'rigrsure' uses the principle of
Stein's Unbiased Risk.
'heursure' is an heuristic variant
of the first option.
'sqtwolog' for the universal threshold
'minimaxi' for minimax thresholding
thselect for more information)
is for soft or hard thresholding (see
SCAL defines multiplicative threshold rescaling:
'one' for no rescaling
'sln' for rescaling using a single estimation
of level noise based on first-level coefficients
'mln' for rescaling done using level-dependent
estimation of level noise
XD = wden(C,L,TPTR,SORH,SCAL,N, returns
the same output arguments, using the same options as above, but obtained
directly from the input wavelet decomposition structure
the signal to be de-noised, at level
N and using
XD = wden(W,'modwtsqtwolog',SORH,'mln',N,WNAME) returns
the denoised signal obtained by operating on the MODWT transform matrix
W is the output of MODWT. You must use the
same wavelet in both
[XD,CXD] = wden(...) returns the denoised
wavelet coefficients. For DWT denoising,
a vector (see
CXD is a matrix with N+1 rows
modwt). The number of
columns is equal to the length of the input signal
[XD,CXD,LXD] = wden(...) returns the number
of coefficients by level for DWT denoising. See
wavedec for details. The
is not supported for MODWT denoising. The additional output arguments
the wavelet decomposition structure (see
more information) of the de-noised signal
% The current extension mode is zero-padding (see dwtmode). % Set signal to noise ratio and set rand seed. snr = 3; init = 2055615866; % Generate original signal and a noisy version adding % a standard Gaussian white noise. [xref,x] = wnoise(3,11,snr,init); % De-noise noisy signal using soft heuristic SURE thresholding % and scaled noise option, on detail coefficients obtained % from the decomposition of x, at level 5 by sym8 wavelet. lev = 5; xd = wden(x,'heursure','s','one',lev,'sym8'); % Plot signals. subplot(611), plot(xref), axis([1 2048 -10 10]); title('Original signal'); subplot(612), plot(x), axis([1 2048 -10 10]); title(['Noisy signal - Signal to noise ratio = ',... num2str(fix(snr))]); subplot(613), plot(xd), axis([1 2048 -10 10]); title('De-noised signal - heuristic SURE'); % De-noise noisy signal using soft SURE thresholding xd = wden(x,'heursure','s','one',lev,'sym8'); % Plot signal. subplot(614), plot(xd), axis([1 2048 -10 10]); title('De-noised signal - SURE'); % De-noise noisy signal using fixed form threshold with % a single level estimation of noise standard deviation. xd = wden(x,'sqtwolog','s','sln',lev,'sym8'); % Plot signal. subplot(615), plot(xd), axis([1 2048 -10 10]); title('De-noised signal - Fixed form threshold'); % De-noise noisy signal using minimax threshold with % a multiple level estimation of noise standard deviation. xd = wden(x,'minimaxi','s','sln',lev,'sym8'); % Plot signal. subplot(616), plot(xd), axis([1 2048 -10 10]); title('De-noised signal - Minimax'); % If many trials are necessary, it is better to perform % decomposition once and threshold it many times: % decomposition. [c,l] = wavedec(x,lev,'sym8'); % threshold the decomposition structure [c,l]. xd = wden(c,l,'minimaxi','s','sln',lev,'sym8'); % Editing some graphical properties, % the following figure is generated.
Denoise a signal consisting of a 2-Hz sine wave with transients at 0.3 and 0.72 seconds. Use Donoho and Johnstone's universal threshold with level-dependent estimation of the noise. Obtain denoised versions using the DWT and MODWT. Compare the results.
N = 1000; t = linspace(0,1,N); x = 4*sin(4*pi*t); x = x - sign(t-.3)-sign(.72 - t); y = x+0.15*randn(size(t)); xdDWT = wden(y,'sqtwolog','s','mln',3,'db2'); xdMODWT = wden(y,'modwtsqtwolog','s','mln',3,'db2'); subplot(2,1,1) plot(xdDWT), title('DWT Denoising'); axis tight; subplot(2,1,2) plot(xdMODWT), title('MODWT Denoising'); axis tight;
Denoise a blocky signal using the Haar wavelet with MODWT and DWT denoising. Compare the L2 and L-infty norms of the difference between the original signal and the denoised versions.
[x,xn] = wnoise('blocks',10,3); xdMODWT = wden(xn,'modwtsqtwolog','s','mln',6,'haar'); xd = wden(xn,'sqtwolog','s','mln',6,'haar'); plot(x) hold on plot(xd,'r--') plot(xdMODWT,'k-.') legend('Original','DWT','MODWT') hold off norm(abs(x-xd),2), norm(abs(x-xd),Inf) norm(abs(x-xdMODWT),2), norm(abs(x-xdMODWT),Inf)
The underlying model for the noisy signal is basically of the following form:
where time n is equally spaced.
In the simplest model, suppose that e(n) is a Gaussian white noise N(0,1) and the noise level σ a is supposed to be equal to 1.
The de-noising objective is to suppress the noise part of the signal s and to recover f.
The de-noising procedure proceeds in three steps:
Decomposition. Choose a wavelet, and
choose a level
N. Compute the wavelet decomposition
of the signal s at level
Detail coefficients thresholding. For
each level from 1 to
N, select a threshold and
apply soft thresholding to the detail coefficients.
Reconstruction. Compute wavelet reconstruction
based on the original approximation coefficients of level
the modified detail coefficients of levels from 1 to
More details about threshold selection rules are in Wavelet Denoising and Nonparametric Function Estimation, in the
User's Guide, and in the help of the
Let us point out that
The detail coefficients vector is the superposition of the coefficients of f and the coefficients of e, and that the decomposition of e leads to detail coefficients that are standard Gaussian white noises.
Minimax and SURE threshold selection rules are more
conservative and are more convenient when small details of function f lie
in the noise range. The two other rules remove the noise more efficiently.
'heursure' is a compromise.
In practice, the basic model cannot be used directly. This section
examines the options available, to deal with model deviations. The
scal has to be specified. It
corresponds to threshold rescaling methods.
to the basic model.
In general, you can ignore the noise level that must be estimated. The detail coefficients CD1 (the finest scale) are essentially noise coefficients with standard deviation equal to σ. The median absolute deviation of the coefficients is a robust estimate of σ. The use of a robust estimate is crucial because if level 1 coefficients contain f details, these details are concentrated in few coefficients to avoid signal end effects, which are pure artifacts due to computations on the edges.
scal = 'sln' handles
threshold rescaling using a single estimation of level noise based
on the first-level coefficients.
When you suspect a nonwhite noise e,
thresholds must be rescaled by a level-dependent estimation of the
level noise. The same kind of strategy is used by estimating σlev level
by level. This estimation is implemented in the file
which handles the wavelet decomposition structure of the original
signal s directly.
threshold rescaling using a level-dependent estimation of the level
Antoniadis, A.; G. Oppenheim, Eds. (1995), Wavelets and statistics, 103, Lecture Notes in Statistics, Springer Verlag.
Donoho, D.L. (1993), "Progress in wavelet analysis and WVD: a ten minute tour," in Progress in wavelet analysis and applications, Y. Meyer, S. Roques, pp. 109–128. Frontières Ed.
Donoho, D.L.; I.M. Johnstone (1994), "Ideal spatial adaptation by wavelet shrinkage," Biometrika, Vol. 81, pp. 425–455.
Donoho, D.L. (1995), "De-noising by soft-thresholding," IEEE Trans. on Inf. Theory, 42 3, pp. 613– 627.
Donoho, D.L.; I.M. Johnstone, G. Kerkyacharian, D. Picard (1995), "Wavelet shrinkage: asymptotia," Jour. Roy. Stat. Soc., series B, Vol. 57, No. 2, pp. 301–369.