adipu/ECE374N
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
b586f0093c3404e19ede1544eb1770b373f585fc
ECE374N/Final Project/Group 2 - Glove/motor_cortex_analysis.ipynb

1210 lines
1.5 MiB
Plaintext
Raw Normal View History

final proj
2026-04-20 14:22:05 -05:00
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Motor Cortex Activity Analysis — Group 2 (Glove)\n",
"**Visualizations:**\n",
"1. **Topoplot of MAV features** — spatial distribution of signal amplitude during MI vs Rest\n",
"2. **Frequency spectrum (color bands)** — power in delta/theta/mu/beta/gamma bands during MI vs Rest"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "578c9128",
"metadata": {},
"outputs": [],
"source": [
"# Install dependencies if needed\n",
"# !pip install pyxdf mne scipy numpy matplotlib"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "857b22c0",
"metadata": {},
"outputs": [],
"source": [
"import os\n",
"import glob\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import matplotlib.patches as mpatches\n",
"import matplotlib.colors as mcolors\n",
"from matplotlib.gridspec import GridSpec\n",
"from matplotlib.colorbar import ColorbarBase\n",
"import mne\n",
"import pyxdf\n",
"from scipy.signal import welch, butter, filtfilt\n",
"from scipy.ndimage import gaussian_filter1d\n",
"\n",
"mne.set_log_level('WARNING')\n",
"plt.rcParams.update({'font.size': 11, 'figure.dpi': 120})"
]
},
{
"cell_type": "markdown",
"id": "fe68bf0e",
"metadata": {},
"source": [
"## Configuration"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "dc4b2c55",
"metadata": {},
"outputs": [],
"source": [
"DATA_DIR = os.path.dirname(os.path.abspath('__file__')) # same folder as notebook\n",
"\n",
"# Use OFFLINE sessions (most complete, hardcoded timings)\n",
"# Change to None to process ALL sessions\n",
"SESSION_FILTER = 'OFFLINE'\n",
"\n",
"# Trigger codes\n",
"MI_BEGIN = 200\n",
"MI_END = 220\n",
"REST_BEGIN = 100\n",
"REST_END = 120\n",
"TARGET_MARKERS = [100, 120, 200, 220]\n",
"\n",
"# Epoch window: t=0 is marker timestamp; baseline [-1, 0], analyze [0, 5]\n",
"T_PRE = -1.0 # seconds before marker (for baseline)\n",
"T_POST = 5.0 # seconds after marker\n",
"BASELINE = (-1.0, 0.0)\n",
"\n",
"# Channels to drop\n",
"NON_EEG = {'AUX1', 'AUX2', 'AUX3', 'AUX7', 'AUX8', 'AUX9', 'TRIGGER'}\n",
"\n",
"# MNE channel name rename (XDF uses uppercase, MNE standard_1020 uses mixed case)\n",
"RENAME = {\n",
" 'FP1': 'Fp1', 'FPZ': 'Fpz', 'FP2': 'Fp2',\n",
" 'FZ': 'Fz', 'CZ': 'Cz', 'PZ': 'Pz',\n",
" 'POZ': 'POz', 'OZ': 'Oz'\n",
"}\n",
"\n",
"# Frequency bands for spectrum plot\n",
"BANDS = {\n",
" 'Delta\\n(0.5–4 Hz)': (0.5, 4, '#4B6FA5'),\n",
" 'Theta\\n(4–8 Hz)': (4, 8, '#5BA85B'),\n",
" 'Mu/Alpha\\n(8–13 Hz)':(8, 13, '#C07B3A'),\n",
" 'Beta\\n(13–30 Hz)': (13, 30, '#A855A8'),\n",
" 'Gamma\\n(30–50 Hz)': (30, 50, '#CC4444'),\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "21a40df3",
"metadata": {},
"source": [
"## Helper Functions"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "e798b039",
"metadata": {},
"outputs": [],
"source": [
"def get_channel_names_from_xdf(eeg_stream):\n",
" \"\"\"Extract channel names from XDF stream metadata.\"\"\"\n",
" ch_desc = eeg_stream['info']['desc'][0]\n",
" channels = ch_desc.get('channels', [{}])[0].get('channel', [])\n",
" return [ch['label'][0] for ch in channels]\n",
"\n",
"\n",
"def load_xdf_file(filepath):\n",
" \"\"\"\n",
" Load one XDF file and return:\n",
" eeg_data : (n_eeg_channels, n_samples) float64\n",
" eeg_timestamps: (n_samples,)\n",
" marker_data : (n_markers,) int — marker codes\n",
" marker_ts : (n_markers,) float — marker timestamps\n",
" channel_names : list[str]\n",
" sfreq : float\n",
" \"\"\"\n",
" streams, _ = pyxdf.load_xdf(filepath)\n",
"\n",
" eeg_stream = marker_stream = None\n",
" for s in streams:\n",
" stype = s['info']['type'][0].lower()\n",
" if stype == 'eeg':\n",
" eeg_stream = s\n",
" elif stype == 'markers':\n",
" marker_stream = s\n",
"\n",
" if eeg_stream is None or marker_stream is None:\n",
" # Fallback: first stream is EEG, second is markers\n",
" eeg_stream = streams[0]\n",
" marker_stream = streams[1] if len(streams) > 1 else None\n",
"\n",
" # EEG\n",
" eeg_timestamps = np.array(eeg_stream['time_stamps'])\n",
" eeg_data = np.array(eeg_stream['time_series']).T # (n_channels, n_samples)\n",
" channel_names = get_channel_names_from_xdf(eeg_stream)\n",
" sfreq = float(eeg_stream['info']['nominal_srate'][0])\n",
"\n",
" # Drop non-EEG channels\n",
" valid_idx = [i for i, ch in enumerate(channel_names) if ch not in NON_EEG]\n",
" channel_names = [channel_names[i] for i in valid_idx]\n",
" eeg_data = eeg_data[valid_idx, :]\n",
"\n",
" # Rename for MNE montage compatibility\n",
" channel_names = [RENAME.get(ch, ch) for ch in channel_names]\n",
"\n",
" # Markers\n",
" marker_data = np.array([int(m[0]) for m in marker_stream['time_series']])\n",
" marker_ts = np.array([float(v[1]) for v in marker_stream['time_series']])\n",
"\n",
" # Keep only target markers\n",
" keep = np.isin(marker_data, TARGET_MARKERS)\n",
" marker_data = marker_data[keep]\n",
" marker_ts = marker_ts[keep]\n",
"\n",
" return eeg_data, eeg_timestamps, marker_data, marker_ts, channel_names, sfreq\n",
"\n",
"\n",
"def extract_epochs(eeg_data, eeg_ts, marker_data, marker_ts, sfreq,\n",
" begin_code, end_code, t_pre=T_PRE, t_post=T_POST):\n",
" \"\"\"\n",
" Epoch EEG around BEGIN markers, ending at the paired END marker.\n",
" Returns list of arrays: each (n_channels, epoch_samples).\n",
" \"\"\"\n",
" epochs = []\n",
" n_pre = int(abs(t_pre) * sfreq)\n",
" n_post = int(t_post * sfreq)\n",
"\n",
" begin_idxs = np.where(marker_data == begin_code)[0]\n",
"\n",
" for bi in begin_idxs:\n",
" t_start = marker_ts[bi]\n",
"\n",
" # Find the paired end marker\n",
" future = np.where((marker_data == end_code) & (marker_ts > t_start))[0]\n",
" if len(future) == 0:\n",
" continue\n",
" t_end = marker_ts[future[0]]\n",
"\n",
" # EEG sample indices\n",
" i0 = np.searchsorted(eeg_ts, t_start + t_pre)\n",
" i1 = np.searchsorted(eeg_ts, t_start + t_post)\n",
"\n",
" if i0 < 0 or i1 > eeg_data.shape[1]:\n",
" continue\n",
"\n",
" epoch = eeg_data[:, i0:i1].copy()\n",
"\n",
" # Baseline correction: subtract mean of baseline window [t_pre, 0]\n",
" n_baseline = int(abs(t_pre) * sfreq)\n",
" if epoch.shape[1] > n_baseline:\n",
" baseline_mean = epoch[:, :n_baseline].mean(axis=1, keepdims=True)\n",
" epoch -= baseline_mean\n",
"\n",
" epochs.append(epoch)\n",
"\n",
" return epochs\n",
"\n",
"\n",
"def bandpass(data, lo, hi, sfreq, order=4):\n",
" nyq = sfreq / 2.0\n",
" lo = max(lo, 0.5) # avoid 0 Hz\n",
" hi = min(hi, nyq - 0.1)\n",
" b, a = butter(order, [lo / nyq, hi / nyq], btype='band')\n",
" return filtfilt(b, a, data, axis=-1)"
]
},
{
"cell_type": "markdown",
"id": "98d225db",
"metadata": {},
"source": [
"## Load Data"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "d266216b",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Found 6 file(s) to process:\n",
" sub-S26CLASS_SUBJ_003_ses-S001OFFLINE_FES_task-Default_run-001_eeg.xdf\n",
" sub-S26CLASS_SUBJ_003_ses-S004OFFLINE_NOFES_task-Default_run-001_eeg.xdf\n",
" sub-S26CLASS_SUBJ_005_ses-S001OFFLINE_FES_task-Default_run-001_eeg.xdf\n",
" sub-S26CLASS_SUBJ_005_ses-S004OFFLINE_NOFES_task-Default_run-001_eeg.xdf\n",
" sub-S26CLASS_SUBJ_009_ses-S001OFFLINE_FES_task-Default_run-001_eeg.xdf\n",
" sub-S26CLASS_SUBJ_009_ses-S004OFFLINE_NOFES_task-Default_run-001_eeg.xdf\n",
"\n",
"Loading sub-S26CLASS_SUBJ_003_ses-S001OFFLINE_FES_task-Default_run-001_eeg.xdf ... 45 MI, 45 REST epochs\n",
"\n",
"Loading sub-S26CLASS_SUBJ_003_ses-S004OFFLINE_NOFES_task-Default_run-001_eeg.xdf ... 45 MI, 45 REST epochs\n",
"\n",
"Loading sub-S26CLASS_SUBJ_005_ses-S001OFFLINE_FES_task-Default_run-001_eeg.xdf ... 45 MI, 45 REST epochs\n",
"\n",
"Loading sub-S26CLASS_SUBJ_005_ses-S004OFFLINE_NOFES_task-Default_run-001_eeg.xdf ... 45 MI, 45 REST epochs\n",
"\n",
"Loading sub-S26CLASS_SUBJ_009_ses-S001OFFLINE_FES_task-Default_run-001_eeg.xdf ... 45 MI, 45 REST epochs\n",
"\n",
"Loading sub-S26CLASS_SUBJ_009_ses-S004OFFLINE_NOFES_task-Default_run-001_eeg.xdf ... 45 MI, 45 REST epochs\n",
"\n",
"Total: 270 MI epochs, 270 REST epochs\n",
"Channels (32): ['Fp1', 'Fpz', 'Fp2', 'F7', 'F3', 'Fz', 'F4', 'F8', 'FC5', 'FC1', 'FC2', 'FC6', 'M1', 'T7', 'C3', 'Cz', 'C4', 'T8', 'M2', 'CP5', 'CP1', 'CP2', 'CP6', 'P7', 'P3', 'Pz', 'P4', 'P8', 'POz', 'O1', 'Oz', 'O2']\n",
"Sampling rate: 512.0 Hz\n"
]
}
],
"source": [
"xdf_files = sorted(glob.glob(os.path.join(DATA_DIR, '*.xdf')))\n",
"\n",
"if SESSION_FILTER:\n",
" xdf_files = [f for f in xdf_files if SESSION_FILTER in os.path.basename(f)]\n",
"\n",
"print(f'Found {len(xdf_files)} file(s) to process:')\n",
"for f in xdf_files:\n",
" print(' ', os.path.basename(f))\n",
"\n",
"all_mi_epochs = [] # list of (n_channels, n_samples) arrays\n",
"all_rest_epochs = []\n",
"channel_names_global = None\n",
"sfreq_global = None\n",
"\n",
"for filepath in xdf_files:\n",
" print(f'\\nLoading {os.path.basename(filepath)} ...', end=' ', flush=True)\n",
" try:\n",
" eeg_data, eeg_ts, marker_data, marker_ts, ch_names, sfreq = load_xdf_file(filepath)\n",
" except Exception as e:\n",
" print(f'ERROR: {e}')\n",
" continue\n",
"\n",
" if channel_names_global is None:\n",
" channel_names_global = ch_names\n",
" sfreq_global = sfreq\n",
"\n",
" mi_epochs = extract_epochs(eeg_data, eeg_ts, marker_data, marker_ts, sfreq,\n",
" MI_BEGIN, MI_END)\n",
" rest_epochs = extract_epochs(eeg_data, eeg_ts, marker_data, marker_ts, sfreq,\n",
" REST_BEGIN, REST_END)\n",
"\n",
" all_mi_epochs.extend(mi_epochs)\n",
" all_rest_epochs.extend(rest_epochs)\n",
" print(f'{len(mi_epochs)} MI, {len(rest_epochs)} REST epochs')\n",
"\n",
"print(f'\\nTotal: {len(all_mi_epochs)} MI epochs, {len(all_rest_epochs)} REST epochs')\n",
"print(f'Channels ({len(channel_names_global)}): {channel_names_global}')\n",
"print(f'Sampling rate: {sfreq_global} Hz')"
]
},
{
"cell_type": "markdown",
"id": "7b8c8bea",
"metadata": {},
"source": [
"## Set Up MNE Montage"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "611baf23",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Channels with known montage positions (32): ['Fp1', 'Fpz', 'Fp2', 'F7', 'F3', 'Fz', 'F4', 'F8', 'FC5', 'FC1', 'FC2', 'FC6', 'M1', 'T7', 'C3', 'Cz', 'C4', 'T8', 'M2', 'CP5', 'CP1', 'CP2', 'CP6', 'P7', 'P3', 'Pz', 'P4', 'P8', 'POz', 'O1', 'Oz', 'O2']\n"
]
},
{
"data": {
"text/html": [
"<script type=\"text/javascript\">\n",
" // must be `var` (not `const`) because this can get embedded multiple times on a page\n",
"var toggleVisibility = (className) => {\n",
"\n",
" const elements = document.querySelectorAll(`.${className}`);\n",
"\n",
" elements.forEach(element => {\n",
" if (element.classList.contains(\"mne-repr-section-header\")) {\n",
" return // Don't collapse the section header row\n",
" }\n",
" element.classList.toggle(\"mne-repr-collapsed\");\n",
" });\n",
"\n",
" // trigger caret to rotate\n",
" var sel = `.mne-repr-section-header.${className} > th.mne-repr-section-toggle > button`;\n",
" const button = document.querySelector(sel);\n",
" button.classList.toggle(\"collapsed\");\n",
"\n",
" // adjust tooltip\n",
" sel = `tr.mne-repr-section-header.${className}`;\n",
" const secHeadRow = document.querySelector(sel);\n",
" secHeadRow.classList.toggle(\"collapsed\");\n",
" secHeadRow.title = secHeadRow.title === \"Hide section\" ? \"Show section\" : \"Hide section\";\n",
"}\n",
"</script>\n",
"\n",
"<style type=\"text/css\">\n",
" /*\n",
"Styles in this section apply both to the sphinx-built website docs and to notebooks\n",
"rendered in an IDE or in Jupyter. In our web docs, styles here are complemented by\n",
"doc/_static/styles.css and other CSS files (e.g. from the sphinx theme, sphinx-gallery,\n",
"or bootstrap). In IDEs/Jupyter, those style files are unavailable, so only the rules in\n",
"this file apply (plus whatever default styling the IDE applies).\n",
"*/\n",
".mne-repr-table {\n",
" display: inline; /* prevent using full container width */\n",
"}\n",
".mne-repr-table tr.mne-repr-section-header > th {\n",
" padding-top: 1rem;\n",
" text-align: left;\n",
" vertical-align: middle;\n",
"}\n",
".mne-repr-section-toggle > button {\n",
" all: unset;\n",
" display: block;\n",
" height: 1rem;\n",
" width: 1rem;\n",
"}\n",
".mne-repr-section-toggle > button > svg {\n",
" height: 60%;\n",
"}\n",
"\n",
"/* transition (rotation) effects on the collapser button */\n",
".mne-repr-section-toggle > button.collapsed > svg {\n",
" transition: 0.1s ease-out;\n",
" transform: rotate(-90deg);\n",
"}\n",
".mne-repr-section-toggle > button:not(.collapsed) > svg {\n",
" transition: 0.1s ease-out;\n",
" transform: rotate(0deg);\n",
"}\n",
"\n",
"/* hide collapsed table rows */\n",
".mne-repr-collapsed {\n",
" display: none;\n",
"}\n",
"\n",
"\n",
"@layer {\n",
" /*\n",
" Selectors in a `@layer` will always be lower-precedence than selectors outside the\n",
" layer. So even though e.g. `div.output_html` is present in the sphinx-rendered\n",
" website docs, the styles here won't take effect there as long as some other rule\n",
" somewhere in the page's CSS targets the same element.\n",
"\n",
" In IDEs or Jupyter notebooks, though, the CSS files from the sphinx theme,\n",
" sphinx-gallery, and bootstrap are unavailable, so these styles will apply.\n",
"\n",
" Notes:\n",
"\n",
" - the selector `.accordion-body` is for MNE Reports\n",
" - the selector `.output_html` is for VSCode's notebook interface\n",
" - the selector `.jp-RenderedHTML` is for Jupyter notebook\n",
" - variables starting with `--theme-` are VSCode-specific.\n",
" - variables starting with `--jp-` are Jupyter styles, *some of which* are also\n",
" available in VSCode. Here we try the `--theme-` variable first, then fall back to\n",
" the `--jp-` ones.\n",
" */\n",
" .mne-repr-table {\n",
" --mne-toggle-color: var(--theme-foreground, var(--jp-ui-font-color1));\n",
" --mne-button-bg-color: var(--theme-button-background, var(--jp-info-color0, var(--jp-content-link-color)));\n",
" --mne-button-fg-color: var(--theme-button-foreground, var(--jp-ui-inverse-font-color0, var(--jp-editor-background)));\n",
" --mne-button-hover-bg-color: var(--theme-button-hover-background, var(--jp-info-color1));\n",
" --mne-button-radius: var(--jp-border-radius, 0.25rem);\n",
" }\n",
" /* chevron position/alignment; in VSCode it looks ok without adjusting */\n",
" .accordion-body .mne-repr-section-toggle > button,\n",
" .jp-RenderedHTML .mne-repr-section-toggle > button {\n",
" padding: 0 0 45% 25% !important;\n",
" }\n",
" /* chevron color; MNE Report doesn't have light/dark mode */\n",
" div.output_html .mne-repr-section-toggle > button > svg > path,\n",
" .jp-RenderedHTML .mne-repr-section-toggle > button > svg > path {\n",
" fill: var(--mne-toggle-color);\n",
" }\n",
" .accordion-body .mne-ch-names-btn,\n",
" div.output_html .mne-ch-names-btn,\n",
" .jp-RenderedHTML .mne-ch-names-btn {\n",
" -webkit-border-radius: var(--mne-button-radius);\n",
" -moz-border-radius: var(--mne-button-radius);\n",
" border-radius: var(--mne-button-radius);\n",
" border: none;\n",
" background-image: none;\n",
" background-color: var(--mne-button-bg-color);\n",
" color: var(--mne-button-fg-color);\n",
" font-size: inherit;\n",
" min-width: 1.5rem;\n",
" padding: 0.25rem;\n",
" text-align: center;\n",
" text-decoration: none;\n",
" }\n",
" .accordion-body .mne-ch-names-btn:hover,\n",
" div.output_html .mne.ch-names-btn:hover,\n",
" .jp-RenderedHTML .mne-ch-names-btn:hover {\n",
" background-color: var(--mne-button-hover-bg-color);\n",
" text-decoration: underline;\n",
" }\n",
" .accordion-body .mne-ch-names-btn:focus-visible,\n",
" div.output_html .mne-ch-names-btn:focus-visible,\n",
" .jp-RenderedHTML .mne-ch-names-btn:focus-visible {\n",
" outline: 0.1875rem solid var(--mne-button-bg-color) !important;\n",
" outline-offset: 0.1875rem !important;\n",
" }\n",
"}\n",
"</style>\n",
"\n",
"\n",
"\n",
"<table class=\"table mne-repr-table\">\n",
" \n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"<tr class=\"mne-repr-section-header general-dfa4643f-fed2-4327-bef1-314ddb73b3f5\"\n",
" title=\"Hide section\" \n",
" onclick=\"toggleVisibility('general-dfa4643f-fed2-4327-bef1-314ddb73b3f5')\">\n",
" <th class=\"mne-repr-section-toggle\">\n",
" <button >\n",
" <svg xmlns=\"http://www.w3.org/2000/svg\" viewBox=\"0 0 512 512\"><!--!Font Awesome Free 6.6.0 by @fontawesome - https://fontawesome.com License - https://fontawesome.com/license/free Copyright 2024 Fonticons, Inc.--><path d=\"M233.4 406.6c12.5 12.5 32.8 12.5 45.3 0l192-192c12.5-12.5 12.5-32.8 0-45.3s-32.8-12.5-45.3 0L256 338.7 86.6 169.4c-12.5-12.5-32.8-12.5-45.3 0s-12.5 32.8 0 45.3l192 192z\"/></svg>\n",
" </button>\n",
" </th>\n",
" <th colspan=\"2\">\n",
" <strong>General</strong>\n",
" </th>\n",
"</tr>\n",
"\n",
"\n",
"<tr class=\"repr-element general-dfa4643f-fed2-4327-bef1-314ddb73b3f5 \">\n",
" <td class=\"mne-repr-section-toggle\"></td>\n",
" <td>MNE object type</td>\n",
" <td>Info</td>\n",
"</tr>\n",
"<tr class=\"repr-element general-dfa4643f-fed2-4327-bef1-314ddb73b3f5 \">\n",
" <td class=\"mne-repr-section-toggle\"></td>\n",
" <td>Measurement date</td>\n",
" \n",
" <td>Unknown</td>\n",
" \n",
"</tr>\n",
"<tr class=\"repr-element general-dfa4643f-fed2-4327-bef1-314ddb73b3f5 \">\n",
" <td class=\"mne-repr-section-toggle\"></td>\n",
" <td>Participant</td>\n",
" \n",
" <td>Unknown</td>\n",
" \n",
"</tr>\n",
"<tr class=\"repr-element general-dfa4643f-fed2-4327-bef1-314ddb73b3f5 \">\n",
" <td class=\"mne-repr-section-toggle\"></td>\n",
" <td>Experimenter</td>\n",
" \n",
" <td>Unknown</td>\n",
" \n",
"</tr>\n",
" \n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"<tr class=\"mne-repr-section-header acquisition-8b291e68-ae7a-4787-990f-89dfe8d1e8f8\"\n",
" title=\"Hide section\" \n",
" onclick=\"toggleVisibility('acquisition-8b291e68-ae7a-4787-990f-89dfe8d1e8f8')\">\n",
" <th class=\"mne-repr-section-toggle\">\n",
" <button >\n",
" <svg xmlns=\"http://www.w3.org/2000/svg\" viewBox=\"0 0 512 512\"><!--!Font Awesome Free 6.6.0 by @fontawesome - https://fontawesome.com License - https://fontawesome.com/license/free Copyright 2024 Fonticons, Inc.--><path d=\"M233.4 406.6c12.5 12.5 32.8 12.5 45.3 0l192-192c12.5-12.5 12.5-32.8 0-45.3s-32.8-12.5-45.3 0L256 338.7 86.6 169.4c-12.5-12.5-32.8-12.5-45.3 0s-12.5 32.8 0 45.3l192 192z\"/></svg>\n",
" </button>\n",
" </th>\n",
" <th colspan=\"2\">\n",
" <strong>Acquisition</strong>\n",
" </th>\n",
"</tr>\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"<tr class=\"repr-element acquisition-8b291e68-ae7a-4787-990f-89dfe8d1e8f8 \">\n",
" <td class=\"mne-repr-section-toggle\"></td>\n",
" <td>Sampling frequency</td>\n",
" <td>512.00 Hz</td>\n",
"</tr>\n",
"\n",
"\n",
"\n",
" \n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"<tr class=\"mne-repr-section-header channels-88352e21-0631-4afd-b4c9-99bbff024342\"\n",
" title=\"Hide section\" \n",
" onclick=\"toggleVisibility('channels-88352e21-0631-4afd-b4c9-99bbff024342')\">\n",
" <th class=\"mne-repr-section-toggle\">\n",
" <button >\n",
" <svg xmlns=\"http://www.w3.org/2000/svg\" viewBox=\"0 0 512 512\"><!--!Font Awesome Free 6.6.0 by @fontawesome - https://fontawesome.com License - https://fontawesome.com/license/free Copyright 2024 Fonticons, Inc.--><path d=\"M233.4 406.6c12.5 12.5 32.8 12.5 45.3 0l192-192c12.5-12.5 12.5-32.8 0-45.3s-32.8-12.5-45.3 0L256 338.7 86.6 169.4c-12.5-12.5-32.8-12.5-45.3 0s-12.5 32.8 0 45.3l192 192z\"/></svg>\n",
" </button>\n",
" </th>\n",
" <th colspan=\"2\">\n",
" <strong>Channels</strong>\n",
" </th>\n",
"</tr>\n",
"\n",
"\n",
" \n",
"<tr class=\"repr-element channels-88352e21-0631-4afd-b4c9-99bbff024342 \">\n",
" <td class=\"mne-repr-section-toggle\"></td>\n",
" <td>EEG</td>\n",
" <td>\n",
" <button class=\"mne-ch-names-btn sd-sphinx-override sd-btn sd-btn-info sd-text-wrap sd-shadow-sm\" onclick=\"alert('Good EEG:\\n\\nFp1, Fpz, Fp2, F7, F3, Fz, F4, F8, FC5, FC1, FC2, FC6, M1, T7, C3, Cz, C4, T8, M2, CP5, CP1, CP2, CP6, P7, P3, Pz, P4, P8, POz, O1, Oz, O2')\" title=\"(Click to open in popup)&#13;&#13;Fp1, Fpz, Fp2, F7, F3, Fz, F4, F8, FC5, FC1, FC2, FC6, M1, T7, C3, Cz, C4, T8, M2, CP5, CP1, CP2, CP6, P7, P3, Pz, P4, P8, POz, O1, Oz, O2\">\n",
" 32\n",
" </button>\n",
"\n",
" \n",
" </td>\n",
"</tr>\n",
"\n",
"\n",
"<tr class=\"repr-element channels-88352e21-0631-4afd-b4c9-99bbff024342 \">\n",
" <td class=\"mne-repr-section-toggle\"></td>\n",
" <td>Head & sensor digitization</td>\n",
" \n",
" <td>35 points</td>\n",
" \n",
"</tr>\n",
" \n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"<tr class=\"mne-repr-section-header filters-efa6aa4a-b77b-44b6-af3a-6feabded0c12\"\n",
" title=\"Hide section\" \n",
" onclick=\"toggleVisibility('filters-efa6aa4a-b77b-44b6-af3a-6feabded0c12')\">\n",
" <th class=\"mne-repr-section-toggle\">\n",
" <button >\n",
" <svg xmlns=\"http://www.w3.org/2000/svg\" viewBox=\"0 0 512 512\"><!--!Font Awesome Free 6.6.0 by @fontawesome - https://fontawesome.com License - https://fontawesome.com/license/free Copyright 2024 Fonticons, Inc.--><path d=\"M233.4 406.6c12.5 12.5 32.8 12.5 45.3 0l192-192c12.5-12.5 12.5-32.8 0-45.3s-32.8-12.5-45.3 0L256 338.7 86.6 169.4c-12.5-12.5-32.8-12.5-45.3 0s-12.5 32.8 0 45.3l192 192z\"/></svg>\n",
" </button>\n",
" </th>\n",
" <th colspan=\"2\">\n",
" <strong>Filters</strong>\n",
" </th>\n",
"</tr>\n",
"\n",
"\n",
"<tr class=\"repr-element filters-efa6aa4a-b77b-44b6-af3a-6feabded0c12 \">\n",
" <td class=\"mne-repr-section-toggle\"></td>\n",
" <td>Highpass</td>\n",
" <td>0.00 Hz</td>\n",
"</tr>\n",
"\n",
"\n",
"<tr class=\"repr-element filters-efa6aa4a-b77b-44b6-af3a-6feabded0c12 \">\n",
" <td class=\"mne-repr-section-toggle\"></td>\n",
" <td>Lowpass</td>\n",
" <td>256.00 Hz</td>\n",
"</tr>\n",
"\n",
"\n",
"</table>"
],
"text/plain": [
"<Info | 8 non-empty values\n",
" bads: []\n",
" ch_names: Fp1, Fpz, Fp2, F7, F3, Fz, F4, F8, FC5, FC1, FC2, FC6, M1, T7, ...\n",
" chs: 32 EEG\n",
" custom_ref_applied: False\n",
" dig: 35 items (3 Cardinal, 32 EEG)\n",
" highpass: 0.0 Hz\n",
" lowpass: 256.0 Hz\n",
" meas_date: unspecified\n",
" nchan: 32\n",
" projs: []\n",
" sfreq: 512.0 Hz\n",
">"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"montage = mne.channels.make_standard_montage('standard_1020')\n",
"montage_ch_names = [ch.lower() for ch in montage.ch_names]\n",
"\n",
"# Keep only channels that exist in the standard montage\n",
"valid_ch = [ch for ch in channel_names_global\n",
" if ch.lower() in montage_ch_names]\n",
"\n",
"# Indices into the epoch arrays\n",
"valid_idx = [channel_names_global.index(ch) for ch in valid_ch]\n",
"\n",
"print(f'Channels with known montage positions ({len(valid_ch)}): {valid_ch}')\n",
"\n",
"# Build MNE Info object for plotting\n",
"info = mne.create_info(ch_names=valid_ch, sfreq=sfreq_global, ch_types='eeg')\n",
"info.set_montage(montage)"
]
},
{
"cell_type": "markdown",
"id": "70922abb",
"metadata": {},
"source": [
"## Compute MAV per Channel"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "f5e80da3",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"MAV (µV) summary:\n",
" MI — mean: 33.034 range: [15.140, 87.473]\n",
" REST — mean: 32.833 range: [15.077, 67.931]\n",
"\n",
"If values are ~1e6x too large, change scale = 1.0 → scale = 1e-6 (data may be in nV).\n",
"If values are ~1e-6x too small, change scale = 1.0 → scale = 1e6 (data may be in V).\n"
]
}
],
"source": [
"def compute_mav(epochs, valid_idx, n_pre_samples):\n",
" \"\"\"\n",
" Mean Absolute Value averaged over the active window [0, T_POST]\n",
" (i.e., after the baseline portion).\n",
" Returns: (n_valid_channels,)\n",
" \"\"\"\n",
" mavs = []\n",
" for ep in epochs:\n",
" active = ep[valid_idx, n_pre_samples:] # drop baseline\n",
" mavs.append(np.mean(np.abs(active), axis=1))\n",
" return np.mean(mavs, axis=0) # average across trials\n",
"\n",
"\n",
"n_pre = int(abs(T_PRE) * sfreq_global)\n",
"\n",
"mav_mi = compute_mav(all_mi_epochs, valid_idx, n_pre)\n",
"mav_rest = compute_mav(all_rest_epochs, valid_idx, n_pre)\n",
"\n",
"# scale = 1.0 → data is already in µV\n",
"# scale = 1e6 → data is in Volts (raw from amp)\n",
"# The MAV values below will reveal which: if mean is ~1–100 µV it's already µV\n",
"scale = 1.0\n",
"mav_mi *= scale\n",
"mav_rest *= scale\n",
"\n",
"# Difference map: MI − Rest (positive = more active during MI)\n",
"mav_diff = mav_mi - mav_rest\n",
"\n",
"print('MAV (µV) summary:')\n",
"print(f' MI — mean: {mav_mi.mean():.3f} range: [{mav_mi.min():.3f}, {mav_mi.max():.3f}]')\n",
"print(f' REST — mean: {mav_rest.mean():.3f} range: [{mav_rest.min():.3f}, {mav_rest.max():.3f}]')\n",
"print()\n",
"print('If values are ~1e6x too large, change scale = 1.0 → scale = 1e-6 (data may be in nV).')\n",
"print('If values are ~1e-6x too small, change scale = 1.0 → scale = 1e6 (data may be in V).')"
]
},
{
"cell_type": "markdown",
"id": "2ab81600",
"metadata": {},
"source": [
"---\n",
"## Figure 1 — Topoplots of MAV Features (Motor Cortex Activity)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "d53e63b9",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAACF4AAAMKCAYAAAC25hJXAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAASdAAAEnQB3mYfeAABAABJREFUeJzs3Qd4FFUXBuCTTug19N47UqUJKIioYAMVRQFBEWygqCgqKHbBLiAqqIBY+W0IiDQLTZDee+81IT3Z//luMpM7s7Ob3U0hwPf6rOzu9Lo3d849N8jlcrmEiIiIiIiIiIiIiIiIiIiIiPwW7P8kRERERERERERERERERERERMTACyIiIiIiIiIiIiIiIiIiIqIsYMYLIiIiIiIiIiIiIiIiIiIiogAx8IKIiIiIiIiIiIiIiIiIiIgoQAy8ICIiIiIiIiIiIiIiIiIiIgoQAy+IiIiIiIiIiIiIiIiIiIiIAsTACyIiIiIiIiIiIiIiIiIiIqIAMfCCiIiIiIiIiIiIiIiIiIiIKEAMvCAiIiIiIiIiIiIiIiIiIiIKEAMviIiIiIiIiIiIiIiIiIiIiALEwAsiIiIiIiIiIiIiIiIiIiKiADHwgoiIiIhIExQU5Pb65JNPHPdRSkqKVKpUyW380aNHX/b7NCEhQT7//HO5/fbbpVq1alK4cGEJCwuTqKgoad++vTz33HOyadOmy34/ZUWtWrXczr2vvvrqgu1TfT06duyYLfPEOaTPF58vd1WqVHE77iNHjvQ4ftu2bd3G79evn1xOpk+f7rYPateu7XUa3Mf18RctWiR51aV6nTz00EOW7RozZozX8UeNGmUZ/7HHHsvS8nEf0+d3sbnuuuvMdX/zzTctw3APsF8TuFd4gnuMfXzci7Ljusuu47xixQrzu9KlS8u5c+fkQvj999/Vvi9VqpSEhoaa6zR06NBMp92zZ49jOdTpZd//9vuAt5e334C4uDj57LPP5I477pAaNWpIkSJF1HagHFe3bl255ZZbZNy4cbJ79+5s2V9ERERERJcSBl4QEREREWViwoQJjt//8ssvsn///hzdf/aHI6iUz+vmzJkjVatWlf79+8t3332nKuejo6MlOTlZjh8/Ln///be88sorUr9+fTl79uwFfXCd2YOjvGrZsmWyfft2t++nTp16UQRUXM7XR07AQ7LExES379esWSNLliyRy53TdbFt2zZZvnz5BVmfi4n9QXBuBu3cc889ls+ZBZbZh997771yufrtt99k7ty56n3x4sVlyJAhmU6De8XatWvdvse95dNPP82x6y67jnPLli3l2muvVe+PHTumyhm5bd68edKtWze170+cOKECdC8m3377rQooHjhwoHq/c+dOFcCC7UA5bsuWLfLjjz/K8OHDVVDt+fPnL/QqExERERHlKaEXegWIiIiIiPK61atXy9KlS6V169aW7z/66KMLtk551RdffKECLlwul/ldSEiINGvWTMqUKaMCLfAw2Ai40Mcj33355ZceH/ocOXJE7evcdtttt5nvEVSTHRAYo8/3Yg2UyWlHjx6VH374QXr37m35nvcokcOHD8sff/zh8cFwq1atcukokb+uvPJKqVmzphlkhoe+//33nzRt2tRtXGQ72LFjh/kZLfPxu3M5wu/qU089ZX5G0EXBggV9mnb8+PHy8ccfW75DACUCGXLqusvO4/z000+rjBPw3nvvySOPPCIVKlSQ3AyCS01NNT83aNDAzE7VpEmTgOap/wbqkEHMm5IlS0qHDh0ch7Vo0cLtuzfeeENGjBhh+S44OFgaNmwoFStWVAE4e/fuVcEzRtmNZTgiIiIiIisGXhARERER+fgwQg+82Lp1q8yfP5/7TrN+/Xq5//77LRXxaH2KlrKotDeg5SSyhbzwwgvcfwHAw49vvvnG/IwuXJKSksx9i9bAjz/+eK7v2++//z7b54kMG3kty0ZevkfpgRdnzpy5oF3P5BXYB3qrc/16+frrr+Wdd95R31HehGwI+m8Fuq9weiCP7+3TXa4Q8LBx40bzc9++fX2eFvvxrbfeUt1KZCWAy9/rLruOc6dOnVR5A9nI0OXZxIkT5eWXX5bcDILTrVy5UiIiIi7IbysCIH2dFkGbzzzzjOW77t27ywcffCCVK1e2fI/MZciG8fbbbwe0XkRERERElzJ2NUJERERE5EGxYsUkX758ZotPpI3WH3IaAQblypXzaR8i2AAtF5HGGfNFC1T0d45gBWTV0I0ePVq1kEQGCR268PDW57w/y/DU9QYeViBFd7169SQyMtLnLAPPP/+8+WAFrrjiCrU+etCFkQHj5ptvVi1a9Yc7Ob0NeICP4WixacB7b11pYD4IHOnatavqMz48PFydF23btpV3331X9YWu++effyx9umN+eiDKvn37VH/pxnD0n4703f6YNWuWnDp1yvyMIIsCBQr43N0I1hktmtEHfdmyZdVDoaJFi6qWw0gvvm7dOss+1S1evNhjtwOe9iOOnfE9rhW9NbDhscces0yPVs1Ofdbjsz/XBx6C4RwyPt94442O+wStr41xcHztxzUvM+4/6MIHwU+GKVOmSGxsrGUcX7qwwTFF63OcU7h20FobLeadurYxlnPfffdJ8+bNVctyTIdzCllXrrnmGvXgLj4+3m26RYsWWY4VjunBgwflwQcfVPcMXGt44Pfkk09m6Xjo2WGwbsOGDTM/nzx5UnXJ4Cs8RO3Ro4eUKFFC8ufPr7bZfg4adu3apVrbo8V7oUKF1EPmUqVKqYehd955p3rwjOAYp0wBzz33nGqRjnMR+wEt27Ev8QA8kH2RWXdBnrpfwmdcUzpsr/242R/KvvjiiyqLgbH+uHfi2ps5c6bf644H6/p9CA/t7fcQPODXg9Ew/t13363eIzvCyJEj5frrr1fnNbIA4FjgmOD3Ab8r+C3yl6d7k7fz2w7bgQfYOKfKly+vrhv8JuK8GjNmjJw+fVqy2j0aziP8zmTGuEeg6wj9nDYyfunj5MR1l9XjrH+H68uA32+9XOIrXGco53Xu3FldfziP8TuJY4PzCdepU7dXOO46/ffHPiwvQaYQvazSpUsX1aWIPegCcB976KGHVACyXvYgIiIiIqK0tHBERERERJQORWTjVblyZVffvn3Nz6+//roaJyYmxlWkSBH1XXBwsOuFF16wTDdq1CjL/oyNjXV1797dMo79hfm8+OKL5jSYh7fxjVdWlmHAdhrjlC1b1tWpUye3/ZCZ6OhoV1hYmGW6n3/+2a/zKqe3oUOHDpnuT4xj2LNnj6tRo0Zex69fv75r7969lnXB+unjvPfee+r71NRUV8eOHc3vsb9WrFjh8tctt9ximf/WrVtdd9xxh+W79evXO067YcMGV82aNb1u0zvvvOO2Tz29cH0YPO1HXDf6sD/++MOyTsnJya7SpUubwxs2bGgOmzJlimVafPb3+ujfv7/l/MFx1e3YscMyzcMPP+zKy+zHRd8XDz74oHmu6cd59OjRHo+bMf7QoUO97st8+fK5pk+f7rY+1atXz/Q44Do6ffq0ZbqFCxdaxunWrZurWLFijtNjWCDWrl1rmU/v3r1dmzdvtnx32223OU5rP8ceffRRV2hoqOP6PfLII5ZpN23aZP5GeHv9+++/lul+/fXXTKerU6eOa/v27ZbpPF0nmV2bTueUfr/35RrTf+9wbZcoUcLr+NjfCQkJfh3H9u3bW+Yxf/58y/C5c+dahuP+b5g6dWqm2xASEuKaNGmS23Ltvxn+7HP7+W0vF5w6dcrye+D0qlChgmvNmjV+7au4uDhXRESEOY9nnnnGcTy9bGO/j9StW9ccb8CAAR7vI57KBoFed1k5zro5c+ZYxlu0aJFf+xDXF64zb8cG1ymuV0/70+mFcyIzu3fvdpvOV/Zz0ulad4IyhH2Zq1ev9nm5RERERESUgRkviIiIiIi8QKs+A7IEoAXmtGnT5OzZs+q7bt26ubUItnvggQdUFgcDWpGj3220nDRaeGK+o0aNksmTJ6vPaImLrA/21oZYHr43XllZhhO04ly4cKFqDYzpr776ap9aNK5atcrSqhR
"text/plain": [
"<Figure size 2160x720 with 6 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Saved: topoplot_MAV.png\n"
]
}
],
"source": [
"MOTOR_CH = ['C3', 'Cz', 'C4'] # channels to highlight\n",
"\n",
"fig, axes = plt.subplots(1, 3, figsize=(18, 6))\n",
"fig.suptitle('Motor Cortex Activity — Mean Absolute Value (MAV) of EEG\\n'\n",
" f'Averaged over {len(all_mi_epochs)} MI / {len(all_rest_epochs)} REST trials | '\n",
" f'Analysis window: [0, {T_POST}] s',\n",
" fontsize=13, fontweight='bold', y=1.02)\n",
"\n",
"vmin_abs = min(mav_mi.min(), mav_rest.min())\n",
"vmax_abs = max(mav_mi.max(), mav_rest.max())\n",
"\n",
"plots = [\n",
" (mav_mi, 'MI (Motor Imagery)', 'Oranges', vmin_abs, vmax_abs),\n",
" (mav_rest, 'REST', 'Blues', vmin_abs, vmax_abs),\n",
" (mav_diff, 'Difference (MI − REST)', 'RdBu_r',\n",
" -np.abs(mav_diff).max(), np.abs(mav_diff).max()),\n",
"]\n",
"\n",
"for ax, (data, title, cmap, vmin, vmax) in zip(axes, plots):\n",
" im, cn = mne.viz.plot_topomap(\n",
" data, info,\n",
" axes=ax,\n",
" cmap=cmap,\n",
" vlim=(vmin, vmax),\n",
" show=False,\n",
" contours=6,\n",
" sensors=True, # draw all sensor dots\n",
" names=valid_ch, # label every electrode\n",
" )\n",
" ax.set_title(title, fontsize=12, fontweight='bold', pad=10)\n",
" plt.colorbar(im, ax=ax, fraction=0.046, pad=0.06,\n",
" label='MAV (µV)' if 'Diff' not in title else 'ΔMAV (µV)')\n",
"\n",
" # ── Shrink the auto-generated name texts so they don't overlap ──────────\n",
" for txt in ax.texts:\n",
" txt.set_fontsize(6.5)\n",
" txt.set_color('#222222')\n",
"\n",
" # ── Overlay larger, coloured markers for motor cortex channels ──────────\n",
" pos_dict = {ch: info['chs'][i]['loc'][:2]\n",
" for i, ch in enumerate(valid_ch)}\n",
"\n",
" for ch in MOTOR_CH:\n",
" if ch not in pos_dict:\n",
" continue\n",
" x_ch, y_ch = pos_dict[ch]\n",
" # bright ring so it stands out against the heatmap\n",
" ax.plot(x_ch, y_ch, 'o',\n",
" markersize=11, markerfacecolor='none',\n",
" markeredgecolor='lime', markeredgewidth=2.0,\n",
" zorder=10)\n",
" ax.text(x_ch, y_ch + 0.055, ch,\n",
" ha='center', va='bottom', fontsize=8,\n",
" fontweight='bold', color='lime',\n",
" bbox=dict(boxstyle='round,pad=0.15', fc='black', alpha=0.55),\n",
" zorder=11)\n",
"\n",
"# Legend for the highlight\n",
"from matplotlib.lines import Line2D\n",
"legend_el = [Line2D([0], [0], marker='o', color='w', markerfacecolor='none',\n",
" markeredgecolor='lime', markeredgewidth=2, markersize=9,\n",
" label='Motor cortex (C3 / Cz / C4)')]\n",
"fig.legend(handles=legend_el, loc='lower center', ncol=1,\n",
" fontsize=10, framealpha=0.8, bbox_to_anchor=(0.5, -0.04))\n",
"\n",
"plt.tight_layout()\n",
"plt.savefig('topoplot_MAV.png', dpi=150, bbox_inches='tight')\n",
"plt.show()\n",
"print('Saved: topoplot_MAV.png')"
]
},
{
"cell_type": "markdown",
"id": "248740bd",
"metadata": {},
"source": [
"---\n",
"## Figure 2 — Frequency Composition During MI vs REST"
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "393042a0",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"PSD channels used: Motor cortex: C3, Cz, C4\n"
]
}
],
"source": [
"# ── Compute PSD (Welch) for each epoch, then average ──────────────────────────\n",
"\n",
"def epochs_to_psd(epochs, valid_idx, n_pre, sfreq):\n",
" \"\"\"\n",
" Compute PSD for the active window of each epoch.\n",
" Returns: freqs (F,), psd_mean (n_channels, F), psd_sem (n_channels, F)\n",
" \"\"\"\n",
" psds = []\n",
" for ep in epochs:\n",
" active = ep[valid_idx, n_pre:] # (n_ch, n_active_samples)\n",
" f, p = welch(active, fs=sfreq,\n",
" nperseg=min(int(sfreq * 2), active.shape[1]),\n",
" noverlap=int(sfreq),\n",
" window='hann', axis=-1)\n",
" psds.append(p) # (n_ch, F)\n",
" psds = np.array(psds) # (n_trials, n_ch, F)\n",
" return f, psds.mean(axis=0), psds.std(axis=0) / np.sqrt(len(psds))\n",
"\n",
"\n",
"freqs, psd_mi_ch, sem_mi_ch = epochs_to_psd(all_mi_epochs, valid_idx, n_pre, sfreq_global)\n",
"freqs, psd_rest_ch, sem_rest_ch = epochs_to_psd(all_rest_epochs, valid_idx, n_pre, sfreq_global)\n",
"\n",
"# Average across all valid channels (or choose motor channels below)\n",
"motor_ch_names = ['C3', 'Cz', 'C4'] # primary motor cortex channels\n",
"motor_idx = [valid_ch.index(ch) for ch in motor_ch_names if ch in valid_ch]\n",
"if len(motor_idx) == 0:\n",
" motor_idx = list(range(len(valid_ch))) # fallback: all channels\n",
" motor_label = 'All channels'\n",
"else:\n",
" motor_label = 'Motor cortex: ' + ', '.join([valid_ch[i] for i in motor_idx])\n",
"\n",
"psd_mi = psd_mi_ch[motor_idx, :].mean(axis=0)\n",
"psd_rest = psd_rest_ch[motor_idx, :].mean(axis=0)\n",
"sem_mi = sem_mi_ch[motor_idx, :].mean(axis=0)\n",
"sem_rest = sem_rest_ch[motor_idx, :].mean(axis=0)\n",
"\n",
"# Restrict to 0.5–50 Hz for display\n",
"freq_mask = (freqs >= 0.5) & (freqs <= 50)\n",
"f_plot = freqs[freq_mask]\n",
"psd_mi_p = psd_mi[freq_mask]\n",
"psd_rest_p = psd_rest[freq_mask]\n",
"sem_mi_p = sem_mi[freq_mask]\n",
"sem_rest_p = sem_rest[freq_mask]\n",
"\n",
"# Convert to dB\n",
"psd_mi_db = 10 * np.log10(psd_mi_p + 1e-30)\n",
"psd_rest_db = 10 * np.log10(psd_rest_p + 1e-30)\n",
"sem_mi_db = 10 * np.log10(np.maximum(psd_mi_p - sem_mi_p, 1e-30) + 1e-30)\n",
"sem_rest_db = 10 * np.log10(np.maximum(psd_rest_p - sem_rest_p, 1e-30) + 1e-30)\n",
"err_mi = psd_mi_db - sem_mi_db\n",
"err_rest = psd_rest_db - sem_rest_db\n",
"\n",
"print(f'PSD channels used: {motor_label}')"
]
},
{
"cell_type": "code",
"execution_count": 22,
"id": "62a254bb",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABXkAAAL1CAYAAACSQOeBAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAASdAAAEnQB3mYfeAABAABJREFUeJzs3QWYG+XaBuAnWXeru7sXh1IKBUpxKO6cgx/k4O6u58fhoMWt+MELRUqRtkDdqHe77i75r+dLJ51ks7vZ3aw/d69c3SSTyWR83nm/93O4XC4XRERERERERERERKRdcrb2BIiIiIiIiIiIiIhI4ynIKyIiIiIiIiIiItKOKcgrIiIiIiIiIiIi0o4pyCsiIiIiIiIiIiLSjinIKyIiIiIiIiIiItKOKcgrIiIiIiIiIiIi0o4pyCsiIiIiIiIiIiLSjinIKyIiIiIiIiIiItKOKcgrIiIiIiIiIiIi0o4pyCsiIiIiIiIiIiLSjinIKyIiIiIiIiIiItKOKcgrIiIiIh7fffcdHA6H58HnDTVgwADP588880zN3Ua69dZbvZaFiIiIiEhtFOQVEZF2bcOGDV5BkLoe0rFkZ2fjwQcfxMEHH4xevXohMjISUVFRGDRoEE4++WTMmTMHJSUlrT2ZHXp7e+mll9Ce2YPRfISEhJh1qFu3bhg3bhxOOOEEvPrqqygtLW3tSW3zfPe33B7T0tL8DnvkkUfWGL4xNxPa4/ppn7799tuvtSdH2tBNHOsRGhqK5ORk7LrrrrjxxhuRkZHRatPIbcg+bdzGRESk7Qpt7QkQERERaaiXX34ZF198MQoKCmq8t379evN444038OKLLyqTtIEGDx6MBx54wOt5Q91www3Iy8szf48ZMwbtRXV1tQno8sHAypIlS/D222/jmmuuwSuvvIL999+/RafnoIMOQmxsLNqjsrIyPPHEE7j99tu9Xl+9ejU+/vjjVpsukbauqqoKOTk5+P33382Dx7HffvvN3MwUERGpi4K8IiLSoeyyyy4m+66xGNxhNl9YWFhQp0uC57HHHsMll1zi9dq0adOw9957myzMLVu2YO7cuVi1apVmeyP07dsXV155ZZPm3TnnnNPu5n1SUhKuv/56VFRUYNu2bfjqq6886xCfM+D6wQcf4LDDDmv2acnPz0d8fDz22msv82ivnn76aTNPmdVreeSRR+ByudDRWMtMgo838+Li4jr8rD3//PPNTbWsrCy8+eabnqxZ7n8eeugh8xAREamTS0REpB1bv349owWexxlnnFHvZ+zD33LLLa558+a5DjjgAFdCQoJ5jeO0bNq0yXXllVe6xo4d64qNjXWFh4e7Bg4c6PrnP//pWrlypd/x5+TkuC699FJXnz59XBEREa5hw4a57rnnHldFRUWN77a8+OKLXu/Zp4H69+9f529s6HTyu+3fV1pa6rr77rtdw4cPN5/t3r2767zzznPl5eX5/Y3r1q1zXX755a7x48e74uPjzWf4e2fMmOF68803zTDXXHONZ/zdunVzlZeXe42jsLDQFR0d7RnmrrvuqnfZrVq1yhUaGur5TFRUlOuzzz7zO+yXX35plq1dWVmZ66mnnnLtt99+rpSUFDOu5ORk17777ut67LHHzHzw5bvMON69997bfHePHj1cl1xyifkt9N5777l22WUXV2RkpJmH5557ris3N9drfN9++63XOPn8tddec+26665mnElJSa5jjz221vUrLS3NdeONN7omTZpk5n1YWJirZ8+eriOOOML10Ucf+f3MF198Yd7v1auXGZ7f07dvXzMfrrrqKjNf65o+33Wwtkeg6yvX78suu8w1evRoV0xMjFl/+vXr5zrhhBNc33//fY3hfbcPrn9cjlz/OK+5LE866STX1q1bXQ1hn07+bVddXe169NFHXQ6HwzNMYmKiKysryzMMf1ttn/fdN/E31PZ71qxZ47rvvvtcI0aMMPNi6tSpfrdTOw5jvc6/t2/f7jr//PPNMuY4hgwZ4rr//vvN7/DF7Zrbb0P2UfWxfy4kJMTz97PPPusZJjMz07PN24exr2d2S5YscZ1zzjmuoUOHmnWWy3rw4MGus88+2/XHH3/UuizrWz8pPz/fzPM999zTLFfuC7p27eo68MADXbNnz3ZVVVXVuzznzJljPs99bqCXVfZxWMu5tmXK9YL7Ak4fj09HHXWUa+3atWZY/v6ZM2e64uLizIN/L1++vMb3Pf/882a7GjVqlPl93P65zXFd4/7J32eauo5wnFwXeTzh8uZy4+e5zW/ZsqXG8L7bUUZGhuuCCy5w9e7d26wnN9xwg5kOaxiOx5d9m3I6na7NmzfXugy4bO3rC+eDLy5b++/866+/zOvZ2dmu66+/3ux7ON85fdz/cP6ecsoprv/+97+uQPlu3/ZtYMWKFV7v8djqz4IFC1ynnXaaOd5zPnN+8zzgpptu8tpXWQKdft/13d/Dd/0VEZHWpyCviIh06iAvL9B9gw1WgPXTTz81F0G1XeDwgopBPd8L4zFjxvgdnkG25gjyNmY6fS8up0yZ4vez06ZNqzH/GJBkwKW27zvyyCM9gWd7QPatt97yGs8bb7zhFRQKJEB34YUXen3XAw884AoUAwcMjNZ10TphwgRXenq61+fs7/Pz9qCf9dh///1dDz/8cEDz0DeIyhsM/j7HwM6ff/5Z44KegZq6fgMv+O0BqldffbXei3V7ALK5g7wMRDPIVNd4GISw890+altfGbjyF6hvTJDXwmCT/TsYGAx2kNf39zQ0yDto0CAT3PU3T2699Vavz/GGBNfzhu6j6mP/HAOO3O/wbwaPrEDzHXfc4RmGwUt/65nlmWeeMQHJ2tYR7luefPJJv8uyvvWTwVPOs7qGnT59uqu4uLjW5bnPPvvUOv5gBHkZtOMNKN/v4A2zDz/80O8+mPsG7ufsJk+eXOfvZADXd943ZR157rnnzE2G2r6PN7F+/PFHr8/Yt6MuXbqY7dj3O3gT0nrO+VJSUuI1joMPPtjz/iGHHFLvcrjtttu8pok3AO3s6+duu+1mXuO+pbbje337kYYGeXkTwv7eqaee6vc3+DseWQ/eOLPfLGzI9CvIKyLSPqlcg4iIdCjLli0znXH5Yl3QGTNm1Hj9559/RnR0tOmoq1+/fqYGJ0s1bNy4EccddxyKi4vNcAMHDsTxxx9vmh1/+OGH+OOPP0xph1NOOcV8Jzv7InaSsnTpUs/4x48fjyOOOALr1q0zNWKDrbHT6euHH37A0UcfjVGjRuG1117zNBP99ttv8csvv2D33Xc3z1kf8IwzzkBlZaV5zo5YDj/8cEycONF0hPb99997NfvnON955x3z/JlnnjHTZmFzVMvMmTMDqjf4zTffeP7md5999tkBz6vTTjsNixYt8jxnh2177LGHqXX4v//9z7zG+cV59eWXX/odBz8/evRoHHPMMfj888/NZ4nlIfgYO3YsjjrqKHzyySdYvHix33no7zftu+++phMmjp+fpdzcXFNP2JpmNgfnumR1wsPOefib+vTpg48++gh//vmneZ21Y0eMGGGayVvlLSzDhw8360t4eLgpa7FixQqzDQRaZ5frxd133+15jaVRWCIlUKyVzM9YHeJx2+NvTEhIwFtvvYW///7bvM7v4Hzmdlnb+nrAAQeYUgYsocDtllauXGmeN6Vki69zzz0XTz31lNfyuvrqqxFM/D0jR440y9fpdDa4w0DON27zF1xwgSlZwum1xvHwww+bdcEqQXPzzTeb9dzCdZYdoa1Zs8Ysg2Do2rWrWTf/+9//Yvny5fjiiy9MPWPW6LVKY5x11ll47733/H5+/vz55rewRjJ16dLFjI+ldGbPno309HSzD7rooovM9O+zzz4Br5+sd8pt1FrXiNsE931ctj/++KN57euvv8all16KZ5991u80criUlBTzHeyoz9oXBAu3FY6f6xqn9d133zWv87dzeXEesy45axxznSfuG55//nlTQ9rC4VhiZMiQIWa+cz3Yvn073n//fWzevNnUTv7Xv/7lddxq7DrC/Ry3F2u5WZ9jbJv7ex4HWWuWxwWOi9u
"text/plain": [
"<Figure size 1680x840 with 3 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Saved: frequency_spectrum.png\n"
]
}
],
"source": [
"# ── Plot ──────────────────────────────────────────────────────────────────────\n",
"\n",
"fig = plt.figure(figsize=(14, 7))\n",
"gs = GridSpec(2, 2, figure=fig, height_ratios=[3, 1], hspace=0.08, wspace=0.35)\n",
"\n",
"ax_main = fig.add_subplot(gs[0, :])\n",
"ax_mi = fig.add_subplot(gs[1, 0])\n",
"ax_rest = fig.add_subplot(gs[1, 1])\n",
"\n",
"# ── Coloured frequency-band background ────────────────────────────────────────\n",
"y_min = min(psd_mi_db.min(), psd_rest_db.min()) - 2\n",
"y_max = max(psd_mi_db.max(), psd_rest_db.max()) + 4\n",
"\n",
"for label, (lo, hi, color) in BANDS.items():\n",
" lo_c = max(lo, f_plot[0])\n",
" hi_c = min(hi, f_plot[-1])\n",
" ax_main.axvspan(lo_c, hi_c, alpha=0.18, color=color, zorder=0)\n",
" mid = (lo_c + hi_c) / 2\n",
" ax_main.text(mid, y_max - 0.5, label, ha='center', va='top',\n",
" fontsize=8, color=color, fontweight='bold')\n",
"\n",
"# ── PSD curves ────────────────────────────────────────────────────────────────\n",
"ax_main.plot(f_plot, psd_mi_db, color='#E05C2A', lw=2.2, label='Motor Imagery', zorder=3)\n",
"ax_main.fill_between(f_plot,\n",
" psd_mi_db - err_mi,\n",
" psd_mi_db + err_mi,\n",
" color='#E05C2A', alpha=0.20, zorder=2)\n",
"\n",
"ax_main.plot(f_plot, psd_rest_db, color='#2A7BE0', lw=2.2, label='Rest', zorder=3)\n",
"ax_main.fill_between(f_plot,\n",
" psd_rest_db - err_rest,\n",
" psd_rest_db + err_rest,\n",
" color='#2A7BE0', alpha=0.20, zorder=2)\n",
"\n",
"ax_main.set_xlim(f_plot[0], f_plot[-1])\n",
"ax_main.set_ylim(y_min, y_max)\n",
"ax_main.set_ylabel('Power Spectral Density (dB)', fontsize=11)\n",
"ax_main.set_xlabel('')\n",
"ax_main.set_xticklabels([])\n",
"ax_main.legend(fontsize=11, loc='upper right')\n",
"ax_main.set_title(\n",
" f'Frequency Composition During Motor Imagery vs Rest\\n'\n",
" f'{motor_label} | Welch PSD | Analysis window [0, {T_POST}] s | ±1 SEM shaded',\n",
" fontsize=12, fontweight='bold')\n",
"ax_main.grid(True, alpha=0.3, lw=0.7)\n",
"\n",
"# ── Band power bar charts ─────────────────────────────────────────────────────\n",
"\n",
"def band_power(f, psd, lo, hi):\n",
" mask = (f >= lo) & (f < hi)\n",
" if mask.sum() == 0:\n",
" return 0.0\n",
" return np.trapezoid(psd[mask], f[mask])\n",
"\n",
"\n",
"band_labels = []\n",
"bp_mi_vals = []\n",
"bp_rest_vals = []\n",
"bar_colors = []\n",
"\n",
"for label, (lo, hi, color) in BANDS.items():\n",
" band_labels.append(label)\n",
" bp_mi_vals.append(band_power(f_plot, psd_mi_p, lo, hi))\n",
" bp_rest_vals.append(band_power(f_plot, psd_rest_p, lo, hi))\n",
" bar_colors.append(color)\n",
"\n",
"# Normalise so bars show relative contribution (% of total power)\n",
"total_mi = sum(bp_mi_vals) + 1e-30\n",
"total_rest = sum(bp_rest_vals) + 1e-30\n",
"bp_mi_pct = [v / total_mi * 100 for v in bp_mi_vals]\n",
"bp_rest_pct = [v / total_rest * 100 for v in bp_rest_vals]\n",
"\n",
"x = np.arange(len(band_labels))\n",
"w = 0.38\n",
"\n",
"for ax, bp, title, line_color in [\n",
" (ax_mi, bp_mi_pct, 'MI — Band Power (%)', '#E05C2A'),\n",
" (ax_rest, bp_rest_pct, 'Rest — Band Power (%)', '#2A7BE0'),\n",
"]:\n",
" bars = ax.bar(x, bp, width=0.65, color=bar_colors, edgecolor='white',\n",
" linewidth=0.8, zorder=2)\n",
" for bar, val in zip(bars, bp):\n",
" ax.text(bar.get_x() + bar.get_width() / 2, bar.get_height() + 0.3,\n",
" f'{val:.1f}%', ha='center', va='bottom', fontsize=7.5)\n",
" ax.set_xticks(x)\n",
" ax.set_xticklabels(band_labels, fontsize=8)\n",
" ax.set_ylabel('% of total power', fontsize=9)\n",
" ax.set_title(title, fontsize=10, fontweight='bold', color=line_color)\n",
" ax.set_ylim(0, max(max(bp_mi_pct), max(bp_rest_pct)) * 1.25)\n",
" ax.grid(axis='y', alpha=0.3, lw=0.7)\n",
" ax.spines[['top', 'right']].set_visible(False)\n",
"\n",
"plt.savefig('frequency_spectrum.png', dpi=150, bbox_inches='tight')\n",
"plt.show()\n",
"print('Saved: frequency_spectrum.png')"
]
},
{
"cell_type": "markdown",
"id": "fcb6d19d",
"metadata": {},
"source": [
"---\n",
"## Figure 3 — Per-Channel Topoplots for Each Frequency Band"
]
},
{
"cell_type": "code",
"execution_count": 23,
"id": "75df404b",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAClkAAAPNCAYAAAAAu57jAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAASdAAAEnQB3mYfeAABAABJREFUeJzs3Qd41PQfBvC3u5S9y94b2VumDEEFBBmCKAgiCqjgQARUUFFE3MpfETeKyEaGyJYte+/RMtoyuuhe93++6eWaS3OrvdIC7+d58txdksvlcndJLnnz/XmYTCYTiIiIiIiIiIiIiIiIiIiIiIiIiIjIiqf1QyIiIiIiIiIiIiIiIiIiIiIiIiIiYsiSiIiIiIiIiIiIiIiIiIiIiIiIiMgGVrIkIiIiIiIiIiIiIiIiIiIiIiIiIjLAkCURERERERERERERERERERERERERkQGGLImIiIiIiIiIiIiIiIiIiIiIiIiIDDBkSURERERERERERERERERERERERERkgCFLIiIiIiIiIiIiIiIiIiIiIiIiIiIDDFkSERERERERERERERERERERERERERlgyJKIiIiIiIiIiIiIiIiIiIiIiIiIyABDlkREREREREREREREREREREREREREBhiyJCIiIiIiIiIiIiIiIiIiIiIiIiIywJAlERERERHRPcbDw8Ow8/HxQdGiRdG4cWOMGTMGBw4cyO1ZdcnUqVOt3s/mzZudfu5PP/1kc7l4e3ujePHiaNeuHT788ENERkbm6Pu4k3Xs2NHmcrTXyfPuJcOGDbN6/xcvXsztWborlqN0999/v83xJ0+enGn8ypUrW40jn4V+nDuVvXVasWLF0Lx5c0yYMAHBwcGGz5d1aHZ/w0uWLEHv3r1Rvnx5+Pn5IX/+/KhQoQKaNGmCIUOG4KOPPsL58+ct48vnkZV1iHwX3OmDDz6wTPuhhx7CnSg720RyDtfleZd+/SW/B7q9cmI/z2QyoX79+pbp/vHHH26ZLhEREREREZEzGLIkIiIiIiIiRUpKihIgPHjwIGbPno2mTZtizpw59/zSSU1NRXh4OLZt24aJEycqJ3cPHz58zy8Xynv0AUF3B8/uBDt27MChQ4cy9U9KSsLcuXNzZZ7y4jotIiICe/fuVUKOderUwb///uv27Un//v3x2GOPYcWKFbhy5YryGcTFxeHy5ctKiP+3335TQp7//PMP8pKwsDAlZKmaMmVKrs4PERGlk30buWBCJfvlCQkJXDxERERERER0W3jfnpchIiIiIiKivKpHjx4ICAhQTlLu3LlTCRSq1WJeeuklJSQjlRzvJSVKlECHDh2U+6Ghodi1a5cSTBISFpLKbCdOnIC/v38uz2neIstMlp3W8ePHlWWlqlSpEpo1a2Y1Tr169W7bPNLdT0Li3377rVW/hQsX4tq1a7iXqet6CRHKul5dp0nwccSIEThz5ozd58tzZRpG9L/hzz//HIsWLbI8zpcvn1I5U6olR0VFKesEmQ89qRqp/5wkDBoUFGR5LKHQunXrWo0j03aXd955B7du3VLut23bFm3atHHbtImIKHsGDhyISZMmKReWyLbhyy+/xGuvvcbFSkRERERERDmOIUsiIiIiIqJ7nASS1CZzJfRSs2ZNREdHK48leCmV4Xr27Il7iQSGtAEhaXayc+fOSEtLUx7Lid1ly5bh8ccfz8W5zHumTZuWqZ800antL01GSvPsRDlFKiRKhcZChQpZ+n399df3/ALXruvXrl2L7t27W5bJ2bNnlZBljRo1bC6nkiVLWq0X7dH+xuU19+zZkymALUHLP//8E6VLl7aaRz2pyPrzzz9bHg8YMCDHmv6VAKh23ocOHZojr0NERFnj6emJJ554AtOnT7ds31955RWlPxEREREREVFO4j9PIiIiIiIispCwS61atayWiL4ZPgkcSoXL9u3bo0qVKihcuDB8fHyUCmUtWrRQqstcvXrVcKlqmzKWsF18fLxyklSqkklVSAnhDBo0SAkxGklMTFTGl3mU8cuVK4dRo0YZVkRzJ5nXdu3aWfWT0JDeX3/9pVT+rFixojJ/BQoUUOZ15MiRSvO4ehJoUpdHkyZNrIZ9/PHHlmFFihSxBDzF4sWLrZbl77//nmk5SdPIDz74oPKZ+vr6Kp/P/fffj88++0xZ7nryuWqnKSEmCV0NHjwYgYGB8PLyyrFgkwgJCVGa5VWr3ck8lypVSgm3ygl0o3mWMJR2nuXxsWPHlCpH8lz5DKR5908//VRpvtiITFeCXV26dFGeI68ry1uqbUqTlDJfzjTLLRVgx48frwTK/Pz8lO+A/E7UyrA5vTxkPuT3qCXBNP1nqpImtYcPH66EqvPnz6/8huW70rBhQyVYJq8hzTvfScqWLavcxsbGWoXy5LcnlRu147ibfB7yvVGXdcuWLQ3Hkyq46jje3t6WdaVUDpZ5lt9smTJllM9bPhep/CrrWgmQbNy40W3zK68j626tmzdvum362qqYDRo0yBSwVCtSvv322+jTpw/yCvkMpLKnkM+gX79+NseVCsfy25f1uHxWUq1Tfk+jR4+2WRX0xx9/VH53sn4pX7688jxZX8g6Vn7bUpHNXtO3sl6Q37FU15QK0/K7lfCrTE+aXne0vjl48CD69u2rfB6yfrzvvvswZ84cuErmQbtuke2HNPvetWtXZX2l9lPJ+nfevHl4+OGHLd9vCUHL648bN04J+WolJycry0advsyz1gsvvGAZJtOwte2UTi4UcYa7foObNm1SpiHrA6n+KvtFS5cuzTSefM+kWXoJDct2Sr4D8l2Q58h2RN6z0fOysvxl/0ECzb169VL2m+R1ZPnL9+bdd99FRESE4evkxHZClvPKlSuV7XTVqlWV6arvWfafli9fnuk5kZGR+PDDD5XKsur3vlixYso+zYwZM5ThWeHq99LW9t/Rvm52Pjt7/v33X+ViH1mOsv6Rz1W2cU2bNsWzzz5r84KWc+fO4dVXX1XGU/cv5Psny/Ott95Sfn/u2Od3RlbWo0L2TVVSzVK+U0REREREREQ5zkRERERERET3FPkrqO0uXLhgGRYaGmoqVKiQZZiHh4fpzJkzVs8fMWJEpmnouyJFipj27t1r97Xr1atnatiwoeHzy5UrZ7p586bVcxMSEkwdOnQwHD8wMNA0ePBgq36bNm1yepn8+OOPVs+V19Hr37+/1TjPPvusZVhcXJypZ8+edpeJp6enadq0aVbTlGloh0dFRVmG9erVy+r5+/fvtwx74YUXrIaFhIRYhl28eNHUoEEDu/Miyz4oKMhqXmR5aceR91OwYEGrfm+//bbJVfIc7TSGDh2aaZyVK1eaChcubHeea9eunem7qP/cBg0aZPLz8zN8/qOPPmpKTU21er5MT6Zr73VlvmT+tOQ3ox2nU6dOpkqVKhk+v0aNGsrvSkuWga3fYFaXh6PfpPbzk8/a19fX4fjXr1835WX65aj9rtWpU8dwnTV16lSr58jnZu+zlc5Zzz33nNXzzp49azVc1mk+Pj5WvzHVyJEjHX4eDz/8sEvLR/98/fdMu66XztE6Qb+s7NGuO7y8vEyTJ09W1mEpKSkmd3zOOUW7jWndurXhOGlpaaZx48bZ/az8/f1Nv/32W6bnVqtWzeHnLOvviIiITM9duHBhpnWyvjtw4IDNde+oUaNM3t7ehs+bMWOGS8tJP2399le7Db527ZqyLO3Nt6y3586da/Ua3bp1swwvWbKk1TDtNk72U7T7C9ptsSyv5ORkp95TVn6D+u/mmDFjbD53/vz5Vs+9dOmSU+vtp556SvnOZXX5h4eHmzp27Gj3NcqXL286ePCg1WvkxHZC9nEefPBBu9Pr3bu31XN2795tKlu2rN3nyHAZTz//9tYbWf1e6rcRRvs09vYnXfns7Jk3b57y3bc3//nz58/0vC+++MLh56pd/7hrn1+/HLKzHlXJvr867rBhwxwuMyIiIiIiIqLsYiVLIiIiIiKie5xUi5FqXY888ojSTLbaVLiQ6jXVq1fP9BypwCaVlzp06KBUZuvWrZtSlUslVYWk+pE9UnFQqiRVq1Z
"text/plain": [
"<Figure size 2700x960 with 12 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Saved: topoplot_band_power.png\n"
]
}
],
"source": [
"# Band power per channel (used for band-specific topoplots)\n",
"\n",
"def band_power_per_channel(epochs, valid_idx, n_pre, sfreq, lo, hi):\n",
" \"\"\"\n",
" Returns mean band power (µV²/Hz integrated) per channel across epochs.\n",
" Shape: (n_valid_channels,)\n",
" \"\"\"\n",
" channel_bp = []\n",
" for ep in epochs:\n",
" active = ep[valid_idx, n_pre:]\n",
" f, p = welch(active, fs=sfreq,\n",
" nperseg=min(int(sfreq * 2), active.shape[1]),\n",
" noverlap=int(sfreq), window='hann', axis=-1)\n",
" mask = (f >= lo) & (f < hi)\n",
" bp = np.trapezoid(p[:, mask], f[mask], axis=-1) # (n_ch,)\n",
" channel_bp.append(bp)\n",
" return np.mean(channel_bp, axis=0)\n",
"\n",
"\n",
"n_bands = len(BANDS)\n",
"fig, axes = plt.subplots(2, n_bands, figsize=(4.5 * n_bands, 8))\n",
"fig.suptitle('Band Power Topoplots — MI vs REST (each row shares colour scale)',\n",
" fontsize=13, fontweight='bold', y=1.01)\n",
"\n",
"row_labels = ['Motor Imagery', 'REST']\n",
"epoch_sets = [all_mi_epochs, all_rest_epochs]\n",
"\n",
"for col, (band_label, (lo, hi, bcolor)) in enumerate(BANDS.items()):\n",
" bp_mi_ch = band_power_per_channel(all_mi_epochs, valid_idx, n_pre, sfreq_global, lo, hi) * scale**2\n",
" bp_rest_ch = band_power_per_channel(all_rest_epochs, valid_idx, n_pre, sfreq_global, lo, hi) * scale**2\n",
"\n",
" vmin = min(bp_mi_ch.min(), bp_rest_ch.min())\n",
" vmax = max(bp_mi_ch.max(), bp_rest_ch.max())\n",
"\n",
" for row, (bp_ch, rlabel) in enumerate([(bp_mi_ch, 'MI'), (bp_rest_ch, 'REST')]):\n",
" ax = axes[row, col]\n",
" im, _ = mne.viz.plot_topomap(\n",
" bp_ch, info, axes=ax,\n",
" cmap='hot_r',\n",
" vlim=(vmin, vmax),\n",
" show=False, contours=4, sensors=True,\n",
" )\n",
" if row == 0:\n",
" ax.set_title(band_label.replace('\\n', ' '), fontsize=9.5,\n",
" fontweight='bold', color=bcolor)\n",
" if col == 0:\n",
" ax.set_ylabel(rlabel, fontsize=10)\n",
" if col == n_bands - 1:\n",
" plt.colorbar(im, ax=ax, fraction=0.046, pad=0.08, label='µV²/Hz')\n",
"\n",
"plt.tight_layout()\n",
"plt.savefig('topoplot_band_power.png', dpi=150, bbox_inches='tight')\n",
"plt.show()\n",
"print('Saved: topoplot_band_power.png')"
]
},
{
"cell_type": "markdown",
"id": "b3db60ba",
"metadata": {},
"source": [
"---\n",
"## Summary Statistics"
]
},
{
"cell_type": "code",
"execution_count": 24,
"id": "cf55268e",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"=== Band Power Summary (motor channels: C3, Cz, C4) ===\n",
"Band MI (µV²/Hz) REST (µV²/Hz) MI/REST ratio\n",
"------------------------------------------------------------------\n",
"Delta (0.5–4 Hz) 31.3986 25.8102 1.217 ← ERD/ERS\n",
"Theta (4–8 Hz) 6.6702 4.4569 1.497 ← ERD/ERS\n",
"Mu/Alpha (8–13 Hz) 7.5159 9.0951 0.826 ← ERD/ERS\n",
"Beta (13–30 Hz) 3.7107 2.7744 1.337 ← ERD/ERS\n",
"Gamma (30–50 Hz) 1.2799 0.8454 1.514 ← ERD/ERS\n",
"\n",
"ERD = Event-Related Desynchronization (ratio < 1 → suppressed during MI)\n",
"ERS = Event-Related Synchronization (ratio > 1 → enhanced during MI)\n"
]
}
],
"source": [
"print('=== Band Power Summary (motor channels: C3, Cz, C4) ===')\n",
"print(f'{\"Band\":<22} {\"MI (µV²/Hz)\":>14} {\"REST (µV²/Hz)\":>14} {\"MI/REST ratio\":>14}')\n",
"print('-' * 66)\n",
"for label, (lo, hi, _) in BANDS.items():\n",
" bp_mi_m = band_power_per_channel(all_mi_epochs, valid_idx, n_pre, sfreq_global, lo, hi)\n",
" bp_rest_m = band_power_per_channel(all_rest_epochs, valid_idx, n_pre, sfreq_global, lo, hi)\n",
" m_mi = bp_mi_m[motor_idx].mean() * scale**2 if motor_idx else bp_mi_m.mean() * scale**2\n",
" m_rest = bp_rest_m[motor_idx].mean() * scale**2 if motor_idx else bp_rest_m.mean() * scale**2\n",
" ratio = m_mi / (m_rest + 1e-30)\n",
" flag = ' ← ERD/ERS' if abs(ratio - 1) > 0.1 else ''\n",
" clean_label = label.replace('\\n', ' ')\n",
" print(f'{clean_label:<22} {m_mi:>14.4f} {m_rest:>14.4f} {ratio:>14.3f}{flag}')\n",
"\n",
"print('\\nERD = Event-Related Desynchronization (ratio < 1 → suppressed during MI)')\n",
"print('ERS = Event-Related Synchronization (ratio > 1 → enhanced during MI)')"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "venv",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.13.12"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Reference in New Issue Copy Permalink