Na Prática: Análise de EEG com MNE-Python

pip install mne mne-bids scipy scikit-learn matplotlib

import mne
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.model_selection import cross_val_score
from sklearn.pipeline import Pipeline
print(f"MNE version: {mne.__version__}")

# The MNE sample dataset (MEG/EEG data from a visual-auditory experiment)
# Download on first run (~1.5 GB)
sample_data_folder = mne.datasets.sample.data_path()
raw_fname = sample_data_folder / "MEG" / "sample" / "sample_audvis_raw.fif"
raw = mne.io.read_raw_fif(raw_fname, preload=True)

print(raw.info)
print(f"Duration: {raw.times[-1]:.1f} s")
print(f"Sampling rate: {raw.info['sfreq']} Hz")
print(f"Channel types: {raw.get_channel_types()}")

# EDF format (common clinical EEG, Biosemi)
# raw = mne.io.read_raw_edf("recording.edf", preload=True)

# BrainVision (.vhdr + .vmrk + .eeg)
# raw = mne.io.read_raw_brainvision("recording.vhdr", preload=True)

# EEGLAB (.set)
# raw = mne.io.read_raw_eeglab("recording.set", preload=True)

# CNT format (Neuroscan)
# raw = mne.io.read_raw_cnt("recording.cnt", preload=True)

# CSV / NumPy array (manual import)
# data: (n_channels, n_samples) float array
# ch_names = ['Fz', 'Cz', 'Pz', ...]
# info = mne.create_info(ch_names=ch_names, sfreq=256, ch_types='eeg')
# raw = mne.io.RawArray(data, info)

# Keep only EEG channels
raw.pick_types(meg=False, eeg=True, eog=True, stim=True)

# Set electrode positions (montage)
# Standard 10-20 system
montage = mne.channels.make_standard_montage("standard_1020")
# raw.set_montage(montage)  # uncomment if EEG channels match standard names

print(f"EEG channels: {raw.ch_names[:10]}...")

# High-pass filter to remove DC drift and slow artifacts
raw.filter(l_freq=1.0, h_freq=None, method="fir", fir_design="firwin")

# Notch filter to remove line noise
raw.notch_filter(freqs=60)  # 60 Hz for North America, 50 Hz for Europe

# Band-pass for specific analyses (e.g., ERP: 0.1–40 Hz)
raw_erp = raw.copy().filter(l_freq=0.1, h_freq=40.0)

# Band-pass for oscillatory analysis (e.g., alpha: 8–12 Hz)
raw_alpha = raw.copy().filter(l_freq=8.0, h_freq=12.0)

# Average reference (recommended for EEG)
raw.set_eeg_reference("average", projection=True)
raw.apply_proj()

# Or reference to specific electrode(s)
# raw.set_eeg_reference(ref_channels=['TP9', 'TP10'])

# Linked mastoids reference
# raw.set_eeg_reference(ref_channels=['M1', 'M2'])

# Fit ICA on band-pass filtered data (1–40 Hz for ICA)
raw_for_ica = raw.copy().filter(1.0, 40.0)

ica = mne.preprocessing.ICA(n_components=20, random_state=42, max_iter="auto")
ica.fit(raw_for_ica)

# Automatic EOG artifact detection
# Find ICA components correlated with EOG channels
eog_indices, eog_scores = ica.find_bads_eog(raw)
print(f"EOG-related components: {eog_indices}")

# Automatic ECG artifact detection (if ECG channel present)
# ecg_indices, ecg_scores = ica.find_bads_ecg(raw)

# Visualize ICA components (in interactive session)
# ica.plot_components()
# ica.plot_scores(eog_scores)

# Exclude artifact components
ica.exclude = eog_indices
ica.apply(raw)  # remove artifacts from raw

print(f"Excluded {len(ica.exclude)} ICA components")

# Read events from stimulus channel (STI channel)
events = mne.find_events(raw, stim_channel="STI 014")
print(f"Found {len(events)} events")
print(f"Event IDs: {np.unique(events[:, 2])}")

# Define event names
event_id = {
    "auditory/left": 1,
    "auditory/right": 2,
    "visual/left": 3,
    "visual/right": 4,
}

# Create epochs around events
epochs = mne.Epochs(
    raw,
    events=events,
    event_id=event_id,
    tmin=-0.2,   # 200 ms before stimulus
    tmax=0.5,    # 500 ms after stimulus
    baseline=(-0.2, 0),   # baseline correction using pre-stimulus period
    reject={"eeg": 150e-6},  # reject epochs with amplitude > 150 μV
    preload=True
)

print(f"Epochs after rejection: {len(epochs)}")
print(f"Epoch duration: {epochs.times[0]:.2f} to {epochs.times[-1]:.2f} s")

# Compute evoked response (average across epochs)
evoked_aud = epochs["auditory/left"].average()
evoked_vis = epochs["visual/left"].average()

# Plot the ERP
fig, axes = plt.subplots(2, 1, figsize=(12, 8))

# Plot single channel ERP
channel = "EEG 060"  # Oz for visual, or use channel_idx
t = evoked_vis.times * 1000  # convert to ms

if channel in evoked_vis.ch_names:
    ch_idx = evoked_vis.ch_names.index(channel)
    axes[0].plot(t, evoked_vis.data[ch_idx] * 1e6, 'b-', label='Visual', linewidth=2)
    axes[0].plot(t, evoked_aud.data[ch_idx] * 1e6, 'r-', label='Auditory', linewidth=2)
    axes[0].axvline(x=0, color='k', linestyle='--', alpha=0.5, label='Stimulus onset')
    axes[0].axhline(y=0, color='k', alpha=0.3)
    axes[0].set_xlabel('Time (ms)')
    axes[0].set_ylabel('Amplitude (μV)')
    axes[0].set_title(f'ERP at {channel}')
    axes[0].legend()
    axes[0].grid(alpha=0.3)

# Plot all channels as butterfly plot
for ch_data in evoked_vis.data * 1e6:
    axes[1].plot(t, ch_data, 'b-', alpha=0.1, linewidth=0.5)
axes[1].axvline(x=0, color='k', linestyle='--', alpha=0.7)
axes[1].set_xlabel('Time (ms)')
axes[1].set_ylabel('Amplitude (μV)')
axes[1].set_title('Butterfly plot — Visual ERP (all channels)')
axes[1].grid(alpha=0.3)

plt.tight_layout()
plt.savefig("erp_comparison.png", dpi=150)
plt.show()

from mne.time_frequency import psd_welch  # or psd_array_welch

# Compute PSD for the continuous signal
# spectrum = raw.compute_psd(method='welch', fmin=1, fmax=50, picks='eeg')
# psd, freqs = spectrum.get_data(return_freqs=True)

# Manual Welch PSD
def compute_channel_psd(data: np.ndarray, sfreq: float, fmin: float = 1, fmax: float = 50) -> tuple:
    freqs, psd = signal.welch(
        data,
        fs=sfreq,
        nperseg=int(sfreq * 4),  # 4-second windows
        noverlap=int(sfreq * 2)  # 50% overlap
    )
    mask = (freqs >= fmin) & (freqs <= fmax)
    return freqs[mask], psd[..., mask]

# Get EEG data for analysis
eeg_picks = mne.pick_types(raw.info, eeg=True)
data = raw.get_data(picks=eeg_picks)
sfreq = raw.info["sfreq"]

freqs, psd = compute_channel_psd(data, sfreq)

# Plot mean PSD across channels
fig, ax = plt.subplots(figsize=(10, 5))
mean_psd = 10 * np.log10(psd.mean(axis=0))  # convert to dB
ax.plot(freqs, mean_psd, 'b-', linewidth=2)
ax.set_xlabel('Frequency (Hz)')
ax.set_ylabel('Power Spectral Density (dB/Hz)')
ax.set_title('Mean EEG PSD')
ax.axvspan(8, 12, alpha=0.2, color='green', label='Alpha (8–12 Hz)')
ax.axvspan(13, 30, alpha=0.2, color='blue', label='Beta (13–30 Hz)')
ax.legend()
ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig("psd.png", dpi=150)
plt.show()

def bandpower(psd: np.ndarray, freqs: np.ndarray, band: tuple) -> np.ndarray:
    """Compute average power in a frequency band. Returns per-channel values."""
    idx = (freqs >= band[0]) & (freqs <= band[1])
    return psd[:, idx].mean(axis=1)

bands = {
    "delta": (1, 4),
    "theta": (4, 8),
    "alpha": (8, 12),
    "beta": (13, 30),
    "gamma": (30, 50),
}

print("Mean bandpower across all EEG channels:")
for band_name, band_range in bands.items():
    power = bandpower(psd, freqs, band_range)
    mean_db = 10 * np.log10(power.mean())
    print(f"  {band_name:8s} ({band_range[0]:2d}–{band_range[1]:2d} Hz): {mean_db:.1f} dB")

from mne.decoding import CSP
from sklearn.pipeline import Pipeline
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis

# Load the motor imagery dataset (PhysioNet EEG Motor Movement/Imagery)
# For this tutorial, we'll simulate motor imagery epochs
np.random.seed(42)

def simulate_motor_imagery_epochs(n_trials_per_class: int = 60, n_channels: int = 64, sfreq: float = 160, tmin: float = 0, tmax: float = 4.0) -> mne.EpochsArray:
    """Simulate left vs. right hand motor imagery epochs with alpha/beta ERD."""
    n_times = int((tmax - tmin) * sfreq)
    n_classes = 2
    n_trials = n_trials_per_class * n_classes

    # Simulate EEG: mostly noise + lateralized alpha/beta suppression
    data = np.random.randn(n_trials, n_channels, n_times) * 5e-6  # ~5 μV background noise

    # Left hand (class 0): ERD over right hemisphere (channels 32-48)
    # Right hand (class 1): ERD over left hemisphere (channels 16-32)
    t = np.linspace(tmin, tmax, n_times)
    alpha_beta = np.sin(2 * np.pi * 10 * t)  # 10 Hz oscillation

    for trial in range(n_trials_per_class):
        # Left hand: suppress right hemisphere alpha/beta
        data[trial, 32:48, :] *= 0.3 * (1 + alpha_beta)
        # Right hand: suppress left hemisphere alpha/beta
        data[n_trials_per_class + trial, 16:32, :] *= 0.3 * (1 + alpha_beta)

    # Create MNE EpochsArray
    ch_names = [f"CH{i:03d}" for i in range(n_channels)]
    info = mne.create_info(ch_names=ch_names, sfreq=sfreq, ch_types="eeg")
    events = np.column_stack([
        np.arange(n_trials) * int(sfreq * 6),  # event times
        np.zeros(n_trials, dtype=int),
        np.array([1] * n_trials_per_class + [2] * n_trials_per_class)  # labels
    ])
    event_id = {"left": 1, "right": 2}
    epochs = mne.EpochsArray(data, info, events=events, tmin=tmin, event_id=event_id)
    return epochs

epochs_mi = simulate_motor_imagery_epochs(n_trials_per_class=60)
print(f"Simulated epochs: {len(epochs_mi)} trials, {epochs_mi.get_data().shape}")

# Build CSP + LDA pipeline
csp = CSP(n_components=6, reg=0.1, log=True, norm_trace=False)
lda = LinearDiscriminantAnalysis()

pipeline = Pipeline([("CSP", csp), ("LDA", lda)])

# Labels: 0 = left, 1 = right
X = epochs_mi.get_data()
y = (epochs_mi.events[:, 2] - 1)  # convert to 0/1

# Cross-validate
scores = cross_val_score(pipeline, X, y, cv=5, scoring="accuracy")
print(f"\nCSP + LDA cross-validation:")
print(f"  Accuracy: {scores.mean():.1%} ± {scores.std():.1%}")
print(f"  (Chance: 50%)")

from mne.connectivity import spectral_connectivity_time

def compute_alpha_coherence(epochs: mne.Epochs, ch_pairs: list = None) -> np.ndarray:
    """
    Compute alpha-band coherence between channel pairs.
    Returns coherence matrix (n_channels × n_channels).
    """
    data = epochs.get_data()  # (n_epochs, n_channels, n_times)
    sfreq = epochs.info["sfreq"]
    n_channels = data.shape[1]
    n_times = data.shape[2]

    # Band-pass filter to alpha
    from scipy.signal import butter, filtfilt
    b, a = butter(4, [8, 12], btype='band', fs=sfreq)
    data_alpha = filtfilt(b, a, data, axis=2)

    # Compute pairwise coherence via cross-spectral density
    coh_matrix = np.zeros((n_channels, n_channels))

    # Simplified: use Pearson correlation in alpha band as proxy
    data_mean = data_alpha.mean(axis=0)  # average over epochs: (n_channels, n_times)
    coh_matrix = np.corrcoef(data_mean)

    return coh_matrix

coh_matrix = compute_alpha_coherence(epochs_mi)

fig, ax = plt.subplots(figsize=(8, 7))
im = ax.imshow(coh_matrix, cmap='RdYlGn', vmin=-1, vmax=1)
ax.set_title('Alpha-band Coherence Matrix')
ax.set_xlabel('Channel index')
ax.set_ylabel('Channel index')
plt.colorbar(im, ax=ax)
plt.tight_layout()
plt.savefig("coherence_matrix.png", dpi=150)
plt.show()

Raw data (mne.io.read_raw_*) 
    → filter (raw.filter)
    → reference (raw.set_eeg_reference)
    → ICA artifact removal (mne.preprocessing.ICA)
    → event detection (mne.find_events)
    → epoching (mne.Epochs)
    → ERP analysis (epochs.average())
    → time-frequency (mne.time_frequency)
    → decoding (mne.decoding.CSP + sklearn)