Iterated Masking Empirical Mode Decomposition

Iterates a uniform-phase masked decomposition, updating each mask frequency to the amplitude-weighted mean instantaneous frequency of its mode until a relative-change test passes or the iteration limit is reached (itEMD).

Usage

decomposeIteratedMaskedEMD(signal, samplingInterval, initialFrequencies,
                           phaseCount, amplitudeRatio, tolerance = 0.1,
                           maxIterations = 15, sifting)

Arguments

signal: Numeric vector of finite real values without dim, in the units of the quantity.
samplingInterval: Sampling interval in seconds, a finite number greater than zero, without default.
initialFrequencies: Numeric vector of $K \ge 1$ initial mask frequencies in Hz, each in $(0, 1/(2\Delta t)]$, without default; its length is the number of masks.
phaseCount: Whole number $P \ge 2$ of mask phases, without default; the package tests use 4.
amplitudeRatio: Dimensionless $\varepsilon_0 > 0$, without default; the mask amplitude of each stage is $\varepsilon_0$ times the standard deviation of the residue with divisor $N$. The package tests use 1.
tolerance: Threshold of the largest relative change of the mask frequencies, a finite number greater than zero; default 0.1.
maxIterations: Whole number $\ge 1$ of frequency updates; default 15. The defaults of tolerance and maxIterations are those of iterated_mask_sift in emd 0.8.1, the code of a coauthor of itEMD.
sifting: Named list of sifting controls, as in decomposeEMD; the rule "consecutive" is rejected.

Details

Iteration $q$ decomposes the signal with the uniform-phase masks of decomposeMaskedEMD, one stage per frequency of $\mathbf{f}^{(q)}$, phases $2 \pi p / P$ and the amplitude above. Each frequency is then updated to the mean instantaneous frequency of its mode weighted by the squared instantaneous amplitude, over the samples where the frequency is defined, $$f_j^{(q+1)} = \frac{\sum_n A_j[n]^2 \, \nu_j[n]}{\sum_n A_j[n]^2},$$ with $A_j$ and $\nu_j$ from computeInstantaneousAttributes. The iteration stops with "tolerance" when every relative change $|f_j^{(q+1)} - f_j^{(q)}| / f_j^{(q)}$ is below tolerance, or with "iterationLimit" after maxIterations updates. The result is the last decomposition together with the frequencies that produced it.

Value

A list with the fields of decomposeMaskedEMD for the last decomposition (modes, residue, iterations, converged) and

stopReason: "tolerance" or "iterationLimit".
maskFrequencies: Numeric vector, the $K$ mask frequencies in Hz that produced the returned decomposition.
maskIterations: Integer, the number of frequency updates.
innerStopReason: The termination of the returned masked decomposition, with the vocabulary of decomposeMaskedEMD.

Conditions

When the masked decomposition has fewer modes than masks the update is undefined and the function stops with an error of class xemdInsufficientModes. A mode with zero squared amplitude where its frequency is defined stops with an error of class xemdZeroEnergy. An updated frequency that is not positive stops with an error when the iteration continues.

References

Fabus, M. S., Quinn, A. J., Warnaby, C. E. and Woolrich, M. W. (2021). Automatic decomposition of electrophysiological data into distinct nonsinusoidal oscillatory modes. Journal of Neurophysiology, 126(5), 1670–1684. doi:10.1152/jn.00315.2021

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 = 3000)
Signal <- sin(pi * Time) + sin(2 * pi * Time) + sin(6 * pi * Time) +
    0.5 * Time
Sifting <- list(stopRule = "doubleThreshold", thresholds = c(0.05, 0.5),
                fraction = 0.05, maxSift = 2000, boundary = "reflection",
                maxModes = 3)
# Three dyadic initial frequencies, with one sample per second.
OUT <- decomposeIteratedMaskedEMD(Signal, samplingInterval = 1,
                                  initialFrequencies = 55 / 3000 / 2 / 2^(0:2),
                                  phaseCount = 4, amplitudeRatio = 1,
                                  tolerance = 0.1, maxIterations = 15,
                                  sifting = Sifting)
OUT$maskFrequencies
#> [1] 0.009136117 0.002925373 0.001369251
OUT$maskIterations
#> [1] 3
OUT$stopReason
#> [1] "tolerance"