Introduction

An observed signal can contain oscillations at different time scales, with amplitudes and frequencies that vary during the record. Empirical mode decomposition (EMD) uses the signal’s local extrema to construct components, rather than prescribing sinusoidal components of constant frequency. xemd provides EMD, its extensions and related decomposition methods, followed by operations for examining and reconstructing their outputs.

This guide follows two examples: successive extraction by EMD and joint fitting by variational mode decomposition (VMD). They illustrate different meanings of a mode, the controls used to obtain it and the diagnostics returned by the calculation. The methods guide covers the other families, and validation and interpretation explains the evidence needed to select components.

Modes and intrinsic mode functions

A mode is a component of a decomposition. In EMD, the target component is an intrinsic mode function (IMF): its numbers of extrema and zero crossings differ by at most one, and the mean of its upper and lower envelopes is zero at each point. These conditions describe oscillation geometry; they allow amplitude and frequency to vary and do not require a sinusoid of constant frequency. The envelope mean is different from the average of all samples.

Sifting repeatedly subtracts the envelope mean to approach those conditions. A numerical stopping rule specifies how that process ends; an iteration limit or a prescribed number of subtractions is not proof that the ideal IMF conditions hold. The returned diagnostics distinguish those outcomes.

The residual signal is what remains after extraction. It can contain a trend or a constant, and remains part of the additive reconstruction. Neither exact reconstruction nor approximate IMF conditions establish that a mode is a separate physical source. Different representations can reconstruct the same signal, and one physical oscillation can be spread over several extracted components.

Data model

decomposeEMD() and decomposeVMD() receive a finite numeric vector signal, in the units of the measured quantity. EMD operates on sample positions. VMD also receives samplingInterval in seconds and center frequencies in Hz. The seismic records guide describes the long-table interface for records, quantities and components.

The returned list depends on the method. For the EMD and VMD examples, the fields have the following meanings:

Field	EMD	VMD
`modes`	Matrix of samples by extracted modes	Matrix of samples by the prescribed modes
Remainder	`residue`, left after successive extractions	`fitResidual`, the input minus the sum of the fitted modes
`iterations`	Number of sifting subtractions for each mode	Number of joint sweeps, repeated for each mode
`converged`	Acceptance of each mode by its sifting rule	Result of the joint stopping test, repeated for each mode
`stopReason`	Reason extraction ended	Reason the joint iteration ended
`centers`	Absent	Final center frequencies in Hz

Empirical mode decomposition

EMD extracts a mode from the current residue and subtracts it before the next extraction. Within an extraction, sifting subtracts the mean of the upper and lower envelopes until the selected rule accepts the mode or the iteration ends. The example signal contains three sines and a trend, sampled at 3000 points. Its controls select the double-threshold rule, reflected boundaries and up to 2000 sifting subtractions per mode; maxModes = Inf lets the residue determine when extraction ends.

library(xemd)

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)
EMD <- decomposeEMD(Signal, sifting = Sifting)
str(EMD)
#> List of 5
#>  $ modes     : num [1:3000, 1:3] -0.815 -0.733 -0.652 -0.572 -0.492 ...
#>  $ residue   : num [1:3000] 0.468 0.467 0.467 0.467 0.466 ...
#>  $ iterations: int [1:3] 2 4 2
#>  $ converged : logi [1:3] TRUE TRUE TRUE
#>  $ stopReason: chr "insufficientExtrema"

EMD$converged describes the individual extractions: TRUE means that the stopping rule accepted the mode; FALSE means that maxSift was reached or that the envelopes could not be built. The fixed-sifting rule, available through decomposeEMD(), instead reports NA because it performs a prescribed number of subtractions.

EMD$stopReason describes the outer extraction process. "terminalResidue" means that the remaining signal has at most one extremum; "modeLimit" means that the requested extraction limit was reached; "insufficientExtrema" means that the boundary rule could not build the next envelopes. These reasons are separate from the acceptance of each extracted mode.

The controls affect the decomposition itself. Envelope construction near the record boundaries affects the extracted components, and additional sifting can alter their amplitude modulation. The example’s controls are not a universal setting for another signal.

The sum of the modes and the residue reconstructs the signal up to rounding. The following expression reports the largest absolute reconstruction difference in the units of the input signal.

max(abs(rowSums(EMD$modes) + EMD$residue - Signal))
#> [1] 1.332268e-15

This difference checks the additive identity. It does not measure whether each extracted component matches a generating oscillation of the signal.

assessIMF() provides numerical diagnostics of the IMF conditions for each column. The count result tests whether the numbers of extrema and zero crossings differ by at most one. The deviation result tests the envelope mean relative to the envelope amplitude with the double threshold, an approximation to the ideal zero envelope mean; it is NA when the boundary rule cannot construct the envelopes. The example also passes the residue as the final column, so its diagnosis remains distinct from those of the extracted modes.

assessIMF(cbind(EMD$modes, EMD$residue))
#> $count
#> [1]  TRUE  TRUE  TRUE FALSE
#> 
#> $deviation
#> [1] TRUE TRUE TRUE   NA

Variational mode decomposition

VMD fits a prescribed number of modes concentrated around center frequencies. It updates the mode spectra, their centers and a multiplier together. The length of initialCenters sets the number of modes; fixedZeroCenter = FALSE allows every center to change. VMD modes are not defined by imposing the geometric IMF conditions of EMD. Their defining property in this formulation is the bandwidth objective. An IMF diagnostic can still be computed for them, but it tests a property additional to that objective.

The signal contains tones of 2, 24 and 288 Hz sampled every millisecond. initialCenters is in Hz and alpha, the bandwidth penalty, is in \mathrm{s}^2. The expression for alpha converts the value 2000 from the cycles-per-sample convention used by sktime and VMDecomp into the units used by xemd. A positive tau selects the constrained problem; tau = 0 selects the penalized problem with a zero multiplier.

SamplingInterval <- 1 / 1000
TimeV <- seq_len(1000) * SamplingInterval
SignalV <- cos(2 * pi * 2 * TimeV) + 0.25 * cos(2 * pi * 24 * TimeV) +
    0.0625 * cos(2 * pi * 288 * TimeV)
VMD <- decomposeVMD(SignalV, samplingInterval = SamplingInterval,
                    initialCenters = (0:2) / (6 * SamplingInterval),
                    alpha = 2000 * SamplingInterval^2 / (8 * pi^2),
                    tau = 0.5, fixedZeroCenter = FALSE, tolerance = 1e-6,
                    maxIterations = 98)
VMD$centers
#> [1]   2.000004  23.999532 287.962079
VMD$stopReason
#> [1] "tolerance"

VMD$centers contains the final center frequencies. A stopReason of "tolerance" means that the relative change of the mode spectra passed the stopping test; "iterationLimit" means that the last permitted sweep was reached. In either case the returned modes come from the last sweep. VMD$fitResidual contains the difference between the signal and their sum; with tau = 0, the formulation permits a nonzero fit residual.

The prescribed mode count and bandwidth penalty affect how oscillations are distributed among the modes. Passing the stopping test indicates small changes between iterations; it does not certify a global optimum or a physically unique decomposition. The returned center frequency summarizes a mode’s spectral concentration and is not its instantaneous frequency at every sample.

Instantaneous attributes and Hilbert spectrum

The Hilbert transform associates each real component with an analytic signal. Its magnitude gives amplitude and its angle gives phase. The time derivative of phase in radians, divided by 2\pi, gives instantaneous frequency in Hz when time is in seconds. Applying this operation to a sum of oscillations produces the phase of that sum, not the individual frequencies of its summands; decomposition therefore precedes this part of Hilbert-Huang analysis.

computeInstantaneousAttributes() returns matrices of amplitude, phase and frequency, with samples in rows and modes in columns. Amplitude has the units of the input modes, phase is in radians and frequency is in Hz. The EMD example uses a time step of 9 / 2999 seconds; the expression below reports the median frequency of each mode.

Attributes <- computeInstantaneousAttributes(EMD$modes,
                                             samplingInterval = 9 / 2999)
apply(Attributes$frequency, 2, stats::median)
#> [1] 2.9992761 0.9952460 0.4857283

The analytic signal is computed by a discrete Fourier transform, which implies periodic extension of the record. Phase is continuous within runs of nonzero analytic amplitude; a zero leaves the phase undefined and starts a new run. Frequency uses neighboring phase values and is NA where the required phase values are undefined.

Small analytic amplitude can make the phase derivative sensitive to the signal, and the periodic extension can affect the record boundaries. The resulting frequency is a property of the analyzed component; its physical interpretation also depends on how that component was separated.

buildHilbertSpectrum() places each instantaneous amplitude into the cell containing its frequency, sample by sample. With exponent = 2, it accumulates squared amplitudes. The edges in this example define the intervals [1,2), [2,3), [3,4) and [4,5) in Hz. Frequencies outside these intervals or equal to NA contribute nothing.

Spectrum <- buildHilbertSpectrum(Attributes$frequency, Attributes$amplitude,
                                 edges = c(1, 2, 3, 4, 5), exponent = 2,
                                 samplingInterval = 9 / 2999)
Spectrum$marginal
#> [1] 3.44528895 4.71656719 4.12037251 0.05831411

Spectrum$spectrum holds samples by frequency cells, in squared signal units for this example. Spectrum$marginal sums each cell over time and multiplies by the sampling interval, giving squared signal units times seconds. The values are cell masses; no division by frequency-cell width is applied.

The marginal spectrum accumulates contributions over time and therefore does not retain when an oscillation occurred. It is constructed from instantaneous amplitude and frequency and is not a Fourier spectrum.

Modes and intrinsic mode functions

Data model

Empirical mode decomposition

Variational mode decomposition

Instantaneous attributes and Hilbert spectrum

Further reading