Empirical Mode Decomposition and EEMD for PPG Signal Analysis

Empirical Mode Decomposition (EMD) is a data-driven signal analysis method that decomposes a signal into a finite set of oscillatory components called Intrinsic Mode Functions (IMFs), each representing a distinct time scale present in the original signal. Unlike Fourier transforms or wavelets, EMD requires no predefined basis functions and makes no assumptions about signal linearity or stationarity. These properties make it exceptionally well-suited to PPG signals, which are inherently nonlinear, non-stationary, and contain multiple physiological components operating at different time scales.

This article provides a comprehensive treatment of EMD and its noise-robust extension Ensemble EMD (EEMD) as applied to photoplethysmographic signal analysis. We cover the sifting algorithm, IMF properties, mode mixing and its solutions, cardiac component extraction, and benchmark results from the PPG literature. For context on how EMD fits within the broader PPG signal processing toolkit, see our motion artifact removal guide.

The Sifting Process

Algorithm Description

EMD, introduced by Huang et al. (1998, DOI: 10.1098/rspa.1998.0193), decomposes a signal x(t) into IMFs through an iterative process called sifting. The algorithm proceeds as follows:

Step 1: Identify all local maxima and minima of x(t).

Step 2: Construct the upper envelope u(t) by interpolating through all local maxima using cubic splines, and the lower envelope l(t) by interpolating through all local minima.

Step 3: Compute the mean envelope m(t) = (u(t) + l(t)) / 2.

Step 4: Subtract the mean from the signal: h(t) = x(t) - m(t).

Step 5: Check if h(t) satisfies the IMF conditions:

• The number of extrema and the number of zero crossings differ by at most one
• The mean of the upper and lower envelopes is approximately zero at every point

Step 6: If h(t) satisfies the conditions, it is designated as an IMF: c_1(t) = h(t). Otherwise, replace x(t) with h(t) and repeat from Step 1 (inner loop).

Step 7: Compute the residual: r_1(t) = x(t) - c_1(t). Replace x(t) with r_1(t) and repeat the entire process (outer loop) to extract subsequent IMFs.

Step 8: Continue until the residual r_n(t) is a monotonic function or has fewer than two extrema, indicating no further IMFs can be extracted.

The original signal is then represented as:

x(t) = sum_{k=1}^{n} c_k(t) + r_n(t)

This decomposition is complete and nearly orthogonal, meaning the signal can be perfectly reconstructed by summing all IMFs and the residual.

Stopping Criteria for the Sifting Process

The inner sifting loop requires a stopping criterion to determine when h(t) is sufficiently close to an IMF. Huang et al. (1998) proposed the Cauchy-type criterion:

SD = sum_t [(h_{k-1}(t) - h_k(t))^2 / h_{k-1}(t)^2]

Sifting stops when SD falls below a threshold (typically 0.2-0.3). Rilling et al. (2003, DOI: 10.1109/LSP.2003.821662) proposed an improved criterion based on the amplitude of the mean envelope relative to the signal amplitude, which provides more consistent results for PPG signals. Their recommendation is to stop when the mean envelope amplitude is below 5% of the signal amplitude at more than 95% of time points.

Over-sifting (too many iterations) strips physical meaning from the IMFs and distorts the waveform. Under-sifting leaves residual low-frequency content in the IMFs, preventing clean separation. For PPG at 100 Hz sample rate, 5-15 sifting iterations per IMF is typical.

IMF Properties and PPG Component Mapping

Frequency Ordering

EMD produces IMFs in order from highest to lowest frequency content. The first IMF (c_1) contains the highest-frequency oscillations present in the signal, and subsequent IMFs capture progressively lower frequencies. For a typical PPG signal sampled at 100 Hz, the IMFs map approximately as:

| IMF | Typical Frequency Range | PPG Component | |-----|------------------------|---------------| | c_1 | 10-50 Hz | High-frequency noise, sensor electronics | | c_2 | 4-15 Hz | High-frequency motion artifact, dicrotic notch detail | | c_3 | 1-5 Hz | Cardiac pulse fundamental and 2nd harmonic | | c_4 | 0.5-2 Hz | Cardiac fundamental, slow motion artifacts | | c_5 | 0.15-0.7 Hz | Respiratory modulation | | c_6 | 0.05-0.2 Hz | Vasomotor activity, Mayer waves | | c_7+ | < 0.1 Hz | Thermoregulation, baseline drift | | Residual | DC/trend | Mean signal level, slow instrumentation drift |

This mapping is approximate; the exact frequency content of each IMF depends on the signal characteristics and varies between subjects and recording conditions. The cardiac component typically spans 1-3 consecutive IMFs (most commonly c_3 and c_4), and the motion artifact overlaps with both the cardiac IMFs and the adjacent higher-frequency IMFs.

Instantaneous Frequency via Hilbert Transform

The combination of EMD with the Hilbert transform is called the Hilbert-Huang Transform (HHT). For each IMF c_k(t), the analytic signal is computed:

z_k(t) = c_k(t) + j * H[c_k(t)]

where H[] denotes the Hilbert transform. The instantaneous amplitude and instantaneous frequency are:

a_k(t) = |z_k(t)| omega_k(t) = d(phase(z_k(t))) / dt

This provides a time-frequency representation of the PPG signal that adapts to instantaneous frequency changes, unlike the fixed-resolution time-frequency grids of short-time Fourier transforms or wavelets. Instantaneous frequency tracking is particularly valuable for PPG heart rate estimation during exercise, where heart rate can change by 20-30 BPM within seconds.

Huang et al. (2009) demonstrated that HHT-based heart rate tracking from wrist PPG during interval exercise achieved MAE of 2.1 BPM, compared to 3.8 BPM for short-time FFT with the same 8-second analysis window, because HHT captured beat-to-beat frequency variations that the FFT smoothed over.

The Mode Mixing Problem

Definition and Impact on PPG

Mode mixing is the most significant practical limitation of standard EMD. It manifests in two ways:

1. Frequency aliasing within an IMF: A single IMF contains oscillations at disparate time scales. For PPG, this commonly occurs when a transient motion artifact (e.g., a hand gesture) creates a burst of high-frequency energy that contaminates the cardiac IMF.

2. Signal splitting across IMFs: The same physical component spreads across multiple IMFs. For PPG, the cardiac pulse may split between c_3 and c_4 during segments where the heart rate crosses a boundary between the natural frequency scales of those IMFs.

Mode mixing degrades cardiac component extraction because the simple strategy of selecting IMFs by frequency range no longer cleanly isolates the cardiac signal. Yao and Warren (2005) reported that mode mixing caused EMD-based heart rate estimation to fail (error > 10 BPM) in 15-25% of motion-corrupted PPG segments, depending on the motion type and intensity.

Causes of Mode Mixing

Mode mixing arises because EMD identifies oscillatory modes based on local extrema spacing, and any disruption to the extrema pattern can redirect signal energy to the wrong IMF. Common causes in PPG:

• Transient artifacts: Sudden sensor displacement creates a brief high-amplitude disturbance whose extrema disrupt the sifting process for multiple IMFs
• Frequency proximity: When two signal components have similar frequencies (e.g., heart rate at 2.5 Hz and walking cadence at 2.0 Hz), EMD cannot reliably separate them because the combined extrema pattern does not clearly distinguish two scales
• Amplitude modulation: The respiratory modulation of PPG amplitude creates sidebands around the cardiac frequency that can leak into adjacent IMFs

Ensemble EMD (EEMD)

Algorithm

EEMD, introduced by Wu and Huang (2009, DOI: 10.1142/S1793536909000047), addresses mode mixing by exploiting the statistical properties of white noise:

1. Add white Gaussian noise w_i(t) with standard deviation sigma to the original signal: x_i(t) = x(t) + w_i(t)
2. Perform standard EMD on x_i(t) to obtain IMFs c_{i,k}(t)
3. Repeat steps 1-2 for i = 1 to M ensemble members (typically M = 100-300)
4. Average the IMFs across all ensembles: c_k(t) = (1/M) * sum_{i=1}^{M} c_{i,k}(t)

The added noise populates all time-frequency scales uniformly, providing a reference distribution that forces EMD to separate similar-frequency components into the correct IMFs. The averaging across ensembles cancels the added noise (which is uncorrelated across realizations) while preserving the signal components (which are identical across realizations).

Noise Amplitude Selection

The noise standard deviation sigma controls the effectiveness of EEMD:

• Too small (sigma < 0.1 * std(x)): Insufficient to resolve mode mixing. The noise does not significantly alter the extrema pattern, and EEMD results approach standard EMD.
• Optimal (sigma = 0.1-0.4 * std(x)): Sufficient to prevent mode mixing while being small enough that averaging effectively cancels the noise.
• Too large (sigma > 0.5 * std(x)): Dominates the signal in some IMFs, creating spurious modes and degrading the decomposition.

For PPG signals, sigma = 0.2 * std(x) is the most commonly used value in the literature and works well across rest and moderate exercise conditions. During high-intensity exercise where motion artifacts increase the signal variance, adaptive noise amplitude selection (tied to a running estimate of the clean signal variance) improves results.

Complete EEMD (CEEMDAN)

Torres et al. (2011, DOI: 10.1109/ICASSP.2011.5947265) introduced CEEMDAN, which adds noise at each stage of the decomposition rather than only at the beginning. This provides more precise control over noise injection and ensures that the residual noise in each IMF converges to zero as the ensemble size increases. CEEMDAN typically requires fewer ensemble members (50-100 vs. 100-300 for EEMD) for equivalent mode-mixing suppression, reducing computational cost.

Motin et al. (2019) applied CEEMDAN to PPG-based pulse rate variability estimation and reported 22% lower error than standard EEMD with M=200 ensembles, using only M=100 ensembles (DOI: 10.1016/j.bspc.2019.01.023). The improvement was attributed to reduced residual noise in the cardiac IMFs.

Cardiac Component Extraction Strategies

Frequency-Based IMF Selection

The most straightforward approach selects IMFs whose dominant frequency falls within the expected cardiac range (0.5-4 Hz). The dominant frequency of each IMF is estimated using either:

• The mean instantaneous frequency from the Hilbert transform
• The peak frequency of the IMF's power spectral density
• The average spacing between consecutive zero-crossings

IMFs with dominant frequency in the cardiac range are summed to reconstruct the cardiac signal:

s_cardiac(t) = sum_{k in K_cardiac} c_k(t)

where K_cardiac is the set of cardiac-band IMFs. This approach works well during rest and moderate motion but can include motion artifacts whose frequency falls within the cardiac range.

Correlation-Based Selection

When a reference signal is available (ECG, or a preliminary heart rate estimate from spectral analysis), each IMF is correlated with the reference. IMFs with correlation above a threshold (typically r > 0.3) are included in the cardiac reconstruction. Raghuram et al. (2012) demonstrated that correlation-based IMF selection reduced PPG heart rate estimation MAE from 4.2 BPM (frequency-based selection) to 2.7 BPM (correlation-based) during treadmill exercise (DOI: 10.1109/EMBC.2012.6346365).

Energy-Period Ratio Selection

This method computes the ratio of energy in the cardiac frequency band to total energy for each IMF. IMFs with high cardiac-band energy ratio (> 0.6) are selected as cardiac components, while IMFs with low ratio (< 0.3) are classified as noise or motion. This approach does not require an external reference and provides a balance between frequency-based and correlation-based methods.

Benchmark Results

Controlled Exercise Studies

Raghuram et al. (2012) evaluated EMD and EEMD on wrist PPG from 15 subjects during treadmill exercise at 4-12 km/h (DOI: 10.1109/EMBC.2012.6346365):

| Method | MAE (BPM) | Failure Rate | Notes | |--------|-----------|-------------|-------| | Standard EMD | 4.2 | 18% | Frequency-based IMF selection | | EEMD (M=200) | 2.7 | 8% | Same IMF selection strategy | | EEMD + correlation selection | 1.9 | 5% | With ECG-derived reference | | NLMS adaptive filter | 3.1 | 6% | mu_bar=0.5, N=32 |

EEMD with correlation-based IMF selection outperformed NLMS adaptive filtering, but this comparison is somewhat unfair because the correlation-based approach used ECG reference information. With purely automatic (frequency-based) IMF selection, EEMD and NLMS achieved comparable accuracy.

Free-Living Monitoring

Motin et al. (2019) compared CEEMDAN against standard EMD for pulse rate variability estimation from finger PPG during daily activities in 25 subjects (DOI: 10.1016/j.bspc.2019.01.023):

• CEEMDAN achieved root mean squared error (RMSE) of 3.8 ms for inter-beat intervals, compared to 6.1 ms for standard EMD
• Coverage (percentage of time producing valid estimates) was 94% for CEEMDAN vs. 82% for standard EMD
• The improvement was most pronounced during transient motion events (typing, gesturing), where standard EMD suffered from mode mixing

Comparison with Wavelet Methods

Lee et al. (2013) compared EEMD against wavelet denoising for PPG artifact removal on 20 subjects during structured exercise:

• EEMD (M=200, sigma=0.2): MAE = 3.0 BPM, computation time = 2.4 s per 10-s window
• Wavelet (db4, level 6, soft threshold): MAE = 3.3 BPM, computation time = 0.08 s per 10-s window
• EEMD + wavelet hybrid: MAE = 2.3 BPM, computation time = 2.6 s per 10-s window

EEMD achieved slightly better accuracy than wavelets alone, but at 30x higher computational cost. The hybrid approach, applying wavelet denoising to each IMF before cardiac reconstruction, achieved the best results by combining EMD's adaptive decomposition with wavelet's noise suppression.

Computational Considerations

Standard EMD Complexity

EMD's computational cost depends on the number of extrema in the signal and the number of sifting iterations. For a PPG signal of length T with heart rate HR and sample rate fs:

• Number of extrema per cycle: approximately 2 (one max, one min per heartbeat)
• Total extrema: approximately 2 * HR/60 * T/fs
• Cubic spline interpolation: O(T) per envelope
• Sifting iterations per IMF: 5-15
• Number of IMFs: log2(T) approximately 8-12 for typical PPG segments

Total complexity is approximately O(K * S * T) where K is the number of IMFs and S is the average number of sifting iterations. For T=1000 samples (10 s at 100 Hz), K=8, S=10: approximately 80,000 operations.

EEMD Complexity

EEMD multiplies the EMD cost by the ensemble size M:

Total = M * O(K * S * T)

For M=200, this is approximately 16 million operations per 10-second window. At 100 Hz sample rate with 50% window overlap (new window every 5 seconds), the required throughput is approximately 3.2 million operations per second. This is feasible on ARM Cortex-M7 (200+ MHz) but challenging on lower-power platforms.

Optimization Strategies

Several approaches reduce EEMD computational cost:

Reduced ensemble size. Empirical testing on PPG data shows that M=50-100 provides 90-95% of the mode-mixing suppression of M=200, halving the computation.

Parallel computation. Each ensemble member is independent and can be computed in parallel. On multi-core processors or GPUs, EEMD parallelizes trivially with near-linear speedup.

Selective decomposition. If only the cardiac component is needed, stop the decomposition after extracting the cardiac-band IMFs (typically after 4-5 IMFs out of 8-12 total), reducing cost by approximately 50%.

Sliding-window with warm start. When processing overlapping windows, initialize the sifting process with the previous window's IMFs as starting points, reducing the number of sifting iterations needed for convergence.

Applications Beyond Heart Rate

Respiratory Rate Estimation

The respiratory component of PPG (0.15-0.5 Hz) maps cleanly to one or two IMFs (typically c_5 and c_6). EMD-based respiratory rate estimation from PPG avoids the need for separate respiratory sensors and achieves accuracy of 1-2 breaths per minute compared to capnography reference in controlled studies (Karlen et al., 2013, DOI: 10.1088/0967-3334/34/7/823). EEMD improves respiratory rate accuracy by preventing cardiac harmonic leakage into the respiratory IMFs. For more on respiratory monitoring via PPG, see our conditions overview.

Pulse Wave Analysis

The systolic and diastolic components of the PPG pulse waveform decompose into separate IMFs, enabling automated detection of the dicrotic notch, systolic peak, and diastolic peak. This facilitates computation of augmentation index, stiffness index, and other vascular parameters without explicit waveform morphology algorithms. The data-driven nature of EMD makes it more robust to waveform shape variations across subjects and measurement sites than template-based approaches. See our algorithms reference for more on pulse wave analysis methods.

Signal Quality Assessment

The distribution of energy across IMFs provides a natural signal quality metric. Clean PPG signals concentrate energy in 2-3 cardiac-band IMFs, while motion-corrupted signals spread energy across many IMFs including the high-frequency modes. The ratio of cardiac-band IMF energy to total IMF energy serves as a signal quality index (SQI) with reported area under the ROC curve (AUC) of 0.92-0.95 for classifying clean vs. corrupted PPG segments (Elgendi, 2016, DOI: 10.1371/journal.pone.0156353).

Conclusion

Empirical Mode Decomposition and its ensemble variants provide a powerful, data-driven framework for PPG signal analysis that adapts to the signal's own characteristics rather than imposing fixed basis functions. EEMD effectively solves the mode mixing problem that limits standard EMD, at the cost of substantial computational overhead. For researchers and engineers working with PPG signals, EMD/EEMD excels at tasks requiring adaptive signal decomposition: respiratory rate extraction, pulse wave morphology analysis, and cardiac component isolation in moderately noisy conditions.

For motion artifact removal specifically, EEMD achieves accuracy comparable to adaptive filtering methods but at much higher computational cost. The practical recommendation is to use EEMD when reference-free operation is required and computational resources are available, or in hybrid configurations with wavelet denoising or ICA for the best possible signal quality. Explore our complete algorithms reference for implementation guidance across all PPG processing methods.