cwt: Continuous Wavelet Transform (CWT)

In zeehio/MassSpecWavelet: Peak Detection for Mass Spectrometry data using wavelet-based algorithms

Continuous Wavelet Transform (CWT)

Description

CWT(Continuous Wavelet Transform) with Mexican Hat wavelet (by default) to match the peaks in Mass Spectrometry spectrum

Usage

cwt(ms, scales = 1, wavelet = "mexh")

Arguments

ms	Mass Spectrometry spectrum (a vector of MS intensities)
scales	a vector represents the scales at which to perform CWT. See the Details section. Additionally, a prepared_wavelets object is also accepted (see prepareWavelets()).
wavelet	The wavelet base, Mexican Hat by default. User can provide wavelet Psi(x) as a form of two row matrix. The first row is the x value, and the second row is Psi(x) corresponding to x.

Details

The default mother wavelet is a Mexican Hat evaluated in the [-8,8] range using 1024 points (a step of roughly 1/64):

\psi(x) = \frac{2}{\sqrt{3}} \pi^{-0.25} ( 1 - x^2 ) \exp{-x^2/2}

The \sigma of the mother Mexican Hat is of one x unit.

The scaled wavelet is a downsampled version of the mother wavelet. The scale determines how many samples per x unit are taken. For instance, with the default Mexican Hat wavelet, a scales = 2 will evaluate the [-8, 8] range sampling twice per x unit, this is with a step of 0.5. This generates a 33 points long scaled wavelet. Choosing this type of scaling is convenient because the scaled wavelet becomes a wavelet of \sigma = `scales` points. Using the default wavelet, a scales smaller than 1 will show sampling issues, while a scales larger than 64 will resample points from the original mother wavelet. If you need to use scales larger than 64, consider providing your own mother wavelet. See the examples.

According to \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1063/1.3505103")}, if your spectrum has a gaussian peak shape of variance m^2 points then the scales range should cover [1, 1.9 m]. If your spectrum has a Lorentzian peak shape of half-width-half-maximum L points then the scales range should cover [1, 2.9 L]. They also suggest using a \log_{1.18} spacing. Take these values as suggestions for your analysis.

Value

The return is the 2-D CWT coefficient matrix, with column names as the scale. Each column is the CWT coefficients at that scale. If the scales are too big for the given signal, the returned matrix may include less columns than the given scales.

Author(s)

Pan Du, Simon Lin

Examples

data(exampleMS)
scales <- seq(1, 64, 3)
wCoefs <- cwt(exampleMS[5000:11000], scales = scales, wavelet = "mexh")

## Plot the 2-D CWT coefficients as image (It may take a while!)
xTickInterval <- 1000
image(5000:11000, scales, wCoefs,
    col = terrain.colors(256), axes = FALSE,
    xlab = "m/z index", ylab = "CWT coefficient scale", main = "CWT coefficients"
)
axis(1, at = seq(5000, 11000, by = xTickInterval))
axis(2, at = c(1, seq(10, 64, by = 10)))
box()

## Provide a larger wavelet:
scales <- c(seq(1, 64, 3), seq(72, 128, 8))
psi_xval <- seq(-8, 8, length.out = 4096)
psi <- (2 / sqrt(3) * pi^(-0.25)) * (1 - psi_xval^2) * exp(-psi_xval^2 / 2)
wCoefs <- cwt(exampleMS[5000:11000], scales = scales, wavelet = rbind(psi_xval, psi))
xTickInterval <- 1000
image(5000:11000, scales, wCoefs,
    col = terrain.colors(256), axes = FALSE,
    xlab = "m/z index", ylab = "CWT coefficient scale", main = "CWT coefficients"
)
axis(1, at = seq(5000, 11000, by = xTickInterval))
axis(2, at = c(1, seq(10, 128, by = 10)))
box()

## Custom log1.18 spaced scales:
get_scales <- function(scale_min, scale_max, num_scales) {
    (seq(0, 1, length.out = num_scales)^1.18) * (scale_max - scale_min) + scale_min
}
scales <- get_scales(scale_min = 1, scale_max = 64, num_scales = 32)
print(scales)
# For detecting a gaussian peak of 10 points:
gaussian_peak_sigma <- 10 # points
scales <- get_scales(scale_min = 1, scale_max = 1.9 * gaussian_peak_sigma, num_scales = 32)
print(scales)
# For detecting a lorentzian peak of 10 points:
lorentzian_peak_gamma <- 10 # points
scales <- get_scales(scale_min = 1, scale_max = 2.9 * lorentzian_peak_gamma, num_scales = 32)
print(scales)

zeehio/MassSpecWavelet documentation built on May 6, 2023, 1:32 a.m.