Ensemble Empirical Mode Decomposition

Decomposes noise-perturbed copies of a signal by scalar sifting and averages their completed modes by extraction order.

Usage

decomposeEEMD(signal, noise, amplitudeRatio = 0.2, divisor, sifting)

Arguments

signal: Numeric vector of at least two finite values without dim, in the units of the quantity.
noise: Finite numeric matrix with one row per sample of signal and one column per realization, $M \ge 1$. Column $m$ is a realization $z^{(m)}$ of unit variance, used as supplied: the function neither generates nor normalizes it, and it draws no random numbers.
amplitudeRatio: Dimensionless ratio $\varepsilon \ge 0$ of the noise standard deviation to the standard deviation of signal (divisor $N - 1$). Default 0.2, the value common to libeemd and emd 0.8.1.
divisor: Categorical, without default: "realizations" divides the sum of the modes of order $j$ by $M$; "available" divides it by the number $M_j$ of realizations that have that mode.
sifting: Named list of sifting controls, as in decomposeEMD; the rule "consecutive" is rejected.

Details

Realization $m$ is $x + w^{(m)}$ with $w^{(m)} = \varepsilon \, \mathrm{sd}(x) \, z^{(m)}$; it is decomposed by scalar sifting with sifting, in the order of the columns of noise. Each realization stops before adding an extraction whose sifting lost its envelopes and keeps the input of that attempt in its residue. Accepted adaptive extractions, completed fixed sifting and extractions that reached the iteration limit are retained, the last ones flagged FALSE. Absent modes are extended with zeros, and mode $j$ is $$D_j = \frac{1}{M} \sum_{m=1}^{M} c_j^{(m)}$$ with "realizations", or the same sum divided by $M_j$ with "available". The residue is the mean of the $M$ final residues and meanNoise the mean of the added perturbations, so that with "realizations" the modes, the residue and signal satisfy rowSums(modes) + residue = signal + meanNoise up to rounding. A realization without modes contributes zeros to the modes and its input to the residue. Matching by extraction order does not establish a physical correspondence between the modes of different realizations.

The method has no outer termination, so the result has no stopReason. The number of realizations has no default; the audited libraries use 100 (PyEMD), 250 (libeemd) and 4 (emd 0.8.1). hht divides by $M$ and PyEMD by $M_j$.

Value

A list with

modes: Double matrix, samples by modes, with as many columns as the largest number of modes among the realizations.
residue: Numeric vector, the mean final residue.
meanNoise: Numeric vector, the mean added perturbation, in the units of signal.
iterations: Integer vector by mode, the sum of the subtractions of the realizations that have the mode.
converged: Logical vector by mode: TRUE only when every realization that has the mode accepted it; NA for completed fixed sifting. Failed attempts are reported in termination.
termination: Data frame with one row per realization: realization, its column index in noise; mode, the next extraction index after the retained modes; iterations, subtractions of the interrupted attempt, zero if it did not start; and stopReason, "terminalResidue", "modeLimit" or "insufficientExtrema". These rows describe individual decompositions, not the averaged mode indices.

References

Wu, Z. and Huang, N. E. (2009). Ensemble empirical mode decomposition: A noise-assisted data analysis method. Advances in Adaptive Data Analysis, 1(1), 1–41. doi:10.1142/S1793536909000047

The selection of the reflected extrema and the spline interpolation of the envelope nodes are adapted from emd.m (G. Rilling, March 2007), distributed with CEEMDAN revision af9a0f83, copyright (c) 2020 Marcelo A. Colominas, MIT License; see the NOTICE file of the installed package.

Examples

Time <- seq(0, 9, length.out = 1000)
Signal <- sin(pi * Time) + sin(2 * pi * Time) + sin(6 * pi * Time) +
    0.5 * Time
set.seed(628)
Noise <- t(matrix(stats::rnorm(10L * length(Signal)), nrow = 10L,
                  byrow = TRUE))
Sifting <- list(stopRule = "doubleThreshold", thresholds = c(0.05, 0.5),
                fraction = 0.05, maxSift = 2000, boundary = "reflection",
                maxModes = Inf)
# Two of the ten realizations keep the example short.
OUT <- decomposeEEMD(Signal, noise = Noise[, 1:2], amplitudeRatio = 0.2,
                     divisor = "realizations", sifting = Sifting)
dim(OUT$modes)
#> [1] 1000    7
OUT$iterations
#> [1] 15 25 42  9 14 13  2
max(abs(rowSums(OUT$modes) + OUT$residue - Signal - OUT$meanNoise))
#> [1] 1.720846e-15