Learn R Programming

wmtsa (version 2.0-3)

wavCWTPeaks: Peak detection in a time series via the CWT

Description

Finds the local maxima in a time series via a CWT tree.

Usage

wavCWTPeaks(x, snr.min=3, scale.range=NULL, length.min=10,
    noise.span=NULL, noise.fun="quantile", noise.min=NULL)

Arguments

x

an object of class wavCWTTree.

length.min

the minimum number of points along a CWT tree branch and within the specified scale.range needed in order for that branches peak to be considered a peak candidate. Default: 10.

noise.fun

a character string defining the function to apply to the local noise estimates in order to sumarize and quantify the local noise level into a scalar value. See the DETAILS section for more information. Supported values are

"quantile"

quantile(x, probs=0.95)

"sd"

sd(x)

"mad"

mad(x, center=0)

where x is a vector of smallest-scale CWT coefficients whose time indices are near that of the branch termination time. Default: "quantile".

noise.min

the minimum allowed estimated local noise level. Default: quantile(attr(x,"noise"), prob=0.05), where x is the input wavCWTTree object.

noise.span

the span in time surrounding each branche's temrination point to use in forming local noise estimates and (ultimately) peak SNR estimates. Default: NULL,max(0.01 * diff(range(times)), 5*sampling.interval), where times and sampling.interval are attributes of the input wavCWTTree object.

scale.range

the range of CWT scales that a peak must fall into in order to be considered a peak candidate. Default: scale[range(which(branch.hist > quantile(branch.hist,prob=0.8)))], where branch.hist is an attribute of the input wavCWTTree object. This default results in isolating the bulk of energetic CWT branches, but the user in encouraged to reduce the scale range in order to attenuate the computational burden of the peak detection scheme.

snr.min

the minimum allowed peak signal-to-noise ratio. Default: 3.

Value

a list of x and y vectors identifying the peaks in the original time series. The pruning criteria (snr.min, scale.range, length.min, noise.span, noise.fun, noise.min) are attached are attached as attributes. In addition, a peaks attribute is attached and corresponds to a data.frame containing the following information for each peak:

branch

index of the associated branch in the CWT tree

itime

index location in time

iscale

index location in scale

time

location in time

scale

location in scale

extrema

CWT value

iendtime

index location of branch termination time, i.e., the index of the point in the time series corresponding to the current peak

Details

The local maxima of the CWT are linked together to form so-called branches, where each branch represents one ridge of the CWT time-scale terrain. The collection of branches forms a tree, as output by the wavCWTTree function. The wavCWTpeaks function prunes the branches of the input CWT tree and records the termination time (i.e., the time associated with point of the branch that is closest to scale zero) as the time index associated with the local peak of the corresponding time series. Information regarding the collection of isolated peaks is returned as a data.frame object.

The tree branches are pruned in the following ways:

peak SNR

an estimate of SNR at peak value is greater than or equal to the specified snr.min. A peak SNR estimate is formed as follows: For each branch of the input CWT tree, a subset of CWT coefficients is collected such that the CWT coefficients are both local to the branch termination time and correspond to the smallest analyzed CWT scale. The user specified noise.span argument is used to define the boundaries of each subset in time ala \([B - \)noise.span, \(B + \)noise.span\(]\), where \(B\) is the branch termination time. Each CWT subset is assumed to be representative of the local noise levels near the corresponding branch termination time and noise.fun is used to quantify (and summarize) each level resulting in a scalar $z$. The minimum value of $z$ is specified by the user ala the noise.min argument. Finally, the ratio \(|P|/|W|\) is used to form an estimate of the local signal-to-noise ration (SNR) for the corresponding branch, where \(P\) is the maximum CWT value along the branch in the CWT time-scale plane.

scale

the scale corresponding to the peak is larger than the minimum of the specified scale.range.

branch length

the length of the branch within the specified scale.range is greater than or equal to the specified minimum length.min.

endpoint

the index of the terminating time of the branch is on the interval \((W, N-W)\), where \(N\) is the length of the series and \(W\) is integer equivalent of \(1/4\) the length of the noise.span or \(3\), whichever is greater.

NOTE: For peak detection, the wavelet filters used to form the CWT must maintain an (approximate) zero phase property so that the CWT coefficients can be meaningfully aligned with the events of the original time series. Currently, only the so-called Mexican hat wavelet maintains this property due to the even-symmetry of the filter's impulse response. Therefore, only the Mexican hat wavelet ("gaussian2") is currently supported for CWT-based peak detection. See the wavCWTFilters and wavCWT function for more information.

References

Pan Du, Warren A. Kibbe, and Simon M. Lin, ``Improved peak detection in mass spectrum by incorporating continuous wavelet transform-based pattern matching", Bioinformatics, 22, 2059--2065 (2006).

J.F. Muzy, E. Bacry, and A. Arneodo., ``The multifractal formalism revisited with wavelets.", International Journal of Bifurcation and Chaos, 4, 245--302 (1994).

See Also

wavCWTTree, wavCWT, wavCWTFilters.

Examples

Run this code
# NOT RUN {
## create linchirp series 
linchirp <- make.signal("linchirp")

## calculate the CWT 
W <- wavCWT(linchirp)

## form CWT tree 
z <- wavCWTTree(W)

## estimate the peak locations using default 
## scale.range 
p <- wavCWTPeaks(z)

## plot an overlay of the original series and the 
## peaks 
x <- as(linchirp@positions,"numeric")
y <- linchirp@data
plot(x, y, type="l", xlab="time", ylab="linchirp")
points(p, pch=16, col="red", cex=1.2)
# }

Run the code above in your browser using DataLab