threshold.wp(wp, levels = 3:(wd$nlevels - 1), dev = madmad, policy = "universal", value = 0, by.level = F, type = "soft", verbose = F, return.threshold = F, cvtol = 0.01, cvnorm = l2norm, add.history = T)
nlevels(wp)-1
inclusive.
Only the levels in this vector contribute to the computation of the
threshold and its application.
var()
function. A popular, useful and
robust alternative is the madmad function.
"universal"
and
"mannum"
.
The policies are described in detail below.
policy="manual"
then value
is the actual
threshold value.
levels
. If TRUE a threshold
is computed and applied separately to each scale level.
"hard"
or "soft"
.
T
then a history statement is added to the object for
displaying.
return.threshold
option is set to TRUE then the threshold
values will be returned rather than the thresholded object.
The basic idea of thresholding is very simple. In a signal plus noise
model the wavelet packet transform of signal is very sparse, the wavelet
packet
transform of noise is not (in particular, if the noise is iid Gaussian
then so if the noise contained in the wavelet packet coefficients [within a
particular orthogonal basis chosen from the packet tree. In general, the
noise is not independent]). Thus since
the signal gets concentrated in the wavelet coefficients and the noise
remains "spread" out it is "easy" to separate the signal from noise
by keeping large coefficients (which correspond to signal) and
delete the small ones (which correspond to noise). However, one has
to have some idea of the noise level (computed using the dev
option in threshold functions). If the noise level is very large then
it is possible, as usual, that no signal "sticks up" above the noise.
There are many components to a successful thresholding procedure. Some components have a larger effect than others but the effect is not the same in all practical data situations. Here we give some rough practical guidance, although you must refer to the papers below when using a particular technique. You cannot expect to get excellent performance on all signals unless you fully understand the rationale and limitations of each method below. I am not in favour of the "black-box" approach. The thresholding functions of WaveThresh3 are not a black box: experience and judgement are required!
Some issues to watch for:
levels = 3:(wd$nlevels - 1)
for the
levels
option most certainly does not work globally for
all data problems and situations. The level at which thresholding begins
(i.e. the given threshold and finer scale wavelets) is called the
primary resolution and is unique to a particular problem.
In some ways choice of the primary resolution is very similar to choosing
the bandwidth in kernel regression albeit on a logarithmic scale.
See Hall and Patil, (1995) and
Hall and Nason (1997) for more information.
For each data problem you need to work out which is the best
primary resolution. This can be done by gaining experience at what works
best, or using prior knowledge. It is possible to "automatically" choose
a "best" primary resolution using cross-validation (but not in WaveThresh).
Secondly the levels argument computes and applies the threshold at the
levels specified in the levels
argument. It does this for
all the levels specified. Sometimes, in wavelet shrinkage, the threshold
is computed using only the finest scale coefficients (or more precisely
the estimate of the overall noise level). If you want your threshold
variance estimate only to use the finest scale coefficients (e.g.
with universal thresholding) then you will have to apply the
threshold.wp
function twice. Once (with levels set equal to
nlevels(wd)-1
and with
return.threshold=TRUE
to return the threshold computed on
the finest scale and then apply the threshold function with the
manual
option supplying the value of the previously computed
threshold as the value
options.
# # Generate some test data # test.data <- example.1()$y tsplot(test.data) # # Generate some noisy data # ynoise <- test.data + rnorm(512, sd=0.1) # # Plot it # tsplot(ynoise) # # Now take the discrete wavelet packet transform # N.b. I have no idea if the default wavelets here are appropriate for # this particular example. # ynwp <- wp(ynoise) # # Now do thresholding. We'll use a universal policy, # and madmad deviance estimate on the finest # coefficients and return the threshold. We'll also get it to be verbose # so we can watch the process. # ynwpT1 <- threshold(ynwp, policy="universal", dev=madmad) # # This is just another wp object. Is it sensible? # Probably not as we have just thresholded the scaling function coefficients # as well. So the threshold might be more sensibly computed on the wavelet # coefficients at the finest scale and then this threshold applied to the # whole wavelet tree??