Empirical Mode Decomposition (EMD)

Initialize NeuroAnalyzer 
using NeuroAnalyzer
using Plots
eeg = load("files/eeg.hdf")
NeuroAnalyzer.filter!(eeg,
                      fprototype = :fir,
                      ftype = :bp,
                      cutoff = (5, 25),
                      order = 91,
                      ch = "all")
e2 = epoch(eeg, ep_len = 2)
keep_epoch!(e2, ep = 11:20)
detrend!(e2, type = :mean, ch = "all")

Empirical: The decomposition is data-driven—the signal itself dictates how it is broken down, rather than using predefined bases (e.g., Fourier or wavelet transforms).

Mode: Refers to a portion of the signal that represents a specific oscillation or trend.

Decomposition: The process of separating a signal into simpler, physically meaningful components.

EMD is particularly effective for analyzing non-linear and non-stationary signals, such as EEG, seismic data, or financial time series.

Intrinsic Mode Functions (IMFs)

Definition: IMFs are the components extracted from the signal during EMD. Each IMF represents a simple oscillatory mode embedded in the data.

Criteria for IMFs:

  • The number of local extrema (maxima and minima) and the number of zero-crossings must either be equal or differ by at most one.
  • At any point, the mean value of the envelope defined by the local maxima and the envelope defined by the local minima should be close to zero (symmetry).

Decomposition Process

EMD decomposes a signal \(s\) into a sum of IMFs and a residual:

\[ s = \sum_i IMF_i + \text{residual} \]

IMFs: Ordered from highest to lowest frequency (smallest to largest scales).

Residual: Represents the trend or mean of the signal after all IMFs are extracted.

EMD Algorithm:

  1. Find Local Extrema: Identify all local maxima and minima of the signal \(s\).
  2. Fit Envelopes: Use cubic splines to fit upper (\(E_{up}\)) and lower (\(E_{low}\)) envelopes to the maxima and minima.
  3. Compute Mean Envelope: Calculate the mean of the upper and lower envelopes:

\[ E_{mean} = \frac{E_{up} + E_{low}}{2} \]

  1. Extract IMF: Subtract the mean envelope from the signal:

\[ IMF_1 = s - E_{mean} \]

If \(IMF_1\) satisfies the IMF criteria, it is designated as the first IMF. Otherwise, repeat steps 1–4 with \(IMF_1\) as the new signal.

  1. Update Residual: Subtract the extracted IMF from the signal to obtain the residual:

\[ \text{residual} = s - IMF_1 \]

  1. Repeat: Use the residual as the new signal and repeat the process to extract subsequent IMFs until the residual is a monotonic function or contains no more than one extremum.

Stopping Criteria

The algorithm stops when the sum of differences (SD) between consecutive iterations is smaller than a predefined threshold (\(\epsilon\)).

This ensures that the extracted IMF meets the criteria and that the decomposition process converges.

Applications of EMD: Decompose brain signals into meaningful components for studying oscillations, artifacts, or epileptic spikes.

Advantages of EMD:

  • Adaptive: No need for predefined basis functions; the decomposition is data-driven.
  • Non-linear and Non-stationary: Effective for signals where traditional methods (e.g., Fourier analysis) fail.
  • Physically Meaningful: IMFs represent intrinsic modes of the signal, making interpretation intuitive.

Limitations of EMD:

  • Mode Mixing: IMFs may contain components of different scales, especially in intermittent signals.
  • End Effects: Artifacts can occur at the boundaries of the signal.
  • Computational Cost: Iterative sifting process can be computationally intensive for long signals.
  • Sensitivity to Noise: EMD can be sensitive to noise, requiring preprocessing (e.g., filtering).

Performing and Visualizing EMD

EMD can be performed on one channel and one epoch at a time:

ch = "O2"
ep = 2
imf = emd(e2,
          ch = ch,
          ep = ep,
          epsilon = 0.3)

imf will contain a set a matrix of IMFs (by rows) and the residual at the last row:

res = imf[end, :]

Summing of all IMFs and the residual gives the original signal:

signal_reconstructed = sum(imf; dims = 1)[:]

Plotting IMFs:

plot_imf(imf, t = t)

Tip: plot_imf() will automatically plot the reconstructed signal.