Applying a Formula

using NeuroAnalyzer
using Statistics
NeuroAnalyzer.progress_bar = false
eeg = import_edf("files/eeg/eeg.edf")
# 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)
# 10-second epochs
e10 = epoch(eeg; ep_len = 10)

apply() evaluates a Julia expression for each selected channel and each epoch. It is a quick way to compute per-channel, per-epoch measures (statistics, amplitudes, custom features) or to transform the signals with any Julia function, without writing loops.

How apply() Works

The formula is a string in which obj stands for the signal of one channel in one epoch (a vector of samples):

r = apply(e10; ch = ["Fp1", "Fp2", "O1"], f = "std(obj)")
println("size: ", size(r))

size: (3, 1, 120)

The formula is compiled once and evaluated for every channel and epoch (in parallel); the result is an array of size channels × (length of the result) × epochs:

• a single number per channel and epoch (e.g. mean(obj), std(obj)) gives channels × 1 × epochs,
• a vector (e.g. obj .^ 2) gives channels × (vector length) × epochs; its length must be the same for all channels and epochs.

The result must be real; complex results (e.g. fft(obj)) must be converted, e.g. with abs.(). Any function available in Julia and in the packages used by NeuroAnalyzer (e.g. Statistics, FFTW) can be used.

apply() vs vch():

	apply()	vch() (Virtual channels)
formula operates on	one channel in one epoch (a whole signal)	all channels at one sample
combines channels	no	yes (e.g. Fp1 - F7)
can reduce over time	yes (e.g. mean(obj), std(obj))	no (element-wise)
result	channels × length × epochs	1 × samples × epochs

Per-Epoch Measures

Standard deviation of each channel in each epoch:

sd = apply(e10; ch = "eeg", f = "std(obj)")
println("channels × 1 × epochs: ", size(sd))
println("mean SD across epochs [μV]:")
for (idx, ch) in enumerate(labels(e10)[get_channel(e10; ch = "eeg")][1:5])
    println("  ", rpad(ch, 4), round(mean(sd[idx, 1, :]); digits = 1))
end

channels × 1 × epochs: (19, 1, 120)
mean SD across epochs [μV]:
  Fp1 6.2
  Fp2 10.2
  F3  6.1
  F4  3.5
  C3  6.2

Example: finding epochs with large amplitudes

The peak-to-peak amplitude of each channel and epoch, and the epochs in which it exceeds 150 μV in any channel (e.g. because of eye blinks or movements):

ptp = apply(e10; ch = "eeg", f = "maximum(obj) - minimum(obj)")
bad_ep = findall(vec(any(ptp[:, 1, :] .> 150; dims = 1)))
println("epochs with peak-to-peak amplitude > 150 μV: ", bad_ep)

epochs with peak-to-peak amplitude > 150 μV: [1, 6, 7, 10, 15, 16, 29, 32, 34, 36, 37, 62, 63, 73, 92, 108]

Such epochs can then be inspected or removed with delete_epoch(); see also Bad channels and epochs for automatic detection methods.

Other examples

Formula	Result (per channel and epoch)
mean(obj)	mean (DC offset)
std(obj)	standard deviation
maximum(abs.(obj))	maximum absolute amplitude
sqrt(mean(obj .^ 2))	root mean square (RMS)
sum(abs.(diff(obj)))	line length (signal complexity, used in seizure detection)
count(abs.(obj) .> 100)	number of samples above 100 μV
quantile(obj, 0.95)	95th percentile

rms = apply(e10; ch = "O1", f = "sqrt(mean(obj .^ 2))")
ll = apply(e10; ch = "O1", f = "sum(abs.(diff(obj)))")
println("O1, first 5 epochs")
println("  RMS [μV]:    ", round.(rms[1, 1, 1:5]; digits = 1))
println("  line length: ", round.(ll[1, 1, 1:5]; digits = 0))

O1, first 5 epochs
  RMS [μV]:    [62.4, 6.2, 5.3, 5.8, 5.6]
  line length: [3272.0, 1981.0, 1947.0, 1939.0, 2037.0]

Transforming Signals

Formulas returning a vector transform the signal:

Formula	Result
obj .- mean(obj)	signal with the mean removed (per epoch)
(obj .- mean(obj)) ./ std(obj)	z-scored signal (per epoch)
abs.(obj)	rectified signal
obj .^ 2	instantaneous power
cumsum(obj)	cumulative sum (integral)
diff(obj)	first derivative (one sample shorter)
abs.(fft(obj))	amplitude spectrum

Z-scoring each channel in each epoch:

z = apply(e10; ch = "eeg", f = "(obj .- mean(obj)) ./ std(obj)")
println("size: ", size(z))
println("mean: ", round(mean(z[1, :, 1]); digits = 6) + 0.0,
        ", SD: ", round(std(z[1, :, 1]); digits = 6))

size: (19, 2560, 120)
mean: 0.0, SD: 1.0

A transformed signal of the same length can be written back into a copy of the object (or added as a new channel with add_channel!()):

e10_z = deepcopy(e10)
eeg_ch = get_channel(e10_z; ch = "eeg")
e10_z.data[eeg_ch, :, :] = z

Tip: NeuroAnalyzer has dedicated functions for many of these operations (e.g. normalize(), remove_dc(), amp(), rms(), psd()), which are faster and record the operation in the processing history; apply() is meant for measures and transformations not available otherwise.