Package {xemd}

Type:	Package
Title:	eXtended Empirical Mode Decomposition
Version:	0.1.5
Description:	Signals are decomposed into oscillatory modes by empirical mode decomposition and its ensemble, complementary ensemble, complete ensemble, masked, bivariate and multivariate variants, by variational mode decomposition and its multivariate and successive variants, and by local mean decomposition, iterative filtering, time varying filter empirical mode decomposition and the empirical wavelet transform; the synchrosqueezed wavelet transform is also provided. Empirical mode decomposition is due to Huang et al. (1998) <doi:10.1098/rspa.1998.0193>, its ensemble variant to Wu and Huang (2009) <doi:10.1142/S1793536909000047> and variational mode decomposition to Dragomiretskiy and Zosso (2014) <doi:10.1109/TSP.2013.2288675>. The modes can be described by their instantaneous amplitude and frequency and their Hilbert spectrum, tested as intrinsic mode functions or against white noise, selected by energy and compared through orthogonality, mode splitting, reconstruction quality, ensemble dispersion and similarity measures. Acceleration, velocity and displacement records stored as a long table of time, value, quantity and component are decomposed for every record and quantity in one call, with a warning and a regularization of the time grid when it is not uniform.
URL:	https://kanameishi.github.io/xemd/
License:	GPL-3
Copyright:	Authors, copyright notices and licenses of the third-party code adapted in xemd are listed in file inst/NOTICE.
Encoding:	UTF-8
Depends:	R (≥ 4.1.0)
Imports:	data.table, pracma, splines, stats, utils
Suggests:	hht, knitr, rmarkdown
VignetteBuilder:	knitr
NeedsCompilation:	no
Packaged:	2026-10-24 20:27:01 UTC; averrik
Author:	Alejandro Verri Kozlowski [aut, cre, cph, trl], G. Rilling [ctb] (emd.m code), Marcelo A. Colominas [cph] (CEEMDAN distribution of emd.m), Pierre Masselot [ctb] (emdr code), Vinícius Rezende Carvalho [ctb, cph] (ewtpy code), Eduardo Mazoni [cph] (ewtpy code), Jerome Gilles [ctb, cph] (pyewt code; MATLAB EWT code translated in ewtpy), Kathryn Heal [ctb] (MATLAB EWT boundary detection translated in ewtpy), Emanuele Papini [ctb, cph] (fifpy code), Antonio Cicone [ctb] (MATLAB iterative filtering code translated in fifpy), Shubhra Prakash [ctb] (LMD code), Stefano Bianchi [ctb, cph] (pytvfemd code), Alessandro Longo [ctb, cph] (pytvfemd code), Eugene Brevdo [ctb] (Synchrosqueezing Toolbox code ported in SynchWave), Hau-Tieng Wu [ctb] (Synchrosqueezing Toolbox code ported in SynchWave), Jianfeng Lu [ctb] (curve_ext of the Synchrosqueezing Toolbox, ported in SynchWave), Dongik Jang [ctb] (SynchWave code), Hee-Seok Oh [ctb] (SynchWave code), Donghoh Kim [ctb] (SynchWave code)
Maintainer:	Alejandro Verri Kozlowski <averri@fi.uba.ar>
Repository:	CRAN
Date/Publication:	2026-11-08 17:10:02 UTC

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

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:

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

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:

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

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:

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

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

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

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:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Multivariate Empirical Mode Decomposition

Description

Decomposes a signal of two or more channels by successive extraction of vector modes, with vector envelopes built at the maxima of its projections onto directions of the channel space.

Usage

decomposeMEMD(signal, directionCount = 64, sifting)

Arguments

Details

The directions map D Hammersley points, built from the radical inverses in the bases of the first C - 1 primes, to unit vectors \mathbf{v}_d through the angles of get.projections of emdr 1.2.3 (Pierre Masselot); their antipodes complete 2D directions. For each direction, \mathbf{e}_d interpolates the vector samples at the maxima of the projection p_d = \mathbf{v}_d^{T} \mathbf{x}, extended by reflection, and the sifting subtracts the mean of the 2D envelopes,

\mathbf{h}_{s+1} = \mathbf{h}_s - \frac{1}{2D} \sum_{d=1}^{2D} \mathbf{e}_d[\mathbf{h}_s].

The adaptive rules act on the projections. "consecutive" accepts when every projection keeps |E - Z| \le 1 in consecutiveCount successive states after a subtraction, the run restarting when one does not. "doubleThreshold" accepts a state, the signal included, whose projections all keep |E - Z| \le 1 and whose ratio |\mathbf{m}| / a meets the double threshold of decomposeEMD; a is the mean distance of the envelopes to their mean with "globalMean", or the mean half distance between the envelopes of antipodal directions with "envelopePairs". Before extraction k the decomposition ends with "terminalResidue" when every projection of \mathbf{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 all the channels, so they share the number of modes.

Value

A list with

References

Rehman, N. and Mandic, D. P. (2010). Multivariate empirical mode decomposition. Proceedings of the Royal Society A, 466, 1291–1302. doi:10.1098/rspa.2009.0502

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

decomposeNAMEMD, decomposeBEMD, decomposeMVMD, decomposeEMD

Examples

Time <- 2 * pi * (0:255) / 256
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 <- decomposeMEMD(Signal, directionCount = 64,
                     sifting = list(stopRule = "fixed", fixedSift = 10,
                                    boundary = "reflection", maxModes = 3))
dim(OUT$modes)
OUT$stopReason

Multivariate Variational Mode Decomposition

Description

Decomposes the channels of a sampled multichannel signal into modes whose channel coordinates share one center frequency.

Usage

decomposeMVMD(signal, samplingInterval, initialCenters, alpha, tau,
              fixedZeroCenter, tolerance, maxIterations = 499)

Arguments

Details

Each channel is extended and transformed as in decomposeVMD. In each sweep the coordinate c of mode k is

\hat u_{k,c} = \frac{\hat x_c - \sum_{i \ne k} \hat u_{i,c} + \hat\lambda_c / 2}{1 + 2 \alpha (\omega - \omega_k)^2},

with a denominator common to the channels; the new center weights |\omega| by the power summed over the channels, and the multiplier of each channel is updated with step \tau. The iteration stops with "tolerance" when \sum_k \sum_c \|\hat u_{k,c}^{(n+1)} - \hat u_{k,c}^{(n)}\|^2 / \|\hat u_{k,c}^{(n)}\|^2 < \epsilon, which cannot happen before the second sweep, or with "iterationLimit" after maxIterations sweeps. A coordinate without energy in the previous sweep leaves its ratio undefined, so the test fails until maxIterations. The channels keep the order received. The joint bandwidth formulation does not impose the geometrical IMF conditions of EMD. A shared spectral center does not require identical amplitudes or instantaneous frequencies across channels.

Value

A list with

Conditions

A mode without energy over all its channels leaves its center undefined and stops the function with an error of class xemdZeroEnergy.

References

ur Rehman, N. and Aftab, H. (2019). Multivariate Variational Mode Decomposition. arXiv:1907.04509.

See Also

decomposeVMD, decomposeSVMD, decomposeTSL

Examples

SamplingInterval <- 1 / 1000
Time <- seq_len(1000) * SamplingInterval
Signal <- cbind(cos(2 * pi * 2 * Time) + 0.25 * cos(2 * pi * 24 * Time) +
                    0.0625 * cos(2 * pi * 288 * Time),
                cos(2 * pi * 2 * Time),
                0.25 * cos(2 * pi * 24 * Time))
OUT <- decomposeMVMD(Signal, samplingInterval = SamplingInterval,
                     initialCenters = (0:2) / (6 * SamplingInterval),
                     alpha = 2000 * SamplingInterval^2 / (8 * pi^2),
                     tau = 0.5, fixedZeroCenter = FALSE, tolerance = 1e-6,
                     maxIterations = 98)
dim(OUT$modes)
OUT$centers
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

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

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

Noise-Assisted Multivariate Empirical Mode Decomposition

Description

Decomposes a signal together with supplied noise channels by one multivariate decomposition and keeps the coordinates of the signal.

Usage

decomposeNAMEMD(signal, noiseChannels, directionCount = 64, sifting)

Arguments

Details

The function runs decomposeMEMD on the augmented matrix \mathbf{y} = [\mathbf{x}, \mathbf{w}] and keeps the first C coordinates of every mode, \mathbf{c}_k = P_x \tilde{\mathbf{c}}_k, and of the residue. Projections, stopping, iterations and termination are those of the augmented decomposition. In decomposeTSL the noise channels are drawn by group with noiseChannelCount and scaled by noiseRatio times the reference deviation of the group.

Value

A list with the fields of decomposeMEMD; modes has the C channels of signal and residue its C columns.

References

ur Rehman, N. and Mandic, D. P. (2011). Filter Bank Property of Multivariate Empirical Mode Decomposition. IEEE Transactions on Signal Processing, 59(5), 2421–2426. doi:10.1109/TSP.2011.2106779

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, decomposeTSL

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))
set.seed(628)
# Two noise channels with 0.02 times the reference deviation of the signal.
Noise <- 0.02 * sqrt(mean(apply(Signal, 2, stats::sd)^2)) *
    matrix(stats::rnorm(512 * 2), nrow = 512)
OUT <- decomposeNAMEMD(Signal, noiseChannels = Noise, directionCount = 64,
                       sifting = list(stopRule = "fixed", fixedSift = 1,
                                      boundary = "reflection",
                                      maxModes = 3))
dim(OUT$modes)
dim(OUT$residue)

Successive Variational Mode Decomposition

Description

Extracts variational modes one at a time from the whole signal, keeping the centers of the modes already extracted fixed, and stops after a known number of modes or when the residual power meets a known noise power.

Usage

decomposeSVMD(signal, samplingInterval, initialCenter = 0, alpha, maxAlpha,
              tau, tolerance, maxIterations = 299, stopping)

Arguments

Details

The signal is extended and transformed as in decomposeVMD. Mode L is extracted from the whole signal with the centers \omega_1, \ldots, \omega_{L-1} of the previous modes fixed. Its penalty starts at alpha and doubles while below maxAlpha; when doubling would exceed maxAlpha, one last computation runs at maxAlpha. Each penalty restarts the inner iteration from the spectrum and center reached, with a zero multiplier and the penalty H_{L-1}(\omega) = \sum_{i<L} 1 / (\alpha^2 (\omega - \omega_i)^4). With A = \alpha^2 (\omega - \omega_L)^4, the inner update is

\hat u_L \leftarrow \frac{\hat x + A \hat u_L + \hat\lambda / 2}{(1 + A)\left(1 + 2 \alpha (\omega - \omega_L)^2 + H_{L-1}\right)},

followed by the new center and the dual ascent, and it stops when the relative change of \hat u_L is below tolerance or after maxIterations. Extraction stops with "modeCount" after modeCount modes; with "noisePower" when |\sigma^2 - P_L| / \sigma^2 < \epsilon_2, P_L the mean square of the residual on the record; or with "modeLimit" after maxModes modes without meeting it. With maxModes = Inf the extraction under "noisePower" continues until the condition holds. The modes keep the extraction order. Their bandwidth formulation does not impose the geometrical IMF conditions of EMD.

Value

A list with

Conditions

A mode without spectral energy leaves its center undefined and stops the function with an error of class xemdZeroEnergy.

References

Nazari, M. and Sakhaei, S. M. (2020). Successive variational mode decomposition. Signal Processing, 107610. doi:10.1016/j.sigpro.2020.107610

See Also

decomposeVMD, decomposeMVMD

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 <- decomposeSVMD(Signal, samplingInterval = SamplingInterval,
                     initialCenter = 0,
                     alpha = 10 * SamplingInterval^2 / (4 * pi^2),
                     maxAlpha = 20000 * SamplingInterval^2 / (4 * pi^2),
                     tau = 0, tolerance = 1e-6, maxIterations = 299,
                     stopping = list(stopRule = "modeCount", modeCount = 5))
OUT$centers
OUT$stopReason

Decomposition of Time-Series Records

Description

Decomposes the series of a long table of records, grouped by record and quantity, with one method of the package, and returns homogeneous tables for all the methods.

Usage

decomposeTSL(TSL, method, control, by)

Arguments

Details

A group is a record and an ID; its OCID are the channels of the group and must share one time grid. Groups follow by ascending and the IDs "AT", "VT" and "DT" first, then the other IDs ascending; the OCID of a group ascend. The order of the rows of TSL does not intervene. A series needs two samples or more, a time cannot repeat within a series, and the metadata other than by must be constant within each record, ID and OCID.

The grid of a group is uniform when diff(t) equals the mean step \Delta t = (t_N - t_1)/(N - 1) under all.equal; it is then used unchanged, with samplingInterval \Delta t. A grid that is not uniform is regularized with one warning per group: the step becomes 1/F_s, with F_s the smallest 2^a 5^b not below 1/\Delta t; the grid is t_1 + k/F_s, and each OCID is interpolated on it by stats::splinefun(method = "monoH.FC"). The table signal holds the series actually decomposed.

Scalar methods decompose each OCID separately with the same controls. The coupled methods decompose the matrix of the OCID of a group in ascending order, without standardization: decomposeBEMD needs exactly two OCID and decomposeMEMD at least two; their iterations, converged, stopReason and centers are common and repeat for every OCID. The ensembles draw realizationCount unit normal realizations for each OCID from the random stream of the caller, in the order group, OCID, realization. decomposeNAMEMD draws noiseChannelCount channels for each group and scales them by noiseRatio times \sqrt{\mathrm{mean}(s_c^2)}, s_c the sample deviations of the OCID of the group; the noise channels give no rows. A control valid for one group and not for another, such as a center above the Nyquist frequency of a group, stops the operation. Units are not inferred, validated or converted; the metadata are carried to the output.

Channels are coupled within a record and an ID. The operation does not impose a kinematic relationship between separate AT, VT and DT groups. Equal mode indices across those groups do not establish corresponding physical contributions.

Value

A named list of data.table tables, each with the by columns, the other metadata, OCID and ID, followed by

Rows follow the by columns, OCID, ID, mode or boundary, and t.

References

The regularization of the grid is adapted from .regularize and .rationalize of gmsp 0.8.0, copyright (c) 2026 Alejandro Verri Kozlowski, MIT License: gmsp starts from the smallest difference of the times and returns the grid from zero, whereas here the step starts from the mean and the grid keeps the origin t_1. See the NOTICE file of the installed package.

See Also

listDecompositions, decomposeEMD, decomposeMVMD

Examples

Time <- seq(0, 9, length.out = 3000)
Signal <- sin(pi * Time) + sin(2 * pi * Time) + sin(6 * pi * Time) +
    0.5 * Time
# The same series under three IDs and three components.
TSL <- data.table::rbindlist(lapply(c("AT", "VT", "DT"), function(id) {
    data.table::rbindlist(lapply(c("H1", "H2", "UP"), function(ocid) {
        data.table::data.table(ID = id, OCID = ocid, t = Time, s = Signal)
    }))
}))
OUT <- decomposeTSL(TSL, method = "decomposeEMD",
                    control = list(sifting = list(stopRule = "doubleThreshold")),
                    by = character())
names(OUT)
OUT$termination
head(OUT$modes)

# A coupled method: three channels of one record decomposed together.
SamplingInterval <- 1 / 1000
TimeV <- seq_len(1000) * SamplingInterval
Series <- cbind(cos(2 * pi * 2 * TimeV) + 0.25 * cos(2 * pi * 24 * TimeV) +
                    0.0625 * cos(2 * pi * 288 * TimeV),
                cos(2 * pi * 2 * TimeV), 0.25 * cos(2 * pi * 24 * TimeV))
Record <- data.table::data.table(ID = "AT",
                                 OCID = rep(c("H1", "H2", "UP"), each = 1000),
                                 t = rep(TimeV, 3), s = as.vector(Series))
Coupled <- decomposeTSL(Record, method = "decomposeMVMD",
                        control = list(initialCenters = (0:2) / (6 * SamplingInterval),
                                       alpha = 2000 * SamplingInterval^2 / (8 * pi^2),
                                       tau = 0.5, fixedZeroCenter = FALSE,
                                       tolerance = 1e-6, maxIterations = 98),
                        by = character())
Coupled$centers

Time-Varying Filter Empirical Mode Decomposition

Description

Decomposes a signal by successive extraction of components, each obtained by sifting with a local mean that is a B-spline approximation on knots placed by a time-varying cutoff frequency (TVF-EMD).

Usage

decomposeTVFEMD(signal, bandwidthThreshold, splineOrder,
                intermittencyThreshold, maxSift, maxModes = Inf)

Arguments

Details

Before extraction k the decomposition ends with "terminalResidue" when r_{k-1} has fewer than four extrema, maxima plus minima, and with "modeLimit" when k - 1 equals maxModes. A component is obtained by sifting,

h_{s+1} = h_s - m_s,

where m_s is the periodic least-squares B-spline approximation of order n of h_s on knots at the extrema of \cos\varphi, \varphi the cumulative realigned cutoff. Each iteration extends h_s by symmetric reflection of \mathrm{round}(N/2) samples per side and restores the record afterwards. The cutoff follows from the amplitude A and the frequency of the analytic signal, interpolated at the minima and at the maxima of A; values above 0.9\pi are lowered to 0.9\pi, undefined or non-positive values are interpolated, and the zones marked by intermittency events, where the relative variation of the cutoff between consecutive maxima of the signal exceeds \rho, are replaced by interpolation.

The stopping rule aggregates the relative instantaneous bandwidth \theta over the samples from \mathrm{round}(0.05N) to \mathrm{round}(0.95N): h_s is accepted when the mean of \theta is below \xi (1 + (s + 1)/4), from s = 1 in the first component and from s = 0 in the later ones, or when the samples with \theta > \xi are fewer than 0.2N. This aggregation of pytvfemd 0.3.14 replaces the local criterion of the monograph. converged is FALSE when maxSift subtractions end without acceptance or when the cutoff or the knots cannot be built. Acceptance refers to the aggregated rule above; it does not establish the pointwise narrow-band condition. Neither that rule nor a returned component guarantees the two classical IMF conditions.

Value

A list with

References

The stopping rule and the reflection of the sifting, the end of the decomposition, the interpolation, the bounds and the realignment of the cutoff, the relative bandwidth, the knots and the periodic B-spline approximation are adapted from pytvfemd 0.3.14 (Stefano Bianchi and Alessandro Longo), GPL-3; see the NOTICE file of the installed package.

See Also

decomposeEMD, computeInstantaneousAttributes

Examples

Time <- seq(0, 9, length.out = 3000)
Signal <- sin(pi * Time) + sin(2 * pi * Time) + sin(6 * pi * Time) +
    0.5 * Time
OUT <- decomposeTVFEMD(Signal, bandwidthThreshold = 0.1, splineOrder = 28,
                       intermittencyThreshold = 0.25, maxSift = 5)
dim(OUT$modes)
OUT$iterations
OUT$stopReason

Variational Mode Decomposition

Description

Decomposes a sampled signal into a given number of modes, each concentrated around a center frequency, by alternating updates of the mode spectra, the centers and a multiplier.

Usage

decomposeVMD(signal, samplingInterval, initialCenters, alpha, tau,
             fixedZeroCenter, tolerance, maxIterations)

Arguments

Details

The record of N samples is extended to T = 2N samples by reflecting \lceil N/2 \rceil samples at the start and the rest at the end, and transformed with the unnormalized discrete Fourier transform; the grid holds |\omega| in rad/s. Sweeps start from zero mode spectra and multiplier. In each sweep mode k, updated after modes 1, \ldots, k - 1, is

\hat u_k = \frac{\hat x - \sum_{i \ne k} \hat u_i + \hat\lambda / 2}{1 + 2 \alpha (\omega - \omega_k)^2},

its center becomes the power-weighted mean of |\omega|, and the sweep closes with \hat\lambda \leftarrow \hat\lambda + \tau (\hat x - \sum_k \hat u_k). The iteration stops with "tolerance" when \sum_k \|\hat u_k^{(n+1)} - \hat u_k^{(n)}\|^2 / \|\hat u_k^{(n)}\|^2 < \epsilon, which cannot happen before the second sweep, or with "iterationLimit" after maxIterations sweeps; the last sweep is returned. The modes are the inverse transforms restricted to the record. With tau = 0 the penalized formulation permits a nonzero fit residual. VMD defines its modes through a bandwidth objective, without imposing the geometrical IMF conditions of EMD.

Value

A list with

Conditions

A mode without spectral energy leaves its center undefined and stops the function with an error of class xemdZeroEnergy.

References

Dragomiretskiy, K. and Zosso, D. (2014). Variational Mode Decomposition. IEEE Transactions on Signal Processing, 62(3), 531–544. doi:10.1109/TSP.2013.2288675

See Also

decomposeMVMD, decomposeSVMD, decomposeTSL

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 <- decomposeVMD(Signal, samplingInterval = SamplingInterval,
                    initialCenters = (0:2) / (6 * SamplingInterval),
                    alpha = 2000 * SamplingInterval^2 / (8 * pi^2),
                    tau = 0.5, fixedZeroCenter = FALSE, tolerance = 1e-6,
                    maxIterations = 98)
OUT$centers
OUT$stopReason
max(abs(OUT$fitResidual))

Identity of the Installed Module

Description

Returns the name and the installed version of the package.

Usage

describeModule()

Value

A list with package, the string "xemd", and version, the installed version as a string.

See Also

listOperations, listDecompositions

Examples

describeModule()

Mask Frequency of the 2005 Masking

Description

Estimates the frequency of the masking signal of 2005 from the instantaneous amplitude and frequency of the first sifted component of a signal.

Usage

estimateMaskFrequency(signal, samplingInterval, sifting)

Arguments

Details

With A_1 and \nu_1 the instantaneous amplitude and frequency, computed by computeInstantaneousAttributes, of the first sifted component y_1 = E_1(x) without mask, the estimate is

\bar\nu = \frac{\sum_n A_1[n] \, \nu_1[n]^2}{\sum_n A_1[n] \, \nu_1[n]},

with the sums over the samples where \nu_1 is defined. The result is the frequency of the rule "pair" of decomposeMaskedEMD.

Value

A number, the mask frequency in Hz.

Conditions

A signal without a first extracted component stops with an error. When the weights A_1 \nu_1 sum to zero the estimate is undefined, and the function stops with an error of class xemdZeroEnergy.

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

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)
# One sample per second, so the result is in cycles per sample.
estimateMaskFrequency(Signal, samplingInterval = 1, sifting = Sifting)
# The time step of the signal, 9 / 2999 s.
estimateMaskFrequency(Signal, samplingInterval = 9 / 2999, sifting = Sifting)

Components Along the Ridges of a Synchrosqueezed Transform

Description

Extracts ridges of a synchrosqueezed transform one after another by dynamic programming and reconstructs a component from the coefficients in a band around each ridge.

Usage

extractSSTModes(transform, frequency, modeCount, penalty = 1000,
                exclusionWidth = 4, reconstructionWidth, morletCenter)

Arguments

Details

With T_{n,j} the coefficients, f_j the frequencies and z_j = \log_2 f_j, the cost of a cell is E_{n,j} = -\log(|T_{n,j}|^2 / Q + 10^{-8}), Q the sum of |T|^2 over the matrix. The accumulated cost is F_{1,j} = E_{1,j} and

F_{n,j} = E_{n,j} + \min_k \left[F_{n-1,k} + \lambda (z_j - z_k)^2\right]

over every cell k; the ridge ends at the first cell of least F_{N,j}, and each earlier cell is the first whose residual in that recurrence is below 10^{-8} in modulus. After each ridge C_n, the cells from C_n - b to C_n + b, clipped to the cells, take the value \sqrt{\epsilon}, \epsilon the machine precision, and the next ridge is extracted with the energy recomputed. Component k sums the original coefficients of its band,

c_k[n] = \frac{1}{C} \mathrm{Re} \sum_{j = C_{k,n} - w}^{C_{k,n} + w} T_{n,j},

over the cells of the band, with C = C_\psi / (\sqrt{2\pi} \, 2 \ln 2) and C_\psi = \int_0^\infty \hat\psi(\xi) / \xi \, d\xi for the Morlet wavelet of center morletCenter. The bands may overlap, the components need not add up to the signal, and no residue is returned. The components are defined by the selected ridge bands, without imposing the IMF conditions of EMD.

Value

Double matrix, samples by components, one column per ridge.

Conditions

A transform of zero energy stops the function with an error of class xemdZeroEnergy. A sample where no cell satisfies the recurrence within the tolerance stops it with an error.

References

Daubechies, I., Lu, J. and Wu, H.-T. (2009). Synchrosqueezed wavelet transforms: a tool for empirical mode decomposition. arXiv:0912.2437.

Brevdo, E. and Thakur, G. Synchrosqueezing Toolbox, revision 7e9fec0c6c9ed478dafac4479c0658d24fc6d8ef, functions curve_ext (original code by Jianfeng Lu), curve_ext_multi and synsq_cwt_iw, https://github.com/ebrevdo/synchrosqueezing.

The extraction of the ridges is adapted from curve_ext (original author Jianfeng Lu, maintainer Eugene Brevdo) and curve_ext_multi of the Synchrosqueezing Toolbox, as ported in SynchWave 1.1.2, revision 9ac80b2c3e60b417d0ee5567a66e8e64faa1a979 (port by Dongik Jang), and the inversion from SynchWave 1.1.2 (MATLAB original by Eugene Brevdo; R port by Dongik Jang, Hee-Seok Oh and Donghoh Kim), LGPL-2.0-or-later, which xemd distributes under GPL-3 as section 3 of the GNU Library General Public License version 2 permits; see the NOTICE file of the installed package.

See Also

transformSST, decomposeEWT

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)
SST <- transformSST(Signal, samplingInterval = SamplingInterval,
                    voiceCount = 16, morletCenter = 2 * pi,
                    threshold = sqrt(.Machine$double.eps))
Modes <- extractSSTModes(SST$transform, frequency = SST$frequency,
                         modeCount = 3, reconstructionWidth = 0,
                         morletCenter = 2 * pi)
dim(Modes)
sqrt(colMeans(Modes^2))

Methods of the Record Operation

Description

Returns the names of the decomposition methods that decomposeTSL accepts.

Usage

listDecompositions()

Details

The scalar decompositions are applied to each OCID of a group and the coupled ones (decomposeBEMD, decomposeMEMD, decomposeNAMEMD and decomposeMVMD) to the OCID of a group together. transformSST, the analysis operations, resiftModes and estimateMaskFrequency have no record operation.

Value

Character vector of function names: the thirteen scalar methods followed by the four coupled ones.

See Also

decomposeTSL, listOperations

Examples

listDecompositions()

Exported Operations and Their Arguments

Description

Lists the exported functions of the package with the names of their formal arguments.

Usage

listOperations()

Details

The list reads the exports of the loaded namespace; it performs no decomposition and reads no files.

Value

A list with one element per exported function, in the order of sort() on their names; each element is a list with name, the function name, and arguments, the names of its formal arguments.

See Also

describeModule, listDecompositions

Examples

Operations <- listOperations()
length(Operations)
Operations[[1]]

Dynamic Time Warping Distance of Modes

Description

Computes the dynamic time warping distance between a signal and each of its modes.

Usage

measureDTW(signal, modes)

Arguments

Details

With the local cost d_{ij} = (x_i - c_j)^2 and the accumulation

G_{ij} = d_{ij} + \min\{G_{i-1,j-1}, G_{i-1,j}, G_{i,j-1}\},

G_{00} = 0 and infinite cost outside the matrix, the distance is \sqrt{G_{NN}}, without normalization by the length of the path. The work grows with the square of the number of samples for each mode.

Value

Numeric vector with one distance per mode, in the units of the signal.

See Also

measureFCF, selectModesEnergy

Examples

Time <- seq_len(1000) / 1000
Tones <- cbind(cos(2 * pi * 2 * Time), 0.25 * cos(2 * pi * 24 * Time),
               0.0625 * cos(2 * pi * 288 * Time))
measureDTW(signal = rowSums(Tones), modes = Tones)

Dispersion of Ensemble Members

Description

Computes, sample by sample and mode by mode, the sample variance of the members of an ensemble around their mean, and that sample variance divided by the number of members. The latter estimates the variance of the mean under the assumptions stated below.

Usage

measureDispersion(realizations)

Arguments

Details

With M realizations C^{(m)} and their mean \bar C, the sample variance is

s^2 = \frac{1}{M - 1} \sum_{m=1}^{M} \left(C^{(m)} - \bar C\right)^2

and the returned estimate for the variance of the mean is s^2 / M. The first describes dispersion between members. Interpreting the second as Monte Carlo precision assumes comparable components and independent realizations with a common mean and a common finite covariance. Opposite perturbations in a complementary pair are dependent; the two signs do not provide independent realizations for this interpretation.

The function does not establish these statistical assumptions or physical correspondence between equal mode indices. Dispersion under perturbations of a fixed record does not quantify all uncertainty in the physical signal, and increasing the number of members does not remove decomposition bias. Attributes of realizations other than dim do not reach the outputs.

Value

A list with two matrices, samples by modes:

See Also

decomposeEEMD

Examples

# Two realizations of four samples: the first has two modes, the second
# one, so its second mode is filled with zeros.
Realizations <- array(c(2, 4, 6, 8, 1, 1, 1, 1,
                        4, 6, 8, 10, 0, 0, 0, 0), dim = c(4, 2, 2))
measureDispersion(Realizations)

Fault Correlation Factor of Modes

Description

Computes the fault correlation factor between a signal and each of its modes, a normalized inner product without removal of means.

Usage

measureFCF(signal, modes)

Arguments

Details

\mathrm{FCF}(x, c) = \frac{\sum_n x[n] c[n]}{\sqrt{\sum_n x[n]^2 \sum_n c[n]^2}}.

The means are not removed, so the factor differs from a correlation of centered vectors when the inputs have a nonzero mean.

Value

Numeric vector with one factor per mode.

Conditions

A signal or a mode of zero norm stops the function with an error of class xemdZeroEnergy.

See Also

measureDTW, selectModesEnergy

Examples

Time <- seq_len(1000) / 1000
Tones <- cbind(cos(2 * pi * 2 * Time), 0.25 * cos(2 * pi * 24 * Time),
               0.0625 * cos(2 * pi * 288 * Time))
measureFCF(signal = rowSums(Tones), modes = Tones)

Metrics of the Modes of Time-Series Records

Description

Measures every mode and every series of the tables of decomposeTSL, for any of its methods: frequency, bandwidth, energy, energy fraction and mean period of each mode, and number of modes and energy fractions of the remainder and the energy difference described by crossFraction for each series.

Usage

measureModesTSL(decomposition)

Arguments

Details

A mode c of a series x with N samples uses the step \Delta t = (t_N - t_1)/(N - 1) of the grid of its group, with its samples ordered by t; the order of the rows does not intervene. Its spectrum is that of the VMD family: the mode is extended to 2N samples by reflecting \lceil N/2 \rceil samples at the start and the others at the end, P_j is the squared modulus of the unnormalized discrete Fourier transform of the extension and |\omega_j| = 2\pi \min(j, 2N - j)/(2N \Delta t), j = 0, \ldots, 2N - 1. Then frequency is the power centroid

f = \frac{\sum_j |\omega_j| P_j}{2\pi \sum_j P_j}

and bandwidth the power deviation around it on the same grid,

b = \left[\frac{\sum_j (|\omega_j|/2\pi - f)^2 P_j}{\sum_j P_j}\right]^{1/2},

both in Hz. frequency equals the centers of decomposeVMD and decomposeSVMD, except a center fixed at zero, which is not a centroid; in decomposeMVMD it is the centroid of each channel, not the common center. A mode of zero energy determines no centroid: its frequency and bandwidth are NA. energy is \Delta t \sum_n c[n]^2, energyFraction is \sum c^2 / \sum x^2, and meanPeriod is N \Delta t over the number of maxima of the mode, or N \Delta t without maxima.

For each series, modeCount is the number of its modes; residueFraction is \sum r^2 / \sum x^2, with r its residue or fitResidual, and zero when the method returns neither, as decomposeEWT; and crossFraction is

1 - \frac{\sum_k \sum c_k^2 + \sum r^2}{\sum x^2},

the input-energy fraction left after subtracting modal and remainder energies. It equals the normalized cross products when the modes and remainder reconstruct that input. Otherwise it also reflects the reconstruction discrepancy, as with EEMD's retained meanNoise. A series without modes has modeCount 0, and its remainder is the series.

Energy here is a quadratic signal measure, with squared signal units times seconds. Its identification with physical energy requires a relation for the measured quantity. A large modal energy fraction does not establish that the component is free of contamination.

Value

A named list of two data.table tables, each with the keys of decomposition followed by

Conditions

A series of signal with zero energy leaves its energy fractions undefined and stops the function with an error of class xemdZeroEnergy.

See Also

decomposeTSL, selectModesTSL, reconstructTSL

Examples

Time <- seq(0, 9, length.out = 3000)
Signal <- sin(pi * Time) + sin(2 * pi * Time) + sin(6 * pi * Time) +
    0.5 * Time
# The same series under three IDs and three components.
TSL <- data.table::rbindlist(lapply(c("AT", "VT", "DT"), function(id) {
    data.table::rbindlist(lapply(c("H1", "H2", "UP"), function(ocid) {
        data.table::data.table(ID = id, OCID = ocid, t = Time, s = Signal)
    }))
}))
Sifting <- list(stopRule = "doubleThreshold")
Decomposition <- decomposeTSL(TSL, method = "decomposeEMD",
                              control = list(sifting = Sifting),
                              by = character())
OUT <- measureModesTSL(Decomposition)
OUT$metrics[OUT$metrics$OCID == "H1" & OUT$metrics$ID == "AT", ]
OUT$metricsSeries

Inner Products and Cosines Between Components

Description

Computes the inner products between the columns of a matrix of components and the cosines of the angles between them.

Usage

measureOrthogonality(components)

Arguments

Details

The inner products are the Euclidean products \sum_n v_i[n] v_j[n], without the sampling interval of the discrete inner product, which cancels in the cosine

\cos\theta_{ij} = \frac{\langle v_i, v_j \rangle}{\|v_i\|_2 \|v_j\|_2}.

The cosine keeps its sign. The cross products are the terms by which the energy of a sum of components differs from the sum of their energies.

Value

A list of two square matrices, components by components:

Conditions

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

See Also

measurePMSI, measureFCF

Examples

Time <- seq_len(1000) / 1000
Tones <- cbind(cos(2 * pi * 2 * Time), 0.25 * cos(2 * pi * 24 * Time),
               0.0625 * cos(2 * pi * 288 * Time))
# The sum of the tones and the tones.
OUT <- measureOrthogonality(cbind(rowSums(Tones), Tones))
round(OUT$cosine, 6)

Pseudo-Mode Splitting Index

Description

Computes the pseudo-mode splitting index of each pair of adjacent modes.

Usage

measurePMSI(modes, denominator)

Arguments

Details

For each pair of adjacent columns,

\mathrm{PMSI}(c_i, c_{i+1}) = \max\left\{\frac{\langle c_i, c_{i+1} \rangle}{\|c_i\|_2^2 + \|c_{i+1}\|_2^2}, 0\right\},

with 10^{-8} added to the denominator under "regularized".

Value

Numeric vector with one value per pair of adjacent modes, of length one less than the number of columns.

Conditions

A pair of adjacent modes with zero energy under "homogeneous" stops the function with an error of class xemdZeroEnergy.

See Also

measureOrthogonality, decomposeMaskedEMD

Examples

Time <- seq_len(1000) / 1000
Tones <- cbind(cos(2 * pi * 2 * Time), 0.25 * cos(2 * pi * 24 * Time),
               0.0625 * cos(2 * pi * 288 * Time))
Mixed <- cbind(rowSums(Tones), Tones)
measurePMSI(Mixed, denominator = "homogeneous")
measurePMSI(Mixed, denominator = "regularized")

Quality of Reconstruction Factor

Description

Computes the quality of reconstruction factor of an estimate against a reference component, in decibels.

Usage

measureQRF(estimate, reference)

Arguments

Details

\mathrm{QRF}(\hat s, s) = 20 \log_{10} \frac{\|s\|_2}{\|s - \hat s\|_2},

Inf when the estimate equals the reference.

Value

A number in dB.

Conditions

A reference of zero norm is outside the domain of the factor and stops the function with an error.

See Also

measureFCF, measureOrthogonality

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)
# The signal without its third tone, against the signal.
measureQRF(estimate = Signal - 0.0625 * cos(2 * pi * 288 * Time),
           reference = Signal)
measureQRF(estimate = Signal, reference = Signal)

Reconstruction of Time-Series Records from Their Modes

Description

Adds up, for each series of the tables of decomposeTSL, the modes that a selection retains and, unless excluded, the remainder of the series.

Usage

reconstructTSL(decomposition, selection, includeResidue = TRUE)

Arguments

Details

The reconstruction of a series is the sum of its retained modes and, with includeResidue, its remainder. A series without retained modes sums zero, and a method without remainder, as decomposeEWT, adds none. With every mode and the remainder, the result follows the reconstruction relation of the selected method. For EMD the sum recovers the input up to rounding. In decomposeEEMD with the divisor "realizations", it gives the signal plus meanNoise. The order of the rows of the tables does not intervene.

A partial sum implements the declared selection; neither the sum nor the choice to include a remainder establishes that the result is free of noise.

Value

A named list with one data.table table, reconstruction: the keys of decomposition, t and s, on the grid and in the order of signal.

See Also

selectModesTSL, measureModesTSL, decomposeTSL

Examples

Time <- seq(0, 9, length.out = 3000)
Signal <- sin(pi * Time) + sin(2 * pi * Time) + sin(6 * pi * Time) +
    0.5 * Time
# The same series under three IDs and three components.
TSL <- data.table::rbindlist(lapply(c("AT", "VT", "DT"), function(id) {
    data.table::rbindlist(lapply(c("H1", "H2", "UP"), function(ocid) {
        data.table::data.table(ID = id, OCID = ocid, t = Time, s = Signal)
    }))
}))
Sifting <- list(stopRule = "doubleThreshold")
Decomposition <- decomposeTSL(TSL, method = "decomposeEMD",
                              control = list(sifting = Sifting),
                              by = character())
Metrics <- measureModesTSL(Decomposition)$metrics
# A minimum energy fraction of zero retains every mode: with the residue,
# the reconstruction is the signal.
Every <- selectModesTSL(Decomposition, metrics = Metrics,
                        criterion = list(rule = "energyFraction", minimum = 0,
                                         action = "keep"))
Complete <- reconstructTSL(Decomposition, selection = Every$selection)
max(abs(Complete$reconstruction$s - Decomposition$signal$s))
# The modes in [2, 15] Hz, without the residue.
Band <- selectModesTSL(Decomposition, metrics = Metrics,
                       criterion = list(rule = "band", fMin = 2, fMax = 15,
                                        action = "keep"))
Filtered <- reconstructTSL(Decomposition, selection = Band$selection,
                           includeResidue = FALSE)
head(Filtered$reconstruction)

Consecutive Post-Processing of Ensemble Modes

Description

Re-sifts a sequence of modes, typically the averaged modes of an ensemble, by extracting the first sifted component of consecutive sums and carrying each remainder to the next mode.

Usage

resiftModes(modes, sifting)

Arguments

Details

With u_1 = D_1 + D_2, stage i extracts a_i = E_1(u_i), keeps the remainder q_i = u_i - a_i and adds the next mode, u_{i+1} = q_i + D_{i+2}. The completed chain returns K - 1 modes and the residue q_{K-1} with "chainComplete", and

\sum_{i=1}^{K-1} a_i + q_{K-1} = \sum_{j=1}^{K} D_j.

When E_1(u_i) does not exist, the chain ends with "insufficientExtrema", the modes a_1, \ldots, a_{i-1} and the residue u_i + \sum_{j > i + 1} D_j, which keeps the same sum.

Value

A list with

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

decomposeEEMD, decomposeEMD

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.
Modes <- decomposeEEMD(Signal, noise = Noise[, 1:2], amplitudeRatio = 0.2,
                       divisor = "realizations", sifting = Sifting)$modes
OUT <- resiftModes(Modes, sifting = Sifting)
c(input = ncol(Modes), output = ncol(OUT$modes))
OUT$stopReason
max(abs(rowSums(OUT$modes) + OUT$residue - rowSums(Modes)))

Selection of Modes by Energy

Description

Orders the modes by decreasing energy and selects those before the largest drop between consecutive energies.

Usage

selectModesEnergy(modes)

Arguments

Details

The energy of mode k is P_k = \sum_n c_k[n]^2. With the permutation \pi that orders the energies decreasingly, the selection keeps \pi(1), \ldots, \pi(v) with

v = \arg\max_{1 \le m < K} \left(P_{\pi(m)} - P_{\pi(m+1)}\right).

The source fixes no tie-break, so a largest drop that is zero or not unique stops the function with an error.

Value

Integer vector of the selected column indices, the most energetic first.

See Also

measureFCF, measureDTW

Examples

Time <- seq_len(1000) / 1000
Tones <- cbind(cos(2 * pi * 2 * Time), 0.25 * cos(2 * pi * 24 * Time),
               0.0625 * cos(2 * pi * 288 * Time))
selectModesEnergy(Tones)
# The same tones in another order.
selectModesEnergy(Tones[, c(3, 1, 2)])

Selection of the Modes of Time-Series Records

Description

Selects the modes of the tables of decomposeTSL with a declared criterion: explicit mode indices, a minimum energy fraction or a frequency band, to keep the modes that meet it or to remove them.

Usage

selectModesTSL(decomposition, metrics, criterion)

Arguments

Details

The criterion is supplied by the caller. Mode indices, energy fractions and centroid frequencies do not by themselves identify uncontaminated physical components. Equal mode indices in independently decomposed series need not represent the same physical contribution.

The band compares frequency alone, whereas the rules keep.Fo and remove.Fo of gmsp 0.8.0 require the interval [F_o - BW, F_o + BW] inside the band.

Value

A named list of two data.table tables:

See Also

measureModesTSL, reconstructTSL, decomposeTSL

Examples

Time <- seq(0, 9, length.out = 3000)
Signal <- sin(pi * Time) + sin(2 * pi * Time) + sin(6 * pi * Time) +
    0.5 * Time
# The same series under three IDs and three components.
TSL <- data.table::rbindlist(lapply(c("AT", "VT", "DT"), function(id) {
    data.table::rbindlist(lapply(c("H1", "H2", "UP"), function(ocid) {
        data.table::data.table(ID = id, OCID = ocid, t = Time, s = Signal)
    }))
}))
Sifting <- list(stopRule = "doubleThreshold")
Decomposition <- decomposeTSL(TSL, method = "decomposeEMD",
                              control = list(sifting = Sifting),
                              by = character())
Metrics <- measureModesTSL(Decomposition)$metrics
# Keep the modes whose frequency lies in [2, 15] Hz.
Band <- selectModesTSL(Decomposition, metrics = Metrics,
                       criterion = list(rule = "band", fMin = 2, fMax = 15,
                                        action = "keep"))
Band$selection[Band$selection$OCID == "H1" & Band$selection$ID == "AT", ]
# Remove the first mode of every series of the ID VT.
Listed <- data.frame(ID = "VT", mode = 1)
First <- selectModesTSL(Decomposition, metrics = Metrics,
                        criterion = list(rule = "modes", modes = Listed,
                                         action = "remove"))
table(First$selection$ID, First$selection$retained)

Synchrosqueezed Wavelet Transform

Description

Computes the Morlet continuous wavelet transform of a sampled signal and reassigns its coefficients to frequency cells by the instantaneous frequency estimated from the transform.

Usage

transformSST(signal, samplingInterval, voiceCount = 16,
             morletCenter = 2 * pi, threshold)

Arguments

Details

The record is extended by symmetric padding to T = 2^{1 + \mathrm{round}(\log_2 N)} samples, and the transform uses the scales a_j = 2^{j / \mathrm{voiceCount}}, j = 1, \ldots, (\log_2 T - 1)\,\mathrm{voiceCount}. For each coefficient above the threshold, the frequency

\frac{\mathrm{Re}\,\Omega_f(a, b)}{2\pi} = \frac{1}{2\pi} \mathrm{Im}\frac{\partial_b W_f(a, b)}{W_f(a, b)}

in Hz goes to the nearest of J logarithmic cells from 1/((N - 1)\Delta t) to 1/(2\Delta t), J the number of scales, where the coefficient adds W a^{-1/2}; the sums are divided by voiceCount. Estimates that are not positive or fall outside the cells are not accumulated, and a transform whose energy is below the machine epsilon gives a zero result. The representation recovers the signal only approximately. The threshold is in the units of the signal, so decomposeTSL has no operation for this transform.

Value

A list with

References

Daubechies, I., Lu, J. and Wu, H.-T. (2009). Synchrosqueezed wavelet transforms: a tool for empirical mode decomposition. arXiv:0912.2437.

The wavelet transform, the Morlet wavelet and the reassignment are adapted from SynchWave 1.1.2, revision 9ac80b2c3e60b417d0ee5567a66e8e64faa1a979 (MATLAB original by Eugene Brevdo; R port by Dongik Jang, Hee-Seok Oh and Donghoh Kim; the Synchrosqueezing toolbox attribution also names Hau-Tieng Wu), LGPL-2.0-or-later, which xemd distributes under GPL-3 as section 3 of the GNU Library General Public License version 2 permits; see the NOTICE file of the installed package.

See Also

decomposeEWT, computeInstantaneousAttributes

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 <- transformSST(Signal, samplingInterval = SamplingInterval,
                    voiceCount = 16, morletCenter = 2 * pi,
                    threshold = sqrt(.Machine$double.eps))
dim(OUT$transform)
range(OUT$frequency)