Documentation Center

  • Trial Software
  • Product Updates

Denoising and Nonparametric Function Estimation


The most general 1-D model for this is








load cuspamax;
plot(y); hold on;
axis tight;




You can summarize these steps as:

  1. Decompose


  2. Threshold detail coefficients


  3. Reconstruct


Threshold Selection Rules

The Wavelet Toolbox supports a number of threshold selection rules. Four threshold selection rules are implemented in the thselect. Each rule corresponds to a tptr option in the command

thr = thselect(y,tptr)

which returns the threshold value.


Threshold Selection Rule


Selection using principle of Stein's Unbiased Risk Estimate (SURE)


Fixed form (universal) threshold equal to

with N the length of the signal.


Selection using a mixture of the first two options


Selection using minimax principle





The following example shows the threshold rules for a 1000-by-1 N(0,1) vector. The signal here is

with f(n) = 0.

rng default;
sig = randn(1e3,1);
thr_rigrsure = thselect(sig,'rigrsure')
thr_univthresh = thselect(sig,'sqtwolog')
thr_heursure = thselect(sig,'heursure')
thr_minimaxi = thselect(sig,'minimaxi')
hold on;




Soft or Hard Thresholding



The soft thresholding is


y = linspace(-1,1,100);
thr = 0.4;

Dealing with Unscaled Noise and Nonwhite Noise

Usually in practice the basic model cannot be used directly. We examine here the options available to deal with model deviations in the main de-noising function wden.

The simplest use of wden is

sd = wden(s,tptr,sorh,scal,n,wav)

which returns the denoised version sd of the original signal s obtained using the tptr threshold selection rule. Other parameters needed are sorh, scal, n, and wav. The parameter sorh specifies the thresholding of details coefficients of the decomposition at level n of s by the wavelet called wav. The remaining parameter scal is to be specified. It corresponds to threshold's rescaling methods.


Corresponding Model


Basic model


Basic model with unscaled noise


Basic model with nonwhite noise

  • Option scal = 'one' corresponds to the basic model.

  • In general, you can ignore the noise level and it 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 for two reasons. The first one is that if level 1 coefficients contain f details, then these details are concentrated in a few coefficients if the function f is sufficiently regular. The second reason is to avoid signal end effects, which are pure artifacts due to computations on the edges.

    Option 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 as in the previous option is used by estimating σlev level by level.

    This estimation is implemented in the file wnoisest, directly handling the wavelet decomposition structure of the original signal s.

    Option scal = 'mln' handles threshold rescaling using a level-dependent estimation of the level noise.

For a more general procedure, the wdencmp function performs wavelet coefficients thresholding for both de-noising and compression purposes, while directly handling one-dimensional and two-dimensional data. It allows you to define your own thresholding strategy selecting in





  • keepapp = 1 to keep approximation coefficients, as previously and

  • keepapp = 0 to allow approximation coefficients thresholding.


Denoising in Action

We begin with examples of one-dimensional de-noising methods with the first example credited to Donoho and Johnstone. You can use the following file to get the first test function using wnoise.

% Set signal to noise ratio and set rand seed. 
sqrt_snr = 4; init = 2055615866;

% Generate original signal xref and a noisy version x adding 
% a standard Gaussian white noise. 

% De-noise noisy signal using soft heuristic SURE thresholding 
% and scaled noise option, on detail coefficients obtained 
% from the decomposition of x, at level 3 by sym8 wavelet. 
xd = wden(x,'heursure','s','one',3,'sym8');

Blocks Signal Denoising


As a second example, let us try the method on the highly perturbed part of the electrical signal studied above.

According to this previous analysis, 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 electrical signal and select part of it. 
load leleccum; indx = 2000:3450; 
x = leleccum(indx);

% Find first value in order to avoid edge effects. 
deb = x(1);

% De-noise signal using soft fixed form thresholding 
% and unknown noise option. 
xd = wden(x-deb,'sqtwolog','s','mln',3,'db3')+deb;

Electrical Signal De-Noising


Extension to Image De-Noising

The de-noising method described for the one-dimensional case applies also to images and applies well to geometrical images. A direct translation of the one-dimensional model is

s(i,j) = f(i,j) + σe(i,j)

where e is a white Gaussian noise with unit variance.

The two-dimensional de-noising procedure has the same three steps and uses two-dimensional wavelet tools instead of one-dimensional ones. For the threshold selection, prod(size(s)) is used instead of length(s) if the fixed form threshold is used.

Note that except for the "automatic" one-dimensional de-noising case, de-noising and compression are performed using wdencmp. As an example, you can use the following file illustrating the de-noising of a real image.

% Load original image. 
load  woman

% Generate noisy image.
x = X + 15*randn(size(X));

% Find default values. 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 is equal to estimated_sigma*sqrt(log(prod(size(X))))

% De-noise image using global thresholding option.
xd = wdencmp('gbl',x,'sym4',2,thr,sorh,keepapp);
% Plots.
colormap(pink(255)), sm = size(map,1);
subplot(221), image(wcodemat(X,sm)), title('Original Image')
subplot(222), image(wcodemat(x,sm)), title('Noisy Image')
subplot(223), image(wcodemat(xd,sm)), title('denoised Image')

The result shown below is acceptable.

Image De-Noising

One-Dimensional Wavelet Variance Adaptive Thresholding

The idea is to define level by level time-dependent thresholds, and then increase the capability of the de-noising strategies to handle nonstationary variance noise models.

More precisely, the model assumes (as previously) that the observation is equal to the interesting signal superimposed on a noise (see Denoising and Nonparametric Function Estimation).





% Generate blocks test signal.
x = wnoise(1,10);       

% Generate noisy blocks with change points.
bb = randn(1,length(x));
cp1 = 200; cp2 = 600;

The aim of this example is to recover the two change points from the signal x. In addition, this example illustrates how the GUI tools locate the change points for interval dependent thresholding.

Step 1. Recover a noisy signal by suppressing an approximation.

% Perform a single-level wavelet decomposition 
% of the signal using db3.
wname = 'db3'; lev = 1;

% Reconstruct detail at level 1.



x = sort(abs(det));
ind = find(abs(det)>v2p100);
det(ind) = mean(det);



  • The maximum number of change points is 5.




Step 4. (Optional) Replace the estimated change points.


cp_est = t_est(kopt+1,1:kopt);

You can replace the estimated change points by computing

% cp_New = t_est(knew+1,1:knew); % where 1 ≤ knew ≤ 5
Was this topic helpful?