任务判别成分分析法（TDCA）\\（TDCA）main-python程序运行

# -*- coding: utf-8 -*- # # Authors: Swolf <swolfforever@gmail.com> # Date: 2021/10/10 # License: MIT License """ Task Decomposition Component Analysis. """ from typing import Optional, List import numpy as np from scipy.linalg import qr from scipy.stats import pearsonr from numpy import ndarray from sklearn.base import BaseEstimator, TransformerMixin, ClassifierMixin from base import FilterBankSSVEP """ -*- coding: utf-8 -*- DSP: Discriminal Spatial Patterns Authors: Swolf <swolfforever@gmail.com> Junyang Wang <2144755928@qq.com> Last update date: 2022-8-11 License: MIT License """ from typing import Optional, List, Tuple from itertools import combinations import numpy as np from scipy.linalg import eigh from numpy import ndarray from sklearn.base import BaseEstimator, TransformerMixin, ClassifierMixin from base import robust_pattern from base import FilterBankSSVEP def xiang_dsp_kernel( X: ndarray, y: ndarray ) -> Tuple[ndarray, ndarray, ndarray, ndarray]: """ DSP: Discriminal Spatial Patterns, only for two classes[1] -Author: Swolf <swolfforever@gmail.com> -Created on: 2021-1-07 -Update log: Parameters ---------- X : ndarray EEG data assuming removing mean, shape (n_trials, n_channels, n_samples) y : ndarray labels of EEG data, shape (n_trials, ) Returns ------- W: ndarray spatial filters, shape (n_channels, n_filters) D: ndarray eigenvalues in descending order M: ndarray template for all classes, shape (n_channel, n_samples) A: ndarray spatial patterns, shape (n_channels, n_filters) Notes ----- the implementation removes regularization on within-class scatter matrix Sw. References ---------- [1] Liao, Xiang, et al. "Combining spatial filters for the classification of single-trial EEG in a finger movement task." IEEE Transactions on Biomedical Engineering 54.5 (2007): 821-831. """ X, y = np.copy(X), np.copy(y) labels = np.unique(y) X = np.reshape(X, (-1, *X.shape[-2:])) X = X - np.mean(X, axis=-1, keepdims=True) # the number of each label n_labels = np.array([np.sum(y == label) for label in labels]) # average template of all trials M = np.mean(X, axis=0) # class conditional template Ms, Ss = zip( *[ ( np.mean(X[y == label], axis=0), np.sum( np.matmul(X[y == label], np.swapaxes(X[y == label], -1, -2)), axis=0 ), ) for label in labels ] ) Ms, Ss = np.stack(Ms), np.stack(Ss) # within-class scatter matrix Sw = np.sum( Ss - n_labels[:, np.newaxis, np.newaxis] * np.matmul(Ms, np.swapaxes(Ms, -1, -2)), axis=0, ) Ms = Ms - M # between-class scatter matrix Sb = np.sum( n_labels[:, np.newaxis, np.newaxis] * np.matmul(Ms, np.swapaxes(Ms, -1, -2)), axis=0, ) D, W = eigh(Sb, Sw) ix = np.argsort(D)[::-1] # in descending order D, W = D[ix], W[:, ix] A = robust_pattern(W, Sb, W.T @ Sb @ W) return W, D, M, A def xiang_dsp_feature( W: ndarray, M: ndarray, X: ndarray, n_components: int = 1 ) -> ndarray: """ Return DSP features in paper [1] -Author: Swolf <swolfforever@gmail.com> -Created on: 2021-1-07 -Update log: Parameters ---------- W : ndarray spatial filters from csp_kernel, shape (n_channels, n_filters) M: ndarray common template for all classes, shape (n_channel, n_samples) X : ndarray eeg test data, shape (n_trials, n_channels, n_samples) n_components : int, optional length of the spatial filters, first k components to use, by default 1 Returns ------- features: ndarray features, shape (n_trials, n_components, n_samples) Raises ------ ValueError n_components should less than half of the number of channels Notes ----- 1. instead of meaning of filtered signals in paper [1]_., we directly return filtered signals. References ---------- [1] Liao, Xiang, et al. "Combining spatial filters for the classification of single-trial EEG in a finger movement task." IEEE Transactions on Biomedical Engineering 54.5 (2007): 821-831. """ W, M, X = np.copy(W), np.copy(M), np.copy(X) max_components = W.shape[1] if n_components > max_components: raise ValueError("n_components should less than the number of channels") X = np.reshape(X, (-1, *X.shape[-2:])) X = X - np.mean(X, axis=-1, keepdims=True) features = np.matmul(W[:, :n_components].T, X - M) return features class DSP(BaseEstimator, TransformerMixin, ClassifierMixin): """ DSP: Discriminal Spatial Patterns -Author: Swolf <swolfforever@gmail.com> -Created on: 2021-1-07 -Update log: Parameters ---------- n_components : int length of the spatial filter, first k components to use, by default 1 transform_method : str method of template matching, by default ’corr‘ (pearson correlation coefficient) classes_ : int number of the EEG classes """ def __init__(self, n_components: int = 1, transform_method: str = "corr"): self.n_components = n_components self.transform_method = transform_method def fit(self, X: ndarray, y: ndarray, Yf: Optional[ndarray] = None): """ import the train data to get a model Parameters ---------- X : ndarray train data, shape(n_trials, n_channels, n_samples) y : ndarray labels of train data, shape (n_trials, ) Yf : ndarray optional parameter Returns ------- W_: ndarray spatial filters, shape (n_channels, n_filters), in which n_channels = n_filters D_: ndarray eigenvalues in descending order, shape (n_filters, ) M_: ndarray template for all classes, shape (n_channel, n_samples) A_: ndarray spatial patterns, shape (n_channels, n_filters) templates_ : ndarray templates of train data, shape (n_channels, n_filters, n_samples) """ X -= np.mean(X, axis=-1, keepdims=True) self.classes_ = np.unique(y) self.W_, self.D_, self.M_, self.A_ = xiang_dsp_kernel(X, y) self.templates_ = np.stack( [ np.mean( xiang_dsp_feature( self.W_, self.M_, X[y == label], n_components=self.W_.shape[1] ), axis=0, ) for label in self.classes_ ] ) return self def transform(self, X: ndarray): """ import the test data to get features Parameters ---------- X : ndarray test data, shape(n_trials, n_channels, n_samples) Returns ------- feature : ndarray features of test data, shape(n_trials, n_classes) """ n_components = self.n_components X -= np.mean(X, axis=-1, keepdims=True) features = xiang_dsp_feature(self.W_, self.M_, X, n_components=n_components) if self.transform_method is None: return features.reshape((features.shape[0], -1)) elif self.transform_method == "mean": return np.mean(features, axis=-1) elif self.transform_method == "corr": return self._pearson_features( features, self.templates_[:, :n_components, :] ) else: raise ValueError("non-supported transform method") def _pearson_features(self, X: ndarray, templates: ndarray): """ pearson correlation coefficient Parameters ---------- X : ndarray features of test data after sp