Empirical Wavelet Transform

Splits the spectrum of a sampled signal at the midpoints between its dominant spectral maxima, or at the boundaries of a detection rule of the EWT code of Gilles, and returns the components of an adapted filter bank, which add up to the signal.

Usage

decomposeEWT(signal, samplingInterval, bandCount, detection = NULL)

Arguments

signal

Numeric vector of at least four finite values without dim, in the units of the quantity.

samplingInterval

Sampling interval $\Delta t$ in seconds, a finite number greater than zero.

bandCount

Whole number $\ge 2$, the largest number of bands, or NULL for the automatic count of detection. No default; ewtpy uses 5, and pyewt, the package of the author of the method, requires the value.

detection

NULL, the default, for the midpoint boundaries, or a named list with these fields:

detector: "locmax", "locmaxmin" or "locmaxminf"; required.
logarithm: TRUE to work on the natural logarithm of the spectral magnitude, which must then be positive; default FALSE.
regularization: "none", "average" or "gaussian"; default "gaussian".
window: Whole number $L \ge 1$ of cells of the regularization filter; default 10.
sigma: Standard deviation $\sigma > 0$ of the gaussian filter, in cells; default 1.5.
thresholdRatio: Number $\alpha$ in $(0, 1)$, given only with bandCount = NULL. No default; Gilles (2013) tries 0.3 and 0.4.

The defaults are the initial values of the EWT code of Gilles; its initial detector, "scalespace", is not available. With bandCount = NULL the list must give detector = "locmax", regularization = "none", logarithm = FALSE and thresholdRatio.

Details

The maxima of the one-sided magnitude of the discrete Fourier transform are ranked, and the bandCount - 1 largest are retained, or all of them when fewer exist, so the number of bands adapts to the spectrum. With $b_0 = 0$ and $b_1 < \cdots < b_{K-1}$ the bins of the retained maxima, the boundaries are the midpoints $\Omega_j = \pi (b_{j-1} + b_j) / N$ rad per sample, and the last band closes at $\pi$. The record is extended by mirroring $\lceil N/2 \rceil - 1$ samples before it and $\lceil N/2 \rceil$ after it. On that extension the bank holds an approximation filter and band filters with cosine transitions of half width $\gamma \Omega_j$, $\gamma$ below the bound that keeps neighbouring transitions apart. The analysis coefficients are $W_j = \mathcal{F}^{-1}[\hat x \, \overline{H_j}]$, and each component applies the dual filter $H_j / \sum_i H_i^2$, $$c_j = \mathcal{F}^{-1}\left[\frac{\hat W_j H_j}{\sum_i H_i^2}\right],$$ restricted to the record, so the components add up to the signal. The components are defined by the filter bank, without imposing the IMF conditions of EMD. The method has no iterative stopping rule.

A detection list works on the one-sided magnitude $s$, or on its natural logarithm, and on its regularization $p$, the central $n$ values of the full convolution of $s$ with the weights $h_r$, $r = 0, \ldots, L-1$, from position $\lfloor L/2 \rfloor + 1$: $h_r = 1/L$ for "average" and $h_r \propto \exp[-u_r^2 / (2\sigma^2)]$, $u_r = r - (L-1)/2$, for "gaussian"; "none" keeps $p = s$. "locmax" places the midpoints above on the maxima of $p$. "locmaxmin" and "locmaxminf" need bandCount - 1 positive maxima of $p$ and retain the largest; on the interval from the first cell to the first retained maximum and on those between consecutive retained maxima, endpoints included, the boundary is the lowest minimum of $p$ for "locmaxmin" or of $s$ for "locmaxminf", the upper central one among equal minima, and a minimum at position $q$ gives $\Omega = 2\pi (q - 1) / N$. No boundary is completed. With bandCount = NULL, the number of bands is the number of maxima of the magnitude whose height exceeds $M_m + \alpha (M_1 - M_m)$, $M_1$ and $M_m$ the largest and smallest heights, and the midpoint boundaries follow; one band returns the signal.

Value

A list with

modes: Double matrix, samples by bands.
boundaries: Numeric vector of the $K - 1$ boundaries in Hz.

Conditions

Equal magnitudes at the last retained maximum and the first discarded one leave the selection undefined and stop with an error, as does a frequency where every filter of the bank vanishes. With a detection list, a zero magnitude under logarithm = TRUE, fewer than bandCount - 1 positive maxima of $p$ or an interval without minimum for "locmaxmin" and "locmaxminf", boundaries that do not increase and an automatic count of zero bands stop with an error.

References

Gilles, J. (2013). Empirical wavelet transform. IEEE Transactions on Signal Processing, 61(16), 3999–4010.

Gilles, J. Empirical Wavelets, MATLAB toolbox, revision 0381208d40352237876b79912013b0f68337133c, https://github.com/jegilles/Empirical-Wavelets.

The detection of the boundaries, the mirror extension and the filter bank are adapted from ewtpy, revision f9ff38bc5cd7859eea58c5b29099f8a503792b63, copyright (c) 2019 Vinícius Carvalho and Eduardo Mazoni, MIT License; the synthesis through the dual bank is adapted from pyewt, revision 9d60a340f34e92356ed4f7a6a8b7044f524455e3, copyright (c) 2024 Jerome Gilles, MIT License. See the NOTICE file of the installed package.

Examples

SamplingInterval <- 1 / 1000
Time <- seq_len(1000) * SamplingInterval
Signal <- cos(2 * pi * 2 * Time) + 0.25 * cos(2 * pi * 24 * Time) +
    0.0625 * cos(2 * pi * 288 * Time)
OUT <- decomposeEWT(Signal, samplingInterval = SamplingInterval,
                    bandCount = 3)
OUT$boundaries
#> [1]  1 13
dim(OUT$modes)
#> [1] 1000    3
max(abs(rowSums(OUT$modes) - Signal))
#> [1] 1.025846e-13
# Automatic count: only the 2 Hz maximum exceeds the threshold of 0.3.
dim(decomposeEWT(Signal, samplingInterval = SamplingInterval,
                 bandCount = NULL,
                 detection = list(detector = "locmax", regularization = "none",
                                  thresholdRatio = 0.3))$modes)
#> [1] 1000    1
# Lowest minima between the maxima of the regularized magnitude.
Time <- seq(0, 9, length.out = 3000)
Signal <- sin(pi * Time) + sin(2 * pi * Time) + sin(6 * pi * Time) +
    0.5 * Time
decomposeEWT(Signal, samplingInterval = 9 / 2999, bandCount = 3,
             detection = list(detector = "locmaxmin"))$boundaries
#> [1] 0.7775185 2.4436296