Common Processing Methods

using NeuroAnalyzer
using CairoMakie
CairoMakie.activate!(type = "png")
using Statistics
NeuroAnalyzer.progress_bar = false

# NEURO object with one channel from a signal vector
# (used in the examples below)
function signal2obj(s::Vector{Float64}, fs::Int64)
    obj = create_object(; data_type = "eeg")
    obj.header.recording[:sampling_rate] = fs
    add_channel!(obj;
                 data = reshape(s, 1, :, 1),
                 label = "s",
                 type = "eeg",
                 unit = "μV")
    create_time!(obj; fs = fs)
    return obj
end
# Makie figure and axis settings matching NeuroAnalyzer plots
fig_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))

# raw (unfiltered) EEG recording, 60-second segment
eeg = import_edf("files/eeg/eeg.edf")
eeg = crop(eeg; seg = (100, 160))

This tutorial covers simple operations on the signal amplitude: removing the DC offset, detrending, normalization and standardization, and the derivative of the signal. All functions work on the selected channels (ch), separately in each epoch, and have in-place (!) versions.

Removing the DC Offset

The DC offset is a constant value added to the signal: the mean of the signal, which is not related to brain activity. It comes from electrode potentials (the skin-electrode interface acts as a small battery) and from the amplifier. DC-coupled amplifiers record it fully, and it can be much larger than the EEG itself:

m = mean(eeg.data[1:4, :, 1]; dims = 2)
for (l, v) in zip(labels(eeg)[1:4], m)
    println(rpad(l, 4), ": mean ", round(v; digits = 1), " μV")
end

Fp1 : mean 562.3 μV
Fp2 : mean 728.4 μV
F3  : mean 689.0 μV
F4  : mean 646.0 μV

remove_dc() subtracts the mean of each channel:

eeg_dc = remove_dc(eeg; ch = "all")
m = mean(eeg_dc.data[1:4, :, 1]; dims = 2)
println("means after remove_dc(): ", round.(vec(m); digits = 6) .+ 0.0)

means after remove_dc(): [0.0, 0.0, 0.0, 0.0]

With n, the mean of a part of the signal is subtracted instead: of the first n samples, or of samples n[1] to n[2] (a tuple). For epochs, this is a baseline correction: the mean of the pre-stimulus interval is subtracted from each epoch:

faces = import_recording("files/eeg/s1_faces.vhdr")
# epochs from 200 ms before to 800 ms after the s71 stimuli
ep71 = epoch(faces; marker = "s71", offset = 0.2, ep_len = 1.0)
# baseline: the first 200 ms
n_bl = round(Int, 0.2 * sr(ep71))
ep71_bl = remove_dc(ep71; ch = "eeg", n = n_bl)

bl = mean(ep71_bl.data[1, 1:n_bl, :]; dims = 1)
println("pre-stimulus mean, first 3 epochs: ",
        round.(vec(bl)[1:3]; digits = 6) .+ 0.0, " μV")

pre-stimulus mean, first 3 epochs: [0.0, 0.0, 0.0] μV

Tip: A high-pass filter also removes the DC offset (see Filtering); remove_dc() removes only the constant and does not change the shape of the signal.

Detrending

A trend is a slow change of the signal over the whole recording or epoch, e.g. an electrode drift caused by sweating, electrode polarization or changes of temperature. detrend() fits a trend to each channel (in each epoch) and subtracts it:

type	Subtracted trend
:linear	a straight line (default)
:ls	a straight line, least-squares fit (closed-form solution, faster)
:poly	a polynomial of degree order
:loess	a Loess smooth (local regression) with span f ∈ (0, 1]
:mean	the mean (as remove_dc())
:constant	the value offset (default 0)

A test signal: a 10 Hz oscillation with a curved drift:

fs = 256
t = collect(0:(1 / fs):(10 - 1 / fs))
drift = 20 .+ 4 .* t .- 0.5 .* t .^ 2
s = 5 .* sin.(2π * 10 .* t) .+ drift
obj = signal2obj(s, fs)
s_lin = detrend(obj; ch = "s", type = :linear).data[1, :, 1]
s_poly = detrend(obj; ch = "s", type = :poly, order = 2).data[1, :, 1]
s_loess = detrend(obj; ch = "s", type = :loess, f = 0.3).data[1, :, 1]
fig = Figure(size = fig_size)
ax = Axis(fig[1, 1];
          ax_style...,
          title = "Detrending",
          xlabel = "Time [s]",
          ylabel = "Amplitude [μV]")
lines!(ax, t, s; color = (:gray, 0.5), label = "signal")
lines!(ax, t, s_lin; label = ":linear")
lines!(ax, t, s_poly; label = ":poly (order 2)")
lines!(ax, t, s_loess; label = ":loess (f = 0.3)")
axislegend(ax; position = :rb)
fig

for (name, x) in [(":linear", s_lin), (":poly", s_poly), (":loess", s_loess)]
    # mean of each 1-s window: the remaining slow trend
    w = mean(reshape(x, fs, :); dims = 1)
    println(rpad(name, 8), ": residual trend ",
            round(maximum(abs, w); digits = 2), " μV")
end

:linear : residual trend 6.04 μV
:poly   : residual trend 0.04 μV
:loess  : residual trend 0.11 μV

The linear trend does not fit the curved drift; the second-order polynomial removes it, and so does the Loess smooth, which can follow drifts of any shape.

Tip: Choose the trend model by the expected shape of the drift. A polynomial of high order or a Loess smooth with a small span also follows slow brain activity and removes it; the span f is a fraction of the signal, so f = 0.3 of a 10-s signal follows changes slower than about 3 s.

Detrending or high-pass filtering? Detrending removes one fitted curve and does not change faster activity; a high-pass filter removes all frequencies below the cutoff, including slow brain activity, and can create filter artifacts. For short epochs, detrending (linear or :mean) is often enough; long continuous recordings with irregular drifts need a high-pass filter.

Normalization

Normalization changes the scale of the signal, e.g. to compare channels, subjects or recordings with different amplitudes, or to prepare features for machine learning. normalize() has many methods (method):

Method	Result
:zscore	\((x - \bar{x}) / SD\): mean 0, standard deviation 1
:mad	\((x - \tilde{x}) / (1.4826 \cdot MAD)\): robust z-score (median, MAD)
:minmax	scaled to \([-n, +n]\) (default \([-1, +1]\))
:n	scaled to \([0, n]\) (default \([0, 1]\))
:perc	scaled to \([0, 1]\)
:pos	shifted so that the minimum is 0
:neg	shifted so that the maximum is 0
:log, :log10	\(\log(1 + x + \vert\min(x)\vert)\) (natural or base 10)
:neglog, :neglog10	negative logarithm
:invroot	\(1 / \sqrt{x}\) (for positive values, e.g. power)
:gauss	rank-based transform to the standard normal distribution
:rank	ranks of the values
:softmax	\(e^{x_i} / \sum e^{x_j}\)
:sigmoid	\(1 / (1 + e^{-x})\)
:fisher	Fisher z-transform (\(\text{atanh}\)) of the signal scaled to \([-1, +1]\)
:none	no change

All channels together or each channel separately

By default (bych = false), the selected channels (and epochs) are normalized together, with one mean and one standard deviation: differences between channels are kept. With bych = true, each channel (in each epoch) is normalized separately:

eeg_f = NeuroAnalyzer.filter(eeg;
                             ch = "eeg",
                             fprototype = :butterworth,
                             ftype = :hp,
                             cutoff = 0.5,
                             order = 4)
z_all = NeuroAnalyzer.normalize(eeg_f; ch = "eeg", method = :zscore)
z_ch = NeuroAnalyzer.normalize(eeg_f;
                               ch = "eeg",
                               method = :zscore,
                               bych = true)

idx = get_channel(eeg_f; ch = ["Fp1", "Fz", "O1"])
println(rpad("", 6), rpad("bych = false", 22), "bych = true")
for i in idx
    a = z_all.data[i, :, 1]
    b = z_ch.data[i, :, 1]
    println(rpad(labels(eeg_f)[i], 6),
            rpad(string("mean ", round(mean(a); digits = 2) + 0.0,
                        ", SD ", round(std(a); digits = 2)), 22),
            "mean ", round(mean(b); digits = 2) + 0.0,
            ", SD ", round(std(b); digits = 2))
end

      bych = false          bych = true
Fp1   mean 0.0, SD 0.72     mean 0.0, SD 1.0
O1    mean 0.0, SD 1.34     mean 0.0, SD 1.0
Fz    mean 0.0, SD 0.06     mean 0.0, SD 1.0

With bych = false, Fz keeps its smaller amplitude relative to the other channels; with bych = true, all channels have SD 1.

Scaling to a range:

x = NeuroAnalyzer.normalize(eeg_f; ch = "Fp1", method = :minmax).data
println(":minmax:        ", round.(extrema(x[1, :, 1]); digits = 2))
x = NeuroAnalyzer.normalize(eeg_f; ch = "Fp1", method = :n, n = 10).data
println(":n with n = 10: ", round.(extrema(x[1, :, 1]); digits = 2))

:minmax:        (-1.0, 1.0)
:n with n = 10: (0.0, 10.0)

Tip: The normalization functions are also available for vectors and arrays (normalize_zscore(), normalize_minmax(), normalize_mad() …), e.g. to normalize computed measures such as band powers.

Standardization

standardize() z-scores each channel in each epoch (as normalize(…; method = :zscore, bych = true)), and also returns the fitted transformations (ZScoreTransform of StatsBase.jl, one per epoch), which store the means and standard deviations, so that the original scale can be restored later:

eeg_std, scaler = standardize(eeg_f; ch = "eeg")
println("transformations: ", length(scaler))
# restore the original values
ch = get_channel(eeg_f; ch = "eeg")
x = NeuroAnalyzer.StatsBase.reconstruct(scaler[1], eeg_std.data[ch, :, 1])
println("maximum difference after restoring: ",
        round(maximum(abs.(x .- eeg_f.data[ch, :, 1])); sigdigits = 2), " μV")

transformations: 1
maximum difference after restoring: 1.1e-13 μV

Power normalization

normpower() divides the signal by its root-mean-square value (RMS), so that the power of the result is 1:

x = normpower(eeg_f; ch = "eeg").data[ch, :, 1]
println("RMS after normpower(): ", round(sqrt(mean(x .^ 2)); digits = 3))

RMS after normpower(): 1.0

Derivative

The derivative measures how fast the signal changes. It emphasizes fast changes (spikes, sharp waves, steps, electrode pops) and attenuates slow activity: it acts as a high-pass filter whose gain increases with frequency. derivative() uses the symmetric difference quotient:

\[ x'[n] = \frac{x[n+1] - x[n-1]}{2} \]

(at the first and last sample, the one-sided difference is used, so the result has the same length as the signal). The result is the change per sample; multiply it by the sampling rate to get the change per second (μV/s).

A slow 2 Hz wave with a short spike:

fs = 256
t = collect(0:(1 / fs):(2 - 1 / fs))
s = 30 .* sin.(2π * 2 .* t)
# spike: Gaussian wave of 20 μV, 20 ms wide
s .+= 20 .* exp.(-(t .- 1.3) .^ 2 ./ (2 * 0.01^2))
obj = signal2obj(s, fs)
ds = derivative(obj; ch = "s").data[1, :, 1] .* fs
fig = Figure(size = fig_size)
ax1 = Axis(fig[1, 1];
           ax_style...,
           title = "Signal",
           ylabel = "[μV]")
lines!(ax1, t, s; color = :black)
ax2 = Axis(fig[2, 1];
           ax_style...,
           title = "Derivative",
           xlabel = "Time [s]",
           ylabel = "[μV/s]")
lines!(ax2, t, ds; color = :black)
linkxaxes!(ax1, ax2)
fig

The spike is smaller than the slow wave in the signal, but dominates the derivative. This is used e.g. to detect spikes and artifacts.