wavethresh (version 4.6.1)

threshold.wd: Threshold (DWT) wavelet decomposition object

Description

This function provides various ways to threshold a wd class object.

Usage

## S3 method for class 'wd':
threshold(wd, levels = 3:(nlevelsWT(wd) - 1), type = "soft", policy = "sure",
        by.level = FALSE, value = 0, dev = madmad, boundary = FALSE, verbose = FALSE,
        return.threshold = FALSE, force.sure = FALSE, cvtol = 0.01, Q = 0.05, OP1alpha = 0.05, 
        alpha = 0.5, beta = 1, C1 = NA, C2 = NA, C1.start = 100, ...)

Arguments

wd

The DWT wavelet decomposition object that you wish to threshold.

levels

a vector of integers which determines which scale levels are thresholded in the decomposition. Each integer in the vector must refer to a valid level in the wd object supplied. This is usually any integer from 0

type

determines the type of thresholding this can be "hard" or "soft".

policy

selects the technique by which the threshold value is selected. Each policy corresponds to a method in the literature. At present the different policies are: "universal", "LSuniversal", "sure

by.level

If FALSE then a global threshold is computed on and applied to all scale levels defined in levels. If TRUE a threshold is computed and applied separately to each scale level.

value

This argument conveys the user supplied threshold. If the policy="manual" then value is the actual threshold value; if the policy="mannum" then value conveys the total number of ordered coefficients kept

dev

this argument supplies the function to be used to compute the spread of the absolute values coefficients. The function supplied must return a value of spread on the variance scale (i.e. not standard deviation) such as the var() function. A po

boundary

If this argument is TRUE then the boundary bookeeping values are included for thresholding, otherwise they are not.

verbose

if TRUE then the function prints out informative messages as it progresses.

return.threshold

If this option is TRUE then the actual value of the threshold is returned. If this option is FALSE then a thresholded version of the input is returned.

force.sure

If TRUE then the sure threshold is computed on a vector even when that vector is very sparse. If FALSE then the normal SUREshrink procedure is followed whereby the universal threshold is used for sparse vector

cvtol

Parameter for the cross-validation "cv" policy.

Q

Parameter for the false discovery rate "fdr" policy.

OP1alpha

Parameter for Ogden and Parzen's first "op1" and "op2" policies.

alpha

Parameter for BayesThresh "BayesThresh" policy.

beta

Parameter for BayesThresh "BayesThresh" policy.

C1

Parameter for BayesThresh "BayesThresh" policy.

C2

Parameter for BayesThresh "BayesThresh" policy.

C1.start

Parameter for BayesThresh "BayesThresh" policy.

...

any other arguments

Value

• An object of class wd. This object contains the thresholded wavelet coefficients. Note that if the return.threshold option is set to TRUE then the threshold values will be returned rather than the thresholded object.

RELEASE

Version 3.6 Copyright Guy Nason and others 1997

Details

This function thresholds or shrinks wavelet coefficients stored in a wd object and returns the coefficients in a modified wd object. See the seminal papers by Donoho and Johnstone for explanations about thresholding. For a gentle introduction to wavelet thresholding (or shrinkage as it is sometimes called) see Nason and Silverman, 1994. For more details on each technique see the descriptions of each method below

The basic idea of thresholding is very simple. In a signal plus noise model the wavelet transform of signal is very sparse, the wavelet transform of noise is not (in particular, if the noise is iid Gaussian then so if the noise contained in the wavelet coefficients). 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:

[object Object],[object Object]

References

Various code segments detailed above were kindly donated by Felix Abramovich, Theofanis Sapatinas and Todd Ogden.

See Also

wd, wd.object, wr, wr.wd, threshold.

Examples

#
# Generate some test data
#
test.data <- example.1()$y
ts.plot(test.data)
#
# Generate some noisy data
#
ynoise <- test.data + rnorm(512, sd=0.1)
#
# Plot it
#
ts.plot(ynoise)
#
# Now take the discrete wavelet transform
# N.b. I have no idea if the default wavelets here are appropriate for
# this particular examples. 
#
ynwd <- wd(ynoise)
plot(ynwd)
#
# 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.
#
ynwdT1 <- threshold(ynwd, policy="universal", dev=madmad,
		levels= nlevelsWT(ynwd)-1, return.threshold=TRUE,
	verbose=TRUE)
# threshold.wd:
# Argument checking
# Universal policy...All levels at once
# Global threshold is:  0.328410967430135 
#
# Why is this the threshold? Well in this case n=512 so sqrt(2*log(n)),
# the universal threshold,
# is equal to 3.53223. Since the noise is about 0.1 (because that's what
# we generated it to be) the threshold is about 0.353.
#
# Now let's apply this threshold to all levels in the noisy wavelet object
#
ynwdT1obj <- threshold(ynwd, policy="manual", value=ynwdT1,
	levels=0:(nlevelsWT(ynwd)-1))
#
# And let's plot it
#
plot(ynwdT1obj)
#
# You'll see that a lot of coefficients have been set to zero, or shrunk.
#
# Let's try a Bayesian examples this time!
#
ynwdT2obj <- threshold(ynwd, policy="BayesThresh")
#
# And plot the coefficients
#
plot(ynwdT2obj)
#
# Let us now see what the actual estimates look like
#
ywr1 <- wr(ynwdT1obj)
ywr2 <- wr(ynwdT2obj)
#
# Here's the estimate using universal thresholding
#
ts.plot(ywr1)
#
# Here's the estimate using BayesThresh 
#
ts.plot(ywr2)