{
  "nbformat": 4,
  "nbformat_minor": 5,
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "name": "python",
      "version": "3.10.0"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "id": "md68364437",
      "metadata": {},
      "source": "# \ud83e\udde0 Tutorial 3 \u2014 EEG Analysis & Classification\n###  \u00b7 Biomedical Signals & Applications \u00b7 EEG Python Workshop\n\n---\n\n**Goal of this tutorial:** Analyse the clean EEG epochs and build a classifier that can distinguish left-hand from right-hand motor imagery.\n\n**By the end you will be able to:**\n- Compute and interpret Event-Related Potentials (ERPs)\n- Analyse the power spectrum and extract band power features\n- Build time-frequency maps and identify ERD/ERS patterns\n- Construct a feature matrix for machine learning\n- Train and evaluate an LDA classifier with cross-validation\n- Apply CSP spatial filtering for improved classification\n\n**Time:** approximately 60 minutes\n\n---\n\n> \u26a0\ufe0f **Run Tutorial 2 first** \u2014 it creates the clean epoch file this notebook loads."
    },
    {
      "cell_type": "markdown",
      "id": "md72880868",
      "metadata": {},
      "source": "---\n## Setup \u2014 Imports\n\nRun this cell first every time you open the notebook."
    },
    {
      "cell_type": "code",
      "id": "cd34736613",
      "metadata": {},
      "source": "import numpy as np\nimport matplotlib.pyplot as plt\nimport matplotlib.gridspec as gridspec\nimport mne\nfrom mne.time_frequency import tfr_morlet\nfrom scipy.signal import welch\nfrom scipy.stats import ttest_ind\n\nmne.set_log_level(\"WARNING\")\nprint(\"All imports successful \u2713\")",
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "id": "md55454484",
      "metadata": {},
      "source": "---\n## Step 0 \u2014 Load Preprocessed Epochs\n\nWe load the clean epochs saved in Tutorial 2 and split them by condition (left hand vs right hand)."
    },
    {
      "cell_type": "code",
      "id": "cd69124891",
      "metadata": {},
      "source": "epochs = mne.read_epochs(\"/tmp/motor_imagery_epochs-epo.fif\", verbose=False)\n\n# Split by condition\nepochs_left  = epochs[\"left_hand\"]\nepochs_right = epochs[\"right_hand\"]\n\nprint(\"Epochs loaded \u2713\")\nprint(f\"  Left hand  : {len(epochs_left)} trials\")\nprint(f\"  Right hand : {len(epochs_right)} trials\")\nprint(f\"  Channels   : {len(epochs.ch_names)}\")\nprint(f\"  Time axis  : {epochs.times[0]:.1f} s  to  {epochs.times[-1]:.1f} s\")\nprint(f\"  Shape      : {epochs.get_data().shape}  (trials \u00d7 channels \u00d7 time points)\")",
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "id": "md43560491",
      "metadata": {},
      "source": "---\n## Step 1 \u2014 Event-Related Potentials (ERPs)\n\n### What is an ERP?\n\nAn **ERP** is the average of many EEG epochs time-locked to the same event. Random background noise averages out (positive and negative cancel each other). Only the brain's consistent, time-locked response survives.\n\n> \ud83d\udfe2 **Analogy:** Photographing a faint star in the night sky. One 1-second exposure shows mostly noise. Average 100 photos and the noise cancels \u2014 the star becomes clearly visible.\n\n### The contralateral rule\n\nThe motor cortex is organised **contralaterally** \u2014 the LEFT hemisphere controls the RIGHT hand:\n\n| Electrode | Brain region | Most active during |\n|---|---|---|\n| **C3** | Left motor cortex | **Right** hand imagery |\n| **Cz** | Central / bilateral | Legs, bilateral |\n| **C4** | Right motor cortex | **Left** hand imagery |\n\nThe C3 vs C4 comparison should show this lateralisation directly."
    },
    {
      "cell_type": "code",
      "id": "cd22752489",
      "metadata": {},
      "source": "# Average all trials within each condition\nevoked_left  = epochs_left.average()\nevoked_right = epochs_right.average()\n\nprint(f\"Evoked left  : {evoked_left.nave} trials averaged\")\nprint(f\"Evoked right : {evoked_right.nave} trials averaged\")",
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "id": "cd64074813",
      "metadata": {},
      "source": "# Butterfly plot \u2014 all channels, both conditions\nfig, axes = plt.subplots(1, 2, figsize=(13, 4))\nevoked_left.plot( axes=axes[0], show=False, titles=\"Evoked \u2014 Left Hand Imagery\")\nevoked_right.plot(axes=axes[1], show=False, titles=\"Evoked \u2014 Right Hand Imagery\")\nplt.tight_layout()\nplt.show()",
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "id": "cd63657293",
      "metadata": {},
      "source": "# Topographic maps at several time points during imagery\ntimes_topo = np.arange(0.5, 3.0, 0.5)   # 0.5, 1.0, 1.5, 2.0, 2.5 s\nfig = evoked_left.plot_topomap(times=times_topo, ch_type=\"eeg\", show=True)\nplt.suptitle(\"Left Hand Imagery \u2014 Scalp Distribution Over Time\")\nplt.show()",
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "id": "cd85237226",
      "metadata": {},
      "source": "# Compare C3 and C4 for both conditions \u2014 this is the lateralisation plot\nidx_c3 = epochs.ch_names.index(\"C3\")\nidx_c4 = epochs.ch_names.index(\"C4\")\n\nfig, axes = plt.subplots(1, 2, figsize=(13, 4), sharey=True)\nfor ax, ch_idx, ch_name in zip(axes, [idx_c3, idx_c4], [\"C3\", \"C4\"]):\n    ax.plot(evoked_left.times,  evoked_left.data[ch_idx]  * 1e6,\n            \"steelblue\", lw=2, label=\"Left hand imagery\")\n    ax.plot(evoked_right.times, evoked_right.data[ch_idx] * 1e6,\n            \"tomato\",    lw=2, label=\"Right hand imagery\")\n    ax.axvline(0, color=\"k\", lw=1, linestyle=\"--\", label=\"Cue onset\")\n    ax.axhline(0, color=\"k\", lw=0.5)\n    ax.set_title(f\"ERP at {ch_name}\")\n    ax.set_xlabel(\"Time (s)\")\n    ax.set_ylabel(\"Amplitude (\u00b5V)\")\n    ax.legend(fontsize=9)\n    ax.grid(True, alpha=0.3)\nplt.suptitle(\"Lateralised Motor Response: C3 (left hemisphere) vs C4 (right hemisphere)\")\nplt.tight_layout()\nplt.show()",
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "id": "md12482242",
      "metadata": {},
      "source": "> \ud83d\udc41 **What to check:**\n> - C3 and C4 should show **opposite patterns** \u2014 when C3 is more negative for one condition, C4 should be more negative for the other\n> - This is the contralateral lateralisation: right-hand imagery activates C3 (left hemisphere), left-hand activates C4 (right hemisphere)\n> - The difference between conditions builds up gradually after the cue (t = 0)"
    },
    {
      "cell_type": "markdown",
      "id": "md67416480",
      "metadata": {},
      "source": "---\n## Step 2 \u2014 Power Spectral Density (PSD) & Band Power\n\n### EEG frequency bands\n\nDifferent mental states produce different dominant frequencies:\n\n| Band | Range | What it means for motor imagery |\n|---|---|---|\n| **Delta** | 0.5\u20134 Hz | Deep sleep (not relevant here) |\n| **Theta** | 4\u20138 Hz | Memory, frontal activity |\n| **Alpha / Mu** | 8\u201313 Hz | **Motor cortex idle rhythm \u2014 decreases during imagery (ERD)** |\n| **Beta** | 13\u201330 Hz | **Active motor processing \u2014 also decreases during imagery** |\n| **Gamma** | 30\u201345 Hz | Higher cognition; muscle noise overlaps here |\n\n> \ud83d\udfe2 **Analogy:** The PSD is like a music equaliser \u2014 it shows the \"volume\" at each frequency. We are looking for whether the \"mu/alpha volume knob\" is turned down during motor imagery.\n\n### Welch's method\n\n`compute_psd(method='welch')` splits the signal into overlapping windows, computes the FFT of each, and averages them. This gives a much smoother spectrum than a single FFT."
    },
    {
      "cell_type": "code",
      "id": "cd14950012",
      "metadata": {},
      "source": "# Plot Welch PSD for a single epoch, channel C3\nfs = epochs.info[\"sfreq\"]\nep_data = epochs_left.get_data()   # (n_epochs, n_channels, n_times)\n\nfreqs_w, psd_w = welch(ep_data[0, idx_c3], fs=fs, nperseg=int(fs * 2))\n\nfig, ax = plt.subplots(figsize=(10, 4))\nax.semilogy(freqs_w, psd_w, \"steelblue\", lw=1.5)\nax.set_xlim(0, 50)\nax.set_xlabel(\"Frequency (Hz)\")\nax.set_ylabel(\"PSD (V\u00b2/Hz)\")\nax.set_title(\"Welch PSD \u2014 Channel C3, Left Hand (single epoch)\")\n\n# Shade the EEG bands\nband_colors = {\"Delta\": \"#c0e8ff\", \"Theta\": \"#a0d8a0\",\n               \"Alpha/Mu\": \"#ffe0a0\", \"Beta\": \"#ffc0c0\", \"Gamma\": \"#e0c0ff\"}\nband_limits = {\"Delta\": (0.5,4), \"Theta\": (4,8), \"Alpha/Mu\": (8,13),\n               \"Beta\": (13,30), \"Gamma\": (30,45)}\nfor band, (lo, hi) in band_limits.items():\n    ax.axvspan(lo, hi, alpha=0.3, color=band_colors[band], label=band)\nax.legend(loc=\"upper right\", fontsize=9)\nax.grid(True, alpha=0.3)\nplt.tight_layout()\nplt.show()",
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "id": "cd26078263",
      "metadata": {},
      "source": "# Compute mean PSD over all trials using MNE \u2014 better estimate\npsd_left  = epochs_left.compute_psd( method=\"welch\", fmin=1, fmax=45)\npsd_right = epochs_right.compute_psd(method=\"welch\", fmin=1, fmax=45)\n\ndef band_power(psd_obj, fmin, fmax):\n    \"\"\"Average PSD power within a frequency band. Returns (n_epochs, n_channels).\"\"\"\n    data, freqs = psd_obj.get_data(return_freqs=True)\n    mask = (freqs >= fmin) & (freqs <= fmax)\n    return data[:, :, mask].mean(axis=-1)\n\n# Compare alpha/mu band power at C3 between conditions\nbp_left_alpha  = band_power(psd_left,  8, 13)   # (n_epochs, 64)\nbp_right_alpha = band_power(psd_right, 8, 13)\n\nprint(\"Alpha/Mu band power at C3:\")\nprint(f\"  Left hand  : {bp_left_alpha[:,  idx_c3].mean():.4e} V\u00b2/Hz\")\nprint(f\"  Right hand : {bp_right_alpha[:, idx_c3].mean():.4e} V\u00b2/Hz\")\n\ndiff_pct = (bp_left_alpha[:, idx_c3].mean() - bp_right_alpha[:, idx_c3].mean()) /             bp_left_alpha[:, idx_c3].mean() * 100\nprint(f\"  Difference : {diff_pct:.1f}%  (positive = more power for left hand)\")",
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "id": "md62269341",
      "metadata": {},
      "source": "> \ud83d\udc41 **What to check:**\n> - The alpha/mu power at C3 should be **lower for right-hand** imagery (ERD is stronger at C3 during contralateral imagery)\n> - If the difference is near zero, this subject may not show strong lateralisation \u2014 try C4"
    },
    {
      "cell_type": "markdown",
      "id": "md19159210",
      "metadata": {},
      "source": "---\n## Step 3 \u2014 Time-Frequency Analysis (TFR) & ERD/ERS\n\n### Why time-frequency?\n\nA power spectrum averages over the whole epoch \u2014 it cannot tell you *when* the frequencies changed. A **time-frequency representation (TFR)** is a 2D image showing how power at each frequency evolves over time.\n\n> \ud83d\udfe2 **Analogy:** A PSD is a recipe card listing all ingredients. A TFR is a video of the cooking process \u2014 you can see when each ingredient was added and how the mixture changed over time.\n\n### ERD and ERS\n\n| Term | Meaning | What it looks like in the TFR |\n|---|---|---|\n| **ERD** (Desynchronisation) | Power **decreases** below baseline | **Blue** region in the TFR |\n| **ERS** (Synchronisation) | Power **increases** above baseline | **Red** region in the TFR |\n\n**Expected pattern for right-hand imagery at C3:**\n- Blue patch at 8\u201313 Hz starting ~0.5 s after cue (mu ERD)\n- Blue patch at 18\u201326 Hz over the same period (beta ERD)\n- Possible red patch at 18\u201326 Hz after 3 s (beta rebound \u2014 motor cortex returning to idle)\n\n### Morlet wavelets\n\n`tfr_morlet()` slides short wavelets at each frequency along the signal and measures how well they match. The parameter `n_cycles` controls the trade-off between time and frequency precision."
    },
    {
      "cell_type": "code",
      "id": "cd35286517",
      "metadata": {},
      "source": "freqs_tfr = np.arange(6, 40, 1)    # analyse 6 to 40 Hz in 1 Hz steps\nn_cycles  = freqs_tfr / 2.0         # wavelet length: frequency/2 cycles (balanced trade-off)\n\nprint(f\"Computing TFR for {len(freqs_tfr)} frequencies \u00d7 {len(epochs.times)} time points...\")\n\ntfr_left  = tfr_morlet(epochs_left,  freqs=freqs_tfr, n_cycles=n_cycles,\n                        use_fft=True, return_itc=False, verbose=False)\ntfr_right = tfr_morlet(epochs_right, freqs=freqs_tfr, n_cycles=n_cycles,\n                        use_fft=True, return_itc=False, verbose=False)\n\n# Express power as percentage change from the pre-stimulus baseline (\u22120.5 to 0 s)\ntfr_left.apply_baseline( (-0.5, 0), mode=\"percent\", verbose=False)\ntfr_right.apply_baseline((-0.5, 0), mode=\"percent\", verbose=False)\n\nprint(\"TFR computed \u2713\")",
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "id": "cd29616153",
      "metadata": {},
      "source": "# Plot TFR at C3 for both conditions side by side\nfig, axes = plt.subplots(1, 2, figsize=(13, 5))\ntfr_left.plot( [idx_c3], axes=axes[0], show=False,\n               title=\"TFR at C3 \u2014 Left Hand Imagery (% change from baseline)\")\ntfr_right.plot([idx_c3], axes=axes[1], show=False,\n               title=\"TFR at C3 \u2014 Right Hand Imagery (% change from baseline)\")\nplt.tight_layout()\nplt.show()",
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "id": "md15321915",
      "metadata": {},
      "source": "> \ud83d\udc41 **What to check:**\n> - For **right-hand imagery** at C3: expect a stronger blue region (bigger ERD) in the mu (8\u201313 Hz) and beta (18\u201326 Hz) bands\n> - For **left-hand imagery** at C3: less ERD (C3 is *ipsilateral* to the left hand \u2014 the opposite hemisphere activates)\n> - This asymmetry is the core signal that makes motor imagery BCI work"
    },
    {
      "cell_type": "code",
      "id": "cd73723867",
      "metadata": {},
      "source": "# Difference map: right hand minus left hand at C3\n# Red = more power for right; Blue = more power for left\ntfr_diff = tfr_right.data[idx_c3] - tfr_left.data[idx_c3]\n\nfig, ax = plt.subplots(figsize=(9, 5))\nim = ax.imshow(\n    tfr_diff,\n    aspect=\"auto\",\n    origin=\"lower\",\n    extent=[epochs.times[0], epochs.times[-1], freqs_tfr[0], freqs_tfr[-1]],\n    cmap=\"RdBu_r\",\n    vmin=-0.5, vmax=0.5,\n)\nax.axvline(0, color=\"k\", lw=1.5, linestyle=\"--\", label=\"Cue onset\")\nax.set_xlabel(\"Time (s)\")\nax.set_ylabel(\"Frequency (Hz)\")\nax.set_title(\"TFR Difference (Right \u2212 Left) at C3\\nRed = more power for right; Blue = more power for left\")\nax.legend()\nplt.colorbar(im, ax=ax, label=\"Relative Power Difference\")\nplt.tight_layout()\nplt.show()",
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "id": "md28780131",
      "metadata": {},
      "source": "---\n## \ud83e\udde0 Interlude: The Motor Cortex Story\n\nBefore we extract features, let's make sure we understand what the TFR is showing us.\n\n### The mu rhythm \u2014 the motor cortex's idle state\n\nThe motor cortex produces a rhythmic oscillation at 8\u201312 Hz when it is at rest \u2014 called the **mu rhythm** (related to the alpha rhythm). When you move your hand \u2014 or even *imagine* moving it \u2014 the mu rhythm disappears. This disappearance is the **ERD** (Event-Related Desynchronisation).\n\n> \ud83d\udfe2 **Analogy:** The mu rhythm is like a car engine idling in neutral. When the driver decides to move \u2014 or even just *thinks* about moving \u2014 the engine engages and the rhythmic idle stops. When the driver parks again, the idle resumes (sometimes briefly overshooting \u2014 that is the beta rebound).\n\n### Why the lateralisation matters for BCI\n\n```\nRight-hand imagery  \u2192  LEFT motor cortex activates  \u2192  C3 shows ERD\nLeft-hand imagery   \u2192  RIGHT motor cortex activates \u2192  C4 shows ERD\n```\n\nBy comparing the mu/beta power at C3 and C4, a computer can predict which hand the subject was imagining \u2014 **without any movement happening**. This is the core of motor imagery BCI."
    },
    {
      "cell_type": "markdown",
      "id": "md13119488",
      "metadata": {},
      "source": "---\n## Step 4 \u2014 Feature Extraction for Machine Learning\n\n### From a matrix to a vector\n\nEach epoch is a 2D matrix: **64 channels \u00d7 561 time points = 35,904 numbers**. A classifier cannot work with this raw matrix. We need to summarise each epoch as a small set of informative numbers \u2014 a **feature vector**.\n\nHere we extract:\n- Alpha/mu band power (8\u201313 Hz) at 5 channels: C3, Cz, C4, FC3, FC4\n- Beta band power (13\u201330 Hz) at the same 5 channels\n\nThis gives **10 features per trial** \u2014 much more manageable than 35,904.\n\n> \ud83d\udfe2 **Analogy:** Summarising a 2-hour football match with statistics: goals, possession %, shots on target. You cannot watch 2 hours to compare 80 matches \u2014 but 5 numbers per match tell you what matters."
    },
    {
      "cell_type": "code",
      "id": "cd78626943",
      "metadata": {},
      "source": "def extract_band_power_features(epochs_obj, channels, bands_dict):\n    \"\"\"\n    Build a feature matrix from band-power across selected channels.\n\n    Returns\n    -------\n    X            : np.ndarray  shape (n_epochs, n_features)\n    feature_names: list of str\n    \"\"\"\n    psd = epochs_obj.compute_psd(method=\"welch\", fmin=1, fmax=45)\n    data_psd, freqs = psd.get_data(return_freqs=True)\n\n    ch_idx = [epochs_obj.ch_names.index(c) for c in channels]\n    features, feature_names = [], []\n\n    for band, (fmin, fmax) in bands_dict.items():\n        mask = (freqs >= fmin) & (freqs <= fmax)\n        bp   = data_psd[:, :, mask].mean(axis=-1)   # (n_epochs, n_channels)\n        for ci, ch_name in zip(ch_idx, channels):\n            features.append(bp[:, ci])\n            feature_names.append(f\"{band}_{ch_name}\")\n\n    return np.column_stack(features), feature_names\n\n\nsel_channels = [\"C3\", \"Cz\", \"C4\", \"FC3\", \"FC4\"]\nsel_bands    = {\"Alpha\": (8, 13), \"Beta\": (13, 30)}\n\nX_left,  feat_names = extract_band_power_features(epochs_left,  sel_channels, sel_bands)\nX_right, _          = extract_band_power_features(epochs_right, sel_channels, sel_bands)\n\nX = np.vstack([X_left,  X_right])\ny = np.hstack([np.zeros(len(X_left)), np.ones(len(X_right))])  # 0=left, 1=right\n\nprint(f\"Feature matrix X shape : {X.shape}  (trials \u00d7 features)\")\nprint(f\"Labels y shape         : {y.shape}\")\nprint(f\"Feature names: {feat_names}\")",
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "id": "cd28809472",
      "metadata": {},
      "source": "# Distribution plot \u2014 how separable are the features?\nfig, axes = plt.subplots(2, len(sel_channels), figsize=(14, 6))\n\nfor row_i, band in enumerate(sel_bands):\n    for col_i, ch in enumerate(sel_channels):\n        feat_label = f\"{band}_{ch}\"\n        fi  = feat_names.index(feat_label)\n        ax  = axes[row_i, col_i]\n\n        ax.hist(X[y == 0, fi], bins=12, alpha=0.6, label=\"Left\",  color=\"steelblue\")\n        ax.hist(X[y == 1, fi], bins=12, alpha=0.6, label=\"Right\", color=\"tomato\")\n\n        # t-test: is this feature significantly different between conditions?\n        t_stat, p_val = ttest_ind(X[y == 0, fi], X[y == 1, fi])\n        title_color = \"darkred\" if p_val < 0.05 else \"black\"\n        marker = \" *\" if p_val < 0.05 else \"\"\n        ax.set_title(f\"{feat_label}{marker}\", fontsize=8, color=title_color)\n        if col_i == 0:\n            ax.set_ylabel(f\"{band} Power\")\n        ax.legend(fontsize=6)\n\nplt.suptitle(\"Feature Distributions: Left (blue) vs Right (red) Hand Imagery\\n* = statistically significant (p < 0.05)\", fontsize=10)\nplt.tight_layout()\nplt.show()",
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "id": "md88901182",
      "metadata": {},
      "source": "> \ud83d\udc41 **What to check:**\n> - Features marked with **\\*** are statistically different between conditions\n> - Features at **C3** and **C4** should show the clearest separation (lateralised ERD)\n> - Overlapping histograms mean the feature alone is not enough \u2014 we need to combine features"
    },
    {
      "cell_type": "markdown",
      "id": "md35092921",
      "metadata": {},
      "source": "---\n## Step 5 \u2014 LDA Classification\n\n### What is LDA?\n\n**Linear Discriminant Analysis (LDA)** finds the straight line (in 2D) or hyperplane (in 10D) that best separates the two classes. It maximises the distance between the class means while minimising the spread within each class.\n\n> \ud83d\udfe2 **Analogy:** Two piles of gold and silver coins mixed on a table. LDA finds the dividing line that puts as many gold coins on one side and silver coins on the other as possible.\n\n### Why log-transform first?\n\nBand power values have a skewed distribution \u2014 a few trials with very high power. LDA assumes Gaussian distributions. The log transform compresses large values and makes the distribution more symmetric.\n\n### Cross-validation \u2014 honest accuracy\n\nNever test on the same data you trained on! **5-fold cross-validation** splits data into 5 groups; trains on 4, tests on 1, rotates 5 times. Every trial is used as a test exactly once. The final accuracy is the average.\n\n| Accuracy | Interpretation |\n|---|---|\n| ~50% | Chance level (guessing randomly) |\n| 55\u201365% | Weak signal \u2014 features need improvement |\n| 65\u201375% | Typical for band power + LDA |\n| 75\u201385% | Good \u2014 clear ERD present |\n| 85%+ | Excellent |"
    },
    {
      "cell_type": "code",
      "id": "cd94863619",
      "metadata": {},
      "source": "from sklearn.discriminant_analysis import LinearDiscriminantAnalysis\nfrom sklearn.model_selection import cross_val_score, StratifiedKFold\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.pipeline import Pipeline\n\n# Log-transform: makes distribution more Gaussian \u2192 better for LDA\nX_log = np.log10(X + 1e-30)\n\npipeline = Pipeline([\n    (\"scaler\", StandardScaler()),              # zero-mean, unit-variance per feature\n    (\"lda\",    LinearDiscriminantAnalysis()),  # find the best separating line\n])\n\ncv     = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)\nscores = cross_val_score(pipeline, X_log, y, cv=cv, scoring=\"accuracy\")\n\nprint(\"LDA cross-validation results:\")\nprint(f\"  Per-fold accuracies : {[f'{s:.1%}' for s in scores]}\")\nprint(f\"  Mean accuracy       : {scores.mean():.1%}  \u00b1  {scores.std():.1%}\")\nprint(f\"  Chance level        : 50.0%\")\n\nif scores.mean() > 0.60:\n    print(\"  \u2713  Above chance \u2014 the classifier is learning the lateralisation pattern\")\nelse:\n    print(\"  \u26a0\ufe0f  Near chance \u2014 check that TFR shows clear ERD difference\")",
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "id": "cd59978693",
      "metadata": {},
      "source": "# Decision boundary visualisation (using the first two features for 2D plot)\npipeline.fit(X_log, y)\n\nf0, f1 = 0, 1  # Alpha_C3 and Alpha_Cz\nx_min, x_max = X_log[:, f0].min() - 0.5, X_log[:, f0].max() + 0.5\ny_min, y_max = X_log[:, f1].min() - 0.5, X_log[:, f1].max() + 0.5\nxx, yy = np.meshgrid(np.linspace(x_min, x_max, 200),\n                     np.linspace(y_min, y_max, 200))\n\n# Fill the remaining dimensions with zeros for visualisation\ngrid_pts = np.c_[xx.ravel(), yy.ravel(),\n                 np.zeros((xx.size, X_log.shape[1] - 2))]\nZ = pipeline.predict(grid_pts).reshape(xx.shape)\n\nfig, ax = plt.subplots(figsize=(8, 6))\nax.contourf(xx, yy, Z, alpha=0.25, cmap=\"bwr\")\nax.scatter(X_log[y==0, f0], X_log[y==0, f1],\n           c=\"steelblue\", edgecolors=\"k\", lw=0.5, label=\"Left hand\",  s=50)\nax.scatter(X_log[y==1, f0], X_log[y==1, f1],\n           c=\"tomato\",    edgecolors=\"k\", lw=0.5, label=\"Right hand\", s=50)\nax.set_xlabel(f\"log\u2081\u2080 [{feat_names[f0]}]\")\nax.set_ylabel(f\"log\u2081\u2080 [{feat_names[f1]}]\")\nax.set_title(f\"LDA Decision Boundary (2D projection)\\nCV accuracy = {scores.mean():.1%}\")\nax.legend()\nax.grid(True, alpha=0.3)\nplt.tight_layout()\nplt.show()",
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "id": "md26602573",
      "metadata": {},
      "source": "---\n## Step 6 \u2014 CSP + LDA: The Motor Imagery Gold Standard\n\n### What is CSP?\n\nIn Step 4 we manually chose which channels to use (C3, Cz, C4, FC3, FC4). **CSP (Common Spatial Patterns)** finds the *optimal* combination of all 64 channels automatically from the data.\n\nCSP creates **virtual channels** (spatial filters) that:\n- **Maximise** variance for one condition (e.g. left hand)\n- **Minimise** variance for the other condition (right hand)\n\nThe features extracted from CSP components are much better at separating the two classes than manually chosen band power.\n\n> \ud83d\udfe2 **Analogy:** Instead of you choosing which microphone to use in a noisy room, CSP finds the ideal *mix* of all microphones that best picks up the signal you want.\n\n### Why CSP + LDA is the standard\n\nCSP + LDA has been the dominant motor imagery BCI method since the mid-2000s. Despite being simple, it consistently outperforms hand-crafted features in BCI competitions."
    },
    {
      "cell_type": "code",
      "id": "cd16533182",
      "metadata": {},
      "source": "from mne.decoding import CSP\nfrom sklearn.pipeline import Pipeline as SKPipeline\n\ncsp = CSP(n_components=4, log=True, reg=None)\n\npipeline_csp = SKPipeline([\n    (\"csp\", csp),\n    (\"lda\", LinearDiscriminantAnalysis()),\n])\n\n# CSP takes raw epoch data (not pre-computed features)\nX_raw = np.vstack([epochs_left.get_data(), epochs_right.get_data()])\n\nscores_csp = cross_val_score(pipeline_csp, X_raw, y, cv=cv, scoring=\"accuracy\")\n\nprint(\"CSP + LDA cross-validation results:\")\nprint(f\"  Per-fold accuracies : {[f'{s:.1%}' for s in scores_csp]}\")\nprint(f\"  Mean accuracy       : {scores_csp.mean():.1%}  \u00b1  {scores_csp.std():.1%}\")\nprint(f\"\\nComparison:\")\nprint(f\"  Band power LDA : {scores.mean():.1%}\")\nprint(f\"  CSP + LDA      : {scores_csp.mean():.1%}\")\nprint(f\"  Improvement    : {(scores_csp.mean() - scores.mean())*100:+.1f} percentage points\")",
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "id": "cd27371282",
      "metadata": {},
      "source": "# Fit CSP on full dataset and plot spatial patterns\npipeline_csp.fit(X_raw, y)\n\nfig = csp.plot_patterns(epochs.info, show=True,\n                        title=\"CSP Spatial Patterns \u2014 where does each component come from?\")\nplt.show()",
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "id": "md69665741",
      "metadata": {},
      "source": "> \ud83d\udc41 **What to check in the CSP patterns:**\n> - **Patterns 1\u20132** (maximise left-hand variance): focal activation near **C4** (right hemisphere)\n> - **Patterns 3\u20134** (maximise right-hand variance): focal activation near **C3** (left hemisphere)\n> - If patterns show no clear spatial structure, the dataset may not have enough trials or the ERD is weak\n>\n> Note: we are plotting `plot_patterns()` not `plot_filters()`. The patterns show *where the signal comes from* on the scalp \u2014 that is what has neurological meaning."
    },
    {
      "cell_type": "markdown",
      "id": "md41658979",
      "metadata": {},
      "source": "---\n## \u2705 Tutorial 3 \u2014 Summary\n\n| Step | MNE / Python Command | What it produces |\n|---|---|---|\n| 0 \u2014 Load | `mne.read_epochs()` | Clean epochs split by condition |\n| 1 \u2014 ERP | `epochs.average()`, `plot_topomap()` | Average waveform + scalp maps |\n| 2 \u2014 PSD | `compute_psd()`, `band_power()` | Power spectrum + band power values |\n| 3 \u2014 TFR | `tfr_morlet()`, `apply_baseline()` | Time \u00d7 frequency ERD/ERS map |\n| 4 \u2014 Features | `extract_band_power_features()` | 10-feature matrix per trial |\n| 5 \u2014 LDA | `Pipeline([StandardScaler, LDA])` | Accuracy via 5-fold CV |\n| 6 \u2014 CSP + LDA | `CSP()` + `Pipeline` | Better accuracy + spatial patterns |\n\n---\n\n### Key findings in motor imagery EEG\n\n- **ERD** (power decrease) in the mu (8\u201312 Hz) and beta (18\u201326 Hz) bands during imagery\n- ERD is **lateralised**: right-hand imagery \u2192 stronger ERD at C3 (left hemisphere)\n- **CSP + LDA** consistently outperforms hand-picked band power features\n- Typical accuracy: 65\u201380% for this single-subject dataset\n\n---\n\n### \ud83c\udfc6 The complete BCI pipeline\n\n```\nRaw EEG (Tutorial 1)\n    \u2193\nPreprocessing: filter, ICA, epoch (Tutorial 2)\n    \u2193\nERD analysis: TFR at C3/C4 (Tutorial 3, Steps 1\u20133)\n    \u2193\nFeature extraction: band power or CSP log-variance (Step 4\u20136)\n    \u2193\nClassification: LDA with cross-validation (Step 5\u20136)\n    \u2193\nPrediction: left or right hand imagery from new EEG\n```\n\nThis pipeline is the foundation of commercial EEG-based BCI systems used clinically for motor rehabilitation and assistive communication."
    }
  ]
}