.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/example_trca.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_example_trca.py: Task-related component analysis for SSVEP detection =================================================== Sample code for the task-related component analysis (TRCA)-based steady -state visual evoked potential (SSVEP) detection method [1]_. The filter bank analysis can also be combined to the TRCA-based algorithm [2]_ [3]_. This code is based on the Matlab implementation from: https://github.com/mnakanishi/TRCA-SSVEP Uses `meegkit.trca.TRCA()`. References ---------- .. [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. .. GENERATED FROM PYTHON SOURCE LINES 28-42 .. code-block:: Python # Authors: Giuseppe Ferraro # Nicolas Barascud import os import time import matplotlib.pyplot as plt import numpy as np import scipy.io from meegkit.trca import TRCA from meegkit.utils.trca import itr, normfit, round_half_up t = time.time() .. GENERATED FROM PYTHON SOURCE LINES 43-45 Parameters ----------------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 45-66 .. code-block:: Python dur_gaze = 0.5 # data length for target identification [s] delay = 0.13 # visual latency being considered in the analysis [s] n_bands = 5 # number of sub-bands in filter bank analysis is_ensemble = True # True = ensemble TRCA method; False = TRCA method alpha_ci = 0.05 # 100*(1-alpha_ci): confidence interval for accuracy sfreq = 250 # sampling rate [Hz] dur_shift = 0.5 # duration for gaze shifting [s] list_freqs = np.array( [[x + 8.0 for x in range(8)], [x + 8.2 for x in range(8)], [x + 8.4 for x in range(8)], [x + 8.6 for x in range(8)], [x + 8.8 for x in range(8)]]).T # list of stimulus frequencies n_targets = list_freqs.size # The number of stimuli # Useful variables (no need to modify) dur_gaze_s = round_half_up(dur_gaze * sfreq) # data length [samples] delay_s = round_half_up(delay * sfreq) # visual latency [samples] dur_sel_s = dur_gaze + dur_shift # selection time [s] ci = 100 * (1 - alpha_ci) # confidence interval .. GENERATED FROM PYTHON SOURCE LINES 67-69 Load data ----------------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 69-82 .. code-block:: Python path = os.path.join("..", "tests", "data", "trcadata.mat") eeg = scipy.io.loadmat(path)["eeg"] n_trials, n_chans, n_samples, n_blocks = eeg.shape # Convert dummy Matlab format to (sample, channels, trials) and construct # vector of labels eeg = np.reshape(eeg.transpose([2, 1, 3, 0]), (n_samples, n_chans, n_trials * n_blocks)) labels = np.array([x for x in range(n_targets)] * n_blocks) crop_data = np.arange(delay_s, delay_s + dur_gaze_s) eeg = eeg[crop_data] .. GENERATED FROM PYTHON SOURCE LINES 83-92 TRCA classification ----------------------------------------------------------------------------- Estimate classification performance with a Leave-One-Block-Out cross-validation approach. To get a sense of the filterbank specification in relation to the stimuli we can plot the individual filterbank sub-bands as well as the target frequencies (with their expected harmonics in the EEG spectrum). We use the filterbank specification described in [2]_. .. GENERATED FROM PYTHON SOURCE LINES 92-126 .. code-block:: Python filterbank = [[(6, 90), (4, 100)], # passband, stopband freqs [(Wp), (Ws)] [(14, 90), (10, 100)], [(22, 90), (16, 100)], [(30, 90), (24, 100)], [(38, 90), (32, 100)], [(46, 90), (40, 100)], [(54, 90), (48, 100)]] f, ax = plt.subplots(1, figsize=(7, 4)) for i, _band in enumerate(filterbank): ax.axvspan(ymin=i / len(filterbank) + .02, ymax=(i + 1) / len(filterbank) - .02, xmin=filterbank[i][1][0], xmax=filterbank[i][1][1], alpha=0.2, facecolor=f"C{i}") ax.axvspan(ymin=i / len(filterbank) + .02, ymax=(i + 1) / len(filterbank) - .02, xmin=filterbank[i][0][0], xmax=filterbank[i][0][1], alpha=0.5, label=f"sub-band{i}", facecolor=f"C{i}") for f in list_freqs.flat: colors = np.ones((9, 4)) colors[:, :3] = np.linspace(0, .5, 9)[:, None] ax.scatter(f * np.arange(1, 10), [f] * 9, c=colors, s=8, zorder=100) ax.set_ylabel("Stimulus frequency (Hz)") ax.set_xlabel("EEG response frequency (Hz)") ax.set_xlim([0, 102]) ax.set_xticks(np.arange(0, 100, 10)) ax.grid(True, ls=":", axis="x") ax.legend(bbox_to_anchor=(1.05, .5), fontsize="small") plt.tight_layout() plt.show() .. image-sg:: /auto_examples/images/sphx_glr_example_trca_001.png :alt: example trca :srcset: /auto_examples/images/sphx_glr_example_trca_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 127-128 Now perform the TRCA-based SSVEP detection algorithm .. GENERATED FROM PYTHON SOURCE LINES 128-168 .. code-block:: Python trca = TRCA(sfreq, filterbank, is_ensemble) print("Results of the ensemble TRCA-based method:\n") accs = np.zeros(n_blocks) itrs = np.zeros(n_blocks) for i in range(n_blocks): # Select all folds except one for training traindata = np.concatenate( (eeg[..., :i * n_trials], eeg[..., (i + 1) * n_trials:]), 2) y_train = np.concatenate( (labels[:i * n_trials], labels[(i + 1) * n_trials:]), 0) # Construction of the spatial filter and the reference signals trca.fit(traindata, y_train) # Test stage testdata = eeg[..., i * n_trials:(i + 1) * n_trials] y_test = labels[i * n_trials:(i + 1) * n_trials] estimated = trca.predict(testdata) # Evaluation of the performance for this fold (accuracy and ITR) is_correct = estimated == y_test accs[i] = np.mean(is_correct) * 100 itrs[i] = itr(n_targets, np.mean(is_correct), dur_sel_s) print(f"Block {i}: accuracy = {accs[i]:.1f}, \tITR = {itrs[i]:.1f}") # Mean accuracy and ITR computation mu, _, muci, _ = normfit(accs, alpha_ci) print(f"\nMean accuracy = {mu:.1f}%\t({ci:.0f}% CI: {muci[0]:.1f}-{muci[1]:.1f}%)") # noqa mu, _, muci, _ = normfit(itrs, alpha_ci) print(f"Mean ITR = {mu:.1f}\t({ci:.0f}% CI: {muci[0]:.1f}-{muci[1]:.1f})") if is_ensemble: ensemble = "ensemble TRCA-based method" else: ensemble = "TRCA-based method" print(f"\nElapsed time: {time.time()-t:.1f} seconds") .. rst-class:: sphx-glr-script-out .. code-block:: none Results of the ensemble TRCA-based method: Block 0: accuracy = 97.5, ITR = 301.3 Block 1: accuracy = 100.0, ITR = 319.3 Block 2: accuracy = 95.0, ITR = 286.3 Block 3: accuracy = 95.0, ITR = 286.3 Block 4: accuracy = 95.0, ITR = 286.3 Block 5: accuracy = 100.0, ITR = 319.3 Mean accuracy = 97.1% (95% CI: 97.0-97.1%) Mean ITR = 299.8 (95% CI: 299.4-300.2) Elapsed time: 12.1 seconds .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 12.104 seconds) .. _sphx_glr_download_auto_examples_example_trca.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: example_trca.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: example_trca.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_