Package 'xemd'

Fractional Gaussian Noise Test of Intrinsic Mode Functions

Description

Compares the energy of each intrinsic mode function, from the second on, with the threshold that fractional Gaussian noise of a tabulated Hurst exponent gives at a tabulated level, with the first function as the noise scale.

Usage

assessFractionalNoise(modes, hurst, probability)

Arguments

modes	Finite numeric matrix of 2048 samples by at least two intrinsic mode functions, in their order of extraction and without the residue.
hurst	Hurst exponent $H$: 0.2, 0.5 or 0.8. No default.
probability	Level $p$ of the threshold: 0.95 or 0.99. No default.

Details

The energy of function $k$ is $W[k] = \sum_n c_k[n]^2$. From $k = 2$ on, the reference energy is

$W_H[k] = \frac{W[1]}{\beta_H} \rho_H^{-2(1-H)k}, \quad \rho_H = 2.01 + 0.2 (H - 0.5) + 0.12 (H - 0.5)^2,$

and the threshold is $T_{H,p}[k] = W_H[k] \, 2^{2^{a_{H,p} k + b_{H,p}}}$, with $\beta_H$, $a_{H,p}$ and $b_{H,p}$ of table 3.1 of Flandrin et al. (2005). A function lies outside when $W[k] > T_{H,p}[k]$. The coefficients come from simulations of 10000 realizations of 2048 samples and are not interpolated to other exponents nor recalibrated for other lengths. The first function is taken to carry mainly noise; the function does not verify it.

The result is conditional on that noise model and calibration. Exceeding the threshold does not prove that a component is free of noise or identify its physical source.

Value

A list of vectors with one element per function:

energy	The energies $W[k]$.
reference	The reference energies $W_H[k]$, NA for the first function.
upper	The thresholds $T_{H,p}[k]$, NA for the first function.
outside	Logical, energy > upper, NA for the first function.

Conditions

A first function of zero energy stops the function with an error of class xemdZeroEnergy.

References

Flandrin, P., Gonçalvès, P. and Rilling, G. (2004). Detrending and denoising with empirical mode decompositions. EUSIPCO 2004, 1581–1584.

Flandrin, P., Gonçalvès, P. and Rilling, G. (2005). EMD equivalent filter banks, from interpretation to applications. In N. E. Huang and S. S. P. Shen (eds.), Hilbert-Huang Transform and Its Applications, World Scientific, 57–74.

See Also

assessWhiteNoise, decomposeEMD

Examples

Time <- seq_len(2048) / 1000
Tones <- cbind(cos(2 * pi * 2 * Time), 0.25 * cos(2 * pi * 24 * Time),
               0.0625 * cos(2 * pi * 288 * Time))
assessFractionalNoise(Tones, hurst = 0.5, probability = 0.95)
# The third tone as the noise scale.
assessFractionalNoise(Tones[, c(3, 1, 2)], hurst = 0.5,
                      probability = 0.95)$outside

---

Intrinsic Mode Function Conditions

Description

Evaluates numerical diagnostics of the intrinsic mode function conditions for each column of a matrix: the balance of extrema and zero crossings, and a double threshold on the envelope mean relative to its amplitude.

Usage

assessIMF(modes, thresholds = c(0.05, 0.5), fraction = 0.05,
          boundary = "reflection")

Arguments

modes	Finite numeric matrix, samples by columns. Attributes other than dim do not reach the outputs.
thresholds	Two dimensionless numbers $0 < \theta_1 < \theta_2$; default c(0.05, 0.5).
fraction	The admitted fraction $\alpha$ in $[0, 1)$; default 0.05.
boundary	"reflection" (default) or "extrapolation", the extension of the envelopes at the borders, as in decomposeEMD. The defaults are those of emd.m, shared with the sifting controls.

Details

The ideal envelope condition is a zero local mean. This diagnostic uses tolerances and a specified boundary rule; passing it does not establish recovery of a physical component. For methods defined by other conditions, such as a bandwidth objective, it is an additional EMD-based diagnostic.

count tests $|N_e - N_z| \le 1$, with the extrema $N_e$ and zero crossings $N_z$ counted as in decomposeEMD. deviation applies to

$\sigma = \frac{|m_e|}{a_e},$

$m_e$ and $a_e$ the half sum and the half difference of the envelopes, the double threshold that accepts a mode under stopRule = "doubleThreshold": the fraction of samples with $\sigma \ge \theta_1$ is at most $\alpha$ and no sample has $\sigma \ge \theta_2$.

Value

A list of two logical vectors with one element per column:

count	The counting condition.
deviation	The double threshold, NA where the boundary rule leaves the column without envelopes.

References

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.

See Also

decomposeEMD, assessWhiteNoise

Examples

Time <- seq(0, 9, length.out = 3000)
Signal <- sin(pi * Time) + sin(2 * pi * Time) + sin(6 * pi * Time) +
    0.5 * Time
OUT <- decomposeEMD(Signal, sifting = list(stopRule = "doubleThreshold"))
# The signal, its modes and its residue as columns.
assessIMF(cbind(Signal, OUT$modes, OUT$residue))

---

White-Noise Significance Test of Modes

Description

Compares the energy of each mode with reference bounds from a white-noise model at the mode's mean period (test of Wu and Huang).

Usage

assessWhiteNoise(modes, probabilities, bound, rescalingMode)

Arguments

modes	Finite numeric matrix, samples by modes. The caller supplies the components of a normalized signal; see Details for the assumptions of the reference distribution.
probabilities	Two numbers $0 < p_- < p_+ < 1$, the probabilities of the lower and upper bounds. No default; PyEMD uses $\alpha = 0.95$, that is $(1 \mp \alpha)/2$.
bound	"gaussian", the line of the logarithmic quantile, or "chiSquare", the quantile of the energy. No default; PyEMD has only the gaussian bound.
rescalingMode	NULL for the test without rescaling, or the index $r$ of the column whose energy is placed on its upper bound. No default; PyEMD rescales on the first component.

Details

The energy of mode $k$ is $E_k = \sum_n c_k[n]^2 / N$ and its mean period $T_k$ is $N$ divided by its number of maxima, or $N$ without maxima. The bound of probability $p$ is, with "gaussian" and $z_p$ the standard normal quantile,

$\ln E_p = -\ln T + z_p \sqrt{2/N}\, T^{1/2},$

and with "chiSquare" $E_p = q(p; r) / N$, the chi-square quantile with $r = N / T$ degrees of freedom. With rescalingMode = r, the energy of component $r$ is placed on its upper bound and every other energy is multiplied by $E_{p_+}(T_r) / E_r$. A mode lies outside when its energy is below the lower bound or above the upper one.

The reference distribution is conditional on the noise model, normalization and EMD procedure. Interpreting the rescaled test assumes that the component selected by rescalingMode is dominated by noise; its mode index does not establish that assumption. Exceeding a bound does not identify a physical source or prove that a component is free of noise. A white-noise interpretation of an acceleration record does not transfer unchanged to modes of its integrated displacement.

Value

A list of vectors with one element per mode:

energy	The energies, rescaled when rescalingMode is given.
period	The mean periods, in samples.
lower, upper	The bounds at those periods.
outside	Logical, energy < lower | energy > upper.

Conditions

A column of zero energy stops the function with an error of class xemdZeroEnergy.

See Also

assessIMF, decomposeEMD

Examples

Time <- seq(0, 9, length.out = 3000)
Signal <- sin(pi * Time) + sin(2 * pi * Time) + sin(6 * pi * Time) +
    0.5 * Time
# The four summands of the signal, divided by its population deviation.
Summands <- cbind(sin(pi * Time), sin(2 * pi * Time), sin(6 * pi * Time),
                  0.5 * Time) / sqrt(mean((Signal - mean(Signal))^2))
assessWhiteNoise(Summands, probabilities = c(0.025, 0.975),
                 bound = "gaussian", rescalingMode = NULL)$outside
assessWhiteNoise(Summands, probabilities = c(0.025, 0.975),
                 bound = "chiSquare", rescalingMode = 1)$outside

---

Hilbert Spectrum by Frequency Cells

Description

Accumulates the instantaneous amplitude, or its square, of every mode in the frequency cell of its instantaneous frequency, sample by sample, and integrates the result over time.

Usage

buildHilbertSpectrum(frequency, amplitude, edges, exponent, samplingInterval)

Arguments

frequency	Numeric matrix, samples by modes, of instantaneous frequencies in Hz, finite or NA.
amplitude	Finite non-negative matrix with the dimensions of frequency, the instantaneous amplitudes in the units of the signal.
edges	At least two finite, strictly increasing cell edges in Hz; cell $b$ is $[e_b, e_{b+1})$. Negative values are admitted. No default.
exponent	The power $p$: 1 for the amplitude spectrum, 2 for the spectrum weighted by squared amplitude. No default.
samplingInterval	Sampling interval $\Delta t$ in seconds, a finite number greater than zero.

Details

The spectrum holds the mass of each cell,

$H[n, b] = \sum_{k:\, f_k[n] \in [e_b, e_{b+1})} a_k[n]^p,$

in units of the signal to the power $p$. A frequency outside the edges or NA adds nothing. The marginal spectrum of a cell is the time integral $\Delta t \sum_n H[n, b]$. No scaling or division by the cell width is applied: the values are masses, not densities. The inputs are usually the frequency and amplitude of computeInstantaneousAttributes.

Value

A list with

spectrum	Double matrix, samples by cells, in units of the signal to the power $p$.
marginal	Numeric vector by cell, in units of the signal to the power $p$ times seconds.

See Also

computeInstantaneousAttributes, decomposeEMD

Examples

Time <- seq(0, 9, length.out = 3000)
Signal <- sin(pi * Time) + sin(2 * pi * Time) + sin(6 * pi * Time) +
    0.5 * Time
Modes <- decomposeEMD(Signal,
                      sifting = list(stopRule = "doubleThreshold"))$modes
Attributes <- computeInstantaneousAttributes(Modes,
                                             samplingInterval = 9 / 2999)
OUT <- buildHilbertSpectrum(Attributes$frequency, Attributes$amplitude,
                            edges = c(1, 2, 3, 4, 5), exponent = 2,
                            samplingInterval = 9 / 2999)
dim(OUT$spectrum)
OUT$marginal

---

Instantaneous Amplitude, Phase and Frequency

Description

Computes, for each column of a matrix of modes, the analytic signal and its instantaneous amplitude, continuous phase and frequency.

Usage

computeInstantaneousAttributes(modes, samplingInterval)

Arguments

modes	Finite double matrix, samples by modes, with at least three rows; it may have zero columns. Attributes other than dim do not reach the outputs.
samplingInterval	Sampling interval $\Delta t$ in seconds, a finite number greater than zero.

Details

The analytic signal $z = c + i \mathcal{H} c$ of each column is computed with the discrete Fourier transform, multiplying each bin by $1 + \mathrm{sgn}(\omega)$: positive frequencies are doubled, negative ones removed, and the zero bin and, for an even number of samples, the Nyquist bin keep the factor 1, so $\mathrm{Re}\, z = c$. The transform implies the periodic extension of the record. The amplitude is $a = |z|$. Within each run of samples with $|z| > 0$, the phase starts at the principal value of $\arg z$ and adds the increments $\arg(z_{n+1} \bar z_n)$ taken in $(-\pi, \pi]$; a zero of $z$ leaves its sample undefined (NA) and the next run starts again at its principal value. The frequency is

$f_n = \frac{\phi_{n+1} - \phi_{n-1}}{4 \pi \Delta t}$

in Hz, the central difference of $\omega = \phi'$; it is NA where one of the three phases is, and each border repeats its neighbour. These local undefined values raise no condition.

Value

A list of three matrices, samples by modes:

amplitude	Instantaneous amplitude, in the units of the modes.
phase	Continuous phase in radians, NA at zeros of the analytic signal.
frequency	Instantaneous frequency in Hz, NA where the phase is undefined.

See Also

buildHilbertSpectrum, decomposeEMD, estimateMaskFrequency

Examples

SamplingInterval <- 1 / 1000
Time <- seq_len(1000) * SamplingInterval
# Three tones of 2, 24 and 288 Hz as columns.
Tones <- cbind(cos(2 * pi * 2 * Time), 0.25 * cos(2 * pi * 24 * Time),
               0.0625 * cos(2 * pi * 288 * Time))
OUT <- computeInstantaneousAttributes(Tones,
                                      samplingInterval = SamplingInterval)
apply(OUT$amplitude, 2, range)
apply(OUT$frequency, 2, range)

---

Bivariate Empirical Mode Decomposition

Description

Decomposes a two-channel signal, read as the complex signal $z = x_1 + i x_2$, by successive extraction of bivariate modes with envelopes built on projections onto equally spaced directions of the plane.

Usage

decomposeBEMD(signal, directionCount = 64, envelope, sifting)

Arguments

signal	Numeric matrix of finite real values with two columns, the real and the imaginary parts, samples in rows.
directionCount	Whole number $D \ge 1$ of directions $\varphi_d = 2 \pi (d - 1) / D$; default 64, the default of Rlibeemd bemd and of emdr memd.
envelope	Mean of the sifting, without default: "trajectory" averages the envelopes of the trajectory itself at the maxima of each projection; "projection" sums the envelopes of the projections, oriented on their directions, with the factor $2/D$.
sifting	Named list of sifting controls, as in decomposeEMD, with stopRule = "fixed" only. The trajectory envelope admits only boundary = "reflection"; the projection envelope also admits "extrapolation".

Details

The projection on direction $\varphi_d$ is $p_d = x_1 \cos\varphi_d + x_2 \sin\varphi_d$. With "trajectory", $e_d$ interpolates the two-dimensional samples at the maxima of $p_d$ and the mean is

$m = \frac{1}{D} \sum_{d=1}^{D} e_d;$

with "projection", the scalar upper envelope of $p_d$ is oriented on $(\cos\varphi_d, \sin\varphi_d)$ and the sum is scaled by $2/D$. Each extraction subtracts the mean fixedSift times. Before extraction $k$ the decomposition ends with "terminalResidue" when every projection of $r_{k-1}$ has at most one extremum, with "modeLimit" when $k - 1$ equals sifting$maxModes, and with "insufficientExtrema" when a direction has no envelope. One extraction serves both channels, so they share the number of modes. If a direction loses its envelope before fixed sifting completes, the decomposition stops before adding that vector candidate and preserves the residue from before the attempt. A candidate that reaches the iteration limit of an adaptive rule is kept and flagged FALSE. Directions are not dropped from the mean.

Value

A list with

modes	Double array, samples by modes by channels.
residue	Numeric matrix, samples by channels.
iterations	Integer vector, the subtractions of each mode.
converged	Logical vector by mode: NA after the fixed number of subtractions. Interrupted attempts are reported in termination.
stopReason	"terminalResidue", "modeLimit" or "insufficientExtrema".
termination	Data frame with one row: mode, the next extraction index after the retained modes; direction, the one-based index of the first direction without an envelope, or NA when none failed; iterations, subtractions of the interrupted attempt, zero if it did not start; and stopReason.

References

Rilling, G., Flandrin, P., Gonçalves, P. and Lilly, J. M. (2007). Bivariate Empirical Mode Decomposition. IEEE Signal Processing Letters, 14, 936–939.

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.

See Also

decomposeMEMD, decomposeEMD

Examples

Time <- 2 * pi * (0:511) / 512
z <- cos(0.3 * Time) * exp(2i * Time) +
    0.3 * abs(sin(2.3 * Time)) * exp(17i * Time)
Signal <- cbind(Re(z), Im(z))
OUT <- decomposeBEMD(Signal, directionCount = 64, envelope = "projection",
                     sifting = list(stopRule = "fixed", fixedSift = 10,
                                    boundary = "extrapolation",
                                    maxModes = 3))
dim(OUT$modes)
OUT$stopReason
max(abs(apply(OUT$modes, c(1, 3), sum) + OUT$residue - Signal))

---

Complementary Ensemble Empirical Mode Decomposition

Description

Decomposes pairs of copies of a signal perturbed with the same noise realization of opposite signs and averages the modes of all the copies by extraction order.

Usage

decomposeCEEMD(signal, noise, amplitudeRatio = 0.2, 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 draws no random numbers.
amplitudeRatio	Dimensionless ratio $\varepsilon \ge 0$; the perturbation is $a z^{(m)}$ with $a = \varepsilon \, \mathrm{sd}(x)$ (divisor $N - 1$). Default 0.2, as in decomposeEEMD.
sifting	Named list of sifting controls, as in decomposeEMD; the rule "consecutive" is rejected.

Details

For each column, the inputs $x + a z^{(m)}$ and $x - a z^{(m)}$ are decomposed by scalar sifting, the positive input first. Each stops before adding an extraction whose sifting lost its envelopes and keeps its input in the residue. Accepted adaptive extractions, completed fixed sifting and extractions that reached the iteration limit are retained. The $2M$ decompositions are extended with zero modes to the largest number of modes and averaged by extraction order,

$d_j = \frac{1}{2M} \sum_{m=1}^{M} \left(c_j^{(m,+)} + c_j^{(m,-)}\right),$

and the residue is the mean of the $2M$ residues, so the perturbations cancel in the reconstruction: rowSums(modes) + residue equals signal up to rounding. The method has no outer termination and returns neither meanNoise nor stopReason.

Value

A list with

modes	Double matrix, samples by modes.
residue	Numeric vector, the mean residue of the $2M$ decompositions.
iterations	Integer vector by mode, the sum of the subtractions of the decompositions that have the mode.
converged	Logical vector by mode: TRUE only when every decomposition that has the mode accepted it; NA for completed fixed sifting. Failed attempts are reported in termination.
termination	Data frame with one row per signed realization: realization, its column index in noise; sign, 1 or -1; and mode, iterations and stopReason as in decomposeEEMD. A stopped realization stays in its pair and in the divisor $2M$; its remaining content stays in the mean residue.

References

Yeh, J.-R., Shieh, J.-S. and Huang, N. E. (2010). Complementary ensemble empirical mode decomposition: A novel noise enhanced data analysis method. Advances in Adaptive Data Analysis, 2(2), 135–156. doi:10.1142/S1793536910000422

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.

See Also

decomposeEEMD, decomposeCEEMDAN, decomposeICEEMDAN

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 <- decomposeCEEMD(Signal, noise = Noise[, 1:2], amplitudeRatio = 0.2,
                      sifting = Sifting)
dim(OUT$modes)
max(abs(rowSums(OUT$modes) + OUT$residue - Signal))

---

Complete Ensemble Empirical Mode Decomposition with Adaptive Noise

Description

Extracts each mode as the average first mode of the common residue perturbed by the modes of supplied noise realizations, with a constant ratio of the noise to the residue at each stage.

Usage

decomposeCEEMDAN(signal, noise, amplitudeRatio = 0.2, 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 $w^{(m)}$ of unit variance, used as supplied. The function draws no random numbers.
amplitudeRatio	Dimensionless ratio $\varepsilon \ge 0$ of the perturbation to the standard deviation of the residue at each stage (divisor $N - 1$). Default 0.2, as in decomposeEEMD.
sifting	Named list of sifting controls, as in decomposeEMD; the rule "consecutive" is rejected. The same list decomposes the noise realizations.

Details

Each column of noise is decomposed once with sifting, $\eta_k^{(m)} = E_k(w^{(m)})$, and $\eta_0^{(m)} = w^{(m)}$. From $r_0 = x$, before stage $k$ the common residue $r_{k-1}$ ends the decomposition, in this order, with "terminalResidue" when it has at most one extremum, with "modeLimit" when $k - 1$ equals sifting$maxModes and with "insufficientExtrema" when its envelopes cannot be built. Stage $k$ perturbs $r_{k-1}$ with $p^{(m)} = \beta_{k-1}^{(m)} \eta_{k-1}^{(m)}$, $\beta_{k-1}^{(m)} = \varepsilon \, \mathrm{sd}(r_{k-1}) / \mathrm{sd}(\eta_{k-1}^{(m)})$, and computes

$d_k = \frac{1}{M} \sum_{m=1}^{M} E_1\left(r_{k-1} + p^{(m)}\right), \qquad r_k = r_{k-1} - d_k.$

A realization without $\eta_{k-1}$ contributes a null perturbation, and a perturbed residue without a first mode contributes zero; both stay in the divisor $M$.

Value

A list with

modes	Double matrix, samples by modes.
residue	Numeric vector, the final common residue.
iterations	Integer vector by mode, the sum of the subtractions of the $M$ extractions of the stage.
converged	Logical vector by mode: TRUE when the $M$ extractions accepted, FALSE when one ended without accepting or could not begin, NA when all completed fixed sifting.
stopReason	"terminalResidue", "modeLimit" or "insufficientExtrema".

Conditions

A noise mode of zero standard deviation leaves its scale undefined and stops the function with an error of class xemdZeroEnergy.

References

Torres, M. E., Colominas, M. A., Schlotthauer, G. and Flandrin, P. (2011). A complete ensemble empirical mode decomposition with adaptive noise. In 2011 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 4144–4147. doi:10.1109/ICASSP.2011.5947265

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.

See Also

decomposeICEEMDAN, decomposeEEMD, decomposeCEEMD

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))
# Two of the ten realizations and three modes keep the example short.
Sifting <- list(stopRule = "doubleThreshold", thresholds = c(0.05, 0.5),
                fraction = 0.05, maxSift = 2000, boundary = "reflection",
                maxModes = 3)
OUT <- decomposeCEEMDAN(Signal, noise = Noise[, 1:2], amplitudeRatio = 0.2,
                        sifting = Sifting)
dim(OUT$modes)
OUT$stopReason
max(abs(rowSums(OUT$modes) + OUT$residue - Signal))

---

Ensemble Empirical Mode Decomposition

Description

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.

See Also

decomposeCEEMD, decomposeCEEMDAN, decomposeICEEMDAN, resiftModes, measureDispersion

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)
OUT$iterations
max(abs(rowSums(OUT$modes) + OUT$residue - Signal - OUT$meanNoise))

---

Empirical Mode Decomposition

Description

Decomposes a signal by successive extraction: each component is obtained by sifting the current residue and is subtracted from that residue. Its numerical acceptance depends on the selected stopping rule.

Usage

decomposeEMD(signal, sifting)

Arguments

signal	Numeric vector of finite real values without dim, in the units of the quantity. Positions are sample indices, so no time vector or sampling interval is needed. Attributes such as names, tsp or class do not reach the outputs.
sifting	Named list of sifting controls: stopRule and the fields of that rule, described in section ‘Sifting controls’. Fields with a default may be omitted; every other field is required, and a field that does not belong to the rule is an error.

Details

With $r_0 = x$, extraction $k$ computes the mode $c_k = E_1(r_{k-1})$ by sifting and updates the residue

$r_k = r_{k-1} - c_k,$

so that $x = \sum_{k=1}^{K} c_k + r_K$ up to rounding. Before extraction $k$ the decomposition ends, in this order, with "terminalResidue" when $r_{k-1}$ has at most one extremum, with "modeLimit" when $k - 1$ equals sifting$maxModes, and with "insufficientExtrema" when the boundary rule cannot build the envelopes of $r_{k-1}$.

Sifting starts from $h_0 = r_{k-1}$ and repeats $h_{j+1} = h_j - m_j$, where $m_j$ is the mean of the upper and lower envelopes of $h_j$. The envelopes interpolate the maxima and the minima, extended at both borders by sifting$boundary: a cubic spline with not-a-knot end conditions through four or more nodes, the parabola through three and the line through two. A run of equal interior samples is one extremum at the exact centre of the run. Each state $h_j$ is examined once, in this order: the stopping rule accepts it (converged TRUE); $j$ equals fixedSift under "fixed" (converged NA); $j$ equals maxSift under an adaptive rule, or the envelopes of $h_j$ cannot be built (converged FALSE); otherwise one more subtraction follows. Non-convergence is reported in converged, without a warning.

Value

A list with

modes	Double matrix, samples by modes, in the units of signal; it has zero columns when no mode is extracted. No names.
residue	Numeric vector $r_K$ of the input length; equal to signal when no mode is extracted.
iterations	Integer vector with the number of subtractions of each mode.
converged	Logical vector by mode: TRUE when the stopping rule accepted the mode, FALSE when maxSift was reached or the envelopes could not be built, NA under "fixed".
stopReason	"terminalResidue", "modeLimit" or "insufficientExtrema".

Sifting controls

The same list configures the extractor of decomposeEEMD, decomposeCEEMD, decomposeCEEMDAN, decomposeICEEMDAN, resiftModes, decomposeMaskedEMD, decomposeIteratedMaskedEMD, estimateMaskFrequency and, with the restrictions stated in their help pages, decomposeBEMD, decomposeMEMD and decomposeNAMEMD. Counts are accepted as doubles with a whole value.

stopRule

Categorical, without default: "SD", "Snumber", "doubleThreshold" or "fixed". The rule "consecutive" belongs to the vector sifting of decomposeMEMD and is rejected here.

tolerance

Under "SD": dimensionless threshold, greater than zero, without default. From $j = 1$ the state is accepted when $\sum_n (h_{j-1}[n] - h_j[n])^2 / h_{j-1}[n]^2$ is finite and below it; the sum runs over all samples, without the factor $1/N$, and a null denominator is not regularized. The 1998 formulation proposes values from 0.2 to 0.3.

stableSift

Under "Snumber": the count $S \ge 1$, without default. From $j = 1$ the state is accepted when the counts of extrema $E$ and zero crossings $Z$ have stayed equal over $S$ consecutive states and $|E - Z| \le 1$. The 2003 formulation reports typical values from three to five.

thresholds

Under "doubleThreshold": two dimensionless numbers $0 < \theta_1 < \theta_2$; default c(0.05, 0.5).

fraction

Under "doubleThreshold": the admitted fraction $\alpha$ in $[0, 1)$; default 0.05. With $\sigma = |m / a|$, $m$ the half sum and $a$ the half difference of the envelopes, the state, the signal included, is accepted when $|E - Z| \le 1$, the fraction of samples with $\sigma \ge \theta_1$ is at most $\alpha$ and no sample has $\sigma \ge \theta_2$. Where $m = a = 0$ the ratio is undefined and counts as no exceedance; where only $a = 0$ it is infinite.

fixedSift

Under "fixed": the count $J \ge 1$ of subtractions of every extraction, without default.

maxSift

Under "SD", "Snumber" and "doubleThreshold": the count $M \ge 1$ of subtractions after which an extraction ends without acceptance; under "Snumber" $M \ge S$. Default 2000 under "doubleThreshold", none under the other rules. The rule "fixed" has no maxSift.

maxModes

Count $\ge 1$ or Inf, the limit of extractions; default Inf, which extracts until the residue ends the decomposition.

boundary

Categorical: "reflection" or "extrapolation"; default "reflection" under "doubleThreshold" and "fixed", none under "SD" and "Snumber". "reflection" mirrors two extrema of each family per side about the first extremum or the border sample and needs three extrema; "extrapolation" anchors the upper envelope at each border at the larger of the border sample and the line through the two nearest maxima, and the lower envelope at the smaller of the border sample and the line through the two nearest minima.

References

Huang, N. E., Shen, Z., Long, S. R., Wu, M. C., Shih, H. H., Zheng, Q., Yen, N.-C., Tung, C. C. and Liu, H. H. (1998). The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proceedings of the Royal Society of London. Series A, 454, 903–995. doi:10.1098/rspa.1998.0193

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 Marcelo A. Colominas, MIT License; see the NOTICE file of the installed package.

See Also

assessIMF, decomposeEEMD, decomposeTSL

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 = Inf)
OUT <- decomposeEMD(Signal, sifting = Sifting)
dim(OUT$modes)
OUT$iterations
OUT$stopReason
max(abs(rowSums(OUT$modes) + OUT$residue - Signal))

# The rule alone: the other fields take their defaults.
identical(decomposeEMD(Signal, sifting = list(stopRule = "doubleThreshold")),
          OUT)

# Ten fixed subtractions per mode, at most three modes.
Fixed <- decomposeEMD(Signal, sifting = list(stopRule = "fixed",
                                             fixedSift = 10, maxModes = 3))
Fixed$converged
Fixed$stopReason

---

Empirical Wavelet Transform

Description

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.

See Also

transformSST, decomposeVMD

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
dim(OUT$modes)
max(abs(rowSums(OUT$modes) - Signal))
# 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)
# 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

---

Improved Complete Ensemble Empirical Mode Decomposition with Adaptive Noise

Description

Computes each residue as the average local mean of the previous residue perturbed by a mode of supplied noise realizations, and each mode as the difference between consecutive residues.

Usage

decomposeICEEMDAN(signal, noise, amplitudeRatio = 0.2, normalizeNoise,
                  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 $w^{(m)}$ of unit variance, used as supplied. The function draws no random numbers.
amplitudeRatio	Dimensionless ratio $\varepsilon_0 \ge 0$. Default 0.2, as in decomposeEEMD.
normalizeNoise	TRUE or FALSE, without default. From the second stage on, TRUE divides each noise mode by its standard deviation before scaling; FALSE keeps its amplitude. normalizeNoise = TRUE corresponds to SNRFlag = 2 of ceemdan.m.
sifting	Named list of sifting controls, as in decomposeEMD; the rule "consecutive" is rejected. The same list decomposes the noise realizations.

Details

Each column of noise is decomposed once with sifting, $\eta_k^{(m)} = E_k(w^{(m)})$. From $r_0 = x$, before stage $k$ the common residue $r_{k-1}$ ends the decomposition, in this order, with "terminalResidue" when it has at most one extremum, with "modeLimit" when $k - 1$ equals sifting$maxModes and with "insufficientExtrema" when its envelopes cannot be built. Stage $k$ perturbs $r_{k-1}$ with $p^{(m)} = \beta_{k-1}^{(m)} \eta_k^{(m)}$, where $\beta_0^{(m)} = \varepsilon_0 \, \mathrm{sd}(x) / \mathrm{sd}(\eta_1^{(m)})$ and, for $k \ge 2$, $\beta_{k-1} = \varepsilon_0 \, \mathrm{sd}(r_{k-1})$, divided by $\mathrm{sd}(\eta_k^{(m)})$ when normalizeNoise is TRUE; standard deviations use the divisor $N - 1$. With $A$ the average first mode of $r_{k-1} + p^{(m)}$ and $\bar p$ the mean perturbation, the new residue is the mean local mean and the mode the difference of residues:

$r_k = r_{k-1} + \bar p - A, \qquad d_k = r_{k-1} - r_k.$

A realization without $\eta_k$ contributes a null perturbation and stays in the divisor $M$.

Value

A list with

modes	Double matrix, samples by modes.
residue	Numeric vector, the final common residue.
iterations	Integer vector by mode, the sum of the subtractions of the $M$ extractions of the stage.
converged	Logical vector by mode: TRUE when the $M$ extractions accepted, FALSE when one ended without accepting or could not begin, NA when all completed fixed sifting.
stopReason	"terminalResidue", "modeLimit" or "insufficientExtrema".

Conditions

A noise mode of zero standard deviation leaves its scale undefined and stops the function with an error of class xemdZeroEnergy.

References

Colominas, M. A., Schlotthauer, G. and Torres, M. E. (2014). Improved complete ensemble EMD: A suitable tool for biomedical signal processing. Biomedical Signal Processing and Control, 14, 19–29. doi:10.1016/j.bspc.2014.06.009

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.

See Also

decomposeCEEMDAN, decomposeEEMD, decomposeCEEMD

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 <- decomposeICEEMDAN(Signal, noise = Noise[, 1:2], amplitudeRatio = 0.2,
                         normalizeNoise = TRUE, sifting = Sifting)
dim(OUT$modes)
OUT$stopReason
max(abs(rowSums(OUT$modes) + OUT$residue - Signal))

---

Iterative Filtering

Description

Decomposes a signal by successive extraction of components, each obtained by iterating the subtraction of a moving average whose filter is a supplied prototype resampled to a support set by the extrema of the residue.

Usage

decomposeIF(signal, filter, boundary, tolerance, maxIterations, stretch,
            maxModes = Inf)

Arguments

signal	Numeric vector of finite values without dim, in the units of the quantity. Positions are sample indices.
filter	Prototype filter, a numeric vector of finite weights of odd length, at least 3. No default. The package tests use the 13 weights that fifpy obtains from its prototype for a half width of 6.
boundary	"zero" extends the record with zeros; "periodic" extends it periodically and evaluates the iteration spectrally (fast iterative filtering). No default.
tolerance	Bound of the relative change between consecutive iterates, greater than zero and less than one. No default; the package tests use $\sqrt{0.001}$.
maxIterations	Whole number $\ge 1$ of iterations of one extraction. No default; the package tests use 200.
stretch	Number greater than zero that scales the support of the filter. No default; the package tests use 1.6.
maxModes	Whole number $\ge 1$ or Inf, the limit of extractions; default Inf.

Details

Before extraction $k$ the decomposition ends with "terminalResidue" when $r_{k-1}$ has at most one extremum and with "modeLimit" when $k - 1$ equals maxModes. Otherwise the half width is $l_k = 2\,\mathrm{round}(\mathrm{stretch}\, N / E)$ for the $E$ extrema of $r_{k-1}$; from the second extraction a value that does not exceed the previous half width becomes $\lceil 1.1\, l_{k-1} \rceil$. The prototype is resampled to $2 l_k + 1$ weights and the iteration

$f^{(q)} = f^{(q-1)} - \mathcal{L}_w f^{(q-1)}, \qquad f^{(0)} = r_{k-1},$

with $\mathcal{L}_w$ the moving average under boundary, stops at the first $q$ with $\|f^{(q)} - f^{(q-1)}\|_2 / \|f^{(q-1)}\|_2 \le$ tolerance (converged TRUE) or after maxIterations (converged FALSE). Then $c_k = f^{(Q)}$ and $r_k = r_{k-1} - c_k$. A shrinking prototype is integrated over rescaled bins; an enlarged one is interpolated linearly and normalized to a unit absolute sum. The components are not shown to satisfy the counting conditions of intrinsic mode functions.

Value

A list with

modes	Double matrix, samples by components.
residue	Numeric vector, the final residue.
iterations	Integer vector, the iterations $Q$ of each component.
converged	Logical vector by component.
stopReason	"terminalResidue" or "modeLimit".

References

Cicone, A. and Zhou, H. (2018). Numerical Analysis for Iterative Filtering with New Efficient Implementations Based on FFT. arXiv:1802.01359.

The support rule, the resampling of the prototype and both iterations are adapted from fifpy, revision 26e51e7d45e551950dea8be37ece447fc4340296, copyright (c) 2022 Emanuele Papini, MIT License; see the NOTICE file of the installed package.

See Also

decomposeLMD, decomposeEMD

Examples

Time <- seq_len(1000) / 1000
Signal <- cos(2 * pi * 2 * Time) + 0.25 * cos(2 * pi * 24 * Time) +
    0.0625 * cos(2 * pi * 288 * Time)
Filter <- c(6.6306270449357867e-06, 0.00053827687799438354,
            0.0082858777882808092, 0.044245208647455882, 0.11927281750086954,
            0.20513332487239222, 0.24503568511284984, 0.20513332487239222,
            0.11927281750086952, 0.044245208647455903, 0.0082858777882808109,
            0.00053827687799438213, 6.6306270449358087e-06)
OUT <- decomposeIF(Signal, filter = Filter, boundary = "periodic",
                   tolerance = sqrt(0.001), maxIterations = 200,
                   stretch = 1.6, maxModes = 3)
OUT$iterations
OUT$stopReason
Zero <- decomposeIF(Signal, filter = Filter, boundary = "zero",
                    tolerance = sqrt(0.001), maxIterations = 200,
                    stretch = 1.6, maxModes = 3)
Zero$iterations

---

Iterated Masking Empirical Mode Decomposition

Description

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.

See Also

decomposeMaskedEMD, computeInstantaneousAttributes

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
OUT$maskIterations
OUT$stopReason

---

Local Mean Decomposition

Description

Decomposes a signal by successive extraction of product functions, each obtained by repeated subtraction of a smoothed local mean and division by a smoothed local magnitude.

Usage

decomposeLMD(signal, maxIterations, envelopeTolerance, changeTolerance,
             maxSmoothPasses, maxModes = Inf)

Arguments

signal	Numeric vector of finite values without dim, in the units of the quantity. Positions are sample indices.
maxIterations	Whole number $\ge 1$ of normalizations allowed in one product function. No default; the package tests use 200.
envelopeTolerance	Number $\ge 0$; an extraction ends when the mean of $|1 - a_q|$ over the samples is at most this value. No default; the package tests use 0.01.
changeTolerance	Number $\ge 0$; an extraction ends when the mean of $|s_{q-1} - s_q|$ is at most this value. No default; the package tests use 0.01.
maxSmoothPasses	Whole number $\ge 1$ of passes of the moving average that smooths the local levels. No default; 12 is the number of passes that LMD 1.2.1 reaches by default.
maxModes	Whole number $\ge 1$ or Inf, the limit of product functions; default Inf.

Details

The local levels of a function use as nodes its interior extrema and its first and last samples: consecutive nodes $n_i, n_{i+1}$ give the mean $m_i = (n_i + n_{i+1}) / 2$ and the magnitude $a_i = |n_i - n_{i+1}| / 2$, held between the two nodes. Both levels are smoothed by a repeated uniform moving average whose half width follows from the largest gap between nodes and shrinks at the borders; passes repeat while two consecutive samples are equal, up to maxSmoothPasses.

From $s_0 = r_{p-1}$, iteration $q$ forms

$s_q = \frac{s_{q-1} - m_q}{a_q}$

with the smoothed levels of $s_{q-1}$. The extraction ends without that normalization, with converged TRUE, when one of the two tolerances is met; it ends with converged FALSE when $s_{q-1}$ has at most one extremum or after maxIterations normalizations. With $Q$ normalizations, the product function is $a_1 \cdots a_Q s_Q$ and $r_p = r_{p-1} - PF_p$. Before extraction $p$ the decomposition ends with "terminalResidue" when $r_{p-1}$ has at most one extremum and with "modeLimit" when $p - 1$ equals maxModes. A product function combines an accumulated amplitude envelope with a normalized oscillation; this construction does not impose the IMF envelope conditions of EMD.

In the first iteration $a_1$ and $s_0 = r_{p-1}$ carry the units of the signal, so both tolerances compare magnitudes in those units with dimensionless quantities; the result depends on the unit only when the initial magnitude approaches 1 in that unit.

Value

A list with

modes	Double matrix, samples by product functions.
residue	Numeric vector, the final residue.
iterations	Integer vector, the number $Q$ of normalizations of each product function.
converged	Logical vector by product function.
stopReason	"terminalResidue" or "modeLimit".

Conditions

A local magnitude that is not positive leaves the normalization undefined and stops the function with an error of class xemdZeroEnvelope.

References

Smith, J. S. (2005). The local mean decomposition and its application to EEG perception data. Journal of the Royal Society Interface, 2(5), 443–454. doi:10.1098/rsif.2005.0058

The local levels and the extraction of a product function are adapted from LMD 1.2.1 (Shubhra Prakash), revision ec11b721ec7b05f0137eeb67a619b28f62ef85a2, Apache License (>= 2); see the NOTICE file of the installed package.

See Also

decomposeIF, decomposeEMD

Examples

Time <- seq(0, 9, length.out = 3000)
Signal <- sin(pi * Time) + sin(2 * pi * Time) + sin(6 * pi * Time) +
    0.5 * Time
OUT <- decomposeLMD(Signal, maxIterations = 200, envelopeTolerance = 0.01,
                    changeTolerance = 0.01, maxSmoothPasses = 12,
                    maxModes = 8)
dim(OUT$modes)
OUT$iterations
OUT$stopReason

---

Masked Empirical Mode Decomposition

Description

Decomposes a signal by successive extraction in which each mode is the average first mode of the residue plus sinusoidal masks of several phases, under one of three mask rules: the opposite sine pair of 2005, dyadic masks (EMD-DMS) and uniformly distributed phases (UPEMD).

Usage

decomposeMaskedEMD(signal, samplingInterval, masking, 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. Mask frequencies are in Hz.
masking	Named list with the field rule and exactly the fields of that rule, without defaults: "pair" frequency, in Hz, in $(0, 1/(2\Delta t)]$, for instance the value of estimateMaskFrequency; amplitude, greater than zero, in the units of the signal. "dyadic" firstFrequency, in Hz, in $(0, 1/(2\Delta t)]$; amplitudeRatio, the dimensionless $V_m > 0$. The example of the EMD-DMS article uses $V_m = 10$. "uniformPhase" amplitudeRatio, the dimensionless $\varepsilon_0 > 0$; phaseCount, the whole number $P \ge 2$ of phases. The package tests use $\varepsilon_0 = 1$ and $P = 4$.
sifting	Named list of sifting controls, as in decomposeEMD; the rule "consecutive" is rejected.

Details

From $r_0 = x$, before stage $k$ the decomposition ends, in this order, with "terminalResidue" when $r_{k-1}$ has at most one extremum, with "modeLimit" when $k - 1$ equals sifting$maxModes or the schedule of masks is exhausted, and with "insufficientExtrema" when the envelopes of $r_{k-1}$ cannot be built. Otherwise the masks of stage $k$ are $a_k \cos(2 \pi f_k n \Delta t + \theta_p)$, $n = 0, \ldots, N - 1$, with the first sample as phase origin, and the mode is the average of the first sifted component of $r_{k-1}$ plus each mask,

$c_k = \frac{1}{P} \sum_{p=1}^{P} E_1\left(r_{k-1} + w_p\right), \qquad r_k = r_{k-1} - c_k.$

A stage is retained unless one of its phases has no first mode or loses its envelopes during sifting; then the decomposition stops before adding the stage mean or subtracting it from the residue. A phase that reaches the iteration limit is kept and flags the stage FALSE. The phases and divisor are unchanged. The rules fix the phases, frequencies and amplitudes:

"pair"

Phases $-\pi/2$ and $\pi/2$, so the masks are $\pm a \sin(2 \pi f n \Delta t)$ and the mode is $(E_1(x + w) + E_1(x - w)) / 2$. The schedule has one stage: the decomposition ends with "modeLimit" after one extraction unless the residue is terminal.

"dyadic"

Phases 0 and $\pi$; $f_k = 2^{1-k} f_1$; amplitude $V_m$ times the mean absolute value of the maxima of the residue.

"uniformPhase"

Phases $2 \pi p / P$, $p = 0, \ldots, P - 1$; $f_k = 1 / (2^k \Delta t)$; amplitude $\varepsilon_0$ times the standard deviation of the residue with divisor $N$.

The masks are deterministic: the function draws no random numbers.

Value

A list with

modes	Double matrix, samples by modes.
residue	Numeric vector, the final residue.
iterations	Integer vector by mode, the sum of the subtractions of the extractions of the stage.
converged	Logical vector by mode: TRUE when every extraction of the stage accepted; NA when all completed fixed sifting. Interrupted stages are reported in termination.
stopReason	"terminalResidue", "modeLimit", "insufficientExtrema" or "incompleteExtraction", the last when at least one phase fails.
maskFrequencies	Numeric vector, the mask frequency of each mode in Hz.
termination	Data frame with mode, the next stage index after the retained modes; phase, its one-based index in the phase sequence; iterations, subtractions of that failed phase; and stopReason, "terminalResidue", "modeLimit", "insufficientExtrema" or "iterationLimit". One row is returned for each phase of the interrupted stage that did not accept, including the phases that reached the iteration limit. If the decomposition stops before a stage starts, one row carries its stopReason, phase = NA and zero iterations.

References

Deering, R. and Kaiser, J. F. (2005). The use of a masking signal to improve empirical mode decomposition. In IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP), 4, 485–488.

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.

See Also

estimateMaskFrequency, decomposeIteratedMaskedEMD, decomposeEMD

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)
# Uniform phases, with one sample per second.
OUT <- decomposeMaskedEMD(Signal, samplingInterval = 1,
                          masking = list(rule = "uniformPhase",
                                         amplitudeRatio = 1, phaseCount = 4),
                          sifting = Sifting)
OUT$maskFrequencies
OUT$stopReason

# The pair of 2005 at the estimated mask frequency.
Frequency <- estimateMaskFrequency(Signal, samplingInterval = 1,
                                   sifting = Sifting)
Pair <- decomposeMaskedEMD(Signal, samplingInterval = 1,
                           masking = list(rule = "pair", frequency = Frequency,
                                          amplitude = sqrt(mean((Signal -
                                              mean(Signal))^2))),
                           sifting = Sifting)
dim(Pair$modes)
Pair$stopReason

---