Electroencephalography is one of the few ways to read brain activity non-invasively in real time, which makes it attractive to machine learners and also makes it easy to fool yourself with. The signal is microvolts riding on millivolts of interference, the "features" that textbooks talk about are spectral, and the usual cross-validation habits leak information in ways that are specific to this kind of data. This post builds a synthetic EEG-like signal in NumPy (no real recordings are used anywhere), runs it through a minimal preprocessing and band-power pipeline, and uses it to show what each step does and where the traps are.
What the signal is
Scalp EEG measures potential differences of a few to a few hundred microvolts between electrodes, produced by the summed post-synaptic currents of large, synchronously active cortical populations. Three facts shape everything downstream:
- Bandwidth is low. Most of the usable power is below 40–50 Hz, so sampling rates of 250–1000 Hz are typical; the Nyquist frequency must stay above the highest frequency of interest and an analogue anti-aliasing filter must act before the ADC.
- The spectrum is roughly
. Power falls with frequency, so a 1 µV oscillation at 30 Hz can be meaningful while a 1 µV fluctuation at 2 Hz is nothing. - Interference dominates. Mains hum at 50 or 60 Hz, electrode drift, eye blinks (hundreds of microvolts, frontal), muscle activity (broadband, above 20 Hz) and heartbeat all appear in the trace at amplitudes larger than the brain signal.
The classical frequency bands are a convention rather than a physiological law, but they organize both the literature and the feature extraction:
| Band | Range (this post) | Textbook association | Relative power, synthetic signal |
|---|---|---|---|
| delta | 1–4 Hz | deep sleep; also slow artefacts and drift | 49.6 % |
| theta | 4–8 Hz | drowsiness, working-memory tasks | 8.7 % |
| alpha | 8–13 Hz | relaxed wakefulness, eyes closed (occipital) | 32.8 % |
| beta | 13–30 Hz | active cognition, motor planning | 6.7 % |
| gamma | 30–45 Hz | attention and binding; often muscle contamination | 1.3 % |
The last column comes from the synthetic signal below after filtering and says nothing about any brain; it is there to show how a single blink inflates the delta band of a ten-second window.
A synthetic signal
The signal is ten seconds at 250 Hz, built from components with known amplitudes: bursting alpha at 10 Hz (20 µV, amplitude-modulated at 0.3 Hz), beta at 20 Hz (5 µV), theta at 6 Hz (6 µV), pink noise (8 µV RMS), 50 Hz mains interference (10 µV), slow electrode drift at 0.2 Hz (30 µV) and one blink-shaped Gaussian bump at
import numpy as np
fs, dur = 250, 10.0
t = np.arange(0, dur, 1 / fs)
rng = np.random.default_rng(0)
def pink(n, rng): # 1/f noise by spectral shaping
X = np.fft.rfft(rng.normal(size=n))
f = np.fft.rfftfreq(n, 1 / fs); f[0] = f[1]
x = np.fft.irfft(X / np.sqrt(f), n)
return x / x.std()
alpha = 20 * (0.6 + 0.4 * np.sin(2 * np.pi * 0.3 * t)) * np.sin(2 * np.pi * 10 * t)
beta = 5 * np.sin(2 * np.pi * 20 * t + 1.0)
theta = 6 * np.sin(2 * np.pi * 6 * t + 2.0)
line = 10 * np.sin(2 * np.pi * 50 * t)
drift = 30 * np.sin(2 * np.pi * 0.2 * t)
blink = 150 * np.exp(-0.5 * ((t - 4.0) / 0.12) ** 2)
raw = alpha + beta + theta + line + drift + blink + 8 * pink(len(t), rng)
Filtering
The minimal offline pipeline is a notch at the mains frequency and a band-pass that removes drift below 1 Hz and everything above the band of interest. Both can be built from a windowed-sinc low-pass:
- A low-pass FIR with cutoff
is the sinc multiplied by a window (Hamming here) and normalized to unit DC gain. - A band-pass is the difference of two low-passes,
. - A band-stop (notch) is an impulse minus a band-pass.
def lowpass(fc, numtaps):
n = np.arange(numtaps) - (numtaps - 1) / 2
h = 2 * fc / fs * np.sinc(2 * fc / fs * n) * np.hamming(numtaps)
return h / h.sum()
def bandpass(lo, hi, numtaps):
return lowpass(hi, numtaps) - lowpass(lo, numtaps)
def bandstop(lo, hi, numtaps):
d = np.zeros(numtaps); d[(numtaps - 1) // 2] = 1.0
return d - bandpass(lo, hi, numtaps)
h_notch, h_bp = bandstop(48.0, 52.0, 751), bandpass(1.0, 45.0, 751)
filtered = np.convolve(np.convolve(raw, h_notch, mode="same"), h_bp, mode="same")
A symmetric FIR applied with mode="same" is zero-phase: it delays nothing, which is what you want offline and cannot have online. The 751-tap filter (3 s long at 250 Hz) has a transition band of roughly 1 Hz; its measured response is 0.00 dB at 10 Hz, −58 dB at the 0.2 Hz drift and −119 dB at 50 Hz.
The middle panel is the first pitfall. A blink is a slow, large deflection with most of its energy between about 0.5 and 5 Hz, inside the delta band, so a 1 Hz high-pass does not remove it; the filtered peak is still 91 µV against an alpha amplitude of 20 µV. That single event accounts for the 49.6 % "delta power" in the table. Artifacts have to be handled by rejection (dropping windows that exceed an amplitude threshold) or by source separation such as ICA, not by frequency filtering, because they overlap the bands of interest.
Spectra and band power
Band power is the integral of the power spectral density over a frequency range. Welch's method estimates the PSD by averaging periodograms of overlapping, windowed segments, which trades frequency resolution for variance:
def welch(x, nperseg=500, noverlap=250): # 2 s Hann segments, 50 % overlap
win = np.hanning(nperseg); step = nperseg - noverlap; acc = []
for s in range(0, len(x) - nperseg + 1, step):
seg = x[s:s + nperseg]
seg = (seg - seg.mean()) * win
acc.append(np.abs(np.fft.rfft(seg)) ** 2)
psd = np.mean(acc, 0) / (fs * np.sum(win ** 2)) # one-sided density in µV²/Hz
psd[1:-1] *= 2
return np.fft.rfftfreq(nperseg, 1 / fs), psd
def bandpower(f, psd, lo, hi):
m = (f >= lo) & (f < hi)
return np.trapezoid(psd[m], f[m])
Why the scaling factor is fs times the window energy
For a windowed segment, the periodogram
With 2 s segments the frequency resolution is 0.5 Hz, which is enough for 4 Hz-wide bands and coarse for separating, say, 9 Hz from 10 Hz. Longer segments resolve better and average less; for a ten-second window there are only a handful of segments either way, which is why single-trial EEG spectra are noisy and why models that consume them need regularization.
Spectral leakage and windows
A finite segment is implicitly multiplied by a rectangular window, whose spectrum has slowly decaying sidelobes. Energy from a strong component then leaks into every other frequency bin, which matters in EEG because the
The tone is at 10.3 Hz rather than 10 Hz on purpose: a frequency that is not an integer multiple of the bin spacing is the worst case for leakage and the typical case for real oscillations.
Pitfalls that inflate accuracy
DANGER
"Leakage" has two meanings in this field and both will hurt you. Spectral leakage is the windowing effect above. Data leakage is fitting any part of the pipeline (filters designed from the data, scalers, ICA components, channel selection, or the classifier itself) on samples that end up in the test set. With EEG the most common form is splitting windows at random when consecutive windows overlap or come from the same trial: adjacent windows share samples and the same slow drift, and a classifier can score far above chance by recognizing the recording rather than the mental state.
The habits that avoid most of the damage:
- Split by session or subject, never by window. Evaluate cross-session and, if the application needs it, cross-subject. Expect a large drop; that drop is the real number.
- Fit every preprocessing step on the training fold only. Filters with fixed coefficients are fine; anything data-dependent (re-referencing to a data-derived component, ICA, CSP spatial filters, normalization statistics) is part of the model.
- Reject or flag artefact windows with a rule that does not look at labels. An amplitude threshold applied uniformly is better than a hand-cleaned dataset whose cleaning correlated with the classes.
- Mind causality when going online. The zero-phase filter above uses future samples. A causal FIR of the same length delays the signal by half its length (1.5 s here); online systems use short IIR filters and accept the phase distortion.
- Report chance level honestly. With unbalanced classes and few trials, chance is not
; permutation tests give the right baseline.
Where machine learning comes in
Once the data are windowed, filtered and honestly split, the decoder is the easy part. Linear discriminants on log band power are a strong baseline for motor-imagery and sleep-staging tasks; common spatial patterns add learned spatial filters; Riemannian methods classify covariance matrices directly; and compact convolutional networks such as EEGNet learn temporal and spatial filters end to end, often matching the hand-crafted pipeline with far fewer assumptions. All of them stand or fall on the preprocessing and the split, which is why this post stops here.
Further reading
- MNE-Python, the standard open-source toolbox for EEG/MEG analysis, whose tutorials cover every step above on real data: https://mne.tools/stable/index.html
- Lotte et al. (2018), A review of classification algorithms for EEG-based brain–computer interfaces: a 10 year update: https://iopscience.iop.org/article/10.1088/1741-2552/aab2f2
- Roy et al. (2019), Deep learning-based electroencephalography analysis: a systematic review: https://arxiv.org/abs/1901.05498
- Lawhern et al. (2016), EEGNet: a compact convolutional network for EEG-based brain–computer interfaces: https://arxiv.org/abs/1611.08024
- Schirrmeister et al. (2017), Deep learning with convolutional neural networks for EEG decoding and visualization: https://arxiv.org/abs/1703.05051
- MOABB, a benchmark framework that enforces proper cross-session and cross-subject evaluation: https://moabb.neurotechx.com/
- BCI Competition IV datasets, the usual public starting point: https://www.bbci.de/competition/iv/
- SciPy
signalmodule reference, for production-grade FIR/IIR design andwelchonce you stop writing them yourself: https://docs.scipy.org/doc/scipy/reference/signal.html