Full EEG/MEG analysis pipeline with MNE-Python: preprocessing, ICA artifact removal, ERP computation, time-frequency analysis, and resting-state power spectral density.
Full EEG/MEG analysis pipeline with MNE-Python: preprocessing, ICA artifact removal, ERP computation, time-frequency analysis, and resting-state power spectral density.
MNE-Python is the standard open-source toolkit for analyzing M/EEG data. This skill provides a
complete pipeline from raw file loading through preprocessing, epoching, ERP computation,
time-frequency analysis, and basic source localization. All functions follow MNE idioms and
produce publication-ready figures.
"""
Load raw EEG data from EDF, FIF, BrainVision, or EEGLab formats.
Parameters
----------
filepath : str
Path to the raw EEG file.
preload : bool
If True, load data into memory (required for most operations).
montage_name : str
Standard montage to apply if the file lacks electrode positions.
Returns
-------
mne.io.BaseRaw
Loaded raw object with standard montage applied.
"""
"""
Apply standard EEG preprocessing: notch filter, bandpass, re-reference, interpolation.
Parameters
----------
raw : mne.io.BaseRaw
Input raw object (must be preloaded).
l_freq : float
High-pass cutoff (Hz). Use 0.1 for ERP, 1.0 for general use.
h_freq : float
Low-pass cutoff (Hz). Use 40 for ERPs, 100+ for oscillations.
notch_freq : float or list of float
Notch filter frequency (50 Hz Europe, 60 Hz North America).
reference : str
Re-referencing scheme: ``"average"``, ``"REST"``, or a channel name.
bad_channels : list of str, optional
Channels to mark as bad before interpolation.
resample_sfreq : float, optional
If provided, resample to this frequency after filtering.
Returns
-------
mne.io.BaseRaw
Preprocessed raw object.
"""
# Mark bad channels
if
"bads"
print
f"Marked bad: {bad_channels}"
# Notch filter (power line noise)
if
isinstance
int
float
else
# Also remove harmonics
for
in
for
in
range
1
int
"sfreq"
2
1
"hamming"
False
# Bandpass filter
filter
"hamming"
False
print
f"Filtered: {l_freq}–{h_freq} Hz, notch at {notch_freqs}"
"""
Compute evoked (ERP) responses for each condition.
Parameters
----------
epochs : mne.Epochs
Epoched data.
conditions : list of str, optional
Subset of conditions to compute. Defaults to all keys in ``epochs.event_id``.
channels : list of str, optional
Channels to pick before averaging (e.g. ``["Cz", "Pz", "Fz"]``).
Returns
-------
dict mapping condition name to mne.Evoked object.
"""
or
list
for
in
if
not
in
print
f"Warning: condition '{cond}' not found in epochs."
"""
Compute time-frequency representation using Morlet wavelets.
Parameters
----------
epochs : mne.Epochs
Input epochs.
freqs : array-like, optional
Frequencies to analyze (Hz). Defaults to 2–40 Hz in log space.
n_cycles : array or int
Number of wavelet cycles. Int = fixed; array = frequency-dependent.
picks : list of str, optional
Channels to include. Defaults to all EEG channels.
return_itc : bool
If True, also return inter-trial coherence (phase locking).
average : bool
If True, return average TFR; if False, return single-trial TFR.
decim : int
Temporal decimation factor to reduce memory.
Returns
-------
tuple : (AverageTFR power, AverageTFR itc) if return_itc else (AverageTFR power,)
"""
"""
Compute PSD and return band-averaged power (delta, theta, alpha, beta, gamma).
Parameters
----------
raw : mne.io.BaseRaw
Input raw (preloaded).
picks : list of str, optional
Channel names; defaults to all EEG channels.
fmin, fmax : float
Frequency range for PSD computation.
method : str
``"welch"`` or ``"multitaper"``.
Returns
-------
pd.DataFrame with columns: channel, delta, theta, alpha, beta, gamma (all in µV²/Hz).
"""
"delta"
1
4
"theta"
4
8
"alpha"
8
13
"beta"
13
30
"gamma"
30
50
or
"eeg"
False
True
1e12
# V²/Hz → µV²/Hz
for
in
enumerate
"channel"
for
in
0
return
Example A: Auditory ERP (N100 / P300) from Oddball Paradigm
An auditory oddball paradigm presents frequent "standard" tones (80%) and rare "deviant" tones
(20%). This example preprocesses the raw data, extracts ERPs, and quantifies the N100 (auditory
response at ~100 ms) and P300 (cognitive response at ~300 ms on Pz).
# ── Example A ─────────────────────────────────────────────────────────────# Uses MNE's built-in sample dataset (auditory oddball data)import mne
import numpy as np
import matplotlib.pyplot as plt
# --- Download MNE sample data (if needed) ------------------------------------
data_path = mne.datasets.sample.data_path()
raw_fif = str(data_path) + "/MEG/sample/sample_audvis_raw.fif"# --- Load (use only EEG channels) -------------------------------------------
raw = mne.io.read_raw_fif(raw_fif, preload=True, verbose=False)
raw.pick_types(eeg=True, eog=True, stim=True)
print(raw.info)
# --- Preprocess --------------------------------------------------------------
raw_prep = preprocess_raw(
raw,
l_freq=0.5,
h_freq=40.0,
notch_freq=60.0, # 60 Hz (North American power line)
reference="average",
bad_channels=["EEG 053"], # Example bad channel
resample_sfreq=None,
)
# --- ICA artifact removal ---------------------------------------------------
raw_clean, ica = run_ica_artifact_removal(
raw_prep,
n_components=20,
method="fastica",
eog_channels=["EOG 061"],
random_state=0,
)
# --- Epoch around auditory events -------------------------------------------# MNE sample dataset: event codes 1=LA, 2=RA, 3=LV, 4=RV, 5=smiley, 32=button
event_id = {"auditory/left": 1, "auditory/right": 2}
epochs = epoch_events(
raw_clean,
event_id=event_id,
tmin=-0.2,
tmax=0.5,
baseline=(-0.2, 0),
reject={"eeg": 100e-6},
)
# --- Compute ERPs -----------------------------------------------------------
evokeds = compute_erp(epochs, conditions=["auditory/left", "auditory/right"])
# --- Plot ERPs at Cz ---------------------------------------------------------
fig = plot_erp_comparison(
evokeds,
channel="EEG 059", # Approximately Cz in the sample data
title="Auditory ERP: Left vs Right (Oddball Paradigm)",
save_path="auditory_erp.png",
)
plt.show()
# --- Quantify N100 and P300 --------------------------------------------------defpeak_amplitude_latency(evoked: mne.Evoked, channel: str, tmin: float, tmax: float) -> dict:
"""Return peak amplitude (µV) and latency (ms) within a time window."""try:
idx = evoked.ch_names.index(channel)
except ValueError:
return {"amplitude_uV": np.nan, "latency_ms": np.nan}
t_mask = (evoked.times >= tmin) & (evoked.times <= tmax)
data = evoked.data[idx, t_mask] * 1e6
times = evoked.times[t_mask] * 1e3
peak_idx = np.argmax(np.abs(data))
return {"amplitude_uV": float(data[peak_idx]), "latency_ms": float(times[peak_idx])}
TARGET_CH = "EEG 059"for cond, evoked in evokeds.items():
n100 = peak_amplitude_latency(evoked, TARGET_CH, 0.070, 0.150)
p300 = peak_amplitude_latency(evoked, TARGET_CH, 0.250, 0.450)
print(f"\n{cond}")
print(f" N100: {n100['amplitude_uV']:.2f} µV @ {n100['latency_ms']:.0f} ms")
print(f" P300: {p300['amplitude_uV']:.2f} µV @ {p300['latency_ms']:.0f} ms")
# --- Difference wave (deviant minus standard) --------------------------------# If you have a full oddball dataset with standard/deviant labels:# diff_wave = mne.combine_evoked([evokeds["deviant"], evokeds["standard"]], weights=[1, -1])# diff_wave.plot(picks=[TARGET_CH], titles={"eeg": "MMN / Difference Wave"})# --- Topographic map at P300 peak (300–400 ms) -------------------------------for cond, evoked in evokeds.items():
fig_topo = evoked.plot_topomap(
times=[0.1, 0.2, 0.3, 0.4],
average=0.05,
show=False,
)
fig_topo.suptitle(f"Topography: {cond}", fontsize=11)
fig_topo.savefig(f"topo_{cond.replace('/', '_')}.png", dpi=120)
Example B: Resting-State Alpha Power Between Eyes-Open and Eyes-Closed
Alpha band power (8–13 Hz) reliably increases during eyes-closed rest. This example compares
band power across conditions and produces a power spectral density plot.
# ── Example B ─────────────────────────────────────────────────────────────# Assumes two raw files: resting_eyes_open.edf and resting_eyes_closed.edf# Replace paths with actual file locations
EYES_OPEN_FILE = os.environ.get("EYES_OPEN_EDF", "resting_eyes_open.edf")
EYES_CLOSED_FILE = os.environ.get("EYES_CLOSED_EDF", "resting_eyes_closed.edf")
OCCIPITAL_CHANNELS = ["O1", "Oz", "O2", "PO3", "PO4", "PO7", "PO8"]
results = {}
psds_dict = {}
spectra_dict = {}
for label, filepath in [("eyes_open", EYES_OPEN_FILE), ("eyes_closed", EYES_CLOSED_FILE)]:
ifnot os.path.exists(filepath):
print(f"Skipping {label}: file not found ({filepath})")
continue
raw = load_raw_eeg(filepath)
raw = preprocess_raw(raw, l_freq=1.0, h_freq=50.0, notch_freq=50.0)
raw_clean, _ = run_ica_artifact_removal(raw, n_components=0.99)
# Keep only available occipital channels
occ_available = [c for c in OCCIPITAL_CHANNELS if c in raw_clean.ch_names]
ifnot occ_available:
occ_available = None# Fall back to all EEG
band_df = compute_band_power(raw_clean, picks=occ_available)
band_df["condition"] = label
results[label] = band_df
# Full PSD for plotting
spectrum = raw_clean.compute_psd(method="welch", fmin=1.0, fmax=50.0,
picks=occ_available or"eeg", verbose=False)
spectra_dict[label] = spectrum
# --- Statistical comparison --------------------------------------------------iflen(results) == 2:
from scipy.stats import ttest_rel, wilcoxon
eo = results["eyes_open"]["alpha"].values
ec = results["eyes_closed"]["alpha"].values
iflen(eo) == len(ec):
t_stat, p_val = ttest_rel(ec, eo)
print(f"\nAlpha power comparison (occipital channels):")
print(f" Eyes open: {eo.mean():.2f} ± {eo.std():.2f} µV²/Hz")
print(f" Eyes closed: {ec.mean():.2f} ± {ec.std():.2f} µV²/Hz")
print(f" Paired t-test: t={t_stat:.3f}, p={p_val:.4f}")
# --- PSD plot comparison -----------------------------------------------------if spectra_dict:
fig, ax = plt.subplots(figsize=(10, 6))
colors = {"eyes_open": "steelblue", "eyes_closed": "firebrick"}
for label, spectrum in spectra_dict.items():
psds, freqs = spectrum.get_data(return_freqs=True)
mean_psd = psds.mean(axis=0) * 1e12# µV²/Hz
se_psd = psds.std(axis=0) / np.sqrt(len(psds)) * 1e12
ax.semilogy(freqs, mean_psd, label=label, color=colors.get(label, "gray"), linewidth=2)
ax.fill_between(
freqs,
mean_psd - se_psd,
mean_psd + se_psd,
alpha=0.2,
color=colors.get(label, "gray"),
)
# Shade alpha band
ax.axvspan(8, 13, alpha=0.12, color="gold", label="Alpha band (8–13 Hz)")
ax.set_xlabel("Frequency (Hz)", fontsize=12)
ax.set_ylabel("Power Spectral Density (µV²/Hz)", fontsize=12)
ax.set_title("Resting-State PSD: Eyes Open vs Eyes Closed (Occipital)", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3, which="both")
fig.tight_layout()
plt.savefig("alpha_power_psd.png", dpi=150)
plt.show()
# --- Time-frequency on resting-state (Eyes Closed) -------------------------# Segment resting state into 2-second pseudo-epochs for TFRif"eyes_closed"in results:
raw_ec = load_raw_eeg(EYES_CLOSED_FILE)
raw_ec = preprocess_raw(raw_ec, l_freq=1.0, h_freq=50.0)
raw_clean_ec, _ = run_ica_artifact_removal(raw_ec, n_components=0.99)
# Create fixed-length epochs (2 s, no overlap)
epochs_rest = mne.make_fixed_length_epochs(
raw_clean_ec, duration=2.0, preload=True, verbose=False
)
epochs_rest.pick_channels(
[c for c in OCCIPITAL_CHANNELS if c in epochs_rest.ch_names] or epochs_rest.ch_names[:4]
)
freqs_tfr = np.arange(4, 30, 1)
n_cycles_tfr = freqs_tfr / 2.0# 0.5 cycles per Hz
power, itc = compute_tfr_morlet(
epochs_rest,
freqs=freqs_tfr,
n_cycles=n_cycles_tfr,
return_itc=True,
decim=4,
)
# Plot average TFR
fig_tfr = power.plot(
picks=[0],
baseline=None,
mode="logratio",
title="Resting-State TFR (Eyes Closed)",
show=False,
)
fig_tfr[0].savefig("resting_tfr_eyes_closed.png", dpi=120)
Notes and Best Practices
File Format Recommendations
Store processed data in .fif format (raw, epochs, evoked) for lossless round-tripping.
Use epochs.save("sub-01_task-oddball-epo.fif", overwrite=True) and reload with
mne.read_epochs("sub-01_task-oddball-epo.fif").
ICA Stability
ICA is sensitive to the high-pass filter. Always high-pass at 1 Hz before fitting ICA, even
if your final analysis uses a lower cutoff (0.1 Hz for ERP). Apply ICA to the original
low-passed raw after fitting on the 1 Hz-filtered version.
Gramfort, A. et al. (2013). MEG and EEG data analysis with MNE-Python. Frontiers in
Neuroscience, 7, 267.
Delorme, A., & Makeig, S. (2004). EEGLAB: an open source toolbox for analysis of single-trial
EEG dynamics. Journal of Neuroscience Methods, 134(1), 9–21.
Makeig, S., et al. (1996). Independent component analysis of electroencephalographic data.
Advances in Neural Information Processing Systems, 8.
