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.

A closed-loop BCI pipeline, from scalp potentials to feedback 1. Acquisition electrodes (µV) amplifier + ADC fs = 250–1000 Hz reference, impedance 2. Preprocessing re-reference band-pass, notch artifact rejection / ICA epoching (windows) 3. Features band power (FFT/Welch) spatial filters (CSP) covariance / Riemann or raw epochs for CNNs 4. Decoder LDA / SVM / logistic EEGNet-style CNN calibrated per subject + rejection threshold 5. Application cursor, speller, prosthesis, game latency budget ~100–500 ms feedback to the user the brain adapts to the decoder while the decoder adapts to the brain Offline (dataset) analysis zero-phase filters, whole-recording statistics cross-validation must split by session / subject leakage risk: fitting filters or scalers on test data Online (closed-loop) constraints causal filters only (group delay adds latency) sliding windows, no future samples non-stationarity: recalibrate or adapt over time
A closed-loop BCI pipeline. Steps 2 and 3 are where most machine-learning effort goes; the feedback arrow is what makes online evaluation different from a benchmark on a saved dataset.

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 1/f. 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 t=4 s (150 µV, width 0.12 s).

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:

  1. A low-pass FIR with cutoff fc is the sinc 2fc/fs⋅sinc(2fcn/fs) multiplied by a window (Hamming here) and normalized to unit DC gain.
  2. A band-pass is the difference of two low-passes, LP(fhi)−LP(flo).
  3. 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.

Synthetic EEG before and after filtering, and its power spectral density
Top: the synthetic raw trace with drift, mains hum and a blink. Middle: after the notch and 1–45 Hz band-pass; the drift and hum are gone, the blink is attenuated but very much present. Bottom: Welch power spectral density of both, with the band boundaries shaded; the 50 Hz line disappears, the 1/f slope and the three oscillation peaks remain.

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 |X[k]|2 has units of (µV)² times samples². Dividing by fs∑nw[n]2 converts it to a density in µV²/Hz that is independent of the window shape and segment length, so that integrating the PSD over all frequencies returns the signal variance (Parseval). The factor of 2 on the interior bins folds the negative frequencies into a one-sided spectrum; the DC and Nyquist bins have no mirror image and are left alone.

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 1/f slope means the delta band is always "strong". Tapering with a Hann window drops the sidelobes by tens of decibels at the price of a wider main lobe:

Spectral leakage of a sinusoid with rectangular versus Hann windows, and the FIR filter response
Left: the spectrum of a pure 10.3 Hz tone in a 2 s window. Five hertz away from the tone the rectangular window's sidelobes are still at about −30 dB, while the Hann window's are below −65 dB. Right: the frequency response of the notch plus band-pass FIR chain used above.

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 1/K; 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