Add Task Related Component Analysis (TRCA)

README.md

@@ -91,3 +91,22 @@ If you use this, you should cite the following article:
 [1] Cohen, M. X., & Gulbinaite, R. (2017). Rhythmic entrainment source separation: Optimizing analyses
     of neural responses to rhythmic sensory stimulation. Neuroimage, 147, 43-56.
 ```
+
+### 4. Task-Related Component Analysis (TRCA)
+
+This code is based on the [Matlab implementation from Masaki Nakanishi](https://github.com/mnakanishi/TRCA-SSVEP), and was adapted to python by [Giuseppe Ferraro](mailto:giuseppe.ferraro@isae-supaero.fr)
+
+If you use this, you should cite the following articles:
+
+```sql
+[1] M. Nakanishi, Y. Wang, X. Chen, Y.-T. Wang, X. Gao, and T.-P. Jung,
+    "Enhancing detection of SSVEPs for a high-speed brain speller using
+    task-related component analysis", IEEE Trans. Biomed. Eng, 65(1): 104-112,
+    2018.
+[2] X. Chen, Y. Wang, S. Gao, T. -P. Jung and X. Gao, "Filter bank
+    canonical correlation analysis for implementing a high-speed SSVEP-based
+    brain-computer interface", J. Neural Eng., 12: 046008, 2015.
+[3] X. Chen, Y. Wang, M. Nakanishi, X. Gao, T. -P. Jung, S. Gao,
+    "High-speed spelling with a noninvasive brain-computer interface",
+    Proc. Int. Natl. Acad. Sci. U. S. A, 112(44): E6058-6067, 2015.
+```

doc/index.rst

@@ -26,6 +26,7 @@ Contents
    ~meegkit.ress
    ~meegkit.sns
    ~meegkit.star
+   ~meegkit.trca
    ~meegkit.tspca
    ~meegkit.utils
 

doc/modules/meegkit.utils.rst

@@ -71,9 +71,19 @@ Signal
 
 ----
 
-
 Statistics
 ----------
 .. automodule:: meegkit.utils.stats
 
    .. autosummary::
+
+
+|
+
+----
+
+TRCA utilities
+--------------
+.. automodule:: meegkit.utils.trca
+
+   .. autosummary::

examples/example_asr.py

@@ -11,7 +11,6 @@
 import matplotlib.pyplot as plt
 
 from meegkit.asr import ASR
-from meegkit.utils.asr import yulewalk_filter
 from meegkit.utils.matrix import sliding_window
 
 # THIS_FOLDER = os.path.dirname(os.path.abspath(__file__))

examples/example_trca.ipynb

@@ -0,0 +1,258 @@
+{
+  "cells": [
+    {
+      "cell_type": "code",
+      "execution_count": 1,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "%matplotlib inline"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "\n",
+        "# Task-related component analysis (TRCA)-based SSVEP detection\n",
+        "\n",
+        "Sample code for the task-related component analysis (TRCA)-based steady\n",
+        "-state visual evoked potential (SSVEP) detection method [1]_. The filter\n",
+        "bank analysis [2, 3]_ can also be combined to the TRCA-based algorithm.\n",
+        "\n",
+        "Uses meegkit.trca.TRCA()\n",
+        "\n",
+        "References:\n",
+        "\n",
+        ".. [1] M. Nakanishi, Y. Wang, X. Chen, Y.-T. Wang, X. Gao, and T.-P. Jung,\n",
+        "   \"Enhancing detection of SSVEPs for a high-speed brain speller using\n",
+        "   task-related component analysis\", IEEE Trans. Biomed. Eng, 65(1): 104-112,\n",
+        "   2018.\n",
+        ".. [2] X. Chen, Y. Wang, S. Gao, T. -P. Jung and X. Gao, \"Filter bank\n",
+        "   canonical correlation analysis for implementing a high-speed SSVEP-based\n",
+        "   brain-computer interface\", J. Neural Eng., 12: 046008, 2015.\n",
+        ".. [3] X. Chen, Y. Wang, M. Nakanishi, X. Gao, T. -P. Jung, S. Gao,\n",
+        "   \"High-speed spelling with a noninvasive brain-computer interface\",\n",
+        "   Proc. Int. Natl. Acad. Sci. U. S. A, 112(44): E6058-6067, 2015.\n",
+        "\n",
+        "This code is based on the Matlab implementation from\n",
+        "https://github.com/mnakanishi/TRCA-SSVEP\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": 2,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "# Author: Giuseppe Ferraro <giuseppe.ferraro@isae-supaero.fr>\n",
+        "import os\n",
+        "import time\n",
+        "\n",
+        "import numpy as np\n",
+        "import scipy.io\n",
+        "from meegkit.trca import TRCA\n",
+        "from meegkit.utils.trca import itr, normfit, round_half_up\n",
+        "\n",
+        "t = time.time()"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Parameters\n",
+        "\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": 3,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [
+        {
+          "output_type": "stream",
+          "name": "stdout",
+          "text": [
+            "Results of the ensemble TRCA-based method:\n\n"
+          ]
+        }
+      ],
+      "source": [
+        "len_gaze_s = 0.5  # data length for target identification [s]\n",
+        "len_delay_s = 0.13  # visual latency being considered in the analysis [s]\n",
+        "n_bands = 5  # number of sub-bands in filter bank analysis\n",
+        "is_ensemble = True  # True = ensemble TRCA method; False = TRCA method\n",
+        "alpha_ci = 0.05   # 100*(1-alpha_ci): confidence interval for accuracy\n",
+        "sfreq = 250  # sampling rate [Hz]\n",
+        "len_shift_s = 0.5  # duration for gaze shifting [s]\n",
+        "list_freqs = np.concatenate(\n",
+        "    [[x + 8 for x in range(8)],\n",
+        "     [x + 8.2 for x in range(8)],\n",
+        "     [x + 8.4 for x in range(8)],\n",
+        "     [x + 8.6 for x in range(8)],\n",
+        "     [x + 8.8 for x in range(8)]])  # list of stimulus frequencies\n",
+        "n_targets = len(list_freqs)  # The number of stimuli\n",
+        "\n",
+        "# Preparing useful variables (DONT'T need to modify)\n",
+        "len_gaze_smpl = round_half_up(len_gaze_s * sfreq)  # data length [samples]\n",
+        "len_delay_smpl = round_half_up(len_delay_s * sfreq)  # visual latency [samples]\n",
+        "len_sel_s = len_gaze_s + len_shift_s  # selection time [s]\n",
+        "ci = 100 * (1 - alpha_ci)  # confidence interval\n",
+        "\n",
+        "# Performing the TRCA-based SSVEP detection algorithm\n",
+        "print('Results of the ensemble TRCA-based method:\\n')"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Load data\n",
+        "\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": 4,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "path = os.path.join('..', 'tests', 'data', 'trcadata.mat')\n",
+        "mat = scipy.io.loadmat(path)\n",
+        "eeg = mat[\"eeg\"]\n",
+        "\n",
+        "n_trials = eeg.shape[0]\n",
+        "n_chans = eeg.shape[1]\n",
+        "n_samples = eeg.shape[2]\n",
+        "n_blocks = eeg.shape[3]\n",
+        "\n",
+        "# Convert dummy Matlab format to (sample, channels, trials) and construct\n",
+        "# vector of labels\n",
+        "eeg = np.reshape(eeg.transpose([2, 1, 3, 0]),\n",
+        "                 (n_samples, n_chans, n_trials * n_blocks))\n",
+        "labels = np.array([x for x in range(n_targets)] * n_blocks)\n",
+        "\n",
+        "crop_data = np.arange(len_delay_smpl, len_delay_smpl + len_gaze_smpl)\n",
+        "eeg = eeg[crop_data]"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## TRCA classification\n",
+        "Estimate classification performance with a Leave-One-Block-Out\n",
+        "cross-validation approach.\n",
+        "\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": 5,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [
+        {
+          "output_type": "stream",
+          "name": "stdout",
+          "text": [
+            "Block 0: accuracy = 97.5, \tITR = 301.3\n",
+            "Block 1: accuracy = 100.0, \tITR = 319.3\n",
+            "Block 2: accuracy = 95.0, \tITR = 286.3\n",
+            "Block 3: accuracy = 95.0, \tITR = 286.3\n",
+            "Block 4: accuracy = 95.0, \tITR = 286.3\n",
+            "Block 5: accuracy = 100.0, \tITR = 319.3\n",
+            "\n",
+            "Mean accuracy = 97.1%\t(95% CI: 97.0-97.1%)\n",
+            "Mean ITR = 299.8\t(95% CI: 299.4-300.2%)\n",
+            "\n",
+            "Elapsed time: 13.8 seconds\n"
+          ]
+        }
+      ],
+      "source": [
+        "# We use the filterbank specification described in [2]_.\n",
+        "filterbank = [[(6, 90), (4, 100)],  # passband freqs, stopband freqs (Wp, Ws)\n",
+        "              [(14, 90), (10, 100)],\n",
+        "              [(22, 90), (16, 100)],\n",
+        "              [(30, 90), (24, 100)],\n",
+        "              [(38, 90), (32, 100)],\n",
+        "              [(46, 90), (40, 100)],\n",
+        "              [(54, 90), (48, 100)]]\n",
+        "trca = TRCA(sfreq, filterbank, is_ensemble)\n",
+        "\n",
+        "accs = np.zeros(n_blocks)\n",
+        "itrs = np.zeros(n_blocks)\n",
+        "for i in range(n_blocks):\n",
+        "\n",
+        "    # Training stage\n",
+        "    traindata = eeg.copy()\n",
+        "\n",
+        "    # Select all folds except one for training\n",
+        "    traindata = np.concatenate(\n",
+        "        (traindata[..., :i * n_trials],\n",
+        "         traindata[..., (i + 1) * n_trials:]), 2)\n",
+        "    y_train = np.concatenate(\n",
+        "        (labels[:i * n_trials], labels[(i + 1) * n_trials:]), 0)\n",
+        "\n",
+        "    # Construction of the spatial filter and the reference signals\n",
+        "    trca.fit(traindata, y_train)\n",
+        "\n",
+        "    # Test stage\n",
+        "    testdata = eeg[..., i * n_trials:(i + 1) * n_trials]\n",
+        "    y_test = labels[i * n_trials:(i + 1) * n_trials]\n",
+        "    estimated = trca.predict(testdata)\n",
+        "\n",
+        "    # Evaluation of the performance for this fold (accuracy and ITR)\n",
+        "    is_correct = estimated == y_test\n",
+        "    accs[i] = np.mean(is_correct) * 100\n",
+        "    itrs[i] = itr(n_targets, np.mean(is_correct), len_sel_s)\n",
+        "    print(f\"Block {i}: accuracy = {accs[i]:.1f}, \\tITR = {itrs[i]:.1f}\")\n",
+        "\n",
+        "# Mean accuracy and ITR computation\n",
+        "mu, _, muci, _ = normfit(accs, alpha_ci)\n",
+        "print()\n",
+        "print(f\"Mean accuracy = {mu:.1f}%\\t({ci:.0f}% CI: {muci[0]:.1f}-{muci[1]:.1f}%)\")  # noqa\n",
+        "\n",
+        "mu, _, muci, _ = normfit(itrs, alpha_ci)\n",
+        "print(f\"Mean ITR = {mu:.1f}\\t({ci:.0f}% CI: {muci[0]:.1f}-{muci[1]:.1f}%)\")\n",
+        "if is_ensemble:\n",
+        "    ensemble = 'ensemble TRCA-based method'\n",
+        "else:\n",
+        "    ensemble = 'TRCA-based method'\n",
+        "\n",
+        "print(f\"\\nElapsed time: {time.time()-t:.1f} seconds\")"
+      ]
+    }
+  ],
+  "metadata": {
+    "kernelspec": {
+      "name": "python388jvsc74a57bd0d64e410d98a0dc7c6b3fb09ececfc32281268599ac952adfc85e199a2f396698",
+      "display_name": "Python 3.8.8 64-bit ('miniconda3': virtualenv)"
+    },
+    "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.8.8-final"
+    }
+  },
+  "nbformat": 4,
+  "nbformat_minor": 0
+}