wmtsa (version 1.1-1)

wavShrink: Nonlinear denoising via wavelet shrinkage

Description

Performs a decimated or undecimated discrete wavelet transform on the input series and "shrinks" (decreases the amplitude towards zero) the wavelet coefficients based on a calculated noise threshold and specified shrinkage function. The resulting shrunken set of wavelet transform coefficients is inverted in a synthesis operation, resulting in a denoised version of the original series.

Usage

wavShrink(x, wavelet="s8",
    n.level=ilogb(length(x), base=2),
    shrink.fun="hard", thresh.fun="universal", threshold=NULL,
    thresh.scale=1, xform="modwt", noise.variance=-1.0,
    reflect=TRUE)

Arguments

x

a vector containing a uniformly-sampled real-valued time series.

n.level

the number of decomposition levels, limited to floor(logb(length(x),2)). Default: floor(logb(length(x),2)).

noise.variance

a numeric scalar representing (an estimate of) the additive Gaussian white noise variance. If unknown, setting this value to 0.0 (or less) will prompt the function to automatically estimate the noise variance based on the median absolute deviation (MAD) o

reflect

a logical value. If TRUE, the last $L_J = (2^{\mbox{n.level}} - 1)(L - 1) + 1$ coefficients of the series are reflected (reversed and appended to the end of the series) in order to attenuate the adverse effect of circular filter operations on

shrink.fun

a character string denoting the shrinkage function. Choices are "hard", "soft", and "mid". Default: "hard".

thresh.fun

a character string denoting the threshold function to use in calculating the waveshrink thresholds. [object Object],[object Object]

Note: if xform == "modwt", then only the "universal" threshold function is (currently) su

thresh.scale

a positive valued numeric scalar which is used to amplify or attenuate the threshold values at each decomposition level. The use of this argument signifies a departure from a model driven estimate of the thresholds and can be used to tweak the levels to o

threshold

explicit setting of the wavelet shrinkage thresholds, one for each level of the decomposition. If a single threshold is given, it is replicated appropriately and (if the chosen transform is additionally a MODWT then) these thresholds are normalized by div

wavelet

a character string denoting the filter type. See wavDaubechies for details. Default: "s8".

xform

a character string denoting the wavelet transform type. Choices are "dwt" and "modwt" for the discrete wavelet transform (DWT) and maximal overlap DWT (MODWT), respectively. The DWT is a decimated transform where (at each level)

Value

• vector containing the denoised series.

concept

nonlinear noise reductionwavelet

Details

Assume that an appropriate model for our time series is $\mathbf{X}=\mathbf{D} + \mathbf{\epsilon}$ where $\mathbf{D}$ represents an unknown deterministic signal of interest and $\mathbf{\epsilon}$ is some undesired stochastic noise that is independent and identically distributed and has a process mean of zero. Waveshrink seeks to eliminate the noise component $\mathbf{\epsilon}$ of $\mathbf{X}$ in hopes of obtaining (a close approximation to) $\mathbf{D}$. The basic algorithm works as follows:

[object Object],[object Object],[object Object]

This function support different shrinkage methods and threshold estimation schemes. Let $W$ represent an arbitrary DWT coefficient and $W^{\mbox{(t)}}$ the correpsonding thresholded coefficient using a threshold of $\delta$. The supported shrinkage methods are

[object Object],[object Object],[object Object]

Hard thresholding reduces to zero all coefficients that do not exceed the threshold. Soft thresholding pushes toward zero any coefficient whose magnitude exceeds the threshold, and zeros the coefficient otherwise. Mid thresholding represents a compromise between hard and soft thresholding such that coefficients whose magnitude exceeds twice the threshold are not adjusted, those between the threshold and twice the trhreshold are shrunk, and those below the threshodl are zeroed.

The threshold is selected based on a model of the noise. The supported techniques for estimating the noise threshold are

[object Object],[object Object],[object Object]

Finally, the user has the choice of using either a decimated (standard) form of the discrete wavelet transform (DWT) or an undecimated version of the DWT (known as the Maximal Overlap DWT (MODWT)). Unlike the DWT, the MODWT is a (circular) shift-invariant transform so that a circular shift in the original time series produces an equivalent shift of the MODWT coefficients. In addition, the MODWT can be interpreted as a cycle-spun version of the DWT, which is achieved by averaging over all non-redundant DWTs of shifted versions of the original series. The z is a smoother version of the DWT at the cost of an increase in computational complexity (for an N-point series, the DWT requires $O(N)$ multiplications while the MODWT requires $O(N\log_2N$ multiplications.

References

Donoho, D. and Johnstone, I. Ideal Spatial Adaptation by Wavelet Shrinkage. Technical report, Department of Statistics, Stanford University, 1992.

Donoho, D. and Johnstone, I. Adapting to Unknown Smoothness via Wavelet Shrinkage. Technical report, Department of Statistics, Stanford University, 1992.

D. B. Percival and A. T. Walden, Wavelet Methods for Time Series Analysis, Cambridge University Press, 2000.

See Also

wavDaubechies, wavDWT, wavMODWT.

Examples

## MODWT waveshrinking using various thresh.scale 
## values on sunspots series 
x <- as.vector(sunspots)
tt <- as.numeric(time(sunspots))
thresh <- seq(0.5,2,length=4)
ws <- lapply(thresh, function(k,x)
    wavShrink(x, wavelet="s8",
       shrink.fun="hard", thresh.fun="universal",
       thresh.scale=k, xform="modwt"), x=x)

stackPlot(x=tt, y=data.frame(x, ws),
    ylab=c("sunspots",thresh),
    xlab="Time")

## DWT waveshrinking using various threshold 
## functions 
threshfuns <- c("universal", "minimax", "adaptive")
ws <- lapply(threshfuns, function(k,x)
    wavShrink(x, wavelet="s8",
       thresh.fun=k, xform="dwt"), x=x)

stackPlot(x=tt, y=data.frame(x, ws),
ylab=c("original", threshfuns),
xlab="Normalized Time")