Full EEG/MEG analysis pipeline with MNE-Python: preprocessing, ICA artifact removal, ERP computation, time-frequency analysis, and resting-state power spectral density.
Standardmäßig ist der Prompt ausgewählt, der zuerst die Quelle prüft. Sie können zu einem direkten Befehl wechseln oder eine lokale Kopie herunterladen.
Quelldateien prüfen
Lesen Sie SKILL.md und alle von SkillsMP angezeigten Begleitdateien, bevor Sie sich für eine Installation entscheiden.
Mit Codex oder Claude installieren Kopieren Sie diesen Prompt, fügen Sie ihn in Codex, Claude oder einen anderen Assistant ein und lassen Sie die Skill-Seite prüfen und installieren.
Ein direkter Befehl überspringt den Prüf-Prompt. Prüfen Sie die Quelle, bevor Sie ihn ausführen.
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.
