Rejecting Bad Channels and Epochs

Initialize NeuroAnalyzer 
using NeuroAnalyzer
using Random
NeuroAnalyzer.progress_bar = false
eeg = import_edf("files/eeg/eeg.edf")
keep_channel!(eeg; ch = "eeg")
# filtered 0.5-35 Hz (removes the DC offsets and the line noise of the raw
# recording)
NeuroAnalyzer.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)
load_locs!(eeg; file_name = "files/locs/standard-10-20-cap19-elmiko.locs")
# 60-second segment
eeg = crop(eeg; seg = (100, 160))

A bad channel records mostly noise for the whole recording, e.g. because the electrode is disconnected, has a high impedance or is bridged with its neighbor. A bad epoch contains a transient artifact (eye blink, movement, muscle burst, electrode pop) in one or more channels. Both must be dealt with before analysis: one noisy channel distorts the average reference, spectra, connectivity and topographic maps, and one large artifact dominates an averaged ERP.

The usual order is:

  1. filter the signal (at least a high-pass filter, so that DC offsets and slow drifts do not mask everything else),
  2. detect bad channels and interpolate or delete them,
  3. re-reference (e.g. to the average reference, which needs clean channels),
  4. split the signal into epochs, detect bad epochs and delete them.

NeuroAnalyzer detects both automatically with channel_reject() and epoch_reject(). The detection should always be checked visually: automatic methods find typical faults, but they cannot tell an artifact from unusual brain activity.

How Detection Works

Both functions compute a statistic (e.g. variance) of each channel in each epoch and flag the values that are abnormal. Most methods do not use fixed thresholds but compare the values with each other using robust z-scores:

\[ z = \frac{x - \mathrm{median}(x)}{1.4826 \cdot \mathrm{MAD}(x)} \]

where MAD is the median absolute deviation. Unlike the mean and the standard deviation, the median and the MAD are not affected by a few extreme values, so a single bad channel cannot hide itself by inflating the spread. For normally distributed data, the factor 1.4826 makes the robust z-score equal to the usual z-score. Values with \(z > 3.5\) (the default z argument) are flagged.

What is compared depends on the function:

  • channel_reject() compares channels with each other, separately in each epoch; a channel is bad if it is flagged in at least a fraction tc (default 0.2) of the epochs,
  • epoch_reject() compares epochs with each other, separately for each channel (so channels with naturally larger signals, e.g. frontal or temporal, are not penalized); an epoch is bad if at least nbad (default 1) channels are flagged in it.

Detection Methods

Method Statistic Detects
:flat fraction of 1-s windows with SD < flat_tol (0.1 μV) disconnected, saturated channels
:amp maximum absolute amplitude > amp_t (400 μV) very large artifacts
:var variance (both too high and too low) noisy, flat, bridged channels
:p2p largest peak-to-peak amplitude in 1-s windows jumps, pops, large artifacts
:tkeo mean Teager-Kaiser energy high-frequency activity (muscle, noise)
:kurt kurtosis spiky signals (pops, spikes)
:z fraction of samples with robust z-score > z pops, short bursts of artifacts
:rmse RMS distance from the median signal of all channels channels different from the others
:rmsd as :rmse, ignoring a constant offset channels different from the others
:euclid as :rmse, Euclidean distance channels different from the others
:ransac correlation with the signal predicted from 3 neighbors channels not matching their neighbors

:flat and :amp use absolute thresholds (in signal units, e.g. μV). :ransac flags channels whose correlation with their neighbors is abnormally low (robust z-score < -z) and below ransac_r (default 0.6); it needs channel locations. All other methods use robust z-scores and need at least 3 channels (channel_reject()) or 3 epochs (epoch_reject()). By default, all methods are used; a list of methods can be given with the method argument:

e2 = epoch(eeg; ep_len = 2)
bad = channel_reject(e2; ch = "eeg", method = [:flat, :var, :p2p])

Tip: The window length used by :flat, :p2p and :ransac is set by w (in samples, default: 1 second). It should be shorter than the epoch, otherwise each epoch is a single window.

Bad Channels

To show what each method detects, four faults are added to a copy of the recording:

  • F3: flat (disconnected electrode),
  • C4: strong noise (high impedance),
  • P3: 15 electrode pops (sudden 300 μV steps lasting 6 samples),
  • O2: slow drift (0 to 400 μV over 60 s).
eeg_bad = deepcopy(eeg)
Random.seed!(1)
n = epoch_len(eeg_bad)
ix(l) = get_channel(eeg_bad; ch = l)[1]
eeg_bad.data[ix("F3"), :, 1] .= 0.0
eeg_bad.data[ix("C4"), :, 1] .+= 40 .* randn(n)
for s in rand(1:(n - 6), 15)
    eeg_bad.data[ix("P3"), s:(s + 5), 1] .+= 300
end
eeg_bad.data[ix("O2"), :, 1] .+= range(0, 400; length = n)
NeuroAnalyzer.plot(eeg_bad;
                   ch = ["Fp1", "F3", "Cz", "C4", "P3", "O2"],
                   seg = (0, 10),
                   gui = false)

Each method applied separately to the clean and to the faulty recording, split into 2-second epochs:

e2 = epoch(eeg; ep_len = 2)
e2_bad = epoch(eeg_bad; ep_len = 2)
methods = [:flat, :amp, :var, :p2p, :tkeo, :kurt, :z,
           :rmse, :rmsd, :euclid, :ransac]
println(rpad("method", 9), rpad("clean", 26), "faulty")
for m in methods
    b0 = labels(e2)[channel_reject(e2; ch = "eeg", method = m)]
    b1 = labels(e2_bad)[channel_reject(e2_bad; ch = "eeg", method = m)]
    println(rpad(m, 9), rpad(join(b0, " "), 26), join(b1, " "))
end
method   clean                     faulty
flat                               F3
amp                                
var      F4 Fz                     F3 F4 C4 P3 Fz
p2p                                C4 P3
tkeo                               C4 P3 O2
kurt                               P3
z                                  P3
rmse                               C4 P3 O2
rmsd     Fp2                       C4 P3
euclid                             C4 P3 O2
ransac   Fp2 F4 F8 Fz              Fp2 F3 C4 Cz

Summary of the faulty recording:

Fault Flagged by
F3 (flat) :flat, :var, :ransac
C4 (noise) :var, :p2p, :tkeo, :rmse, :rmsd, :euclid, :ransac
P3 (pops) :var, :p2p, :tkeo, :kurt, :z, :rmse, :rmsd, :euclid
O2 (drift) :tkeo, :rmse, :euclid

Notes:

  • :amp finds nothing: no channel exceeds 400 μV (the pops are 300 μV); the threshold must be adjusted to the data (e.g. amp_t = 200).
  • :var also flags F4 and Fz in the clean recording: their variance is genuinely much lower than in the other channels; low variance is suspicious (e.g. bridged electrodes), but here it may as well be normal.
  • :ransac flags Fp2, F4, F8 and Fz in the clean recording, whose correlation with their neighbors is indeed low; with faults, it also flags Cz, a neighbor of the bad channels F3 and C4.
  • :rmsd ignores slow offsets, so it misses the drift in O2.

No method detects all faults, and some flag clean channels; the methods complement each other and their results should be combined and checked visually.

Marking Bad Channels

channel_reject() returns a mask with one entry per selected channel:

bad = channel_reject(e2_bad; ch = "eeg", method = [:flat, :p2p, :tkeo])
println("bad channels: ", labels(e2_bad)[bad])
bad channels: ["F3", "C4", "P3", "O2"]

channel_reject!() writes the result into the object header (obj.header.recording[:bad_channel], one entry per channel), where it is kept for the following steps:

channel_reject!(e2_bad; ch = "eeg", method = [:flat, :p2p, :tkeo])
bc = e2_bad.header.recording[:bad_channel]
println("marked as bad: ", labels(e2_bad)[bc])
marked as bad: ["F3", "C4", "P3", "O2"]

Channels can also be marked manually, e.g. after visual inspection:

bc[ix("O2")] = true
println("marked as bad: ", labels(e2_bad)[bc])
marked as bad: ["F3", "C4", "P3", "O2"]

Tip: In the interactive signal plot (plot() with gui = true), bad channels are drawn in light gray; a right-click on the channel label marks or unmarks it as bad.

What to Do with Bad Channels

Marking a channel as bad does not change its signal. The bad channels can be:

  • interpolated from the neighboring channels (plinterpolate_channel(), lrinterpolate_channel(), mlinterpolate_channel()), see Interpolating bad channels; this keeps the number of channels (needed e.g. for averaging across subjects or for topographic maps),
  • deleted with delete_channel(),
  • excluded from analyses: with the exclude_bads preference set to true, channels marked as bad are left out whenever channels are selected with the ch argument (e.g. ch = "eeg"); see Preferences.
e2_clean = delete_channel(e2_bad; ch = labels(e2_bad)[bc])

Bad Epochs

Four artifacts are added to the clean recording split into 2-second epochs:

  • epoch 5: eye blink in Fp1 and Fp2 (150 μV),
  • epoch 12: muscle burst in T3 (0.8 s of noise),
  • epoch 20: jump of 200 μV in O1,
  • epoch 25: flat C3.
e2_art = deepcopy(e2)
Random.seed!(2)
t = e2_art.epoch_time
# eye blink: Gaussian wave of 150 μV, centered at 1 s
for ch in ["Fp1", "Fp2"]
    e2_art.data[ix(ch), :, 5] .+= 150 .* exp.(-(t .- 1) .^ 2 ./ (2 * 0.1^2))
end
e2_art.data[ix("T3"), 200:400, 12] .+= 30 .* randn(201)
e2_art.data[ix("O1"), 256:end, 20] .+= 200
e2_art.data[ix("C3"), :, 25] .= 0.0
println(rpad("method", 9), "bad epochs")
for m in methods
    println(rpad(m, 9), epoch_reject(e2_art; ch = "eeg", method = m))
end
println(rpad("amp 100", 9),
        epoch_reject(e2_art; ch = "eeg", method = :amp, amp_t = 100))
method   bad epochs
flat     [25]
amp      Int64[]
var      [20, 22, 25, 26]
p2p      [5, 12, 20, 22, 25, 26]
tkeo     [12, 20, 22, 24, 25, 26]
kurt     [2, 3, 4, 5, 10, 12, 16, 21, 23, 26, 27, 30]
z        [1, 4, 5, 12, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27]
rmse     [4, 12, 18, 20, 21, 22, 23, 24, 25, 26, 27]
rmsd     [4, 12, 18, 20, 21, 22, 23, 24, 25, 26, 27]
euclid   [4, 12, 18, 20, 21, 22, 23, 24, 25, 26, 27]
ransac   [4, 5, 12, 14, 19, 20, 21, 22, 23, 25, 27, 28]
amp 100  [5, 20, 25, 26]

Besides the injected artifacts, the recording contains real ones: a large artifact in almost all channels in epoch 26 (about 200 μV), smaller ones in F8 (epoch 22) and T3 (epoch 25), and the whole segment of epochs 15-27 is noisier than the rest. Because of the T3 artifact, epoch 25 is flagged by several methods regardless of the flat C3; the flat channel itself is found by :flat and :var.

Method Injected artifacts Other epochs flagged
:flat 25 -
:amp (400 μV) - -
:amp (100 μV) 5, 20, 25 26
:var 20, 25 22, 26
:p2p 5, 12, 20, 25 22, 26
:tkeo 12, 20, 25 22, 24, 26
:kurt 5, 12 2, 3, 4, 10, 16, 21, 23, 26, 27, 30
:z 5, 12, 20, 25 1, 4, 14-19, 21-24, 26, 27
:rmse, :rmsd, :euclid 12, 20, 25 4, 18, 21-24, 26, 27
:ransac 5, 12, 20, 25 4, 14, 19, 21-23, 27, 28

:p2p, :var and :tkeo find the injected artifacts and the real large-amplitude epochs with few other flags. :kurt is sensitive to any spiky activity, :z flags all epochs with amplitudes larger than usual for the channel (the noisier segment), and :rmse and :ransac flag epochs in which a channel differs from the others.

Strictness of Epoch Rejection

With many channels, a strict criterion (nbad = 1: a single flagged channel is enough) flags many epochs, because each of the 19 channels is a separate chance for a false alarm. The detection is made less sensitive by increasing:

  • z: the robust z-score threshold,
  • nbad: the number of channels that must be flagged.
for zz in [3.5, 5.0], nb in [1, 2]
    be = epoch_reject(e2_art;
                      ch = "eeg",
                      method = [:var, :p2p, :tkeo],
                      z = zz,
                      nbad = nb)
    println("z = ", zz, ", nbad = ", nb, ": ", be)
end
z = 3.5, nbad = 1: [5, 12, 20, 22, 24, 25, 26]
z = 3.5, nbad = 2: [20, 25, 26]
z = 5.0, nbad = 1: [12, 20, 25, 26]
z = 5.0, nbad = 2: [26]

Tip: With nbad = 2, artifacts flagged in a single channel are no longer detected: here the muscle burst (epoch 12, T3) and the eye blink (epoch 5, flagged in Fp2 only). Increasing z keeps such artifacts as long as they are large.

Removing Bad Epochs

be = epoch_reject(e2_art; ch = "eeg", method = [:flat, :var, :p2p, :tkeo])
e2_ok = delete_epoch(e2_art; ep = be)
println("epochs: ", nepochs(e2_art), " → ", nepochs(e2_ok))
epochs: 30 → 23

Different Thresholds for Different Channel Types

Absolute thresholds depend on the signal units; for recordings with several channel types (e.g. MEG magnetometers and gradiometers with EEG and EOG), :amp is applied separately to each type, and the epochs are removed together:

# rejection thresholds
thresholds = ["mag" => 4000,   # fT
              "grad" => 4000,  # fT/cm
              "eeg" => 150,    # μV
              "eog" => 250]    # μV
be = Int64[]
for (type, amp_t) in thresholds
    append!(be, epoch_reject(meg; ch = type, method = :amp, amp_t = amp_t))
end
delete_epoch!(meg; ep = unique(be))

Limitations

  • :var, :rmse and the other comparisons between channels assume that most channels are good; if many channels are bad, the median itself is affected.
  • :var flags channels with genuinely low variance (e.g. channels close to the reference) as well as flat or bridged channels.
  • :ransac compares each channel with its neighbors, so the neighbors of a very noisy channel may also be flagged; it also needs channel locations.
  • :kurt and :z are sensitive in epoch_reject(): kurtosis of short epochs varies a lot even without artifacts.
  • Detection finds abnormal signals, not their cause; a short burst of epileptiform activity or a high-amplitude sleep wave may be flagged as an artifact. Check the results visually before deleting data.