-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathSSVEPClassification-SVM-RMDM-BasicCode
180 lines (143 loc) · 4.89 KB
/
SSVEPClassification-SVM-RMDM-BasicCode
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Feb 4 18:10:25 2024
Checkout for more detailed information (MNE)
https://mne.tools/stable/auto_tutorials/time-freq/50_ssvep.html
Machine Learning
https://github.com/pyRiemann/pyRiemann
@author: abinjacob
"""
#%% libraries
import mne
import numpy as np
import matplotlib.pyplot as plt
import os.path as op
import pyriemann
from pyriemann.estimation import Covariances
from pyriemann.tangentspace import TangentSpace
from sklearn.model_selection import cross_val_score
from sklearn.pipeline import make_pipeline
from sklearn.svm import SVC
#%% params
# filter
l_freq = 0.1
h_freq = None
# epoching
tmin = -0.2
tmax = 4
# Events
event_id = {'stim_L15': 13, 'stim_L20': 14, 'stim_R15': 15, 'stim_R20': 16}
event_names = list(event_id.keys())
foi = [15, 20, 15, 20] # Freqs of interest
# PSD computation
fmin = 1.0
fmax = 100.0
# Show filter
show_filter = False
#%% load data
rootpath = r'/Users/abinjacob/Documents/01. Calypso/EEG Coding/SSVEP/DATA'
# EEGLab file to load (.set)
filename = 'SSVEP_pilot2_rawdata.set'
filepath = op.join(rootpath,filename)
# load file in mne
raw = mne.io.read_raw_eeglab(filepath, eog= 'auto', preload= True)
a = raw.info
#Preprocess the data
# extracting events
events, _ = mne.events_from_annotations(raw, verbose= False)
epochs = mne.Epochs(
raw,
events= events,
event_id= [event_id['stim_L15'], event_id['stim_L20'], event_id['stim_R15'], event_id['stim_R20']],
tmin=tmin, tmax=tmax,
baseline= None,
preload= True,
reject={'eeg': 3.0}) # Reject epochs based on maximum peak-to-peak signal amplitude (PTP)
# No rejection due to very high value
#%% Frequency analysis - Calculate power spectral density (PSD)
sfreq = epochs.info["sfreq"]
spectrum = epochs.compute_psd(
method="welch",
n_fft=int(sfreq * (tmax - tmin)),
n_overlap=0,
n_per_seg=None,
tmin=tmin,
tmax=tmax,
fmin=fmin,
fmax=fmax,
window="boxcar",
verbose=False,
)
psds, freqs = spectrum.get_data(return_freqs=True)
print(f"Shape of psds (NumEpochs, NumChans, NumFreqBins): {psds.shape}")
# Plot PSD - all 4 conditions
# fig, axs = plt.subplots(2, 2, sharex="all", sharey="none", figsize=(8, 5))
# axs = axs.reshape(-1)
# freq_range = range(np.where(np.floor(freqs) == fmin)[0][0], np.where(np.ceil(freqs) == fmax - 1)[0][0])
# for i in range(0,4):
# # Get all events for specific event ID (condition)
# condition = event_names[i]
# idx = epochs.events[:,2]==event_id[condition]
# # Select only corresponding epochs
# psds_plot = 10 * np.log10(psds[idx,:,:])
# psds_mean = psds_plot.mean(axis=(0, 1))[freq_range] # mean over trials and channels
# psds_std = psds_plot.std(axis=(0, 1))[freq_range]
# # Plot
# axs[i].plot(freqs[freq_range], psds_mean, color="b")
# axs[i].fill_between(
# freqs[freq_range], psds_mean - psds_std, psds_mean + psds_std, color="b", alpha=0.2
# )
# axs[i].set(title=f"PSD spectrum {condition}", ylabel="Power Spectral Density [dB]")
# fig.show()
## plot Topoplot - all 4 Conditions
# fig, axs = plt.subplots(2, 2, figsize=(8, 5))
# axs = axs.reshape(-1)
# for i in range(0,4):
# # Get all events for specific event ID (condition)
# condition = event_names[i]
# idx = epochs.events[:,2]==event_id[condition]
# # Select only corresponding epochs
# # psds_plot = 10 * np.log10(psds[idx,:,:])
# psds_plot = psds[idx,:,:]
# psds_mean = psds_plot.mean(axis=(0))
# # Get condition-specific frequency
# freq_range = range(np.where(np.floor(freqs) == foi[i]-1)[0][0], np.where(np.ceil(freqs) == foi[i]+1)[0][0])
# # Mean over freq bins
# psds_mean = psds_mean[:,freq_range].mean(axis=1)
# # Plot
# mne.viz.plot_topomap(psds_mean, epochs.info, axes=axs[i])
# axs[i].set(title=f"PSD spectrum {condition}")
#%% SVM Classification
# Create a label vector
labels = epochs.events[:,2]
for i in range(0,len(labels)):
if labels[i]==13 or labels[i]==15:
labels[i] = 15
else:
labels[i] = 20
# Refine psds to frequency range around [12,25]
freq_range = range(np.where(np.floor(freqs) == 12)[0][0], np.where(np.ceil(freqs) == 25)[0][0])
# Mean over freq bins
X = psds[:,:,freq_range]
# X = psds
# X = epochs.get_data() # regularization necessary
y = labels # labels
# estimate covariance matrices
cov = pyriemann.estimation.Covariances().fit_transform(X)
# build your pipeline
covest = Covariances()
ts = TangentSpace()
svc = SVC(kernel='linear')
clf = make_pipeline(covest, ts, svc)
# cross validation
accuracy = cross_val_score(clf, X, y)
# print accuracy
print(accuracy.mean())
#%% Riemannian - MDM Classification
# estimate covariance matrices
cov = pyriemann.estimation.Covariances().fit_transform(X)
# cross validation
mdm = pyriemann.classification.MDM()
accuracy = cross_val_score(mdm, cov, y)
print(accuracy.mean())