using NeuroAnalyzer
using Statistics
NeuroAnalyzer.verbose = false
NeuroAnalyzer.progress_bar = false
const CH = ["Fp1", "Fp2", "F3", "F4", "C3", "C4", "P3", "P4", "O1", "O2",
"F7", "F8", "T3", "T4", "T5", "T6", "Fz", "Cz", "Pz"]
function workflow()
t = Dict{String, Float64}()
t["import"] = @elapsed eeg = import_edf("files/eeg/eeg.edf")
t["channels"] = @elapsed begin
keep_channel!(eeg; ch = CH)
load_locs!(eeg;
file_name = "files/locs/standard-10-20-cap19-elmiko.locs")
end
t["filter"] = @elapsed begin
NeuroAnalyzer.filter!(eeg; ch = "eeg", fprototype = :fir,
ftype = :hp, cutoff = 0.5, bw = 0.5)
NeuroAnalyzer.filter!(eeg; ch = "eeg", fprototype = :fir,
ftype = :lp, cutoff = 45, bw = 11.25)
end
t["reference"] = @elapsed reference_avg!(eeg)
t["resample"] = @elapsed resample!(eeg; new_sr = 128)
t["epoch"] = @elapsed ep = epoch(eeg; ep_len = 2)
t["psd"] = @elapsed bp = band_power(ep; ch = "eeg", flim = (8, 12),
method = :welch,
wlen = 128, woverlap = 64)
# mean alpha power over epochs, μV²
return t, vec(mean(bp; dims = 2))
end
steps = ["import", "channels", "filter", "reference",
"resample", "epoch", "psd"]
t_first = @elapsed workflow() # includes compilation
times = Dict(s => Float64[] for s in steps)
for _ in 1:10
t, _ = workflow()
for s in steps
push!(times[s], t[s])
end
end
for s in steps
println(rpad(s, 10), round(median(times[s]); digits = 3))
end
println(rpad("total", 10),
round(sum(median(times[s]) for s in steps); digits = 3))NeuroAnalyzer vs MNE
Resting-state EEG workflow
The same EDF recording processed with the same steps in NeuroAnalyzer and MNE-Python, with filter and spectral parameters matched so that the two pipelines give the same numbers.
Data: eeg.edf from the tutorial files: a 20-minute resting-state recording, 24 channels at 256 Hz, of which the 19 EEG channels of the 10–20 system are used.
Workflow:
| Step | Operation |
|---|---|
| import | read the EDF file |
| channels | keep 19 EEG channels, load the electrode locations |
| filter | FIR (Hamming-windowed sinc) high-pass 0.5 Hz and low-pass 45 Hz |
| reference | common average reference |
| resample | 256 → 128 Hz |
| epoch | 2-second epochs (602 epochs) |
| psd | Welch PSD per epoch (128-sample Hann windows, 50% overlap), alpha (8–12 Hz) band power |
The filter transition bands are set to MNE’s defaults: 0.5 Hz for the 0.5 Hz high-pass and 11.25 Hz for the 45 Hz low-pass (both designs give 1691 and 77 taps). Band power is the sum of the PSD bins in 8–12 Hz × the frequency resolution (1 Hz) in both libraries.
Each step is timed separately; the table shows medians of 10 runs, after a first (warm-up) run.
Run with julia -t auto (multi-threaded) or julia -t 1 (single thread).
import time
import numpy as np
import mne
mne.set_log_level("WARNING")
CH = ["Fp1", "Fp2", "F3", "F4", "C3", "C4", "P3", "P4", "O1", "O2",
"F7", "F8", "T3", "T4", "T5", "T6", "Fz", "Cz", "Pz"]
def workflow():
t = {}
t0 = time.perf_counter()
raw = mne.io.read_raw_edf("files/eeg/eeg.edf", preload=True)
t["import"] = time.perf_counter() - t0
t0 = time.perf_counter()
raw.rename_channels(lambda name: name.replace("EEG ", ""))
raw.pick(CH)
raw.set_montage(mne.channels.read_custom_montage(
"files/locs/standard-10-20-cap19-elmiko.locs"))
t["channels"] = time.perf_counter() - t0
t0 = time.perf_counter()
raw.filter(0.5, None, fir_design="firwin") # transition 0.5 Hz
raw.filter(None, 45.0, fir_design="firwin") # transition 11.25 Hz
t["filter"] = time.perf_counter() - t0
t0 = time.perf_counter()
raw.set_eeg_reference("average")
t["reference"] = time.perf_counter() - t0
t0 = time.perf_counter()
raw.resample(128)
t["resample"] = time.perf_counter() - t0
t0 = time.perf_counter()
epochs = mne.make_fixed_length_epochs(raw, duration=2.0, preload=True)
t["epoch"] = time.perf_counter() - t0
t0 = time.perf_counter()
psd = epochs.compute_psd(method="welch", n_fft=128, n_per_seg=128,
n_overlap=64, window="hann")
p, f = psd.get_data(return_freqs=True)
band = (f >= 8) & (f <= 12)
bp = p[:, :, band].sum(axis=2) * (f[1] - f[0])
t["psd"] = time.perf_counter() - t0
# mean alpha power over epochs, V² → μV²
return t, bp.mean(axis=0) * 1e12
steps = ["import", "channels", "filter", "reference",
"resample", "epoch", "psd"]
t0 = time.perf_counter()
workflow()
t_first = time.perf_counter() - t0
times = {s: [] for s in steps}
for _ in range(10):
t, _ = workflow()
for s in steps:
times[s].append(t[s])
for s in steps:
print(f"{s:<10}{np.median(times[s]):.3f}")
print(f"{'total':<10}{sum(np.median(times[s]) for s in steps):.3f}")Test system: AMD Ryzen Threadripper 3960X (24 cores, 48 threads), Debian Linux (7.2.11 kernel); Julia 1.13.2; Python 3.14, MNE 1.13.2, NumPy 2.4.6, SciPy 1.17.1.
Results (seconds, median of 10 runs):
| Step | NeuroAnalyzer (48 threads) | NeuroAnalyzer (1 thread) | MNE |
|---|---|---|---|
| import | 0.119 | 0.099 | 0.028 |
| channels | 0.045 | 0.045 | 0.011 |
| filter | 0.425 | 0.655 | 0.490 |
| reference | 0.109 | 0.104 | 0.030 |
| resample | 0.034 | 0.165 | 0.149 |
| epoch | 0.010 | 0.011 | 0.027 |
| psd | 0.288 | 0.222 | 1.862 |
| total | 1.030 | 1.301 | 2.598 |
| first run (incl. compilation) | 14.6 | 18.7 | 3.3 |
- Once compiled, the NeuroAnalyzer workflow is ≈ 2.5× faster than MNE with multiple threads and ≈ 2× faster on a single thread. Most of the difference comes from the PSD (MNE: ≈ 1.9 s, NeuroAnalyzer: ≈ 0.25 s); MNE is faster at importing, channel selection and referencing.
- Filtering and resampling are run per channel in parallel by NeuroAnalyzer; with 19 channels the gain from multithreading is modest.
- The first run in Julia includes just-in-time compilation of all the functions used (≈ 15 s), which is paid once per session. For a single file MNE is faster; for batch processing of many files, NeuroAnalyzer is faster.
Validation
Mean alpha power (μV²) per channel, averaged over epochs:
| Channel | NeuroAnalyzer | MNE | Ratio |
|---|---|---|---|
| Fp1 | 1.903 | 1.894 | 1.004 |
| Fp2 | 94.789 | 93.042 | 1.019 |
| F3 | 1.305 | 1.298 | 1.005 |
| F4 | 2.028 | 2.022 | 1.003 |
| C3 | 0.974 | 0.966 | 1.009 |
| C4 | 0.790 | 0.783 | 1.009 |
| P3 | 0.994 | 0.985 | 1.009 |
| P4 | 1.439 | 1.429 | 1.007 |
| O1 | 1.329 | 1.318 | 1.008 |
| O2 | 1.593 | 1.583 | 1.006 |
| F7 | 1.607 | 1.604 | 1.002 |
| F8 | 1.490 | 1.480 | 1.007 |
| T3 | 1.194 | 1.187 | 1.005 |
| T4 | 1.556 | 1.543 | 1.008 |
| T5 | 1.394 | 1.384 | 1.007 |
| T6 | 1.576 | 1.565 | 1.007 |
| Fz | 1.377 | 1.368 | 1.007 |
| Cz | 0.629 | 0.623 | 1.009 |
| Pz | 0.980 | 0.973 | 1.008 |
The results agree within 1% (median ratio 1.007, Pearson’s r = 1.000; Fp2, which contains a large artifact, within 2%). The remaining differences come from how the FIR filters are applied:
- NeuroAnalyzer filters in two passes (forward and backward,
dir = :twopass), so the magnitude response is squared; MNE applies the filter once and compensates its delay. Both are zero-phase and have the same pass band, but the two-pass filter attenuates more in the transition bands. - The recording contains strong 50 Hz line noise. 50 Hz lies inside the 45–50.6 Hz transition band of the low-pass filter: MNE keeps ≈ 35% of its power, NeuroAnalyzer < 0.01%. This does not affect the alpha band, but it should be considered when comparing broadband measures; a notch filter (see Power line noise) removes it in both libraries.