Project F: The software pipeline
From a raw file to cleaned epochs to a topographic map, in MNE-Python, on public data. Build the pipeline before you have your own data so you can trust it when you do.
You are skimming: the title, the first figure, and the short version. Switch to Read in the header for the full page, or Deep to open every deep dive.
Every EEG analysis is the same six steps: load, filter, re-reference, cut into epochs around events, remove artifacts, and plot something across the scalp. MNE-PythonMNE-PythonThe standard open-source Python library for EEG and MEG analysis: filtering, epoching, ICA, topographies, source localization. Glossary entry does each step in a line or two and is what nearly every lab uses. In this project you build the pipeline on the PhysioNet motor dataset, where the events are known, so that by the time your own board produces data the pipeline is already tested and you know what every step does to the signal.
You cut the recording into one-second epochs and then apply a 0.5 Hz high-pass filter to each epoch. What goes wrong?
The filter cannot do its job on one-second pieces. A 0.5 Hz high-pass has a time constant of a fraction of a second and needs a long stretch of data to settle; applied to a short epoch it smears the edges and can invent a slope across the whole epoch. Filter the continuous recording first, then cut. This is the single most common ordering mistake and the reason the pipeline’s steps come in the order they do.
The dataset
PhysioNet EEG Motor Movement/Imagery, runs 3, 7, and 11 for a subject: the subject opens and closes the left or right fist when a target appears. Events are annotated in the file. 64 channels at 160 Hz.
The pipeline
import mne
from mne.datasets import eegbci
raws = []
for run in (3, 7, 11):
path = eegbci.load_data(subject=1, runs=[run])[0]
raw = mne.io.read_raw_edf(path, preload=True)
eegbci.standardize(raw); raw.set_montage('standard_1005')
raws.append(raw)
raw = mne.concatenate_raws(raws)
# 1. filter the continuous data
raw.filter(1., 40., fir_design='firwin')
# 2. re-reference to the average of all channels
raw.set_eeg_reference('average', projection=False)
# 3. events and epochs: T1 = left fist, T2 = right fist
events, event_id = mne.events_from_annotations(raw, event_id=dict(T1=1, T2=2))
epochs = mne.Epochs(raw, events, event_id, tmin=-1., tmax=4., baseline=None, preload=True)
# 4. artifact removal with ICA
ica = mne.preprocessing.ICA(n_components=20, random_state=0)
ica.fit(epochs)
ica.plot_components() # find the blink component by its frontal topography
ica.exclude = [0] # whichever index looked like a blink
epochs_clean = ica.apply(epochs.copy())
# 5. something across the scalp: mu-band power during movement vs. rest
mu = epochs_clean.copy().filter(8., 13.).apply_hilbert(envelope=True)
move = mu.copy().crop(0.5, 3.5).get_data().mean(axis=(0, 2))
rest = mu.copy().crop(-1., 0.).get_data().mean(axis=(0, 2))
mne.viz.plot_topomap((move - rest) / rest, epochs.info)
The final map should show a decrease in mu power over the motor cortex, opposite the moving hand: the event-related desynchronizationEvent-related desynchronization (ERD)The drop in rhythmic power over a cortical area when it becomes active, such as mu weakening when you imagine moving your hand. Glossary entry that motor imagery BCIs decode. If it does, your pipeline works end to end on real data.
What each step is for
Filter first, on continuous data. The high-pass removes drift; the low-pass removes muscle and hum. Filters need long stretches of data; see the prediction above. MNE’s default FIR filter is zero-phase, so nothing shifts in time.
Re-reference. The file’s original reference is one electrode, and everything is measured against it. The average reference subtracts the mean of all channels, which approximates a neutral reference when electrodes cover the head. The montages explainer is about why this changes every plot.
Epoch. Cut a window around each event. Now you can compare what happened after events to what happened before, and average across events to reduce noise.
Clean. ICAIndependent component analysis (ICA)A method that unmixes multichannel data into statistically independent sources, so blinks and heartbeat can be identified and removed. Glossary entry unmixes the channels into components; blinks and heartbeat come out as their own components with characteristic topographies, and you subtract them. Look at every component’s map before you exclude it. Never exclude blindly.
Plot across the scalp. The topographic map is EEG’s native picture. Anything you compute per channel can be drawn on the head.
- Run the pipeline as written and get a topographic map. Check that the mu decrease is over the correct hemisphere (left hand moves, right hemisphere decreases).
- Break it on purpose: filter after epoching and compare the epoch edges. Skip re-referencing and compare the map. Skip ICA and see whether blinks show up in the frontal channels of the average.
- Wrap the steps in a function that takes a subject number and returns the epochs. Run it for five subjects. Note which ones show clean ERD and which do not.
- Save the cleaned epochs to disk in MNE’s format. Your Phase 3 projects start from that file.