For a broad class of functions (signals, images) that possess
certain smoothness properties, wavelet techniques are optimal or near
optimal for function recovery.

Specifically, the method is efficient for families of functions f that
have only a few nonzero wavelet coefficients. These functions have
a sparse wavelet representation. For example, a smooth function almost
everywhere, with only a few abrupt changes, has such a property.

The general wavelet–based method for denoising and nonparametric
function estimation is to transform the data into the wavelet domain,
threshold the wavelet coefficients, and invert the transform.

You can summarize these steps as:

Decompose

Choose a wavelet and a level N. Compute the
wavelet decomposition of the signal s down to level N.

Threshold detail coefficients

For each level from 1 to N, threshold 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.

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.

Option

Threshold
Selection Rule

'rigrsure'

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

'sqtwolog'

Fixed form (universal) threshold equal to

2ln(N)

with N the length of the signal.

'heursure'

Selection using a mixture of the first two options

'minimaxi'

Selection using minimax principle

Option 'rigrsure' uses for the
soft threshold estimator a threshold selection rule based on Stein's
Unbiased Estimate of Risk (quadratic loss function). You get an estimate
of the risk for a particular threshold value t.
Minimizing the risks in t gives a selection of
the threshold value.

Option 'sqtwolog' uses a fixed form threshold yielding minimax performance multiplied
by a small factor proportional to log(length(s)).

Option 'heursure' is a mixture
of the two previous options. As a result, if the signal-to-noise ratio
is very small, the SURE estimate is very noisy. So if
such a situation is detected, the fixed form threshold is used.

Option 'minimaxi' uses a fixed
threshold chosen to yield minimax performance
for mean square error against an ideal procedure. The minimax principle
is used in statistics to design estimators. Since the denoised signal
can be assimilated to the estimator of the unknown regression function,
the minimax estimator is the option that realizes the minimum, over
a given set of functions, of the maximum mean square error.

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

For Stein's Unbiased Risk Estimate (SURE) and minimax thresholds,
approximately 3% of coefficients are retained. In the case of the
universal threshold, all values are rejected.

We know 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, which are standard Gaussian white noises.

After you use thselect to
determine a threshold, you can threshold each level of a . This second
step can be done using wthcoef, directly handling
the wavelet decomposition structure of the original signal s.

Soft or Hard Thresholding

Hard and soft thresholding are examples of shrinkage rules.
After you have determined your threshold, you have to decide how to
apply that threshold to your data.

The simplest scheme is hard thresholding.
Let T denote the threshold and x your
data. The hard thresholding is

η(x)={x|x|≥T0|x|<T

The soft thresholding is

η(x)={x−Tx>T0|x|≤Tx+Tx<−T

You can apply your threshold using the hard or soft rule with wthresh.

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.

Option

Corresponding
Model

'one'

Basic model

'sln'

Basic model with unscaled noise

'mln'

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.

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

xd = wdencmp(opt,x,wav,n,thr,sorh,keepapp);

where

opt = 'gbl' and thr is
a positive real number for uniform threshold.

opt = 'lvd' and thr is
a vector for level dependent threshold.

keepapp = 1 to keep approximation
coefficients, as previously and

keepapp = 0 to allow approximation
coefficients thresholding.

x is the signal to be denoised
and wav, n, sorh are
the same as above.

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.

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 2450.

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.

But the noise variance can vary with time. There are several
different variance values on several time intervals. The values as
well as the intervals are unknown.

Let us focus on the problem of estimating the change points
or equivalently the intervals. The algorithm used is based on an original
work of Marc Lavielle about detection of change points using dynamic
programming (see [Lav99] in References).

Let us generate a signal from a fixed-design regression model
with two noise variance change points located at positions 200 and
600.

The reconstructed detail at level 1 recovered at this stage
is almost signal free. It captures the main features of the noise
from a change points detection viewpoint if the interesting part of
the signal has a sparse wavelet representation. To remove almost all
the signal, we replace the biggest values by the mean.

Step 2. To remove almost all
the signal, replace 2% of biggest values by the mean.