usingNeuroAnalyzerusingCairoMakieCairoMakie.activate!(type="png")usingFFTWusingStatisticsNeuroAnalyzer.progress_bar =false# NEURO object with one channel from a signal vector# (used in the examples below)functionsignal2obj(s::Vector{Float64}, fs::Int64) obj =create_object(; data_type ="eeg") obj.header.recording[:sampling_rate] = fsadd_channel!(obj; data =reshape(s, 1, :, 1), label ="s",type="eeg", unit ="μV")create_time!(obj; fs = fs)return objend# frequency of the largest peak of the amplitude spectrumfunctionpeak_frq(x::AbstractVector, fs::Int64) a =abs.(rfft(x))return (argmax(a) -1) * fs /length(x)end# Makie figure and axis settings matching NeuroAnalyzer plotsfig_size = (900, 450)ax_style = (titlesize =18, xlabelsize =18, ylabelsize =18, xticklabelsize =12, yticklabelsize =12, xticks =LinearTicks(10), xminorticksvisible =true, xminorticks =IntervalsBetween(10), xautolimitmargin = (0, 0), yautolimitmargin = (0.1, 0.1))eeg =import_edf("files/eeg/eeg.edf")keep_channel!(eeg; ch ="eeg")# filtered 0.5-35 HzNeuroAnalyzer.filter!(eeg; ch ="all", fprototype =:butterworth, ftype =:hp, cutoff =0.5, order =4)NeuroAnalyzer.filter!(eeg; ch ="all", fprototype =:fir, ftype =:lp, cutoff =35, bw =10)# 2-second epochse2 =epoch(crop(eeg; seg = (100, 160)); ep_len =2)
Empirical Mode Decomposition (EMD, Huang et al. 1998) splits a signal into a small number of oscillatory components, Intrinsic Mode Functions (IMFs), and a residual (trend):
Unlike the Fourier or wavelet transforms, EMD uses no predefined basis functions (sine waves, wavelets): the components are derived from the local maxima and minima of the signal itself (“empirical”). The IMFs need not be stationary or sinusoidal: their amplitude and frequency can change over time, which makes EMD suitable for non-stationary and non-linear signals such as EEG.
An IMF is an oscillation that:
has numbers of local extrema and zero crossings that are equal or differ by at most one,
is symmetric about zero: the mean of the envelope through its maxima and the envelope through its minima is (close to) zero.
Such a function has a well-defined instantaneous amplitude and frequency, which can be obtained with the Hilbert transform (the Hilbert-Huang transform).
Algorithm
The IMFs are extracted one by one, from the fastest to the slowest oscillation, by sifting:
Find the local maxima and minima of the signal \(h\) (at the start, the signal itself).
Interpolate the maxima and the minima with cubic splines: the upper and lower envelopes.
Subtract the mean of the envelopes \(m\) from \(h\): \(h \leftarrow h - m\).
Repeat steps 1-3 until the change is small: \(SD = \sum m^2 / \sum h^2 < \epsilon\). The result is an IMF.
Subtract the IMF from the signal; the remainder (residual) is decomposed in the same way.
The decomposition stops when the residual has fewer than three maxima or three minima (too few for a reliable envelope: it is treated as the trend), or after max_imf IMFs. The IMFs and the residual sum exactly to the signal.
In emd(), epsilon (default 0.3) sets the stopping criterion of the sifting: smaller values give more sifting iterations and IMFs closer to the IMF definition (narrower-band), but too much sifting removes the amplitude variation of the IMFs. The ends of the signal are a problem for the spline envelopes (there are no extrema beyond them); to reduce end effects, the two outermost extrema at each end are mirrored about the ends of the signal.
IMFs: 2
IMF 1: peak at 20.0 Hz, SD 3.52
IMF 2: peak at 4.0 Hz, SD 7.11
residual - trend: max 0.14
reconstruction error: 3.552713678800501e-15
The two sine waves (SD of a sine wave = amplitude / √2: 3.54 and 7.07) are separated into two IMFs and the trend is the residual.
Mode mixing
EMD separates oscillations by their local time scale, not by frequency bands. When a fast oscillation appears only intermittently, the first IMF contains the fast oscillation where it is present and the slower one elsewhere (mode mixing):
s2 =10.*sin.(2π *4.* t)# 30 Hz bursts at 1-1.5 s and 2.5-3 sburst = (1.<= t .<1.5) .| (2.5.<= t .<3)s2[burst] .+=3.*sin.(2π *30.* t[burst])imf2 =emd(signal2obj(s2, fs); ch ="s", ep =1)fig =Figure(size = fig_size)ax =Axis(fig[1, 1]; ax_style..., title ="IMF 1 of a 4 Hz wave with 30 Hz bursts", xlabel ="Time [s]", ylabel ="Amplitude [μV]")lines!(ax, t, imf2[1, :]; color =:black)fig
Between the bursts, IMF 1 contains the whole 4 Hz wave: the same IMF holds 30 Hz activity in some parts of the signal and 4 Hz activity in others. Variants of EMD reduce mode mixing by adding noise (Ensemble EMD, CEEMDAN); they are not available in NeuroAnalyzer.
EEG
Decomposition of one 2-s epoch of O1:
imf_eeg =emd(e2; ch ="O1", ep =5)plot_imf(imf_eeg; t = e2.epoch_time)
signal: SD 5.9 μV
IMF 1: peak at 21.0 Hz, SD 1.8 μV
IMF 2: peak at 5.0 Hz, SD 2.3 μV
IMF 3: peak at 1.0 Hz, SD 4.6 μV
residual: SD 1.5 μV
The first IMFs contain the fastest activity, the last ones the slowest. Unlike Fourier components, IMFs are not orthogonal: neighboring IMFs can be partly correlated (or anti-correlated, e.g. when a large transient such as an eye blink is split between two IMFs), so their variances need not add up to the variance of the signal. In short epochs, slow oscillations have only a few cycles, so the number of IMFs is small and the slow IMFs are poorly defined (end effects); decompose longer segments when slow activity is of interest.
Tip: IMFs are not frequency bands: their frequency changes over time and differs between epochs and channels, so the “same” IMF number may contain different activity in different epochs. Compare IMFs by their frequency content, not by their number.
Hilbert-Huang Spectrum
The instantaneous frequency and amplitude of each IMF (from its analytic signal) give a time-frequency representation, the Hilbert-Huang spectrum: hhtspectrogram() (or spectrogram(obj; method = :hht)):
