threshold.wd(wd, levels=3:(wd$nlevels-1), type="hard",
policy="universal", by.level=FALSE, value=0, dev=var, boundary=FALSE,
verbose=FALSE, return.threshold=FALSE)"wd", typically from a wavelet decomposition
using the wd function.FALSE then a global threshold is
applied to all the levels specified by levels,
otherwise a threshold is computed and applied separately to each level."manual"
policy, or the user-supplied quantile level for the
"probability" policy.var is to use the
variance, but you can insert your own measure of deviance.threshold() spurts informative
messages at you.Hard thresholding simply compares a coefficient with a threshold. If it is larger in absolute magnitude it is left alone, if it is smaller it is set to zero. The "soft threshold" formula is
soft(w) = sgn(w)*max(|w| - t, 0)
where "w" is the wavelet coefficient, "t" is the threshold and "sgn" is the sign of w. Soft thresholding causes w to be replaced by soft(w).
There are many ways that the threshold can be computed, we term this the "policy". A universal policy computes a threshold based on Donoho and Johnstone's "universal thresholds". The threshold is sqrt(2*log(n))*noise, where noise is computed as sqrt(dev(w)), i.e. some measure of the variability of the coefficients, and n is the number of data points (or number of wavelet coefficients). By default "dev" is "var", so the noise level is estimated using the sample standard deviation. You can use any other such estimate by writing your own function and supplying it as the "dev" argument. For example you could create the function "myvar" by
myvar <- function(d) mad(d)^2
This computes the square of the mean absolute deviation of the data. It is squared because "dev" should be a measure on the variance scale, not the standard deviation (you know what I mean). If you make the "by.levels" argument T, then a separate threshold is computed for each level in the "levels" vector. This means that the variance is estimated separately for each level.
The "manual" policy is simple. You supply a threshold value ("value") and hard or soft thresholding is performed using that value. The "value" argument is a vector. If it is of length 1 then it is replicated to be the same length as the "levels" vector, otherwise it is repeated as many times as necessary to be the length vector's length. In this way, different thresholds can be given for different levels. Note that the "by.level" argument has no effect with this policy.
The "probability" policy works as follows. All coefficients that are smaller than the "value"th quantile of the coefficients are set to zero. If "by.level" is false, then the quantile is computed for all coefficients in the levels specified by the "levels" vector; if "by.level" is true, then each level's quantile is estimated separately.
This function is a method for the generic function
threshold() for class wd.
It can be invoked by calling threshold(x) for an
object x of the appropriate class, or directly by
calling threshold.wd(x) regardless of the
class of the object.
wd for a list.wr, wd# "Standard" example:
example(wd)
#
# Threshold it
#
(thH.wds <- threshold(wds))
thS.wds <- threshold(wds, type="soft")
#
# Reconstruct from the thresholded coefficients
#
str( trecH <- wr(thH.wds))
summary(trecS <- wr(thS.wds))
#
# Plot the reconstructed functions
#
plot(x, trecH, type="l", col=2)
lines(x,trecS, col=3)
lines(x, fx, col="gray", lty=2, lwd=1.5)# true functionRun the code above in your browser using DataLab