OSCR

Body surface potential mapping of the cortico-muscular axis using smart textile electrode arrays.

Code ↔ Paper

19 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 19 matches
  1. [1] § Methods › Feature extraction ↔ Wireless_BSPM.ipynb, lines 1290–1370 · score 0.95 · overlapping windows, power spectral density, 0.5–4 Hz, 30–100 Hz, 8–13 Hz, 0.5–100 Hz
  2. [2] § Methods › EEG signal processing and analytical procedures › Spectral decomposition and periodic vs. aperiodic components ↔ Additional_figures.ipynb, lines 560–637 · score 0.93 · peak width limits, 1–100 Hz, 1–12 Hz, peak height, Power spectral densities, aperiodic
  3. [3] § Methods › Preprocessing of recorded electrophysiological signals ↔ Wireless_BSPM.ipynb, lines 1290–1370 · score 0.85 · power spectral density, 8–13 Hz, 13–30 Hz, band power, EEG signals, frequency bands
  4. [4] § Methods › Classification and performance evaluation ↔ Wireless_BSPM.ipynb, lines 288–344 · score 0.83 · Confusion matrices, Logistic regression, F1 score, cross validation, single channel, precision
  5. [5] § Methods › PLS regression (EEG → EMG mapping) ↔ Wireless_BSPM.ipynb, lines 3107–3136 · score 0.80 · PLSRegression, EMG maps, EEG features, RMSE, squares, components
  6. [6] § Methods › EEG signal processing and analytical procedures › Evoked vs. induced activity ↔ Additional_figures.ipynb, lines 825–954 · score 0.80 · flattened ERP, confusion matrices, logistic regression, lbfgs, multinomial, solver
  7. [7] § Results and Discussion › Spatiotemporal features for enhanced classification ↔ Wireless_BSPM.ipynb, lines 1434–1505 · score 0.77 · logistic regression, F1 score, Event related, motor activity, precision, recall
  8. [8] § Methods › Feature importance analysis using correlation circles ↔ Wireless_BSPM.ipynb, lines 399–516 · score 0.76 · correlation circle coordinates, unit circle, standardised feature matrix, PCA, variance, component
  9. [9] § Results and Discussion › Spatiotemporal features for enhanced classification ↔ Wireless_BSPM.ipynb, lines 288–344 · score 0.76 · repeated cross validation, logistic regression, F1 score, single channel, precision, recall
  10. [10] § Methods › EEG signal processing and analytical procedures › Event-Related Desynchronisation/Synchronisation (ERD/ERS) ↔ Additional_figures.ipynb, lines 1030–1148 · score 0.75 · movement onset, 8–12 Hz, 13–30 Hz, Hilbert, ERD, ERS
  11. [11] § Methods › Classification and performance evaluation ↔ Wireless_BSPM.ipynb, lines 3138–3185 · score 0.69 · Spearman rank correlation, Pearson correlation, standard deviation, cosine, predictions
  12. [12] § Methods › Machine learning model descriptions ↔ Wireless_BSPM.ipynb, lines 953–1009 · score 0.69 · repeated stratified, F1 score, cross validation, precision, recall, fold
  13. [13] § Methods › EEG signal processing and analytical procedures › Event-Related Desynchronisation/Synchronisation (ERD/ERS) ↔ Additional_figures.ipynb, lines 1030–1148 · score 0.69 · broadband envelope, ERD, ERS, Gaussian, Epochs, onset
  14. [14] § Methods › Machine learning model training ↔ Wireless_BSPM.ipynb, lines 3138–3185 · score 0.67 · Spearman rank correlation, Pearson correlation, cosine, Predictive
  15. [15] § Methods › Feature importance analysis using correlation circles ↔ Wireless_BSPM.ipynb, lines 399–516 · score 0.62 · correlation circles, feature matrices, PCA, variance, component, SNR
  16. [16] § Results and Discussion › Mapping and prediction of cortico-muscular axis dynamics ↔ Wireless_BSPM.ipynb, lines 2691–2744 · score 0.61 · Gaussian fits, hand flexion, wrist extension, Histograms, Spatio temporal, movement
  17. [17] § Methods › EEG signal processing and analytical procedures › Evoked Response Potentials (ERPs) ↔ Additional_figures.ipynb, lines 439–557 · score 0.60 · Anterior Frontal, SEM, Parietal, Occipital, ERPs, Evoked
  18. [18] § Results and Discussion › Mapping and prediction of cortico-muscular axis dynamics ↔ Wireless_BSPM.ipynb, lines 2804–2809 · score 0.53 · 15–30 Hz, beta band, cortico muscular, components, Event, 15 Hz
  19. [19] § Results and Discussion › Mapping and prediction of cortico-muscular axis dynamics ↔ Wireless_BSPM.ipynb, lines 2691–2744 · score 0.52 · hand flexion, wrist extension, Gaussian, 600 ms, movements, reaction

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Jupyter notebook · 3,675 lines · 148 KB · no license · 14 matches

  1. # %% [markdown]
  2. # # Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology
  4. # %% [markdown]
  5. # Ruben Ruiz-Mateos Serrano <br>
  6. # Start date : 10/01/2025
  8. # %% [markdown]
  9. # ## BSPM & ML
  11. # %% [markdown]
  12. # ### Test 16 channel BSPM recordings (with INTAN)
  14. # %% [markdown]
  15. # #### Data import
  17. # %%
  18. import os
  19. import re
  21. from biolab.bspm import BSPM,image2layout
  22. from biolab.utils import apply2all
  24. # Save new recordings into .bspm files
  25. # Define parameters
  26. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Recordings'
  27. # Convert electrode array image into layout
  28. ECG_lyt = image2layout(filename=r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\Masks\Layouts\ECG_lyt.jpg',seq='cw')
  29. EMG_lyt = image2layout(filename=r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\Masks\Layouts\EMG_lyt.jpg',seq='cw')
  30. #EEG_lyt = image2layout(filename=r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\Masks\Layouts\EEG_lyt.jpg',seq='cw')
  31. # Define operation to be applied to each file
  32. def func(filepath):
  33. # Extract filename from filepath
  34. filename = os.path.splitext(os.path.basename(filepath))[0]
  35. # Define the regex pattern to match the required values
  36. pattern = r'^(([^\s]+)_M(\d+))_\d{6}_\d{6}$'
  37. # Search for the pattern in the given filename
  38. match = re.search(pattern,filename)
  39. if match:
  40. filename = match.group(1)
  41. mode = match.group(2) # The type of signal (i.e. ECG, EMG, EEG)
  42. measurement = match.group(3) # The number following 'M'
  43. else:
  44. raise ValueError("Filename \"{}\" does not match the expected pattern.".format(filename))
  46. meta = {
  47. 'mode': mode.split('_')[0],
  48. 'measurement': measurement
  49. }
  50. # Load BSPM data from .rhs file
  51. if mode == 'ECG_BSPM':
  52. MP = BSPM.from_file(filepath=filepath,layout=ECG_lyt,metadata=meta,refchs=[])
  53. elif mode == 'EMG_BSPM':
  54. MP = BSPM.from_file(filepath=filepath,layout=EMG_lyt,metadata=meta,refchs=[])
  55. #elif mode == 'EEG_BSPM':
  56. #MP = BSPM.from_file(filepath=filepath,layout=EEG_lyt,metadata=meta,refchs=[])
  58. MP.save(savepath=r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\\'+filename+'.bspm')
  60. # Call DataSaver function
  61. apply2all(read_directory=readpath,extension='.rhs',func=func)
  63. # %% [markdown]
  64. # #### Electrode array impedance analysis
  66. # %% [markdown]
  67. # Average electrode impedance across all measurements and channels for each signal mode
  69. # %%
  70. import pandas as pd
  71. import matplotlib.pyplot as plt
  73. # Import .bspm files
  74. ECG = BSPM.load(loadpath=r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\ECG_BSPM_M2.bspm')
  75. df_ECG = ECG.channel_data.loc[:,['electrode_impedance_magnitude','custom_order']]
  76. EMG = BSPM.load(loadpath=r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\EMG_BSPM_M2.bspm')
  77. df_EMG = EMG.channel_data.loc[:,['electrode_impedance_magnitude','custom_order']]
  78. # Add a 'mode' column to each DataFrame
  79. df_ECG['mode'] = 'ECG'
  80. df_EMG['mode'] = 'EMG'
  81. # Concatenate the DataFrames
  82. df = pd.concat([df_ECG,df_EMG],ignore_index=True)
  84. # %% [markdown]
  85. # Violin plot of impedance distribution for all modes
  87. # %%
  88. import seaborn as sns
  90. sns.violinplot(x='mode',y='electrode_impedance_magnitude',data=df,cut=0)
  91. plt.title('')
  92. plt.xlabel('Mode')
  93. plt.ylabel('Impedance [$\Omega$]')
  94. plt.savefig(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Results\Impedance_distribution_by_mode.pdf')
  95. plt.show()
  97. # %% [markdown]
  98. # #### Body Surface Potential Mapping
  100. # %%
  101. from biolab.bspm import BSPM
  103. # Import .bspm files
  104. ECG = BSPM.load(loadpath=r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\ECG_BSPM_M1.bspm')
  105. EMG = BSPM.load(loadpath=r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\EMG_BSPM_M2.bspm')
  107. # %%
  108. import matplotlib.pyplot as plt
  110. from biolab.bspm import interpolate
  112. vmap = ECG.potential(time=10)
  113. frame = interpolate(valmap=vmap)
  114. plt.imshow(frame,vmin=-2000,vmax=2000)
  115. cbar = plt.colorbar()
  116. cbar.set_label('Voltage [$\mu$V]') # Setting the label for the colorbar
  117. # Remove the axes
  118. plt.axis('off')
  120. # %% [markdown]
  121. # ### Shape classification from muscular BSPM
  123. # %% [markdown]
  124. # #### Data import
  126. # %%
  127. import os
  128. import re
  130. from biolab.bspm import BSPM,image2layout
  131. from biolab.utils import apply2all
  133. # Save new recordings into .bspm files
  134. # Define parameters
  135. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Recordings\EMG_classification'
  136. # Convert electrode array image into layout
  137. EMG_lyt = image2layout(filename=r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\Masks\Layouts\EMG_lyt.jpg',seq='cw')
  138. # Define operation to be applied to each file
  139. def func(filepath):
  140. # Extract filename from filepath
  141. filename = os.path.splitext(os.path.basename(filepath))[0]
  142. # Define the regex pattern to match the required values
  143. pattern = r'^([A-Za-z]+)_([A-Za-z]+)_(\d+)_\d{6}_\d{6}$'
  144. # Search for the pattern in the given filename
  145. match = re.search(pattern,filename)
  146. if match:
  147. filename = match.group(1)+'_'+match.group(2)+'_'+match.group(3)
  148. mode = match.group(1) # type of signal (i.e. ECG, EMG, EEG)
  149. shape = match.group(2) # shape of the object
  150. measurement = match.group(3) # measurement
  151. else:
  152. raise ValueError("Filename \"{}\" does not match the expected pattern.".format(filename))
  154. meta = {
  155. 'mode': mode,
  156. 'measurement': measurement,
  157. 'shape': {'ball':'sphere','cylinder':'cylinder','emptyhand':'no object'}.get(shape,None)
  158. }
  159. # Load BSPM data from .rhs file
  160. MP = BSPM.from_file(filepath=filepath,layout=EMG_lyt,metadata=meta,refchs=[])
  161. # Save BSPM object
  162. MP.save(savepath=r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\EMG_classification\\'+filename+'.bspm')
  164. # Call DataSaver function
  165. apply2all(read_directory=readpath,extension='.rhs',func=func)
  167. # %% [markdown]
  168. # #### SNR computation
  170. # %%
  171. import pandas as pd
  172. import numpy as np
  174. from biolab.utils import apply2all,snr
  175. from biolab.bspm import BSPM
  177. # Define read directory
  178. rdir = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\EMG_classification'
  179. # Define `apply2all` function
  180. def fun(filepath):
  181. # Load BSPM files
  182. MP = BSPM.load(loadpath=filepath)
  183. # Resample data
  184. newfs = 800
  185. rMP = MP.resample(newfs=newfs)
  186. # Preprocess data
  187. rMP.preprocess(bw=(5,399))
  188. # Read DataFrame with existing data
  189. excelpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\EMG_classification\SNR.xlsx'
  190. with pd.ExcelWriter(excelpath,engine='openpyxl',mode='a',if_sheet_exists='overlay') as writer:
  191. existing_data = pd.read_excel(excelpath)
  192. updated_data = existing_data
  193. # Compute SNR for each channel (5s baseline + 5s contraction)
  194. for col in rMP.preprocessed_data.columns:
  195. acc = 0
  196. # Get the slices as lists of absolute values
  197. arrays = [
  198. np.abs(rMP.preprocessed_data[col][0.5:4.5].values),
  199. np.abs(rMP.preprocessed_data[col][10.5:14.5].values),
  200. np.abs(rMP.preprocessed_data[col][20.5:24.5].values),
  201. np.abs(rMP.preprocessed_data[col][30.5:34.5].values),
  202. np.abs(rMP.preprocessed_data[col][40.5:44.5].values),
  203. np.abs(rMP.preprocessed_data[col][50.5:54.5].values),
  204. np.abs(rMP.preprocessed_data[col][60.5:64.5].values),
  205. np.abs(rMP.preprocessed_data[col][70.5:74.5].values),
  206. np.abs(rMP.preprocessed_data[col][80.5:84.5].values),
  207. np.abs(rMP.preprocessed_data[col][90.5:94.5].values)
  208. ]
  209. # Find the minimum length of all arrays
  210. min_length = min(len(arr) for arr in arrays)
  211. # Trim each array to the minimum length
  212. arrays_trimmed = [arr[:min_length] for arr in arrays]
  213. # Compute the mean
  214. baseline = np.mean(arrays_trimmed,axis=0)
  215. while (acc+10)*newfs<len(rMP.preprocessed_data[col].values):
  216. data = np.abs(rMP.preprocessed_data[col][5+acc:10+acc].values)
  217. snrval = snr(s=data,n=baseline)
  218. new_data = {'SNR':snrval,'shape':rMP.metadata['shape'],'measurement':rMP.metadata['measurement'],'envelope':acc/10+1,'channel':col}
  219. updated_data = updated_data.append(new_data,ignore_index=True)
  220. acc = acc+10
  221. # Write the updated DataFrame back to the Excel file, overwriting the old data
  222. updated_data.to_excel(writer,index=False)
  225. apply2all(read_directory=rdir,extension='.bspm',func=fun)
  227. # %% [markdown]
  228. # #### SNR values inspection
  230. # %%
  231. import pandas as pd
  232. import seaborn as sns
  233. import matplotlib.pyplot as plt
  235. # Import SNR data from Excel
  236. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\EMG_classification\SNR.xlsx'
  237. df = pd.read_excel(io=readpath)
  238. # Display SNR as a function of shape, channel and measurement
  239. # Create a figure with subplots
  240. fig, axes = plt.subplots(1,3,figsize=(18,6))
  241. # Plot 1: Index by shape
  242. df = df.sort_values(by='shape')
  243. sns.violinplot(x='shape',y='SNR',data=df,order=df['shape'].unique(),ax=axes[0])
  244. axes[0].set_title('')
  245. axes[0].set_xlabel('Shape')
  246. axes[0].set_ylabel('SNR')
  247. # Plot 2: Index by channel
  248. df = df.sort_values(by='channel')
  249. df['channel'] = df['channel'].astype(int)
  250. sns.violinplot(x='channel',y='SNR',data=df,order=df['channel'].unique(),ax=axes[1])
  251. axes[1].set_title('')
  252. axes[1].set_xlabel('Channel')
  253. axes[1].set_ylabel('SNR')
  254. # Plot 3: Index by measurement
  255. df = df.sort_values(by='measurement')
  256. sns.violinplot(x='measurement',y='SNR',data=df,order=df['measurement'].unique(),ax=axes[2])
  257. axes[2].set_title('')
  258. axes[2].set_xlabel('Measurement')
  259. axes[2].set_ylabel('SNR')
  260. # Adjust layout to prevent overlap of titles and labels
  261. plt.tight_layout()
  262. # Display SNR as a function of channel for each shape
  263. # Create a figure with subplots
  264. fig, axes = plt.subplots(1,3,figsize=(18,6))
  265. # Plot 1: Index by shape
  266. df_sphere = df[df['shape']=='sphere'].sort_values(by='channel')
  267. sns.violinplot(x='channel',y='SNR',data=df_sphere,order=df_sphere['channel'].unique(),ax=axes[0])
  268. axes[0].set_title('Sphere')
  269. axes[0].set_xlabel('Channel')
  270. axes[0].set_ylabel('SNR')
  271. # Plot 2: Index by channel
  272. df_cylinder = df[df['shape']=='cylinder'].sort_values(by='channel')
  273. df_cylinder['channel'] = df_cylinder['channel'].astype(int)
  274. sns.violinplot(x='channel',y='SNR',data=df_cylinder,order=df_cylinder['channel'].unique(),ax=axes[1])
  275. axes[1].set_title('Cylinder')
  276. axes[1].set_xlabel('Channel')
  277. axes[1].set_ylabel('SNR')
  278. # Plot 3: Index by measurement
  279. df_noobject = df[df['shape']=='no object'].sort_values(by='channel')
  280. sns.violinplot(x='channel',y='SNR',data=df_noobject,order=df_noobject['channel'].unique(),ax=axes[2])
  281. axes[2].set_title('No object')
  282. axes[2].set_xlabel('Channel')
  283. axes[2].set_ylabel('SNR')
  284. # Adjust layout to prevent overlap of titles and labels
  285. plt.tight_layout()
  288. # %% [markdown]
  289. # *There doesn't seem to be a significant difference in SNR across different shapes. This probably indicates inability for a single electrode to perform accurate classification. When displaying SNR as a function of channel for different shapes, a difference in the distribution of SNR can be observed.*
  291. # %% [markdown]
  292. # #### Machine learning classification of shapes
  294. # %% [markdown]
  295. # Single electrode (control - middle electrode)
  297. # %%
  298. import pandas as pd
  299. from sklearn.model_selection import train_test_split
  300. from sklearn.linear_model import LogisticRegression
  301. from sklearn.model_selection import RepeatedStratifiedKFold, cross_val_score
  302. from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay, make_scorer,precision_score,recall_score,f1_score
  303. import matplotlib.pyplot as plt
  305. # Import SNR data from Excel
  306. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Cortico-muscular axis\Experiments\INTAN\Data\EMG_classification\SNR.xlsx'
  307. df = pd.read_excel(io=readpath)
  308. # Generate feature matrix from data
  309. fmatrix = df[df['channel']==9].drop(['measurement','envelope','channel'],axis=1)
  310. # Separate features and labels
  311. X = fmatrix.drop('shape',axis=1) # feature matrix
  312. y = fmatrix['shape'] # labels
  313. # Split the data into training and testing sets
  314. X_train,X_test,y_train,y_test = train_test_split(X,y,test_size=0.2,stratify=fmatrix['shape'])
  315. # Initialise and fit the logistic regression
  316. clf = LogisticRegression()
  317. clf.fit(X_train,y_train)
  318. # Predict the labels for the test set
  319. y_pred = clf.predict(X_test)
  320. # Calculate the confusion matrix
  321. cm = confusion_matrix(y_test,y_pred)
  322. # Display the confusion matrix
  323. disp = ConfusionMatrixDisplay(confusion_matrix=cm, display_labels=clf.classes_)
  324. fig,ax = plt.subplots(figsize=(6,6)) # Adjust figure size if needed
  325. # Set the color range limits
  326. vmin,vmax = 0,20 # Adjust these values as needed
  327. disp.plot(cmap="Reds",ax=ax,colorbar=False,values_format="d")
  328. disp.im_.set_clim(vmin,vmax) # Set the colorbar limits
  329. # Adjust the colorbar size
  330. cbar = ax.figure.colorbar(disp.im_,fraction=0.046,pad=0.04)
  331. cbar.set_label("Count") # Label for better readability
  332. #plt.savefig(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Cortico-muscular axis\Experiments\INTAN\Results\Confusion_matrices\Single_channel_EMG',dpi=1000)
  333. # Repeated Cross-Validation Accuracy
  334. rkf = RepeatedStratifiedKFold(n_splits=10,n_repeats=5)
  335. accuracy_scores = cross_val_score(clf,X,y,cv=rkf)
  336. precision_scores = cross_val_score(clf, X, y, cv=rkf, scoring=make_scorer(precision_score, average='macro'))
  337. recall_scores = cross_val_score(clf, X, y, cv=rkf, scoring=make_scorer(recall_score, average='macro'))
  338. f1_scores = cross_val_score(clf, X, y, cv=rkf, scoring=make_scorer(f1_score, average='macro'))
  340. # Print results
  341. print(f'Repeated Cross-Validation Accuracy: {accuracy_scores.mean():.2f} ± {accuracy_scores.std():.2f}')
  342. print(f"Repeated Cross-Validation Precision: {precision_scores.mean():.2f} ± {precision_scores.std():.2f}")
  343. print(f"Repeated Cross-Validation Recall: {recall_scores.mean():.2f} ± {recall_scores.std():.2f}")
  344. print(f"Repeated Cross-Validation F1 Score: {f1_scores.mean():.2f} ± {f1_scores.std():.2f}")
  347. # %% [markdown]
  348. # Electrode array
  350. # %%
  351. import pandas as pd
  352. from sklearn.model_selection import train_test_split
  353. from sklearn.linear_model import LogisticRegression
  354. from sklearn.model_selection import RepeatedStratifiedKFold, cross_val_score
  355. from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay
  356. import matplotlib.pyplot as plt
  358. # Import SNR data from Excel
  359. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Cortico-muscular axis\Experiments\INTAN\Data\EMG_classification\SNR.xlsx'
  360. df = pd.read_excel(io=readpath)
  361. # Generate feature matrix from data
  362. # Pivot the DataFrame so each channel becomes its own column
  363. fmatrix = df.pivot(index=['shape','envelope','measurement'],columns='channel',values='SNR')
  364. # Rename the columns to match the desired format
  365. fmatrix.columns = [f'snr_channel_{col}' for col in fmatrix.columns]
  366. # Reset index for a tidy DataFrame
  367. fmatrix.reset_index(inplace=True)
  368. # Remove unnecessary columns
  369. fmatrix = fmatrix.drop(columns=['envelope','measurement'])
  370. # Separate features and labels
  371. X = fmatrix.drop('shape',axis=1) # feature matrix
  372. y = fmatrix['shape'] # labels
  373. # Split the data into training and testing sets
  374. X_train,X_test,y_train,y_test = train_test_split(X,y,test_size=0.2,stratify=fmatrix['shape'])
  375. # Initialise and fit the logistic regression
  376. clf = LogisticRegression()
  377. clf.fit(X_train,y_train)
  378. # Predict the labels for the test set
  379. y_pred = clf.predict(X_test)
  380. # Calculate the confusion matrix
  381. cm = confusion_matrix(y_test,y_pred)
  382. # Display the confusion matrix
  383. disp = ConfusionMatrixDisplay(confusion_matrix=cm,display_labels=clf.classes_)
  384. fig,ax = plt.subplots(figsize=(6,6)) # Adjust figure size if needed
  385. # Set the color range limits
  386. vmin,vmax = 0,20 # Adjust these values as needed
  387. disp.plot(cmap="Reds",ax=ax,colorbar=False,values_format="d")
  388. disp.im_.set_clim(vmin,vmax) # Set the colorbar limits
  389. # Adjust the colorbar size
  390. cbar = ax.figure.colorbar(disp.im_,fraction=0.046,pad=0.04)
  391. cbar.set_label("Count") # Label for better readability
  392. #plt.savefig(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Cortico-muscular axis\Experiments\INTAN\Results\Confusion_matrices\BSPM_EMG',dpi=1000)
  393. # Repeated Cross-Validation Accuracy
  394. rkf = RepeatedStratifiedKFold(n_splits=10,n_repeats=5)
  395. scores = cross_val_score(clf,X,y,cv=rkf)
  396. print(f'Repeated Cross-Validation Accuracy: {scores.mean():.2f} ± {scores.std():.2f}')
  399. # %% [markdown]
  400. # Correlation circle for feature importance
  402. # %%
  403. import pandas as pd
  404. import numpy as np
  405. import matplotlib.pyplot as plt
  406. import matplotlib.cm as cm
  408. from sklearn.preprocessing import StandardScaler
  409. from sklearn.decomposition import PCA
  411. # ------------------------------------------------
  412. # Load data
  413. # ------------------------------------------------
  415. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Cortico-muscular axis\Experiments\INTAN\Data\EMG_classification\SNR.xlsx'
  417. df = pd.read_excel(io=readpath)
  419. # ------------------------------------------------
  420. # Build feature matrix (same as your ML pipeline)
  421. # ------------------------------------------------
  423. fmatrix = df.pivot(
  424. index=['shape','envelope','measurement'],
  425. columns='channel',
  426. values='SNR'
  427. )
  429. # Rename columns
  430. fmatrix.columns = [f'snr_channel_{col}' for col in fmatrix.columns]
  432. # Reset index
  433. fmatrix.reset_index(inplace=True)
  435. # Remove unnecessary columns
  436. fmatrix = fmatrix.drop(columns=['envelope','measurement'])
  438. # Separate features and labels
  439. X = fmatrix.drop('shape', axis=1)
  440. y = fmatrix['shape']
  442. # ------------------------------------------------
  443. # Standardise features (IMPORTANT for PCA)
  444. # ------------------------------------------------
  446. scaler = StandardScaler()
  447. X_scaled = scaler.fit_transform(X)
  449. # ------------------------------------------------
  450. # PCA
  451. # ------------------------------------------------
  453. pca = PCA(n_components=2)
  454. X_pca = pca.fit_transform(X_scaled)
  456. # Compute loadings (correlation circle coordinates)
  457. loadings = pca.components_.T * np.sqrt(pca.explained_variance_)
  459. # ------------------------------------------------
  460. # Plot correlation circle
  461. # ------------------------------------------------
  463. fig, ax = plt.subplots(figsize=(8,8), constrained_layout=True)
  465. # Draw unit circle
  466. circle = plt.Circle((0,0),1,fill=False,linestyle='--', color='black')
  467. ax.add_artist(circle)
  469. # Create colour mapping by channel index
  470. num_features = len(X.columns)
  471. colors = cm.viridis(np.linspace(0,1,num_features))
  473. channel_numbers = [int(col.split('_')[-1]) for col in X.columns]
  475. # Plot arrows (THICKER)
  476. for i, feature in enumerate(X.columns):
  478. ax.arrow(
  479. 0, 0,
  480. loadings[i,0],
  481. loadings[i,1],
  482. color=colors[i],
  483. linewidth=2.0, # <-- thicker arrows
  484. head_width=0.03, # slightly larger arrowhead
  485. length_includes_head=True
  486. )
  488. # Add colorbar
  489. norm = plt.Normalize(min(channel_numbers), max(channel_numbers))
  490. sm = plt.cm.ScalarMappable(cmap='viridis', norm=norm)
  491. sm.set_array([])
  493. cbar = fig.colorbar(sm, ax=ax)
  494. cbar.set_label('Channel number')
  496. # Axis formatting (BLACK axes)
  497. ax.axhline(0, color='black', linewidth=1.5)
  498. ax.axvline(0, color='black', linewidth=1.5)
  500. # Make frame/spines black
  501. for spine in ax.spines.values():
  502. spine.set_color('black')
  503. spine.set_linewidth(1.5)
  505. # Tick colors
  506. ax.tick_params(colors='black')
  508. ax.set_xlim(-1.1,1.1)
  509. ax.set_ylim(-1.1,1.1)
  510. ax.set_aspect('equal')
  512. ax.set_xlabel(f'PC1 ({pca.explained_variance_ratio_[0]*100:.1f}% variance)')
  513. ax.set_ylabel(f'PC2 ({pca.explained_variance_ratio_[1]*100:.1f}% variance)')
  514. ax.set_title('Correlation Circle (Muscular BSPM)')
  516. plt.show()
  519. # %% [markdown]
  520. # #### Muscular BSPM
  522. # %%
  523. import matplotlib.pyplot as plt
  524. import numpy as np
  526. from biolab.bspm import interpolate,BSPM
  528. # Read .bspm data for sphere
  529. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\EMG_classification\EMG_ball_1.bspm'
  530. EMG = BSPM.load(loadpath=readpath)
  531. rEMG_sphere = EMG.resample(newfs=800)
  532. # Display muscular BSPM
  533. vmap = rEMG_sphere.potential(time=78.70023333333333)
  534. abmap = np.abs(vmap)
  535. frame = interpolate(valmap=abmap,method='cubic')
  536. plt.imshow(frame,vmin=0,vmax=1500,cmap='viridis')
  537. cbar = plt.colorbar()
  538. cbar.set_label('Voltage [$\mu$V]') # Setting the label for the colorbar
  539. # Remove the axes
  540. plt.axis('off')
  541. plt.savefig(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Results\Figure_3\Sphere_viridis.pdf',dpi=1000)
  542. plt.show()
  544. # Read .bspm data for cylinder
  545. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\EMG_classification\EMG_cylinder_1.bspm'
  546. EMG = BSPM.load(loadpath=readpath)
  547. rEMG_cylinder = EMG.resample(newfs=800)
  548. # Display muscular BSPM
  549. vmap = rEMG_cylinder.potential(time=50.11526666666666)
  550. abmap = np.abs(vmap)
  551. frame = interpolate(valmap=abmap,method='cubic')
  552. plt.imshow(frame,vmin=0,vmax=1500,cmap='viridis')
  553. cbar = plt.colorbar()
  554. cbar.set_label('Voltage [$\mu$V]') # Setting the label for the colorbar
  555. # Remove the axes
  556. plt.axis('off')
  557. plt.savefig(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Results\Figure_3\Cylinder_viridis.pdf',dpi=1000)
  558. plt.show()
  560. # Read .bspm data for no object
  561. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\EMG_classification\EMG_emptyhand_1.bspm'
  562. EMG = BSPM.load(loadpath=readpath)
  563. rEMG_noobject = EMG.resample(newfs=800)
  564. # Display muscular BSPM
  565. vmap = rEMG_noobject.potential(time=35.5237)
  566. abmap = np.abs(vmap)
  567. frame = interpolate(valmap=abmap,method='cubic')
  568. plt.imshow(frame,vmin=0,vmax=1500,cmap='viridis')
  569. cbar = plt.colorbar()
  570. cbar.set_label('Voltage [$\mu$V]') # Setting the label for the colorbar
  571. # Remove the axes
  572. plt.axis('off')
  573. plt.savefig(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Results\Figure_3\No_object_viridis.pdf',dpi=1000)
  574. plt.show()
  577. # %% [markdown]
  578. # ### Position classification from cardiac BSPM
  580. # %% [markdown]
  581. # #### Data import
  583. # %%
  584. import os
  585. import re
  587. from biolab.bspm import BSPM,image2layout
  588. from biolab.utils import apply2all
  590. # Save new recordings into .bspm files
  591. # Define parameters
  592. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Recordings\ECG_classification'
  593. # Convert electrode array image into layout
  594. ECG_lyt = image2layout(filename=r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\Masks\Layouts\ECG_lyt.jpg',seq='cw')
  595. # Define operation to be applied to each file
  596. def func(filepath):
  597. # Extract filename from filepath
  598. filename = os.path.splitext(os.path.basename(filepath))[0]
  599. # Define the regex pattern to match the required values
  600. pattern = r'^([A-Za-z]+)_([A-Za-z-]+)_M(\d+)_\d{6}_\d{6}$'
  601. # Search for the pattern in the given filename
  602. match = re.search(pattern,filename)
  603. if match:
  604. filename = match.group(1)+'_'+match.group(2)
  605. mode = match.group(1) # type of signal (i.e. ECG, EMG, EEG)
  606. position = match.group(2) # position of the subject
  607. measurement = match.group(3) # measurement
  608. else:
  609. raise ValueError("Filename \"{}\" does not match the expected pattern.".format(filename))
  611. meta = {
  612. 'mode': mode,
  613. 'measurement': measurement,
  614. 'position': {'standing':'standing','sitting':'sitting','lying-side':'sideways','lying-up':'lying down'}.get(position,None)
  615. }
  616. # Load BSPM data from .rhs file
  617. MP = BSPM.from_file(filepath=filepath,layout=ECG_lyt,metadata=meta,refchs=[])
  618. # Save BSPM object
  619. MP.save(savepath=r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\ECG_classification\\'+filename+'.bspm')
  621. # Call DataSaver function
  622. apply2all(read_directory=readpath,extension='.rhs',func=func)
  624. # %% [markdown]
  625. # #### Feature extraction
  627. # %% [markdown]
  628. # RR intervals from all channels
  630. # %%
  631. import wfdb.processing
  632. import numpy as np
  633. import pandas as pd
  635. from biolab.bspm import BSPM
  636. from biolab.utils import apply2all
  638. # Define read directory
  639. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\ECG_classification'
  640. # Define apply2all function
  641. def fun(filepath):
  642. # Define new dataframe
  643. df = pd.DataFrame()
  645. # Load .bspm data
  646. ECG = BSPM.load(loadpath=filepath)
  648. # Resample and preprocess the data
  649. rECG = ECG.resample(newfs=1000)
  650. rECG.preprocess(bw=(0.5,20))
  652. for col in rECG.preprocessed_data.columns:
  653. # Find R-R intervals
  654. rpeaks = wfdb.processing.xqrs_detect(rECG.preprocessed_data[col][0:119.75].values,fs=1000,verbose=False)
  655. rrintv = wfdb.processing.calc_rr(qrs_locs=rpeaks)
  656. # Check the current number of rows in the dataframe
  657. current_len = len(df)
  658. new_len = len(rrintv)
  659. # If the DataFrame has more rows than rrintv and the DataFrame is not empty
  660. if (current_len<new_len) & (current_len!=0):
  661. # Calculate the mean of rrintv
  662. rrintv_mean = np.mean(rrintv)
  663. # Calculate the absolute deviation from the mean
  664. deviations = np.abs(rrintv-rrintv_mean)
  665. # Get the indices of the most different (largest deviations) values
  666. indices_to_remove = np.argsort(deviations)[:(new_len-current_len)] # Select top largest deviations
  667. # Remove excess values by excluding them
  668. rrintv = np.delete(rrintv,indices_to_remove)
  669. # If rrintv is shorter than the DataFrame, truncate the DataFrame to match the length of rrintv
  670. elif current_len>new_len:
  671. # Remove rows from the dataframe to match rrintv length
  672. df = df.iloc[:new_len]
  674. df['RRinterval_ch{}'.format(col)] = rrintv
  676. df['cycle'] = np.arange(1,len(rrintv)+1)
  677. df['position'] = rECG.metadata['position']
  678. # Append the DataFrame to an existing Excel file
  679. excel_path = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\ECG_classification\RR_interval_feature_matrix.xlsx'
  680. # Load the existing data from the Excel file
  681. try:
  682. existing_df = pd.read_excel(excel_path)
  683. # Append the new data to the existing data
  684. combined_df = pd.concat([existing_df,df],ignore_index=True)
  685. except FileNotFoundError:
  686. # If the file doesn't exist, just use the new data
  687. combined_df = df
  689. # Write the combined data back to the Excel file
  690. with pd.ExcelWriter(excel_path,mode='w',engine='openpyxl') as writer:
  691. combined_df.to_excel(writer,index=False)
  693. # Apply function to all files in read directory
  694. apply2all(read_directory=readpath,extension='.bspm',func=fun)
  697. # %% [markdown]
  698. # Delays from all channels
  700. # %%
  701. import wfdb.processing
  702. import numpy as np
  703. import pandas as pd
  705. from biolab.bspm import BSPM
  706. from biolab.utils import apply2all
  708. # Define read directory
  709. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\ECG_classification'
  710. # Define apply2all function
  711. def fun(filepath):
  712. # Load .bspm data
  713. ECG = BSPM.load(loadpath=filepath)
  714. # Resample and preprocess the data
  715. rECG = ECG.resample(newfs=2000)
  716. rECG.preprocess(bw=(0.5,20))
  717. # Find R peaks
  718. peaks_per_channel = {}
  719. for col in rECG.preprocessed_data.columns:
  720. peaks_per_channel[col] = wfdb.processing.xqrs_detect(rECG.preprocessed_data[col][0:119.75].values,fs=2000,verbose=False)
  722. # Get the peak indices for channel 1
  723. channel_1_peaks = peaks_per_channel[rECG.preprocessed_data.columns[0]]
  724. # Initialise list to store delays for valid cycles
  725. delays = []
  726. # Define ECG peak delay tolerance
  727. tolerance = 0.01
  728. # Calculate delays for each cycle
  729. for peak_1 in channel_1_peaks:
  730. cycle_delays = []
  731. valid_cycle = True
  732. for col in rECG.preprocessed_data.columns[1:]: # Skip channel 1, compare all others to it
  733. channel_peaks = peaks_per_channel[col]
  734. # Find the closest peak in the current channel within the tolerance window
  735. valid_peaks = np.abs(channel_peaks-peak_1)<=tolerance*rECG.metadata['amplifier_sample_rate']
  736. if np.any(valid_peaks): # If there's a valid peak
  737. closest_peak = channel_peaks[valid_peaks][0] # Get the first valid peak
  738. delay = (closest_peak-peak_1)/rECG.metadata['amplifier_sample_rate'] # Convert to seconds
  739. cycle_delays.append(delay)
  740. else:
  741. valid_cycle = False
  742. break # If any channel does not have a valid peak for this cycle, skip it
  744. # If the cycle was valid, store the delays
  745. if valid_cycle:
  746. delays.append(cycle_delays)
  748. # Convert delays to DataFrame
  749. if delays: # Only create DataFrame if there are valid cycles
  750. delay_df = pd.DataFrame(delays)
  751. delay_df.columns = [f'Delay_ch{col}' for col in rECG.preprocessed_data.columns[1:]]
  752. else:
  753. delay_df = pd.DataFrame() # Return an empty DataFrame if no valid cycles are found
  755. delay_df['position'] = rECG.metadata['position']
  757. # Append the DataFrame to an existing Excel file
  758. excel_path = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\ECG_classification\Delay_feature_matrix.xlsx'
  759. # Load the existing data from the Excel file
  760. try:
  761. existing_df = pd.read_excel(excel_path)
  762. # Append the new data to the existing data
  763. combined_df = pd.concat([existing_df,delay_df],ignore_index=True)
  764. except FileNotFoundError:
  765. # If the file doesn't exist, just use the new data
  766. combined_df = delay_df
  768. # Write the combined data back to the Excel file
  769. with pd.ExcelWriter(excel_path,mode='w',engine='openpyxl') as writer:
  770. combined_df.to_excel(writer,index=False)
  772. # Apply function to all files in read directory
  773. apply2all(read_directory=readpath,extension='.bspm',func=fun)
  776. # %% [markdown]
  777. # #### Feature values inspection
  779. # %% [markdown]
  780. # RR intervals
  782. # %%
  783. import pandas as pd
  784. import seaborn as sns
  785. import matplotlib.pyplot as plt
  787. # Import SNR data from Excel
  788. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\ECG_classification\RR_interval_feature_matrix.xlsx'
  789. df = pd.read_excel(io=readpath)
  790. # Display RR interval as a function of position
  791. df = df.sort_values(by='position')
  792. sns.violinplot(x='position',y='RRinterval_ch9',data=df,order=df['position'].unique())
  793. plt.title('')
  794. plt.xlabel('Position')
  795. plt.ylabel('RR Interval [$\mu$s]')
  796. # Adjust layout to prevent overlap of titles and labels
  797. plt.tight_layout()
  798. plt.show()
  801. # %% [markdown]
  802. # *Mild differences between lying and upward positions can be observed. No clear distinction can be observed between sitting and standing or lying down or sideways.*
  804. # %% [markdown]
  805. # Delays
  807. # %%
  808. import pandas as pd
  809. import seaborn as sns
  810. import numpy as np
  811. import matplotlib.pyplot as plt
  813. # Import SNR data from Excel
  814. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\ECG_classification\Delay_feature_matrix.xlsx'
  815. df = pd.read_excel(io=readpath)
  816. # Display RR interval as a function of position
  817. df = df.sort_values(by='position')
  818. sns.violinplot(x='position',y='Delay_ch16',data=df,order=df['position'].unique())
  819. plt.title('')
  820. plt.xlabel('Position')
  821. plt.ylabel('Delay CH1 - CH9 [s]')
  822. # Adjust layout to prevent overlap of titles and labels
  823. plt.tight_layout()
  824. plt.show()
  825. # Generate delay plots
  826. # Group by position and compute average delay and standard deviation
  827. avg = df.groupby(by='position').mean()
  828. stdv = df.groupby(by='position').std()
  829. avg = avg[['Delay_ch2','Delay_ch3','Delay_ch4','Delay_ch5','Delay_ch6','Delay_ch7','Delay_ch8','Delay_ch9','Delay_ch10','Delay_ch11','Delay_ch12','Delay_ch13','Delay_ch14','Delay_ch15','Delay_ch16']]
  830. stdv = stdv[['Delay_ch2','Delay_ch3','Delay_ch4','Delay_ch5','Delay_ch6','Delay_ch7','Delay_ch8','Delay_ch9','Delay_ch10','Delay_ch11','Delay_ch12','Delay_ch13','Delay_ch14','Delay_ch15','Delay_ch16']]
  831. plt.plot(np.arange(2,17),avg.iloc[0,:].values)
  832. plt.fill_between(np.arange(2,17),avg.iloc[0,:].values-stdv.iloc[0,:].values,avg.iloc[0,:].values+stdv.iloc[0,:].values,color='blue',alpha=0.2)
  833. plt.plot(np.arange(2,17),avg.iloc[1,:].values)
  834. plt.fill_between(np.arange(2,17),avg.iloc[1,:].values-stdv.iloc[1,:].values,avg.iloc[1,:].values+stdv.iloc[1,:].values,color='red',alpha=0.2)
  835. plt.plot(np.arange(2,17),avg.iloc[2,:].values)
  836. plt.fill_between(np.arange(2,17),avg.iloc[2,:].values-stdv.iloc[2,:].values,avg.iloc[2,:].values+stdv.iloc[2,:].values,color='green',alpha=0.2)
  837. plt.plot(np.arange(2,17),avg.iloc[3,:].values)
  838. plt.fill_between(np.arange(2,17),avg.iloc[3,:].values-stdv.iloc[3,:].values,avg.iloc[3,:].values+stdv.iloc[3,:].values,color='orange',alpha=0.2)
  839. plt.show()
  841. # %% [markdown]
  842. # *Clear differences can be observed between delays for all four positions and across channels.*
  844. # %%
  845. import numpy as np
  846. import pandas as pd
  847. import matplotlib.pyplot as plt
  849. # Import SNR data from Excel
  850. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\ECG_classification\Delay_feature_matrix.xlsx'
  851. df = pd.read_excel(io=readpath)
  852. # Set 'Position' as the index
  853. df.set_index('position',inplace=True)
  854. # Group by 'Position' and calculate mean and std for each channel
  855. grouped = df.groupby('position').agg(['mean','std'])
  856. # Channels as columns (without multi-indexing)
  857. channels = df.columns
  858. num_channels = len(channels)
  859. # Prepare the angle array for the radar chart (evenly distributed for each channel)
  860. angles = np.linspace(0,2*np.pi,num_channels,endpoint=False).tolist()
  861. # Close the loop by repeating the first angle
  862. angles += angles[:1]
  863. # Create a figure for the radar chart
  864. fig, ax = plt.subplots(figsize=(6,6),subplot_kw=dict(polar=True))
  865. # Loop through each position and plot the radar chart
  866. color = {'lying down':'blue','sideways':'red','sitting':'green','standing':'orange'}
  867. for position in grouped.index:
  868. # Calculate the mean and standard deviation for each channel
  869. mean_delays = grouped.loc[position,(channels,'mean')].values
  870. std_delays = grouped.loc[position,(channels,'std')].values
  871. # Close the loop by appending the first value of the mean and std
  872. mean_delays = np.append(mean_delays,mean_delays[0])
  873. std_delays = np.append(std_delays,std_delays[0])
  874. # Plot the shaded standard deviation area
  875. ax.fill(angles,mean_delays+std_delays,color=color[position],alpha=0.2) # Shaded area
  876. ax.fill(angles, mean_delays - std_delays, color=color[position], alpha=0.2) # Shaded area
  878. # Plot the mean delay
  879. ax.plot(angles, mean_delays, label=f'Position {position}', linewidth=2, marker='o')
  881. # Remove radial ticks
  882. ax.set_yticklabels([])
  884. # Set axis labels for the channels
  885. ax.set_xticks(angles[:-1]) # Set the angles for the channels
  886. ax.set_xticklabels(channels)
  888. # Set the title of the radar chart
  889. ax.set_title("Radar Chart: Average Delay Between ECG Peaks with Std Dev", size=16, color='black', y=1.1)
  891. # Add a legend
  892. ax.legend(title="Positions", loc='upper right')
  894. # Show the chart
  895. plt.show()
  898. # %% [markdown]
  899. # #### Machine learning classification of positions
  901. # %% [markdown]
  902. # Single electrode (middle - control electrode)
  904. # %%
  905. import pandas as pd
  906. from sklearn.model_selection import train_test_split
  907. from sklearn.linear_model import LogisticRegression
  908. from sklearn.model_selection import RepeatedStratifiedKFold, cross_val_score
  909. from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay,make_scorer,precision_score,recall_score,f1_score
  910. import matplotlib.pyplot as plt
  912. # Import SNR data from Excel
  913. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\ECG_classification\RR_interval_feature_matrix.xlsx'
  914. df = pd.read_excel(io=readpath)
  915. # Generate feature matrix from data
  916. fmatrix = df[['RRinterval_ch9', 'position']]
  917. # Separate features and labels
  918. X = fmatrix.drop('position', axis=1) # feature matrix
  919. y = fmatrix['position'] # labels
  920. # Split the data into training and testing sets
  921. X_train,X_test,y_train,y_test = train_test_split(X,y,test_size=0.2,stratify=fmatrix['position'])
  922. # Initialise and fit the logistic regression
  923. clf = LogisticRegression()
  924. clf.fit(X_train,y_train)
  925. # Predict the labels for the test set
  926. y_pred = clf.predict(X_test)
  927. # Calculate the confusion matrix
  928. cm = confusion_matrix(y_test,y_pred)
  929. # Display the confusion matrix
  930. disp = ConfusionMatrixDisplay(confusion_matrix=cm,display_labels=clf.classes_)
  931. fig,ax = plt.subplots(figsize=(6,6)) # Adjust figure size if needed
  932. # Set the color range limits
  933. vmin,vmax = 0,40 # Adjust these values as needed
  934. disp.plot(cmap="Blues",ax=ax,colorbar=False,values_format="d")
  935. disp.im_.set_clim(vmin,vmax) # Set the colorbar limits
  936. # Adjust the colorbar size
  937. cbar = ax.figure.colorbar(disp.im_,fraction=0.046,pad=0.04)
  938. cbar.set_label("Count") # Label for better readability
  939. plt.savefig(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Results\Confusion_matrices\Single_channel_ECG',dpi=1000)
  940. # Repeated Cross-Validation Accuracy
  941. rkf = RepeatedStratifiedKFold(n_splits=10,n_repeats=5)
  942. accuracy_scores = cross_val_score(clf,X,y,cv=rkf)
  943. precision_scores = cross_val_score(clf, X, y, cv=rkf, scoring=make_scorer(precision_score, average='macro'))
  944. recall_scores = cross_val_score(clf, X, y, cv=rkf, scoring=make_scorer(recall_score, average='macro'))
  945. f1_scores = cross_val_score(clf, X, y, cv=rkf, scoring=make_scorer(f1_score, average='macro'))
  947. # Print results
  948. print(f'Repeated Cross-Validation Accuracy: {accuracy_scores.mean():.2f} ± {accuracy_scores.std():.2f}')
  949. print(f"Repeated Cross-Validation Precision: {precision_scores.mean():.2f} ± {precision_scores.std():.2f}")
  950. print(f"Repeated Cross-Validation Recall: {recall_scores.mean():.2f} ± {recall_scores.std():.2f}")
  951. print(f"Repeated Cross-Validation F1 Score: {f1_scores.mean():.2f} ± {f1_scores.std():.2f}")
  953. # %% [markdown]
  954. # Electrode array
  956. # %%
  957. import pandas as pd
  958. from sklearn.model_selection import train_test_split
  959. from sklearn.linear_model import LogisticRegression
  960. from sklearn.model_selection import RepeatedStratifiedKFold, cross_val_score
  961. from sklearn.preprocessing import StandardScaler
  962. from sklearn.pipeline import Pipeline
  963. from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay,make_scorer,precision_score,recall_score,f1_score
  964. import matplotlib.pyplot as plt
  966. # Import SNR data from Excel
  967. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\ECG_classification\Delay_feature_matrix.xlsx'
  968. df = pd.read_excel(io=readpath)
  969. # Generate feature matrix from data
  970. fmatrix = df
  971. # Separate features and labels
  972. X = fmatrix.drop('position',axis=1) # feature matrix
  973. y = fmatrix['position'] # labels
  974. # Split the data into training and testing sets
  975. X_train,X_test,y_train,y_test = train_test_split(X,y,test_size=0.2,stratify=fmatrix['position'])
  976. # Define a pipeline with scaling and logistic regression
  977. pipeline = Pipeline([
  978. ('scaler',StandardScaler()), # Standardizes the data
  979. ('classifier',LogisticRegression()) # Logistic Regression model
  980. ])
  981. # Fit the model using only training data
  982. pipeline.fit(X_train,y_train)
  983. # Predict the labels for the test set
  984. y_pred = pipeline.predict(X_test)
  985. # Calculate the confusion matrix
  986. cm = confusion_matrix(y_test,y_pred)
  987. # Display the confusion matrix
  988. disp = ConfusionMatrixDisplay(confusion_matrix=cm,display_labels=pipeline.classes_)
  989. fig,ax = plt.subplots(figsize=(6,6)) # Adjust figure size if needed
  990. # Set the color range limits
  991. vmin,vmax = 0,40 # Adjust these values as needed
  992. disp.plot(cmap="Blues",ax=ax,colorbar=False,values_format="d")
  993. disp.im_.set_clim(vmin,vmax) # Set the colorbar limits
  994. # Adjust the colorbar size
  995. cbar = ax.figure.colorbar(disp.im_,fraction=0.046,pad=0.04)
  996. cbar.set_label("Count") # Label for better readability
  997. plt.savefig(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Results\Confusion_matrices\BSPM_ECG',dpi=1000)
  998. # Repeated Cross-Validation Accuracy
  999. rkf = RepeatedStratifiedKFold(n_splits=10,n_repeats=5)
  1000. accuracy_scores = cross_val_score(pipeline,X,y,cv=rkf)
  1001. precision_scores = cross_val_score(pipeline, X, y, cv=rkf, scoring=make_scorer(precision_score, average='macro'))
  1002. recall_scores = cross_val_score(pipeline, X, y, cv=rkf, scoring=make_scorer(recall_score, average='macro'))
  1003. f1_scores = cross_val_score(pipeline, X, y, cv=rkf, scoring=make_scorer(f1_score, average='macro'))
  1005. # Print results
  1006. print(f'Repeated Cross-Validation Accuracy: {accuracy_scores.mean():.2f} ± {accuracy_scores.std():.2f}')
  1007. print(f"Repeated Cross-Validation Precision: {precision_scores.mean():.2f} ± {precision_scores.std():.2f}")
  1008. print(f"Repeated Cross-Validation Recall: {recall_scores.mean():.2f} ± {recall_scores.std():.2f}")
  1009. print(f"Repeated Cross-Validation F1 Score: {f1_scores.mean():.2f} ± {f1_scores.std():.2f}")
  1012. # %% [markdown]
  1013. # #### Cardiac BSPM
  1015. # %%
  1016. import matplotlib.pyplot as plt
  1017. import numpy as np
  1019. from biolab.bspm import interpolate,BSPM
  1021. # Read .bspm data for sphere
  1022. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\ECG_classification\ECG_lying-side.bspm'
  1023. ECG = BSPM.load(loadpath=readpath)
  1024. rECG_sideways = ECG.resample(newfs=200)
  1025. rECG_sideways.preprocess(bw=(0.5,99))
  1026. # Display cardiac BSPM
  1027. vmap = rECG_sideways.potential(time=10.385)
  1028. frame = interpolate(valmap=np.where(np.isnan(-vmap),-vmap,-vmap/np.nanmax(-vmap)),method='cubic')
  1029. plt.imshow(frame,vmin=0.5,vmax=1)
  1030. cbar = plt.colorbar()
  1031. cbar.set_label('Voltage [$\mu$V]') # Setting the label for the colorbar
  1032. # Remove the axes
  1033. plt.axis('off')
  1034. plt.savefig(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Results\Sideways.pdf',dpi=1000)
  1035. plt.show()
  1037. # Read .bspm data for sphere
  1038. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\ECG_classification\ECG_lying-up.bspm'
  1039. ECG = BSPM.load(loadpath=readpath)
  1040. rECG_lying_down = ECG.resample(newfs=200)
  1041. rECG_lying_down.preprocess(bw=(0.5,99))
  1042. # Display cardiac BSPM
  1043. vmap = rECG_lying_down.potential(time=10.47)
  1044. frame = interpolate(valmap=np.where(np.isnan(-vmap),-vmap,-vmap/np.nanmax(-vmap)),method='cubic')
  1045. plt.imshow(frame,vmin=0.5,vmax=1)
  1046. cbar = plt.colorbar()
  1047. cbar.set_label('Voltage [$\mu$V]') # Setting the label for the colorbar
  1048. # Remove the axes
  1049. plt.axis('off')
  1050. plt.savefig(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Results\Lying down.pdf',dpi=1000)
  1051. plt.show()
  1053. # Read .bspm data for sphere
  1054. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\ECG_classification\ECG_sitting.bspm'
  1055. ECG = BSPM.load(loadpath=readpath)
  1056. rECG_sitting = ECG.resample(newfs=200)
  1057. rECG_sitting.preprocess(bw=(0.5,99))
  1058. # Display cardiac BSPM
  1059. vmap = rECG_sitting.potential(time=10.235)
  1060. frame = interpolate(valmap=np.where(np.isnan(-vmap),-vmap,-vmap/np.nanmax(-vmap)),method='cubic')
  1061. plt.imshow(frame,vmin=0.5,vmax=1)
  1062. cbar = plt.colorbar()
  1063. cbar.set_label('Voltage [$\mu$V]') # Setting the label for the colorbar
  1064. # Remove the axes
  1065. plt.axis('off')
  1066. plt.savefig(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Results\Sitting.pdf',dpi=1000)
  1067. plt.show()
  1069. # Read .bspm data for sphere
  1070. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\ECG_classification\ECG_standing.bspm'
  1071. ECG = BSPM.load(loadpath=readpath)
  1072. rECG_standing = ECG.resample(newfs=200)
  1073. rECG_standing.preprocess(bw=(0.5,99))
  1074. # Display cardiac BSPM
  1075. vmap = rECG_standing.potential(time=10.515)
  1076. frame = interpolate(valmap=np.where(np.isnan(-vmap),-vmap,-vmap/np.nanmax(-vmap)),method='cubic')
  1077. plt.imshow(frame,vmin=0.5,vmax=1)
  1078. cbar = plt.colorbar()
  1079. cbar.set_label('Voltage [$\mu$V]') # Setting the label for the colorbar
  1080. # Remove the axes
  1081. plt.axis('off')
  1082. plt.savefig(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Results\Standing.pdf',dpi=1000)
  1083. plt.show()
  1085. # %%
  1086. import pandas as pd
  1088. def delay(bspm,delays) -> np.ndarray:
  1090. # Create copy of mapped layout
  1091. dmap = np.copy(bspm.layout).astype(float)
  1092. dmap[dmap == 0] = np.NaN
  1093. # Convert channel_data to numpy arrays
  1094. locations = bspm.channel_data['location'].values
  1096. # Change value of label in layout to potential taking into account reference channels
  1097. # Create a mapping from location to voltage channel
  1098. location_to_voltage = dict(zip(locations,delays))
  1099. # Replace values in zmap based on location_to_impedance mapping
  1100. mask = np.isin(bspm.layout,locations)
  1101. dmap[mask] = np.vectorize(location_to_voltage.get)(bspm.layout[mask])
  1103. return dmap
  1105. data = pd.read_excel(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\ECG_classification\Delay_feature_matrix.xlsx')
  1106. data_avg = data.groupby('position',as_index=False).mean().drop(columns='position')
  1108. # Display cardiac delay maps
  1109. dmap = delay(rECG_sideways,np.insert(data_avg.iloc[0,:].values,13,0.0)*1000)
  1110. frame = interpolate(valmap=dmap,method='cubic')
  1111. plt.imshow(frame,vmin=0,vmax=10)
  1112. cbar = plt.colorbar()
  1113. cbar.set_label('Delay [ms]') # Setting the label for the colorbar
  1114. # Remove the axes
  1115. plt.axis('off')
  1116. plt.savefig(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Results\Sideways_delay.pdf',dpi=1000)
  1117. plt.show()
  1119. dmap = delay(rECG_lying_down,np.insert(data_avg.iloc[1,:].values,13,0.0)*1000)
  1120. frame = interpolate(valmap=dmap,method='cubic')
  1121. plt.imshow(frame,vmin=0,vmax=10)
  1122. cbar = plt.colorbar()
  1123. cbar.set_label('Delay [ms]') # Setting the label for the colorbar
  1124. # Remove the axes
  1125. plt.axis('off')
  1126. plt.savefig(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Results\Lying_down_delay.pdf',dpi=1000)
  1127. plt.show()
  1129. dmap = delay(rECG_sitting,np.insert(data_avg.iloc[2,:].values,13,0.0)*1000)
  1130. frame = interpolate(valmap=dmap,method='cubic')
  1131. plt.imshow(frame,vmin=0,vmax=10)
  1132. cbar = plt.colorbar()
  1133. cbar.set_label('Delay [ms]') # Setting the label for the colorbar
  1134. # Remove the axes
  1135. plt.axis('off')
  1136. plt.savefig(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Results\Sitting_delay.pdf',dpi=1000)
  1137. plt.show()
  1139. dmap = delay(rECG_standing,np.insert(data_avg.iloc[3,:].values,13,0.0)*1000)
  1140. frame = interpolate(valmap=dmap,method='cubic')
  1141. plt.imshow(frame,vmin=0,vmax=10)
  1142. cbar = plt.colorbar()
  1143. cbar.set_label('Delay [ms]') # Setting the label for the colorbar
  1144. # Remove the axes
  1145. plt.axis('off')
  1146. plt.savefig(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Results\Standing_delay.pdf',dpi=1000)
  1147. plt.show()
  1150. # %% [markdown]
  1151. # ### Sensorimotor classification from cerebral BSPM
  1153. # %% [markdown]
  1154. # #### Data import
  1156. # %%
  1157. import os
  1158. import re
  1160. from biolab.bspm import BSPM,image2layout
  1161. from biolab.utils import apply2all
  1163. # Save new recordings into .bspm files
  1164. # Define parameters
  1165. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Recordings\EEG_classification'
  1166. # Convert electrode array image into layout
  1167. ECG_lyt = image2layout(filename=r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\Masks\Layouts\EEG_lyt.jpg',seq='cw')
  1168. # Define operation to be applied to each file
  1169. def func(filepath):
  1170. # Extract filename from filepath
  1171. filename = os.path.splitext(os.path.basename(filepath))[0]
  1172. # Define the regex pattern to match the required values
  1173. pattern = r'^([A-Za-z_]+)_\d{6}_\d{6}$'
  1174. # Search for the pattern in the given filename
  1175. match = re.search(pattern,filename)
  1176. if match:
  1177. filename = match.group(1)
  1178. if filename == 'muscle activity':
  1179. mode = 'EEG/EMG'
  1180. else:
  1181. mode = 'EEG'
  1182. sensorimotor = filename
  1183. measurement = 'M1'
  1184. else:
  1185. raise ValueError("Filename \"{}\" does not match the expected pattern.".format(filename))
  1187. meta = {
  1188. 'mode': mode,
  1189. 'measurement': measurement,
  1190. 'sensorimotor': {'auditory_stimulus':'auditory','eyes_closed_open':'eye closing','mental_calculus':'mental calculus','muscle_activity':'motor','visual_stimulus':'visual'}.get(sensorimotor,None)
  1191. }
  1192. # Load BSPM data from .rhs file
  1193. MP = BSPM.from_file(filepath=filepath,layout=ECG_lyt,metadata=meta,refchs=[])
  1194. # Save BSPM object
  1195. MP.save(savepath=r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\EEG_classification\\'+filename+'.bspm')
  1197. # Call DataSaver function
  1198. apply2all(read_directory=readpath,extension='.rhs',func=func)
  1201. # %% [markdown]
  1202. # #### Data visualisation
  1204. # %% [markdown]
  1205. # Changes in alpha wave activity when opening and closing eyes
  1207. # %%
  1208. import pandas as pd
  1209. import plotly.express as px
  1210. from scipy.signal import butter,filtfilt
  1212. from biolab.bspm import BSPM
  1214. # Define read directory
  1215. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\EEG_classification\Testing\eyes_closed_open.bspm'
  1216. # Load .bspm data
  1217. EEG = BSPM.load(loadpath=readpath)
  1218. # Resample and preprocess the data
  1219. rEEG = EEG.resample(newfs=1000)
  1220. rEEG.preprocess(bw=(0.5,100))
  1221. # Define bandpass filter function
  1222. def bandpass_filter(data,lowcut,highcut,sf,order=4):
  1223. return filtfilt(*butter(order,[lowcut,highcut],fs=sf,btype='band'),data)
  1225. # Define EEG bands
  1226. bands = {'Delta (0.5-4 Hz)': (1,4), 'Theta (4-8 Hz)': (4,8),
  1227. 'Alpha (8-13 Hz)': (8,13), 'Beta (13-30 Hz)': (13,30),
  1228. 'Gamma (30-100 Hz)': (30,100)}
  1230. # Select one channel (e.g., 'A-000') and apply filters
  1231. channel = 'A-000'
  1232. filtered_signals = {band: bandpass_filter(rEEG.preprocessed_data[channel],low,high,1000) for band,(low,high) in bands.items()}
  1234. # Create DataFrame for Plotly
  1235. df_plot = pd.DataFrame({'Time (s)': rEEG.preprocessed_data.index.values})
  1236. for band, signal in filtered_signals.items():
  1237. df_plot[band] = signal
  1239. # Melt DataFrame for Plotly Express
  1240. df_melted = df_plot.melt(id_vars=['Time (s)'],var_name='Band',value_name='Amplitude')
  1242. # Plot with Plotly Express
  1243. fig = px.line(df_melted,x='Time (s)',y='Amplitude',color='Band',title=f'EEG Trace of {channel} Across Frequency Bands')
  1244. fig.show()
  1247. # %% [markdown]
  1248. # Changes in theta wave activity when performing mental calculus
  1250. # %%
  1251. import pandas as pd
  1252. import plotly.express as px
  1253. from scipy.signal import butter,filtfilt
  1255. from biolab.bspm import BSPM
  1257. # Define read directory
  1258. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\EEG_classification\Testing\mental_calculus.bspm'
  1259. # Load .bspm data
  1260. EEG = BSPM.load(loadpath=readpath)
  1261. # Resample and preprocess the data
  1262. rEEG = EEG.resample(newfs=1000)
  1263. rEEG.preprocess(bw=(0.5,100))
  1264. # Define bandpass filter function
  1265. def bandpass_filter(data,lowcut,highcut,sf,order=4):
  1266. return filtfilt(*butter(order,[lowcut,highcut],fs=sf,btype='band'),data)
  1268. # Define EEG bands
  1269. bands = {'Delta (0.5-4 Hz)': (1,4), 'Theta (4-8 Hz)': (4,8),
  1270. 'Alpha (8-13 Hz)': (8,13), 'Beta (13-30 Hz)': (13,30),
  1271. 'Gamma (30-100 Hz)': (30,100)}
  1273. # Select one channel (e.g., 'A-000') and apply filters
  1274. channel = 'A-000'
  1275. filtered_signals = {band: bandpass_filter(rEEG.preprocessed_data[channel],low,high,1000) for band,(low,high) in bands.items()}
  1277. # Create DataFrame for Plotly
  1278. df_plot = pd.DataFrame({'Time (s)': rEEG.preprocessed_data.index.values})
  1279. for band, signal in filtered_signals.items():
  1280. df_plot[band] = signal
  1282. # Melt DataFrame for Plotly Express
  1283. df_melted = df_plot.melt(id_vars=['Time (s)'],var_name='Band',value_name='Amplitude')
  1285. # Plot with Plotly Express
  1286. fig = px.line(df_melted,x='Time (s)',y='Amplitude',color='Band',title=f'EEG Trace of {channel} Across Frequency Bands')
  1287. fig.show()
  1290. # %% [markdown]
  1291. # #### Feature extraction
  1293. # %%
  1294. import numpy as np
  1295. import pandas as pd
  1296. import scipy.signal as signal
  1298. from biolab.utils import apply2all
  1299. from biolab.bspm import BSPM
  1301. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\EEG_classification'
  1303. def compute_band_power(eeg_signal,sf,band,method='welch'):
  1305. low,high = band
  1306. # Compute Power Spectral Density (PSD)
  1307. if method == 'welch':
  1308. freqs,psd = signal.welch(eeg_signal,sf,nperseg=sf*2) # 2-sec windows
  1309. else: # FFT method
  1310. freqs = np.fft.rfftfreq(len(eeg_signal),d=1/sf)
  1311. psd = np.abs(np.fft.rfft(eeg_signal))**2/len(eeg_signal)
  1313. # Integrate PSD within the frequency band
  1314. band_power = np.trapz(psd[(freqs>=low)&(freqs<=high)],freqs[(freqs>=low)&(freqs<=high)])
  1316. return band_power
  1318. # Define bands
  1319. bands = {'Delta':(0.5,4),'Theta':(4,8),'Alpha':(8,13),'Beta':(13,30),'Gamma':(30,100)}
  1321. def f(filepath):
  1322. EEG = BSPM.load(loadpath=filepath)
  1323. rEEG = EEG.resample(newfs=1000)
  1324. rEEG.preprocess(bw=(0.5,100))
  1326. sf = 1000 # Sampling frequency (Hz)
  1327. window_size = 2*sf # 2-second window
  1328. gap_size = 2*sf # 2-second window
  1329. start_idx = 10*sf # Start after 10 seconds
  1330. num_samples = len(rEEG.preprocessed_data)
  1331. # Create an empty list to store computed band power values
  1332. band_power_list = []
  1333. # Iterate over 2-second non-overlapping windows after the 10-second mark
  1334. for start in range(start_idx,num_samples,window_size+gap_size):
  1335. end = start+window_size
  1336. if end>num_samples: # Ensure the window does not exceed data length
  1337. break
  1339. # Extract windowed EEG data (2-sec segment)
  1340. segment = rEEG.preprocessed_data.iloc[start:end]
  1341. # Compute band power for each channel-band combination
  1342. band_power_dict = {}
  1343. for band_name,band_range in bands.items():
  1344. band_power_values = segment.apply(lambda x: compute_band_power(x.values,sf,band_range),axis=0)
  1345. # Rename columns to match channel + band (e.g., Ch1_Delta, Ch1_Theta)
  1346. for ch in band_power_values.index:
  1347. band_power_dict[f"{ch}_{band_name}"] = band_power_values[ch]
  1349. # Append result as a row in the list
  1350. band_power_list.append(band_power_dict)
  1352. # Convert list to DataFrame
  1353. df_band_power = pd.DataFrame(band_power_list)
  1354. df_band_power['action'] = EEG.metadata['sensorimotor']
  1355. excel_path = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\EEG_classification\EEG_feature_matrix.xlsx'
  1356. # Load the existing data from the Excel file
  1357. try:
  1358. existing_df = pd.read_excel(excel_path)
  1359. # Append the new data to the existing data
  1360. combined_df = pd.concat([existing_df,df_band_power],ignore_index=True)
  1361. except FileNotFoundError:
  1362. # If the file doesn't exist, just use the new data
  1363. combined_df = df_band_power
  1365. # Write the combined data back to the Excel file
  1366. with pd.ExcelWriter(excel_path,mode='w',engine='openpyxl') as writer:
  1367. combined_df.to_excel(writer,index=False)
  1370. apply2all(read_directory=readpath,extension='.bspm',exclude=['Testing'],func=f)
  1373. # %% [markdown]
  1374. # #### Feature values inspection
  1376. # %%
  1377. import matplotlib.pyplot as plt
  1378. import pandas as pd
  1379. import seaborn as sns
  1381. # Import EEG waves data
  1382. df = pd.read_excel(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\EEG_classification\EEG_feature_matrix.xlsx')
  1384. # Reshaping DataFrame: Convert from Wide to Long Format
  1385. df_long = df.melt(id_vars=["action"], var_name="Channel_Band", value_name="Power")
  1387. # Splitting "Channel_Band" into "Channel" and "Band"
  1388. df_long[['Channel', 'Band']] = df_long['Channel_Band'].str.extract(r"(A-\d{3})_(\w+)")
  1390. # Function to Remove Outliers using IQR
  1391. def remove_outliers(group):
  1392. Q1 = group["Power"].quantile(0.25)
  1393. Q3 = group["Power"].quantile(0.75)
  1394. IQR = Q3 - Q1
  1395. lower_bound = Q1 - 1.5 * IQR
  1396. upper_bound = Q3 + 1.5 * IQR
  1397. return group[(group["Power"] >= lower_bound) & (group["Power"] <= upper_bound)]
  1399. # Apply outlier removal
  1400. df_filtered = df_long.groupby(["Channel", "Band"]).apply(remove_outliers).reset_index(drop=True)
  1402. # Drop high-impedance channels
  1403. channels_to_remove = ["A-004", "A-025", "A-030", "A-031"]
  1404. df_filtered = df_filtered[~df_filtered["Channel"].isin(channels_to_remove)]
  1406. # Get unique EEG bands
  1407. bands = df_filtered["Band"].unique()
  1409. # Create subplots (one for each band)
  1410. fig, axes = plt.subplots(nrows=len(bands), figsize=(15, 3 * len(bands)), sharex=True)
  1412. for i, band in enumerate(bands):
  1413. ax = axes[i] if len(bands) > 1 else axes # Adjust for single-band case
  1414. band_data = df_filtered[df_filtered["Band"] == band] # Filter for current band
  1416. sns.violinplot(x="Channel", y="Power", hue="action", data=band_data, dodge=True,palette="muted", ax=ax)
  1418. ax.set_title(f"{band} Band Power Distribution Across Channels")
  1419. ax.set_ylabel("Power")
  1420. ax.set_xlabel("")
  1421. ax.tick_params(axis='x', rotation=90) # Rotate x-axis labels for readability
  1423. # Add x-axis label
  1424. plt.xlabel("Channel")
  1426. # Adjust layout & legend placement
  1427. plt.tight_layout()
  1428. plt.legend(title="Action",loc='upper left') # Move legend outside
  1430. # Show the plot
  1431. plt.show()
  1434. # %% [markdown]
  1435. # 1. Auditory Activity:
  1436. # - Theta (4-8 Hz): Involved in auditory processing, especially during tasks requiring attention or memory.
  1437. # - Gamma (30-100 Hz): Associated with auditory perception and integration, especially in response to complex sounds or language.
  1438. #
  1439. # 2. Visual Activity:
  1440. # - Alpha (8-12 Hz): Typically decreases (desynchronizes) in occipital regions during visual processing, indicating engagement.
  1441. # - Gamma (30-100 Hz): Related to visual perception, feature binding, and conscious awareness.
  1442. #
  1443. # 3. Motor Activity:
  1444. # - Beta (13-30 Hz): In sensorimotor cortex, beta activity often decreases (event-related desynchronization, ERD) during movement or motor preparation and increases (event-related synchronization, ERS) after movement.
  1445. # - Mu Rhythm (8-13 Hz): Similar to alpha but localized over motor cortex; suppressed during actual or imagined movement.
  1446. #
  1447. # Observation
  1448. # - Auditory vs. Visual: Focus on Alpha (visual engagement) and Gamma (auditory processing).
  1449. # - Auditory vs. Motor: Compare Theta/Gamma (auditory) with Beta/Mu (motor).
  1450. # - Visual vs. Motor: Contrast Alpha (visual) with Beta/Mu (motor).
  1452. # %% [markdown]
  1453. # #### Machine learning classification of brain activity
  1455. # %% [markdown]
  1456. # Single electrode (middle - control electrode)
  1458. # %%
  1459. import pandas as pd
  1460. from sklearn.model_selection import train_test_split
  1461. from sklearn.linear_model import LogisticRegression
  1462. from sklearn.model_selection import RepeatedStratifiedKFold, cross_val_score
  1463. from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay, make_scorer, precision_score, recall_score, f1_score
  1464. import matplotlib.pyplot as plt
  1466. # Import SNR data from Excel
  1467. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\EEG_classification\EEG_feature_matrix.xlsx'
  1468. df = pd.read_excel(io=readpath)
  1469. # Generate feature matrix from data
  1470. fmatrix = df[['A-003_Delta','A-003_Theta','A-003_Alpha','A-003_Beta','A-003_Gamma','action']]
  1471. # Separate features and labels
  1472. X = fmatrix.drop('action', axis=1) # feature matrix
  1473. y = fmatrix['action'] # labels
  1474. # Split the data into training and testing sets
  1475. X_train,X_test,y_train,y_test = train_test_split(X,y,test_size=0.2,stratify=fmatrix['action'])
  1476. # Initialise and fit the logistic regression
  1477. clf = LogisticRegression()
  1478. clf.fit(X_train,y_train)
  1479. # Predict the labels for the test set
  1480. y_pred = clf.predict(X_test)
  1481. # Calculate the confusion matrix
  1482. cm = confusion_matrix(y_test,y_pred)
  1483. # Display the confusion matrix
  1484. disp = ConfusionMatrixDisplay(confusion_matrix=cm,display_labels=clf.classes_)
  1485. fig,ax = plt.subplots(figsize=(6,6)) # Adjust figure size if needed
  1486. # Set the color range limits
  1487. vmin,vmax = 0,25 # Adjust these values as needed
  1488. disp.plot(cmap="Greens",ax=ax,colorbar=False,values_format="d")
  1489. disp.im_.set_clim(vmin,vmax) # Set the colorbar limits
  1490. # Adjust the colorbar size
  1491. cbar = ax.figure.colorbar(disp.im_,fraction=0.046,pad=0.04)
  1492. cbar.set_label("Count") # Label for better readability
  1493. plt.savefig(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Results\Confusion_matrices\Single_channel_EEG',dpi=1000)
  1494. # Repeated Cross-Validation Accuracy
  1495. rkf = RepeatedStratifiedKFold(n_splits=10,n_repeats=5)
  1496. accuracy_scores = cross_val_score(clf,X,y,cv=rkf)
  1497. precision_scores = cross_val_score(clf, X, y, cv=rkf, scoring=make_scorer(precision_score, average='macro'))
  1498. recall_scores = cross_val_score(clf, X, y, cv=rkf, scoring=make_scorer(recall_score, average='macro'))
  1499. f1_scores = cross_val_score(clf, X, y, cv=rkf, scoring=make_scorer(f1_score, average='macro'))
  1501. # Print results
  1502. print(f'Repeated Cross-Validation Accuracy: {accuracy_scores.mean():.2f} ± {accuracy_scores.std():.2f}')
  1503. print(f"Repeated Cross-Validation Precision: {precision_scores.mean():.2f} ± {precision_scores.std():.2f}")
  1504. print(f"Repeated Cross-Validation Recall: {recall_scores.mean():.2f} ± {recall_scores.std():.2f}")
  1505. print(f"Repeated Cross-Validation F1 Score: {f1_scores.mean():.2f} ± {f1_scores.std():.2f}")
  1508. # %% [markdown]
  1509. # Electrode array
  1511. # %%
  1512. import pandas as pd
  1513. from sklearn.model_selection import train_test_split
  1514. from sklearn.linear_model import LogisticRegression
  1515. from sklearn.model_selection import RepeatedStratifiedKFold, cross_val_score
  1516. from sklearn.preprocessing import StandardScaler
  1517. from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay, make_scorer, precision_score, recall_score, f1_score
  1518. import matplotlib.pyplot as plt
  1520. # Import SNR data from Excel
  1521. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Cortico-muscular axis\Experiments\INTAN\Data\EEG_classification\EEG_feature_matrix.xlsx'
  1522. df = pd.read_excel(io=readpath)
  1523. # Generate feature matrix from data
  1524. fmatrix = df
  1525. # Separate features and labels
  1526. X = fmatrix.drop('action',axis=1) # feature matrix
  1527. y = fmatrix['action'] # labels
  1528. # Split the data into training and testing sets
  1529. X_train,X_test,y_train,y_test = train_test_split(X,y,test_size=0.2,stratify=fmatrix['action'])
  1530. # Initialise and fit the logistic regression
  1531. clf = LogisticRegression()
  1532. clf.fit(X_train,y_train)
  1533. # Predict the labels for the test set
  1534. y_pred = clf.predict(X_test)
  1535. # Calculate the confusion matrix
  1536. cm = confusion_matrix(y_test,y_pred)
  1537. # Display the confusion matrix
  1538. disp = ConfusionMatrixDisplay(confusion_matrix=cm,display_labels=clf.classes_)
  1539. fig,ax = plt.subplots(figsize=(6,6)) # Adjust figure size if needed
  1540. # Set the color range limits
  1541. vmin,vmax = 0,25 # Adjust these values as needed
  1542. disp.plot(cmap="Greens",ax=ax,colorbar=False,values_format="d")
  1543. disp.im_.set_clim(vmin,vmax) # Set the colorbar limits
  1544. # Adjust the colorbar size
  1545. cbar = ax.figure.colorbar(disp.im_,fraction=0.046,pad=0.04)
  1546. cbar.set_label("Count") # Label for better readability
  1547. #plt.savefig(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Cortico-muscular axis\Experiments\INTAN\Results\Confusion_matrices\BSPM_EEG',dpi=1000)
  1548. # Repeated Cross-Validation Accuracy
  1549. rkf = RepeatedStratifiedKFold(n_splits=10,n_repeats=5)
  1550. accuracy_scores = cross_val_score(clf,X,y,cv=rkf)
  1551. precision_scores = cross_val_score(clf, X, y, cv=rkf, scoring=make_scorer(precision_score, average='macro'))
  1552. recall_scores = cross_val_score(clf, X, y, cv=rkf, scoring=make_scorer(recall_score, average='macro'))
  1553. f1_scores = cross_val_score(clf, X, y, cv=rkf, scoring=make_scorer(f1_score, average='macro'))
  1555. # Print results
  1556. print(f'Repeated Cross-Validation Accuracy: {accuracy_scores.mean():.2f} ± {accuracy_scores.std():.2f}')
  1557. print(f"Repeated Cross-Validation Precision: {precision_scores.mean():.2f} ± {precision_scores.std():.2f}")
  1558. print(f"Repeated Cross-Validation Recall: {recall_scores.mean():.2f} ± {recall_scores.std():.2f}")
  1559. print(f"Repeated Cross-Validation F1 Score: {f1_scores.mean():.2f} ± {f1_scores.std():.2f}")
  1562. # %% [markdown]
  1563. # Correlation circle for feature importance
  1565. # %%
  1566. import numpy as np
  1567. import matplotlib.pyplot as plt
  1568. from matplotlib.colors import to_rgba, Normalize
  1569. from matplotlib.cm import ScalarMappable
  1570. from sklearn.preprocessing import StandardScaler
  1571. from sklearn.decomposition import PCA
  1573. # ------------------------------------------------
  1574. # Extract frequency and channel info from feature names
  1575. # ------------------------------------------------
  1577. feature_names = X.columns
  1579. # Map canonical frequency bands to base colors
  1580. freq_colors = {
  1581. 'delta': 'blue',
  1582. 'theta': 'green',
  1583. 'alpha': 'red',
  1584. 'beta': 'orange',
  1585. 'gamma': 'purple'
  1586. }
  1588. frequencies = []
  1589. channels = []
  1591. for feat in feature_names:
  1592. channel_part, freq_part = feat.split('_')
  1593. channel_num = int(channel_part.split('-')[1])
  1594. frequencies.append(freq_part.lower())
  1595. channels.append(channel_num)
  1597. # Normalize channels for alpha mapping
  1598. chan_min = min(channels)
  1599. chan_max = max(channels)
  1600. channels_norm = [(c - chan_min)/(chan_max - chan_min) for c in channels] # 0 → 1
  1602. # ------------------------------------------------
  1603. # Standardize features and run PCA
  1604. # ------------------------------------------------
  1606. scaler = StandardScaler()
  1607. X_scaled = scaler.fit_transform(X)
  1609. pca = PCA(n_components=2)
  1610. X_pca = pca.fit_transform(X_scaled)
  1611. loadings = pca.components_.T * np.sqrt(pca.explained_variance_)
  1613. # ------------------------------------------------
  1614. # Plot correlation circle
  1615. # ------------------------------------------------
  1617. fig, ax = plt.subplots(figsize=(10,10), constrained_layout=True)
  1619. # Unit circle
  1620. circle = plt.Circle((0,0), 1, fill=False, linestyle='--', color='black')
  1621. ax.add_artist(circle)
  1623. # Plot arrows with frequency color + channel brightness
  1624. for i, feat in enumerate(feature_names):
  1625. base_color = freq_colors.get(frequencies[i], 'gray')
  1627. # Use alpha for channel: 0 = min channel → white overlay, 1 = max channel → base color
  1628. rgba = to_rgba(base_color, alpha=0.4 + 0.6*channels_norm[i])
  1630. ax.annotate(
  1631. "",
  1632. xy=(loadings[i,0], loadings[i,1]),
  1633. xytext=(0,0),
  1634. arrowprops=dict(
  1635. arrowstyle="->",
  1636. lw=2.5,
  1637. color=rgba
  1638. )
  1639. )
  1641. # ------------------------------------------------
  1642. # Legend for frequency bands
  1643. # ------------------------------------------------
  1644. for freq, color in freq_colors.items():
  1645. ax.plot([], [], color=color, lw=2, label=freq.capitalize())
  1646. ax.legend(title='Frequency band', loc='upper right')
  1648. # ------------------------------------------------
  1649. # Colorbar for channel brightness (white → black)
  1650. # ------------------------------------------------
  1651. # We'll create a grayscale colormap from white → black
  1652. from matplotlib.colors import LinearSegmentedColormap
  1654. grayscale = LinearSegmentedColormap.from_list('gray', ['white', 'black'])
  1655. norm = Normalize(vmin=chan_min, vmax=chan_max)
  1657. sm = ScalarMappable(cmap=grayscale, norm=norm)
  1658. sm.set_array([])
  1660. cbar = fig.colorbar(sm, ax=ax)
  1661. cbar.set_label('Channel number')
  1663. # ------------------------------------------------
  1664. # Axis formatting
  1665. # ------------------------------------------------
  1666. ax.axhline(0, color='black', linewidth=1.5)
  1667. ax.axvline(0, color='black', linewidth=1.5)
  1669. for spine in ax.spines.values():
  1670. spine.set_color('black')
  1671. spine.set_linewidth(1.5)
  1673. ax.tick_params(colors='black')
  1675. ax.set_xlim(-1.1,1.1)
  1676. ax.set_ylim(-1.1,1.1)
  1677. ax.set_aspect('equal')
  1679. ax.set_xlabel(f'PC1 ({pca.explained_variance_ratio_[0]*100:.1f}% variance)')
  1680. ax.set_ylabel(f'PC2 ({pca.explained_variance_ratio_[1]*100:.1f}% variance)')
  1681. ax.set_title('Correlation Circle (Cerebral BSPM)')
  1683. plt.show()
  1686. # %% [markdown]
  1687. # #### Cerebral BSPM
  1689. # %% [markdown]
  1690. # ##### BSPM data import
  1692. # %%
  1693. import os
  1694. import re
  1696. from biolab.bspm import BSPM,image2layout
  1697. from biolab.utils import apply2all
  1699. # Save new recordings into .bspm files
  1700. # Define parameters
  1701. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Recordings\EEG_classification\BSPM'
  1702. # Convert electrode array image into layout
  1703. ECG_lyt = image2layout(filename=r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\Masks\Layouts\EEG_lyt.jpg',seq='cw')
  1704. # Define operation to be applied to each file
  1705. def func(filepath):
  1706. # Extract filename from filepath
  1707. filename = os.path.splitext(os.path.basename(filepath))[0]
  1708. # Define the regex pattern to match the required values
  1709. pattern = r'^([A-Za-z_]+)_\d{6}_\d{6}$'
  1710. # Search for the pattern in the given filename
  1711. match = re.search(pattern,filename)
  1712. if match:
  1713. filename = match.group(1)
  1714. if filename == 'muscle activity':
  1715. mode = 'EEG/EMG'
  1716. else:
  1717. mode = 'EEG'
  1718. sensorimotor = filename
  1719. measurement = 'M1'
  1720. else:
  1721. raise ValueError("Filename \"{}\" does not match the expected pattern.".format(filename))
  1723. meta = {
  1724. 'mode': mode,
  1725. 'measurement': measurement,
  1726. 'sensorimotor': {'auditory':'auditory','motor':'motor','visual':'visual'}.get(sensorimotor,None)
  1727. }
  1728. # Load BSPM data from .rhs file
  1729. MP = BSPM.from_file(filepath=filepath,layout=ECG_lyt,metadata=meta,refchs=[])
  1730. # Save BSPM object
  1731. MP.save(savepath=r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\EEG_classification\BSPM\\'+filename+'.bspm')
  1733. # Call DataSaver function
  1734. apply2all(read_directory=readpath,extension='.rhs',func=func)
  1737. # %% [markdown]
  1738. # ##### Extract EEG features
  1740. # %%
  1741. import numpy as np
  1742. import pandas as pd
  1743. import scipy.signal as signal
  1745. from biolab.utils import apply2all
  1746. from biolab.bspm import BSPM
  1748. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\EEG_classification\BSPM'
  1750. # Function to normalize EEG signals (Z-score normalization)
  1751. def normalize_signals(eeg_data):
  1752. return (eeg_data - eeg_data.mean()) / eeg_data.std()
  1754. def compute_band_power(eeg_signal, sf, band, method='welch'):
  1755. """Compute power spectral density (PSD) or average magnitude for a given EEG signal."""
  1756. low, high = band
  1758. if method == 'welch':
  1759. freqs, psd = signal.welch(eeg_signal, sf, nperseg=sf*2) # 2-sec windows
  1760. band_power = np.trapz(psd[(freqs >= low) & (freqs <= high)], freqs[(freqs >= low) & (freqs <= high)])
  1761. return band_power
  1763. elif method == 'fft':
  1764. freqs = np.fft.rfftfreq(len(eeg_signal), d=1/sf)
  1765. psd = np.abs(np.fft.rfft(eeg_signal))**2 / len(eeg_signal)
  1766. band_power = np.trapz(psd[(freqs >= low) & (freqs <= high)], freqs[(freqs >= low) & (freqs <= high)])
  1767. return band_power
  1769. elif method == 'magnitude': # Compute average magnitude using RMS
  1770. filtered_signal = bandpass_filter(eeg_signal, low, high, sf)
  1771. return np.sqrt(np.mean(filtered_signal**2)) # RMS magnitude
  1773. else:
  1774. raise ValueError("Invalid method. Choose 'welch', 'fft', or 'magnitude'.")
  1776. def bandpass_filter(data, lowcut, highcut, sf, order=4):
  1777. """Apply a bandpass filter to the EEG signal."""
  1778. b, a = signal.butter(order, [lowcut, highcut], fs=sf, btype='band')
  1779. return signal.filtfilt(b, a, data)
  1781. # Define bands
  1782. bands = {'Theta': (4, 8), 'Alpha': (8, 13), 'Beta': (13, 30), 'Gamma': (30, 100)}
  1784. def f(filepath, method='power'): # Choose 'power' or 'magnitude'
  1785. EEG = BSPM.load(loadpath=filepath)
  1786. rEEG = EEG.resample(newfs=1000)
  1787. rEEG.preprocess(bw=(0.5, 100))
  1789. sf = 1000 # Sampling frequency (Hz)
  1791. # Define start & end times based on the file type
  1792. segment_map = {
  1793. 'auditory.bspm': (8.5,9.5),
  1794. 'motor.bspm': (10.2,11.2),
  1795. 'visual.bspm': (6.5,7.5)
  1796. }
  1798. for key, (start, end) in segment_map.items():
  1799. if key in filepath:
  1800. break
  1802. # Extract the EEG segment
  1803. segment = rEEG.preprocessed_data.iloc[int(start*sf):int(end*sf)]
  1805. # Normalize EEG data
  1806. segment_normalized = segment.apply(normalize_signals, axis=0)
  1808. # Compute feature (power or magnitude) for each band & channel
  1809. band_feature_list = []
  1810. band_feature_dict = {}
  1812. for band_name, band_range in bands.items():
  1813. feature_values = segment_normalized.apply(lambda x: compute_band_power(x.values, sf, band_range, method=method), axis=0)
  1815. for ch in feature_values.index:
  1816. band_feature_dict[f"{ch}_{band_name}"] = feature_values[ch]
  1818. band_feature_list.append(band_feature_dict)
  1820. # Convert list to DataFrame
  1821. df_band_feature = pd.DataFrame(band_feature_list)
  1822. df_band_feature['action'] = EEG.metadata['sensorimotor']
  1824. # Save to Excel
  1825. excel_path = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\EEG_classification\BSPM\EEG_bands_features.xlsx'
  1827. try:
  1828. existing_df = pd.read_excel(excel_path)
  1829. combined_df = pd.concat([existing_df, df_band_feature], ignore_index=True)
  1830. except FileNotFoundError:
  1831. combined_df = df_band_feature
  1833. with pd.ExcelWriter(excel_path, mode='w', engine='openpyxl') as writer:
  1834. combined_df.to_excel(writer, index=False)
  1836. # Apply to all files
  1837. apply2all(read_directory=readpath, extension='.bspm', exclude=['Testing'], func=lambda fp: f(fp, method='fft')) # Change method here
  1840. # %% [markdown]
  1841. # ##### Generate BSPM
  1843. # %%
  1844. import matplotlib.pyplot as plt
  1845. import numpy as np
  1847. from biolab.bspm import interpolate,BSPM
  1849. limit = 50
  1851. # Read .bspm data for sphere
  1852. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\EEG_classification\BSPM\auditory.bspm'
  1853. EEG = BSPM.load(loadpath=readpath)
  1854. rEEG_auditory = EEG.resample(newfs=1000)
  1855. rEEG_auditory.preprocess(bw=(0.5,100))
  1856. # Display cardiac BSPM
  1857. vmap = rEEG_auditory.potential(time=16.4)
  1858. frame = interpolate(valmap=vmap,method='cubic')
  1859. plt.imshow(frame,vmin=-limit,vmax=limit)
  1860. cbar = plt.colorbar()
  1861. cbar.set_label('Voltage [$\mu$V]') # Setting the label for the colorbar
  1862. # Remove the axes
  1863. plt.axis('off')
  1864. plt.savefig(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Results\Auditory.pdf',dpi=1000)
  1865. plt.show()
  1867. # Read .bspm data for sphere
  1868. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\EEG_classification\BSPM\motor.bspm'
  1869. EEG = BSPM.load(loadpath=readpath)
  1870. rEEG_motor = EEG.resample(newfs=1000)
  1871. rEEG_motor.preprocess(bw=(0.5,100))
  1872. # Display cardiac BSPM
  1873. vmap = rEEG_motor.potential(time=8.2)
  1874. frame = interpolate(valmap=vmap,method='cubic')
  1875. plt.imshow(frame,vmin=-limit,vmax=limit)
  1876. cbar = plt.colorbar()
  1877. cbar.set_label('Voltage [$\mu$V]') # Setting the label for the colorbar
  1878. # Remove the axes
  1879. plt.axis('off')
  1880. plt.savefig(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Results\Motor.pdf',dpi=1000)
  1881. plt.show()
  1883. # Read .bspm data for sphere
  1884. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\EEG_classification\BSPM\visual.bspm'
  1885. EEG = BSPM.load(loadpath=readpath)
  1886. rEEG_visual = EEG.resample(newfs=200)
  1887. rEEG_visual.preprocess(bw=(0.5,99))
  1888. # Display cardiac BSPM
  1889. vmap = rEEG_visual.potential(time=18.2)
  1890. frame = interpolate(valmap=vmap,method='cubic')
  1891. plt.imshow(frame,vmin=-limit,vmax=limit)
  1892. cbar = plt.colorbar()
  1893. cbar.set_label('Voltage [$\mu$V]') # Setting the label for the colorbar
  1894. # Remove the axes
  1895. plt.axis('off')
  1896. plt.savefig(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Results\Visual.pdf',dpi=1000)
  1897. plt.show()
  1899. # %%
  1900. import pandas as pd
  1901. import plotly.express as px
  1902. from scipy.signal import butter, filtfilt
  1904. from biolab.bspm import BSPM
  1906. # Define read directory
  1907. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\EEG_classification\BSPM\motor.bspm'
  1909. # Load .bspm data
  1910. EEG = BSPM.load(loadpath=readpath)
  1912. # Resample and preprocess the data
  1913. rEEG = EEG.resample(newfs=1000)
  1914. rEEG.preprocess(bw=(0.5, 100))
  1916. # Define bandpass filter function
  1917. def bandpass_filter(data, lowcut, highcut, sf, order=4):
  1918. return filtfilt(*butter(order, [lowcut, highcut], fs=sf, btype='band'), data)
  1920. # Define EEG bands
  1921. bands = {
  1922. #'Delta (0.5-4 Hz)': (0.5, 4),
  1923. 'Theta (4-8 Hz)': (4, 8),
  1924. 'Alpha (8-13 Hz)': (8, 13),
  1925. 'Beta (13-30 Hz)': (13, 30),
  1926. 'Gamma (30-100 Hz)': (30, 100)
  1927. }
  1929. # Initialize DataFrame for plotting
  1930. filtered_data = []
  1931. for channel in rEEG.preprocessed_data.columns:
  1932. for band, (low, high) in bands.items():
  1933. filtered_signal = bandpass_filter(rEEG.preprocessed_data[channel], low, high, 1000)
  1934. temp_df = pd.DataFrame({
  1935. 'Time (s)': rEEG.preprocessed_data.index.values,
  1936. 'Amplitude': filtered_signal,
  1937. 'Band': band,
  1938. 'Channel': channel
  1939. })
  1940. filtered_data.append(temp_df)
  1942. # Combine all data into a single DataFrame
  1943. df_melted = pd.concat(filtered_data, ignore_index=True)
  1945. # Plot with Plotly Express (all in one figure)
  1946. fig = px.line(
  1947. df_melted,
  1948. x='Time (s)',
  1949. y='Amplitude',
  1950. color='Band', # Different colors for bands
  1951. line_dash='Channel', # Different line styles for channels
  1952. title='EEG Traces Across Channels and Frequency Bands'
  1953. )
  1955. fig.show()
  1958. # %% [markdown]
  1959. # ##### Generate maps of anything but BSP
  1961. # %%
  1962. import pandas as pd
  1963. import numpy as np
  1964. import matplotlib.pyplot as plt
  1966. from biolab.bspm import BSPM,interpolate
  1968. def mapping(bspm,values) -> np.ndarray:
  1970. # Create copy of mapped layout
  1971. pmap = np.copy(bspm.layout).astype(float)
  1972. pmap[pmap == 0] = np.NaN
  1973. # Convert channel_data to numpy arrays
  1974. locations = bspm.channel_data['location'].values
  1976. # Change value of label in layout to potential taking into account reference channels
  1977. # Create a mapping from location to voltage channel
  1978. location_to_voltage = dict(zip(locations,values))
  1979. # Replace values in zmap based on location_to_impedance mapping
  1980. mask = np.isin(bspm.layout,locations)
  1981. pmap[mask] = np.vectorize(location_to_voltage.get)(bspm.layout[mask])
  1983. return pmap
  1985. # Load data
  1986. file_path = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\EEG_classification\BSPM\EEG_bands_features.xlsx'
  1987. data = pd.read_excel(file_path)
  1988. # Get the lowest value for each action
  1989. data_min = data.groupby('action',as_index=False).min()
  1990. # Select only Gamma band columns
  1991. band = 'Theta'
  1992. signal = 'Power_reg'
  1993. units = {'Power':'$\mu$V$^2$','RMS':'$\mu$V','Power_reg':'$\mu$V$^2$'}
  1994. unit = units[signal]
  1995. limits = {'Theta':500,'Alpha':200,'Beta':150,'Gamma':15} # RMS -> {'Theta':0.75,'Alpha':0.5,'Beta':0.5,'Gamma':0.1}
  1996. limit = limits[band]
  1997. data_min = data_min.loc[:,data_min.columns.str.contains(band,case=False,na=False)]
  1999. # Display power maps
  2000. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\EEG_classification\BSPM\auditory.bspm'
  2001. EEG = BSPM.load(loadpath=readpath)
  2002. rEEG_auditory = EEG.resample(newfs=1000)
  2003. pmap = mapping(rEEG_auditory,data_min.iloc[0,:].values)
  2004. frame = interpolate(valmap=pmap,method='cubic')
  2005. plt.imshow(frame,vmin=0,vmax=limit,cmap='viridis')
  2006. cbar = plt.colorbar()
  2007. cbar.set_label('{} [{}]'.format(signal,unit)) # Setting the label for the colorbar
  2008. # Remove the axes
  2009. plt.axis('off')
  2010. plt.savefig(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Results\\'+'{}_auditory_{}'.format(band,signal)+'_viridis.pdf',dpi=1000)
  2011. plt.show()
  2013. # Display power maps
  2014. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\EEG_classification\BSPM\motor.bspm'
  2015. EEG = BSPM.load(loadpath=readpath)
  2016. rEEG_motor = EEG.resample(newfs=1000)
  2017. pmap = mapping(rEEG_motor,data_min.iloc[1,:].values)
  2018. frame = interpolate(valmap=pmap,method='cubic')
  2019. plt.imshow(frame,vmin=0,vmax=limit,cmap='viridis')
  2020. cbar = plt.colorbar()
  2021. cbar.set_label('{} [{}]'.format(signal,unit)) # Setting the label for the colorbar
  2022. # Remove the axes
  2023. plt.axis('off')
  2024. plt.savefig(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Results\\'+'{}_motor_{}'.format(band,signal)+'_viridis.pdf',dpi=1000)
  2025. plt.show()
  2027. # Display power maps
  2028. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\EEG_classification\BSPM\visual.bspm'
  2029. EEG = BSPM.load(loadpath=readpath)
  2030. rEEG_visual = EEG.resample(newfs=1000)
  2031. pmap = mapping(rEEG_visual,data_min.iloc[2,:].values)
  2032. frame = interpolate(valmap=pmap,method='cubic')
  2033. plt.imshow(frame,vmin=0,vmax=limit,cmap='viridis')
  2034. cbar = plt.colorbar()
  2035. cbar.set_label('{} [{}]'.format(signal,unit)) # Setting the label for the colorbar
  2036. # Remove the axes
  2037. plt.axis('off')
  2038. plt.savefig(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Results\\'+'{}_visual_{}'.format(band,signal)+'_viridis.pdf',dpi=1000)
  2039. plt.show()
  2042. # %% [markdown]
  2043. # ## Multi-modal BSPM
  2045. # %% [markdown]
  2046. # ### Data import
  2048. # %%
  2049. import os
  2050. import re
  2052. from biolab.bspm import BSPM,image2layout
  2053. from biolab.utils import apply2all
  2055. # Save new recordings into .bspm files
  2056. # Define parameters
  2057. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Recordings\Multi-modal_recording'
  2058. # Convert electrode array image into layout
  2059. EMG_EEG_lyt = image2layout(image_path=r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\Masks\Layouts\EEG-EMG_lyt.jpg',seq='cw')
  2060. # Define operation to be applied to each file
  2061. def func(filepath):
  2062. # Extract filename from filepath
  2063. filename = os.path.splitext(os.path.basename(filepath))[0]
  2064. # Define the regex pattern to match the required values
  2065. pattern = r'^([A-Za-z_]+)_\d{6}_\d{6}$'
  2066. # Search for the pattern in the given filename
  2067. match = re.search(pattern,filename)
  2068. if match:
  2069. filename = match.group(1)
  2070. mode = 'EMG/ECG'
  2071. movement = filename
  2072. measurement = 'M1'
  2073. else:
  2074. raise ValueError("Filename \"{}\" does not match the expected pattern.".format(filename))
  2076. meta = {
  2077. 'mode': mode,
  2078. 'measurement': measurement,
  2079. 'movement': movement
  2080. }
  2082. savepath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\Multi-modal\\'+filename+'.bspm'
  2083. if not os.path.isfile(savepath):
  2084. # Load BSPM data from .rhs file
  2085. MP = BSPM.from_file(filepath=filepath,layout=EMG_EEG_lyt,metadata=meta,refchs=[])
  2086. # Save BSPM object
  2087. MP.save(savepath=savepath)
  2089. # Call DataSaver function
  2090. apply2all(read_directory=readpath,extension='.rhs',exclude=['Individual_recordings'],func=func)
  2093. # %% [markdown]
  2094. # ### Joint cerebral and muscular BSPM
  2096. # %% [markdown]
  2097. # Import multimodal data
  2099. # %%
  2100. from biolab.bspm import BSPM
  2102. # Read .bspm data for fist
  2103. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Cortico-muscular axis\Experiments\INTAN\Data\Multi-modal\fist.bspm'
  2104. FIST = BSPM.load(loadpath=readpath)
  2105. rEMGFIST = FIST.resample(newfs=1000)
  2106. rEEGFIST = FIST.resample(newfs=1000)
  2107. rEMGFIST.preprocess(bw=(5,400))
  2108. rEEGFIST.preprocess(bw=(0.5,100))
  2110. # Read .bspm data for up
  2111. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Cortico-muscular axis\Experiments\INTAN\Data\Multi-modal\up.bspm'
  2112. UP = BSPM.load(loadpath=readpath)
  2113. rEMGUP = UP.resample(newfs=1000)
  2114. rEEGUP = UP.resample(newfs=1000)
  2115. rEMGUP.preprocess(bw=(5,400))
  2116. rEEGUP.preprocess(bw=(0.5,100))
  2118. # %% [markdown]
  2119. # Display EEG power bands
  2121. # %%
  2122. import pandas as pd
  2123. import plotly.express as px
  2124. from scipy.signal import butter, filtfilt
  2127. # Define data to be displayed
  2128. data = rEEGUP.preprocessed_data
  2130. # Define bandpass filter function
  2131. def bandpass_filter(data, lowcut, highcut, sf, order=4):
  2132. return filtfilt(*butter(order, [lowcut, highcut], fs=sf, btype='band'), data)
  2134. # Define EEG bands
  2135. bands = {
  2136. #'Delta (0.5-4 Hz)': (0.5, 4),
  2137. 'Theta (4-8 Hz)': (4, 8),
  2138. 'Alpha (8-13 Hz)': (8, 13),
  2139. 'Beta (13-30 Hz)': (13, 30),
  2140. 'Gamma (30-100 Hz)': (30, 100)
  2141. }
  2143. # Initialize DataFrame for plotting
  2144. filtered_data = []
  2145. for channel in data.filter(like="A").columns:
  2146. for band, (low, high) in bands.items():
  2147. filtered_signal = bandpass_filter(data[channel], low, high, 1000)
  2148. temp_df = pd.DataFrame({
  2149. 'Time (s)': data.index.values,
  2150. 'Amplitude': filtered_signal,
  2151. 'Band': band,
  2152. 'Channel': channel
  2153. })
  2154. filtered_data.append(temp_df)
  2156. # Combine all data into a single DataFrame
  2157. df_melted = pd.concat(filtered_data, ignore_index=True)
  2159. # Plot with Plotly Express (all in one figure)
  2160. fig = px.line(
  2161. df_melted,
  2162. x='Time (s)',
  2163. y='Amplitude',
  2164. color='Band', # Different colors for bands
  2165. line_dash='Channel', # Different line styles for channels
  2166. title='EEG Traces Across Channels and Frequency Bands'
  2167. )
  2169. fig.show()
  2171. # %% [markdown]
  2172. # Display EMG envelopes
  2174. # %%
  2175. import plotly.express as px
  2177. data = rEMGFIST.preprocessed_data
  2179. # Ensure the DataFrame is properly referenced
  2180. fig = px.line(
  2181. data.abs(),#.rolling(window=500).mean(),
  2182. x=data.index,
  2183. y=data.columns # Specify the column(s) to plot
  2184. )
  2186. fig.show()
  2188. # %% [markdown]
  2189. # Generate EMG BSPM
  2191. # %%
  2192. import copy
  2194. from biolab.bspm import interpolate,layout2map
  2196. data = copy.deepcopy(rEMGFIST.preprocessed_data.head(30000).abs().rolling(window=500).mean())
  2197. chdata = copy.deepcopy(rEMGFIST.channel_data)
  2198. # Display cardiac BSPM
  2199. chdata['custom_channel_name'] = ['A','A','A','A','A','A','A','A','A','A','A','A','A','A','A','A','B','B','B','B','B','B','B','B','B','B','B','B','B','B','B','B']+chdata['custom_channel_name']
  2200. # Step 1: Create a mapping from 'locations' to 'custom_channel_names'
  2201. location_to_channel = dict(zip(chdata['location'],chdata['custom_channel_name']))
  2202. # Step 2: Rearrange the recdata columns based on the locations order
  2203. # Match the locations to the recdata columns
  2204. reordered_columns = [location_to_channel[loc] for loc in chdata['location'] if location_to_channel[loc] in rEMGFIST.preprocessed_data.columns]
  2205. # Step 3: Get the values at the specific index from recdata with the correct column order
  2206. values_at_index = data.loc[15.045,reordered_columns] # FIST = PRE (14) PERI (15.045) # UP = PRE (9.68) PERI (10.218)
  2207. vmap = layout2map(layout=rEMGFIST.layout,mapping=dict(zip(chdata['location'].values,values_at_index)))
  2208. frame = interpolate(valmap=vmap,method='cubic')
  2211. # %% [markdown]
  2212. # Generate EEG BSPM
  2214. # %%
  2215. import copy
  2217. from biolab.bspm import interpolate,layout2map
  2219. import numpy as np
  2220. import pandas as pd
  2221. import scipy.signal as signal
  2224. # Data to be used
  2225. data = copy.deepcopy(rEMGUP.preprocessed_data.head(30000))
  2226. method = 'fft'
  2228. # Function to normalize EEG signals (Z-score normalization)
  2229. def normalize_signals(eeg_data):
  2230. return (eeg_data - eeg_data.mean()) / eeg_data.std()
  2232. def compute_band_power(eeg_signal, sf, band, method='welch'):
  2233. """Compute power spectral density (PSD) or average magnitude for a given EEG signal."""
  2234. low, high = band
  2236. if method == 'welch':
  2237. freqs, psd = signal.welch(eeg_signal, sf, nperseg=sf*2) # 2-sec windows
  2238. band_power = np.trapz(psd[(freqs >= low) & (freqs <= high)], freqs[(freqs >= low) & (freqs <= high)])
  2239. return band_power
  2241. elif method == 'fft':
  2242. freqs = np.fft.rfftfreq(len(eeg_signal), d=1/sf)
  2243. psd = np.abs(np.fft.rfft(eeg_signal))**2 / len(eeg_signal)
  2244. band_power = np.trapz(psd[(freqs >= low) & (freqs <= high)], freqs[(freqs >= low) & (freqs <= high)])
  2245. return band_power
  2247. elif method == 'magnitude': # Compute average magnitude using RMS
  2248. filtered_signal = bandpass_filter(eeg_signal, low, high, sf)
  2249. return np.sqrt(np.mean(filtered_signal**2)) # RMS magnitude
  2251. else:
  2252. raise ValueError("Invalid method. Choose 'welch', 'fft', or 'magnitude'.")
  2254. def bandpass_filter(data, lowcut, highcut, sf, order=4):
  2255. """Apply a bandpass filter to the EEG signal."""
  2256. b, a = signal.butter(order, [lowcut, highcut], fs=sf, btype='band')
  2257. return signal.filtfilt(b, a, data)
  2259. # Define bands
  2260. bands = {'Theta': (4, 8), 'Alpha': (8, 13), 'Beta': (13, 30), 'Gamma': (30, 100)}
  2262. sf = 1000 # Sampling frequency (Hz)
  2264. segment_map = (14.4,15.3) # FIST = PRE (12.2,14.2) PERI (14.2,16.2) # UP = PRE (13.5,14.4) PERI (14.4,15.3)
  2266. # Extract the EEG segment
  2267. segment = data.iloc[int(segment_map[0]*sf):int(segment_map[1]*sf)]
  2269. # Normalize EEG data
  2270. segment_normalized = segment.apply(normalize_signals, axis=0)
  2272. # Compute feature (power or magnitude) for each band & channel
  2273. band_feature_list = []
  2274. band_feature_dict = {}
  2276. for band_name, band_range in bands.items():
  2277. feature_values = segment_normalized.apply(lambda x: compute_band_power(x.values, sf, band_range, method=method), axis=0)
  2279. for ch in feature_values.index:
  2280. band_feature_dict[f"{ch}_{band_name}"] = feature_values[ch]
  2282. band_feature_list.append(band_feature_dict)
  2284. # Convert list to DataFrame
  2285. df_band_feature = pd.DataFrame(band_feature_list)
  2288. # %%
  2289. chdata = copy.deepcopy(rEEGUP.channel_data)
  2290. wave = 'Beta'
  2291. # Display BSPM
  2292. chdata['custom_channel_name'] = ['A','A','A','A','A','A','A','A','A','A','A','A','A','A','A','A','B','B','B','B','B','B','B','B','B','B','B','B','B','B','B','B']+chdata['custom_channel_name']
  2293. # Step 1: Create a mapping from 'locations' to 'custom_channel_names'
  2294. location_to_channel = dict(zip(chdata['location'],chdata['custom_channel_name']+'_{}'.format(wave)))
  2295. # Step 2: Rearrange the recdata columns based on the locations order
  2296. # Match the locations to the recdata columns
  2297. reordered_columns = [location_to_channel[loc] for loc in chdata['location'] if location_to_channel[loc] in df_band_feature.filter(like=wave)]
  2298. # Step 3: Get the values at the specific index from recdata with the correct column order
  2299. values_at_index = df_band_feature.loc[0,reordered_columns]
  2300. vmap = layout2map(layout=rEEGUP.layout,mapping=dict(zip(chdata['location'].values,values_at_index)))
  2301. frame = interpolate(valmap=vmap,method='cubic')
  2303. # %%
  2304. import matplotlib.pyplot as plt
  2306. limit = 100
  2308. plt.imshow(frame,vmin=0,vmax=limit,cmap='viridis')
  2309. cbar = plt.colorbar()
  2310. cbar.set_label('Voltage [$\mu$V]') # Setting the label for the colorbar
  2311. # Remove the axes
  2312. plt.axis('off')
  2313. plt.savefig(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Results\Figure_4\Multimodal_EMG_peri_fist_viridis.pdf',dpi=1000)
  2314. plt.show()
  2316. # %% [markdown]
  2317. # ### Reaction time
  2319. # %% [markdown]
  2320. # Import multimodal data
  2322. # %%
  2323. from biolab.bspm import BSPM
  2325. # Read .bspm data for fist
  2326. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\Multi-modal\fist.bspm'
  2327. FIST = BSPM.load(loadpath=readpath)
  2328. rEMGFIST = FIST.resample(newfs=1000)
  2329. rEEGFIST = FIST.resample(newfs=1000)
  2330. rEMGFIST.preprocess(bw=(5,400))
  2331. rEEGFIST.preprocess(bw=(0.5,100))
  2333. # Read .bspm data for up
  2334. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\Multi-modal\up.bspm'
  2335. UP = BSPM.load(loadpath=readpath)
  2336. rEMGUP = UP.resample(newfs=1000)
  2337. rEEGUP = UP.resample(newfs=1000)
  2338. rEMGUP.preprocess(bw=(5,400))
  2339. rEEGUP.preprocess(bw=(0.5,100))
  2341. # %% [markdown]
  2342. # Display EEG power bands
  2344. # %%
  2345. import pandas as pd
  2346. import plotly.express as px
  2347. from scipy.signal import butter, filtfilt
  2348. import copy
  2350. # Define data to be displayed
  2351. data = copy.deepcopy(rEEGUP.preprocessed_data.head(30000))
  2353. # Define bandpass filter function
  2354. def bandpass_filter(data, lowcut, highcut, sf, order=4):
  2355. return filtfilt(*butter(order, [lowcut, highcut], fs=sf, btype='band'), data)
  2357. # Define EEG bands
  2358. bands = {
  2359. #'Delta (0.5-4 Hz)': (0.5, 4),
  2360. 'Theta (4-8 Hz)': (4, 8),
  2361. 'Alpha (8-13 Hz)': (8, 13),
  2362. 'Beta (13-30 Hz)': (13, 30),
  2363. 'Gamma (30-100 Hz)': (30, 100)
  2364. }
  2366. # Initialize DataFrame for plotting
  2367. filtered_data = []
  2368. for channel in data.filter(like="A").columns:
  2369. for band, (low, high) in bands.items():
  2370. filtered_signal = bandpass_filter(data[channel], low, high, 1000)
  2371. temp_df = pd.DataFrame({
  2372. 'Time (s)': np.lib.stride_tricks.sliding_window_view(np.abs(data.index.values),500).mean(axis=1),
  2373. 'Amplitude': np.lib.stride_tricks.sliding_window_view(np.abs(filtered_signal),500).mean(axis=1),
  2374. 'Band': band,
  2375. 'Channel': channel
  2376. })
  2377. filtered_data.append(temp_df)
  2379. # Combine all data into a single DataFrame
  2380. df_melted = pd.concat(filtered_data, ignore_index=True)
  2382. # Plot with Plotly Express (all in one figure)
  2383. fig = px.line(
  2384. df_melted,
  2385. x='Time (s)',
  2386. y='Amplitude',
  2387. color='Band', # Different colors for bands
  2388. line_dash='Channel', # Different line styles for channels
  2389. title='EEG Traces Across Channels and Frequency Bands'
  2390. )
  2392. fig.show()
  2394. # %% [markdown]
  2395. # Display EMG envelopes
  2397. # %%
  2398. import plotly.express as px
  2400. data = rEMGUP.preprocessed_data.head(30000)
  2402. # Ensure the DataFrame is properly referenced
  2403. fig = px.line(
  2404. data.abs().rolling(window=500).mean(), #.sub(data.iloc[0],axis=1),
  2405. x=data.index,
  2406. y=data.columns[data.columns.str.contains("B")] # Specify the column(s) to plot
  2407. )
  2409. fig.show()
  2411. # %% [markdown]
  2412. # Compute peaks for EMG data for display
  2414. # %%
  2415. import plotly.graph_objects as go
  2416. from scipy.signal import find_peaks
  2417. import copy
  2419. data = copy.deepcopy(rEMGUP.preprocessed_data.head(90000).abs().rolling(window=500).mean().dropna())
  2420. df_emg = data[[col for col in data.columns if "B" in col]]
  2422. # Function to find peaks
  2423. def get_peaks(df):
  2424. peaks_dict = {}
  2425. for col in df.columns:
  2426. peaks, _ = find_peaks(df[col],height=20, distance=1000) # Adjust height threshold as needed
  2427. peaks_dict[col] = peaks
  2428. return peaks_dict
  2430. # Compute the peaks
  2431. peaks_dict_emg = get_peaks(df_emg)
  2433. # Plot with Plotly
  2434. fig = go.Figure()
  2436. # Add the channel signals
  2437. for col in df_emg.columns:
  2438. fig.add_trace(go.Scatter(
  2439. x=df_emg.index,
  2440. y=df_emg[col],
  2441. mode='lines',
  2442. name=col
  2443. ))
  2445. # Add the peaks as scatter points
  2446. for col, peaks in peaks_dict_emg.items():
  2447. fig.add_trace(go.Scatter(
  2448. x=df_emg.index[peaks], # Use time for x-axis
  2449. y=df_emg[col].iloc[peaks],
  2450. mode='markers',
  2451. name=f"{col} Peaks",
  2452. marker=dict(size=8, symbol='circle', color='red')
  2453. ))
  2455. # Layout adjustments
  2456. fig.update_layout(
  2457. title="Interactive Peaks Plot with Time Axis",
  2458. xaxis_title="Time (seconds)",
  2459. yaxis_title="Amplitude",
  2460. legend=dict(title="Channels & Peaks"),
  2461. height=600,
  2462. width=1000
  2463. )
  2465. # Show the interactive plot
  2466. fig.show()
  2469. # %% [markdown]
  2470. # Compute peaks for EEG data (theta wave, motor planning) for display
  2472. # %%
  2473. import pandas as pd
  2474. import numpy as np
  2475. import plotly.graph_objects as go
  2476. from scipy.signal import butter, filtfilt, find_peaks
  2477. import copy
  2479. # Define bandpass filter function
  2480. def bandpass_filter(data, lowcut, highcut, sf, order=4):
  2481. """Applies a bandpass filter to the data."""
  2482. return filtfilt(*butter(order, [lowcut, highcut], fs=sf, btype='band'), data)
  2484. # Sampling frequency
  2485. sf = 1000
  2487. # Define EEG bands
  2488. bands = {
  2489. 'Theta (4-8 Hz)': (4, 8),
  2490. 'Alpha (8-13 Hz)': (8, 13),
  2491. 'Beta (13-30 Hz)': (13, 30),
  2492. 'Gamma (30-100 Hz)': (30, 100)
  2493. }
  2495. # Load and preprocess the data
  2496. data = copy.deepcopy(rEEGFIST.preprocessed_data.head(90000))
  2498. # Filter for theta waves across all channels
  2499. filtered_data = []
  2500. for channel in data.filter(like="A").columns:
  2501. filtered_signal = bandpass_filter(data[channel], bands['Theta (4-8 Hz)'][0], bands['Theta (4-8 Hz)'][1], sf)
  2503. # Apply rolling amplitude extraction
  2504. amplitude = np.lib.stride_tricks.sliding_window_view(np.abs(filtered_signal), 500).mean(axis=1)
  2506. # Align time index with the amplitude array
  2507. time_aligned = data.index.values[499:]
  2509. # Store filtered data
  2510. temp_df = pd.DataFrame({
  2511. 'Time (s)': time_aligned,
  2512. 'Amplitude': amplitude,
  2513. 'Channel': channel
  2514. })
  2515. filtered_data.append(temp_df)
  2517. # Combine all theta-filtered data into a single DataFrame
  2518. df_eeg = pd.concat(filtered_data, ignore_index=True)
  2520. # Function to detect peaks
  2521. def get_peaks(df, height=0.5, distance=1000):
  2522. """Detect peaks across all channels."""
  2523. peaks_dict = {}
  2524. for channel in df['Channel'].unique():
  2525. channel_data = df[df['Channel'] == channel]
  2526. peaks, _ = find_peaks(channel_data['Amplitude'], height=height, distance=distance)
  2527. peaks_dict[channel] = peaks
  2528. return peaks_dict
  2530. # Detect peaks
  2531. peaks_dict_eeg = get_peaks(df_eeg, height=5, distance=1000) # Adjust height and distance as needed
  2533. # Create interactive plot
  2534. fig = go.Figure()
  2536. # Plot the theta amplitude traces
  2537. for channel in df_eeg['Channel'].unique():
  2538. channel_data = df_eeg[df_eeg['Channel'] == channel]
  2540. fig.add_trace(go.Scatter(
  2541. x=channel_data['Time (s)'],
  2542. y=channel_data['Amplitude'],
  2543. mode='lines',
  2544. name=f"{channel} Amplitude"
  2545. ))
  2547. # Plot the peaks as scatter points
  2548. for channel, peaks in peaks_dict_eeg.items():
  2549. channel_data = df_eeg[df_eeg['Channel'] == channel]
  2551. fig.add_trace(go.Scatter(
  2552. x=channel_data['Time (s)'].iloc[peaks], # Use correct time index
  2553. y=channel_data['Amplitude'].iloc[peaks],
  2554. mode='markers',
  2555. name=f"{channel} Peaks",
  2556. marker=dict(size=8, symbol='circle', color='red')
  2557. ))
  2559. # Layout adjustments
  2560. fig.update_layout(
  2561. title="Theta Wave Peaks Across Channels",
  2562. xaxis_title="Time (seconds)",
  2563. yaxis_title="Amplitude",
  2564. legend=dict(title="Channels & Peaks"),
  2565. height=600,
  2566. width=1000
  2567. )
  2569. # Show the interactive plot
  2570. fig.show()
  2574. # %% [markdown]
  2575. # Apply to all data without visualising
  2577. # %%
  2578. import pandas as pd
  2579. import numpy as np
  2580. import copy
  2581. from scipy.signal import find_peaks, butter, filtfilt
  2583. DATA = (rEMGUP,rEEGUP)
  2585. # Load and preprocess the EMG data
  2586. data = copy.deepcopy(DATA[0].preprocessed_data.abs().rolling(window=900).mean().dropna())
  2587. df_emg = data[[col for col in data.columns if "B" in col]]
  2589. # Define bandpass filter function
  2590. def bandpass_filter(data, lowcut, highcut, sf, order=4):
  2591. """Applies a bandpass filter to the data."""
  2592. return filtfilt(*butter(order, [lowcut, highcut], fs=sf, btype='band'), data)
  2594. # Function to find peaks in EMG
  2595. def get_peaks(df):
  2596. peaks_dict = {}
  2597. for col in df.columns:
  2598. peaks, _ = find_peaks(df[col], height=20, distance=1000) # Adjust height threshold as needed
  2599. peaks_dict[col] = peaks
  2600. return peaks_dict
  2602. # Get peaks for EMG channels
  2603. peaks_dict_emg = get_peaks(df_emg)
  2605. # Load and preprocess the EEG data (filtered for theta band)
  2606. data = copy.deepcopy(DATA[1].preprocessed_data)
  2607. sf = 1000 # Sampling frequency
  2608. filtered_data = []
  2610. # Filter for theta waves across all channels (only A4 in EEG for delay calculation)
  2611. for channel in data.filter(like="A").columns:
  2612. filtered_signal = bandpass_filter(data[channel], 4, 8, sf) # Filtering for theta band (4-8 Hz)
  2613. amplitude = np.lib.stride_tricks.sliding_window_view(np.abs(filtered_signal), 500).mean(axis=1)
  2614. time_aligned = data.index.values[499:] # Align time with amplitude
  2615. temp_df = pd.DataFrame({'Time (s)': time_aligned, 'Amplitude': amplitude, 'Channel': channel})
  2616. filtered_data.append(temp_df)
  2618. df_eeg = pd.concat(filtered_data, ignore_index=True)
  2620. # Function to detect peaks in EEG
  2621. def get_peaks_eeg(df, height=0.5, distance=1000):
  2622. peaks_dict = {}
  2623. for channel in df['Channel'].unique():
  2624. channel_data = df[df['Channel'] == channel]
  2625. peaks, _ = find_peaks(channel_data['Amplitude'], height=height, distance=distance)
  2626. peaks_dict[channel] = peaks
  2627. return peaks_dict
  2629. # Get peaks for EEG (across all channels)
  2630. peaks_dict_eeg = get_peaks_eeg(df_eeg, height=5, distance=1000)
  2632. # Extract the time values for peaks from each EEG channel
  2633. eeg_peaks_times = {}
  2634. for channel in df_eeg['Channel'].unique():
  2635. eeg_peaks_times[channel] = df_eeg[df_eeg['Channel'] == channel]['Time (s)'].iloc[peaks_dict_eeg[channel]].values
  2637. # Extract the time values for peaks in channel A4 specifically
  2638. a4_peaks_times = eeg_peaks_times.get('A4', [])
  2640. # Function to calculate delay between EEG and EMG peaks with tolerance check and no reuse of EEG peaks
  2641. def calculate_delays_with_tolerance(a4_peaks_times, emg_peaks_dict, df_emg, eeg_peaks_times, min_tolerance=0.1, max_tolerance=0.6):
  2642. delays = {}
  2643. # Initialize a dictionary to track the matched peaks for each EEG channel
  2644. matched_eeg_peaks = {channel: set() for channel in eeg_peaks_times}
  2646. # Calculate the maximum number of peaks across all EMG channels
  2647. max_emg_peaks = max(len(peaks) for peaks in emg_peaks_dict.values())
  2649. for emg_channel, emg_peaks in emg_peaks_dict.items():
  2650. emg_peaks_times = df_emg[emg_channel].index[emg_peaks].values # Get time indices for EMG peaks
  2651. emg_delays = []
  2653. for emg_time in emg_peaks_times:
  2654. # For each EMG peak, try to match it with the EEG peaks of each channel
  2655. matched_delay = False
  2656. for eeg_channel, a4_peaks_times in eeg_peaks_times.items():
  2657. # Find the closest EEG peak in the current EEG channel within the tolerance range, excluding already matched peaks
  2658. time_diffs = emg_time - a4_peaks_times # Time differences
  2659. valid_peaks_times = a4_peaks_times[np.where((time_diffs >= min_tolerance) & (time_diffs <= max_tolerance))[0]]
  2661. # Exclude matched EEG peaks for the current EEG channel
  2662. valid_peaks_times = [t for t in valid_peaks_times if t not in matched_eeg_peaks[eeg_channel]]
  2664. if valid_peaks_times:
  2665. # Find the closest peak within the valid range
  2666. closest_eeg_peak_time = valid_peaks_times[np.argmin(np.abs(valid_peaks_times - emg_time))]
  2667. delay = emg_time - closest_eeg_peak_time
  2668. emg_delays.append(delay)
  2670. # Mark this EEG peak as matched (add to the set for this EEG channel)
  2671. matched_eeg_peaks[eeg_channel].add(closest_eeg_peak_time)
  2672. matched_delay = True
  2673. break
  2675. if not matched_delay:
  2676. # If no valid peaks are found within the tolerance range, append NaN to indicate no match
  2677. emg_delays.append(np.nan)
  2679. # Pad the emg_delays array with NaN if its length is less than the maximum number of EMG peaks
  2680. if len(emg_delays) < max_emg_peaks:
  2681. emg_delays.extend([np.nan] * (max_emg_peaks - len(emg_delays)))
  2683. delays[emg_channel] = emg_delays
  2685. return pd.DataFrame(delays)
  2687. # Calculate delays between EEG and all EMG channels with tolerance check and no reuse of EEG peaks for each channel
  2688. delay_df = calculate_delays_with_tolerance(a4_peaks_times, peaks_dict_emg, df_emg, eeg_peaks_times, min_tolerance=0.1, max_tolerance=0.6)
  2691. # %% [markdown]
  2692. # Display histogram of reaction times
  2694. # %%
  2695. import matplotlib.pyplot as plt
  2696. import pandas as pd
  2697. import numpy as np
  2698. from scipy.stats import norm
  2700. # Import data
  2701. delaymatrix = pd.read_excel(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Results\Figure_4\Reaction_time.xlsx')
  2703. # Flatten all numerical data into a single series
  2704. all_data_up = delaymatrix[delaymatrix['movement'] == 'UP'].select_dtypes(include='number').values.flatten()
  2705. all_data_fist = delaymatrix[delaymatrix['movement'] == 'FIST'].select_dtypes(include='number').values.flatten()
  2707. # Remove non-finite values (NaN and inf)
  2708. all_data_up = all_data_up[np.isfinite(all_data_up)]
  2709. all_data_fist = all_data_fist[np.isfinite(all_data_fist)]
  2711. # Calculate the Gaussian fitting parameters (mean and std) for both datasets
  2712. mu_up, std_up = norm.fit(all_data_up)
  2713. mu_fist, std_fist = norm.fit(all_data_fist)
  2715. # Plot the histogram
  2716. plt.figure(figsize=(6,6))
  2718. # Plot histogram for 'UP' data
  2719. count_up, bins_up, patches_up = plt.hist(all_data_up, bins=10, color='blue', edgecolor='black', alpha=0.6, label='Wrist extension')
  2721. # Plot Gaussian fit for 'UP' data
  2722. xmin_up, xmax_up = plt.xlim()
  2723. x_up = np.linspace(xmin_up, xmax_up, 100)
  2724. p_up = norm.pdf(x_up, mu_up, std_up)
  2725. plt.plot(x_up, p_up * len(all_data_up) * (bins_up[1] - bins_up[0]), 'blue')
  2727. # Plot histogram for 'FIST' data
  2728. count_fist, bins_fist, patches_fist = plt.hist(all_data_fist, bins=10, color='green', edgecolor='black', alpha=0.6, label='Hand flexion')
  2730. # Plot Gaussian fit for 'FIST' data
  2731. xmin_fist, xmax_fist = plt.xlim()
  2732. x_fist = np.linspace(xmin_fist, xmax_fist, 100)
  2733. p_fist = norm.pdf(x_fist, mu_fist, std_fist)
  2734. plt.plot(x_fist, p_fist * len(all_data_fist) * (bins_fist[1] - bins_fist[0]), 'green')
  2736. # Title and labels
  2737. plt.title('Histogram of All Numerical Data with Gaussian Fit')
  2738. plt.xlabel('Time [s]')
  2739. plt.ylabel('Frequency')
  2740. plt.legend(frameon=False)
  2742. # Show the plot
  2743. plt.savefig(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Results\Figure_4\Reaction_time_histogram.pdf',dpi=1000)
  2744. plt.show()
  2747. # %% [markdown]
  2748. # Display spatially the difference between reaction time across channels
  2750. # %%
  2751. from biolab.bspm import image2layout
  2753. lyt = image2layout(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\Masks\Layouts\EMG_lyt.jpg',seq='cw')
  2755. # %%
  2756. from biolab.bspm import layout2map
  2758. delaymatrix = pd.read_excel(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Results\Figure_4\Reaction_time.xlsx')
  2759. delaymatrixFIST = delaymatrix[delaymatrix['movement']=='FIST'].drop(columns=['movement'])
  2760. delaymatrixUP = delaymatrix[delaymatrix['movement']=='UP'].drop(columns=['movement'])
  2761. column_meansFIST = delaymatrixFIST.mean()
  2762. mean_dictFIST = {int(col[1:]): column_meansFIST[col] for col in delaymatrixFIST.columns}
  2763. vmapFIST = layout2map(layout=lyt,mapping=mean_dictFIST)
  2764. column_meansUP = delaymatrixUP.mean()
  2765. mean_dictUP = {int(col[1:]): column_meansUP[col] for col in delaymatrixUP.columns}
  2766. vmapUP = layout2map(layout=lyt,mapping=mean_dictUP)
  2769. # %%
  2770. from biolab.bspm import interpolate
  2772. intmapFIST = interpolate(valmap=vmapFIST,method='cubic')
  2773. intmapUP = interpolate(valmap=vmapUP,method='cubic')
  2775. # %%
  2776. plt.imshow(intmapFIST,vmin=0.275,vmax=0.34, cmap='Greens')
  2777. plt.colorbar()
  2778. plt.savefig(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Results\Figure_4\Reaction_time_mapping_FIST_greens.pdf')
  2779. plt.show()
  2781. # %% [markdown]
  2782. # ### CMC computation
  2784. # %% [markdown]
  2785. # Import multimodal data
  2787. # %%
  2788. # =====================================================
  2789. # COMPLETE CORTICO-MUSCULAR COHERENCE PIPELINE
  2790. # USING EVENTS.XLSX INTERVALS (FIST ONLY)
  2791. # REMOVES LOW-FREQUENCY PEAKS
  2792. # =====================================================
  2794. import numpy as np
  2795. import pandas as pd
  2796. import matplotlib.pyplot as plt
  2797. from scipy.signal import butter, filtfilt, coherence
  2798. from biolab.bspm import BSPM
  2800. # =====================================================
  2801. # PARAMETERS
  2802. # =====================================================
  2804. FS = 1000
  2805. BETA_BAND = (15,30)
  2806. LOW_FREQ_CUTOFF = 8 # remove <8 Hz slow components
  2808. FIST_PATH = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Cortico-muscular axis\Experiments\INTAN\Data\Multi-modal\fist.bspm'
  2809. EVENTS_PATH = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Cortico-muscular axis\Experiments\INTAN\Data\Multi-modal\events.xlsx'
  2811. # =====================================================
  2812. # LOAD BSPM
  2813. # =====================================================
  2815. def load_bspm(path):
  2816. obj = BSPM.load(loadpath=path)
  2817. r = obj.resample(newfs=FS)
  2818. r.preprocess(bw=(0.5,100))
  2819. return r
  2821. # =====================================================
  2822. # SPLIT EEG / EMG
  2823. # =====================================================
  2825. def split_eeg_emg(bspm_obj):
  2826. ch_info = bspm_obj.channel_data
  2827. data = bspm_obj.preprocessed_data
  2829. eeg_cols, emg_cols = [], []
  2831. for idx, row in ch_info.iterrows():
  2832. col_name = data.columns[idx]
  2833. if row["port_prefix"] == "A":
  2834. eeg_cols.append(col_name)
  2835. elif row["port_prefix"] == "B":
  2836. emg_cols.append(col_name)
  2838. eeg = data[eeg_cols].values
  2839. emg = data[emg_cols].values
  2841. return eeg, emg
  2843. # =====================================================
  2844. # PREPROCESS EMG (BANDPASS + RECTIFY + HIGH-PASS)
  2845. # =====================================================
  2847. def preprocess_emg(emg):
  2848. # Bandpass 5–400 Hz
  2849. b,a = butter(4, [5/(FS/2),400/(FS/2)], btype='band')
  2850. emg_filt = filtfilt(b,a, emg, axis=0)
  2851. emg_rect = np.abs(emg_filt)
  2853. # High-pass to remove <LOW_FREQ_CUTOFF Hz drift
  2854. b,a = butter(4, LOW_FREQ_CUTOFF/(FS/2), btype='high')
  2855. emg_hp = filtfilt(b,a, emg_rect, axis=0)
  2857. return emg_hp
  2859. # =====================================================
  2860. # LOAD EVENTS + EXTRACT SEGMENTS
  2861. # =====================================================
  2863. def extract_event_segments(eeg, emg, events_path):
  2864. events = pd.read_excel(events_path)
  2865. events = events[events["movement"]=="FIST"]
  2867. eeg_segments = []
  2868. emg_segments = []
  2870. for _, row in events.iterrows():
  2871. start_idx = int(row["start"] * FS)
  2872. end_idx = int(row["end"] * FS)
  2874. if end_idx > len(eeg):
  2875. continue
  2877. eeg_segments.append(eeg[start_idx:end_idx])
  2878. emg_segments.append(emg[start_idx:end_idx])
  2880. eeg_segments = np.array(eeg_segments, dtype=object)
  2881. emg_segments = np.array(emg_segments, dtype=object)
  2883. print("Number of event segments:", len(eeg_segments))
  2884. return eeg_segments, emg_segments
  2886. # =====================================================
  2887. # COMPUTE CMC
  2888. # =====================================================
  2890. def compute_cmc(eeg_segments, emg_segments, fmin=5, fmax=100, fstep=1):
  2891. # Define common frequency grid
  2892. common_freqs = np.arange(fmin, fmax+fstep, fstep)
  2893. coh_trials = []
  2895. NPERSEG = 512
  2897. for eeg_seg, emg_seg in zip(eeg_segments, emg_segments):
  2898. if len(eeg_seg) < NPERSEG:
  2899. continue
  2901. emg_mean = emg_seg.mean(axis=1)
  2903. for e in range(eeg_seg.shape[1]):
  2904. eeg_trial = eeg_seg[:,e] - np.mean(eeg_seg[:,e])
  2905. emg_trial = emg_mean - np.mean(emg_mean)
  2907. f, cxy = coherence(eeg_trial, emg_trial, fs=FS, nperseg=NPERSEG)
  2909. # Keep only fmin–fmax
  2910. mask = (f >= fmin) & (f <= fmax)
  2911. f_sel = f[mask]
  2912. cxy_sel = cxy[mask]
  2914. # Interpolate onto common grid
  2915. cxy_interp = np.interp(common_freqs, f_sel, cxy_sel)
  2916. coh_trials.append(cxy_interp)
  2918. coh_trials = np.vstack(coh_trials)
  2919. return common_freqs, coh_trials
  2921. # =====================================================
  2922. # MAIN PIPELINE
  2923. # =====================================================
  2925. print("\nProcessing FIST using event intervals")
  2927. obj = load_bspm(FIST_PATH)
  2928. eeg, emg = split_eeg_emg(obj)
  2929. emg = preprocess_emg(emg)
  2931. eeg_segments, emg_segments = extract_event_segments(
  2932. eeg,
  2933. emg,
  2934. EVENTS_PATH
  2935. )
  2937. freqs, coh_all = compute_cmc(eeg_segments, emg_segments)
  2939. # =====================================================
  2940. # BETA STRENGTH
  2941. # =====================================================
  2943. idx = (freqs>=BETA_BAND[0]) & (freqs<=BETA_BAND[1])
  2944. beta_strength = coh_all[:,idx].mean(axis=1)
  2946. # =====================================================
  2947. # DEBUG PLOT
  2948. # =====================================================
  2950. plt.figure()
  2951. plt.plot(freqs, np.mean(coh_all, axis=0))
  2952. plt.axvspan(*BETA_BAND, alpha=0.2)
  2953. plt.title("Mean CMC Spectrum (FIST — Event based, low-freq removed)")
  2954. plt.xlabel("Frequency (Hz)")
  2955. plt.ylabel("Coherence")
  2956. plt.show()
  2959. # %% [markdown]
  2960. # ### Muscular BSPM prediction from cerebral BSPM
  2962. # %% [markdown]
  2963. # Import multimodal data
  2965. # %%
  2966. from biolab.bspm import BSPM
  2968. # Read .bspm data for fist
  2969. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\Multi-modal\fist.bspm'
  2970. FIST = BSPM.load(loadpath=readpath)
  2971. rEMGFIST = FIST.resample(newfs=1000)
  2972. rEEGFIST = FIST.resample(newfs=1000)
  2973. rEMGFIST.preprocess(bw=(5,400))
  2974. rEEGFIST.preprocess(bw=(0.5,100))
  2976. # Read .bspm data for up
  2977. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\Multi-modal\up.bspm'
  2978. UP = BSPM.load(loadpath=readpath)
  2979. rEMGUP = UP.resample(newfs=1000)
  2980. rEEGUP = UP.resample(newfs=1000)
  2981. rEMGUP.preprocess(bw=(5,400))
  2982. rEEGUP.preprocess(bw=(0.5,100))
  2984. # %% [markdown]
  2985. # Define activity locations from EMG envelopes and reaction time
  2987. # %%
  2988. import copy
  2989. from scipy.signal import find_peaks
  2992. DATA = (rEMGFIST,rEEGFIST)
  2994. window_size = 500
  2995. # Load and preprocess EMG data
  2996. emgdata = copy.deepcopy(DATA[0].preprocessed_data.abs().rolling(window=window_size).mean().dropna())
  2997. emgdata = emgdata[[col for col in emgdata.columns if "B" in col]]
  2999. # Function to find peaks in EMG
  3000. def emgpeaks(df,window_size):
  3001. peaks_dict = {}
  3002. for col in df.columns:
  3003. peaks, _ = find_peaks(df[col], height=20, distance=1000) # Adjust height threshold as needed
  3004. # Adjust peak indices by the rolling window offset
  3005. adjusted_peaks = peaks + (window_size // 2)
  3006. peaks_dict[col] = adjusted_peaks
  3007. return peaks_dict
  3009. # Get peaks for EMG channels
  3010. emg_peak_locs = emgpeaks(emgdata,window_size)
  3013. # %% [markdown]
  3014. # - DEFINE WINDOW LENGTH 1s FOR BOTH EMG AND EEG ACTIVITY
  3015. # - DEFINE DELAY BETWEEN EMG AND EEG SIGNAL TO BE 300ms
  3017. # %% [markdown]
  3018. # Compute feature matrix from EEG bands power
  3020. # %%
  3021. import numpy as np
  3022. import pandas as pd
  3023. import scipy.signal as signal
  3026. eegdata = copy.deepcopy(DATA[0].preprocessed_data.abs().rolling(window=window_size).mean().dropna())
  3027. eegdata = eegdata[[col for col in eegdata.columns if "A" in col]]
  3028. sf = 1000 # Sampling frequency
  3030. eegbands = {
  3031. 'Delta':(0.5,4),
  3032. 'Theta':(4,8),
  3033. 'Alpha':(8,13),
  3034. 'Beta':(13,30),
  3035. 'Gamma':(30,100)
  3036. }
  3038. # Function to normalize EEG signals (Z-score normalization)
  3039. def normalize_signals(eeg_data):
  3040. return (eeg_data - eeg_data.mean()) / eeg_data.std()
  3042. def compute_band_power(eeg_signal, sf, band):
  3043. """Compute power spectral density (PSD) or average magnitude for a given EEG signal."""
  3044. low, high = band
  3046. freqs = np.fft.rfftfreq(len(eeg_signal), d=1/sf)
  3047. psd = np.abs(np.fft.rfft(eeg_signal))**2 / len(eeg_signal)
  3048. band_power = np.trapz(psd[(freqs >= low) & (freqs <= high)], freqs[(freqs >= low) & (freqs <= high)])
  3049. return band_power
  3051. def bandpass_filter(data, lowcut, highcut, sf, order=4):
  3052. """Apply a bandpass filter to the EEG signal."""
  3053. b, a = signal.butter(order, [lowcut, highcut], fs=sf, btype='band')
  3054. return signal.filtfilt(b, a, data)
  3056. emg_segments = [(int(peak - 0.5*sf), int(peak + 0.5*sf)) for peak in emg_peak_locs['B1']] # 1 second
  3057. eeg_segments = [(int(peak - 0.8*sf), int(peak + 0.2*sf)) for peak in emg_peak_locs['B1']] # 1 second, with 300ms delay
  3059. # Use filter with zip and unpacking
  3060. filtered = [(t1, t2) for t1, t2 in zip(emg_segments, eeg_segments) if all(x >= 0 for x in t1 + t2)]
  3061. # Unzip the filtered pairs
  3062. emg_segments, eeg_segments = zip(*filtered) if filtered else ((), ())
  3064. band_feature_list = []
  3065. for seg in eeg_segments:
  3066. segment = eegdata.iloc[seg[0]:seg[1]]
  3068. # Normalize EEG data
  3069. segment_normalized = segment.apply(normalize_signals, axis=0)
  3071. # Compute feature (power or magnitude) for each band & channel
  3072. band_feature_dict = {}
  3074. for band_name, band_range in eegbands.items():
  3075. feature_values = segment_normalized.apply(lambda x: compute_band_power(x.values, sf, band_range), axis=0)
  3077. for ch in feature_values.index:
  3078. band_feature_dict[f"{ch}_{band_name}"] = feature_values[ch]
  3080. band_feature_list.append(band_feature_dict)
  3082. # Convert list to DataFrame
  3083. feature_matrix = pd.DataFrame(band_feature_list)
  3085. # %% [markdown]
  3086. # Compute output matrix from EMG data
  3088. # %%
  3089. # Deadjust peaks for rolling averaged EMG data
  3090. window_size = 500
  3091. emg_segments = [(int(peak - 0.5*sf), int(peak + 0.5*sf)) for peak in emg_peak_locs['B1']-(window_size // 2)] # 1 second
  3092. eeg_segments = [(int(peak - 0.8*sf), int(peak + 0.2*sf)) for peak in emg_peak_locs['B1']] # 1 second, with 300ms delay
  3094. # Use filter with zip and unpacking
  3095. filtered = [(t1, t2) for t1, t2 in zip(emg_segments, eeg_segments) if all(x >= 0 for x in t1 + t2)]
  3096. # Unzip the filtered pairs
  3097. emg_segments, eeg_segments = zip(*filtered) if filtered else ((), ())
  3099. emg_feature_list = []
  3100. for seg in emg_segments:
  3101. segment = emgdata.iloc[seg[0]:seg[1]].mean()
  3102. emg_feature_list.append(segment)
  3104. # Convert list to DataFrame
  3105. output_matrix = pd.DataFrame(emg_feature_list)
  3107. # %% [markdown]
  3108. # Perform multimodal prediction
  3110. # %%
  3111. import numpy as np
  3112. import copy
  3113. from sklearn.cross_decomposition import PLSRegression
  3114. from sklearn.model_selection import train_test_split
  3115. from sklearn.metrics import r2_score, mean_squared_error
  3117. # Example data (replace with your EEG & EMG features)
  3118. X = copy.deepcopy(feature_matrix.to_numpy()) # EEG features
  3119. y = copy.deepcopy(output_matrix.to_numpy()) # EMG maps
  3121. # Train-test split
  3122. X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
  3124. # PLSR model
  3125. n_components = 10 # Tune this hyperparameter
  3126. plsr = PLSRegression(n_components=n_components)
  3127. plsr.fit(X_train, y_train)
  3129. # Predictions
  3130. y_pred = plsr.predict(X_test)
  3132. # Evaluation
  3133. r2 = r2_score(y_test, y_pred)
  3134. rmse = np.sqrt(mean_squared_error(y_test, y_pred))
  3136. print(f'R²: {r2:.4f}, RMSE: {rmse:.4f}')
  3138. # %%
  3139. import numpy as np
  3140. from scipy.stats import pearsonr, spearmanr
  3141. from sklearn.metrics.pairwise import cosine_similarity
  3143. # Assuming y_test and y_pred are already defined
  3145. # Initialize an array to store Pearson correlation for each test sample
  3146. pearson_corr = np.zeros(y_test.shape[0])
  3148. # Compute Pearson correlation for each test sample (across 16 channels)
  3149. for i in range(y_test.shape[0]):
  3150. corr, _ = pearsonr(y_test[i], y_pred[i]) # Pearson correlation per sample
  3151. pearson_corr[i] = corr
  3153. # Average correlation across all test samples (optional)
  3154. average_pearson_corr = np.mean(pearson_corr)
  3155. print(f"Average Pearson Correlation: {average_pearson_corr}")
  3156. std_pearson_corr = np.std(pearson_corr)
  3157. print(f"Standard deviation Pearson Correlation: {std_pearson_corr}")
  3159. # Initialize an array to store Spearman correlation for each test sample
  3160. spearman_corr = np.zeros(y_test.shape[0])
  3162. # Compute Spearman rank correlation for each test sample (across 16 channels)
  3163. for i in range(y_test.shape[0]):
  3164. corr, _ = spearmanr(y_test[i], y_pred[i]) # Spearman correlation per sample
  3165. spearman_corr[i] = corr
  3167. # Average rank correlation across all test samples (optional)
  3168. average_spearman_corr = np.mean(spearman_corr)
  3169. print(f"Average Spearman Rank Correlation: {average_spearman_corr}")
  3170. std_spearman_corr = np.std(spearman_corr)
  3171. print(f"Standard deviation Spearman Rank Correlation: {std_spearman_corr}")
  3173. # Initialize an array to store Cosine similarity for each test sample
  3174. cosine_sim = np.zeros(y_test.shape[0])
  3176. # Compute Cosine similarity for each test sample (across 16 channels)
  3177. for i in range(y_test.shape[0]):
  3178. cos_sim = cosine_similarity([y_test[i]], [y_pred[i]]) # Cosine similarity per sample
  3179. cosine_sim[i] = cos_sim[0][0]
  3181. # Average Cosine similarity across all test samples (optional)
  3182. average_cosine_sim = np.mean(cosine_sim)
  3183. print(f"Average Cosine Similarity: {average_cosine_sim}")
  3184. std_cosine_sim = np.std(cosine_sim)
  3185. print(f"Standard deviation Cosine Similarity: {std_cosine_sim}")
  3188. # %% [markdown]
  3189. # Display predicted vs expected maps
  3191. # %%
  3192. from biolab.bspm import image2layout
  3194. lyt = image2layout(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\Masks\Layouts\EMG_lyt.jpg',seq='cw')
  3196. # %%
  3197. from biolab.bspm import layout2map
  3200. test = pd.DataFrame(y_test,columns=output_matrix.columns)
  3201. pred = pd.DataFrame(y_pred,columns=output_matrix.columns)
  3202. mapstest = {int(col[1:]): test[col][2] for col in output_matrix.columns}
  3203. mapspred = {int(col[1:]): pred[col][2] for col in output_matrix.columns}
  3204. vmaptest = layout2map(layout=lyt,mapping=mapstest)
  3205. vmappred = layout2map(layout=lyt,mapping=mapspred)
  3207. # %%
  3208. plt.imshow(intmappred,vmin=28,vmax=36,cmap='jet')
  3209. plt.colorbar()
  3210. #plt.savefig(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Results\Figure_4\EMG_BSPM_regression_pred_FIST.pdf')
  3211. plt.show()
  3213. # %% [markdown]
  3214. # ## Scalability
  3216. # %% [markdown]
  3217. # ### Import Matlab format data from 128CH wireless system
  3219. # %% [markdown]
  3220. # Generate BSPM
  3222. # %%
  3223. import os
  3224. import h5py
  3225. import pandas as pd
  3227. from biolab.bspm import BSPM,image2layout
  3228. from biolab.utils import apply2all
  3230. fs = 1954
  3232. # Save new recordings into .bspm files
  3233. # Define parameters
  3234. readpath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Recordings\Scalability'
  3235. # Convert electrode array image into layout
  3236. EMG_126CH_lyt = image2layout(image_path=r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\Masks\Layouts\EMG_126CH_lyt.jpg',seq='cw')
  3237. # Define operation to be applied to each file
  3238. def func(filepath):
  3239. # Extract filename from filepath
  3240. filename = os.path.splitext(os.path.basename(filepath))[0]
  3241. # Split filename into sections
  3242. _,_,muscle,_ = filename.split('_',3)
  3244. meta = {
  3245. 'mode': 'EMG',
  3246. 'measurement': 'M1',
  3247. 'muscle': muscle,
  3248. 'amplifier_sample_rate': fs
  3249. }
  3251. with h5py.File(filepath, 'r') as f:
  3252. # Access the dataset inside the file
  3253. emgdata = np.delete(f['EMG_DATA'][:], [0, 64], axis=0) # remove not-connected channels
  3254. t = np.linspace(0,emgdata.shape[1]/fs,emgdata.shape[1],endpoint=False)
  3255. # Create channel names
  3256. channel_names = [f'A{i+1}' for i in range(emgdata.shape[0])]
  3257. # Create the recorded data
  3258. recdata = pd.DataFrame(emgdata.T, columns=channel_names, index=t)
  3259. # Store channel information
  3260. chns = {}
  3261. chns['custom_channel_name'] = channel_names
  3262. chns['custom_order'] = range(126)
  3263. chns['electrode_impedance_magnitude'] = np.zeros(shape=126,dtype='float')
  3264. chns['electrode_impedance_phase'] = np.zeros(shape=126,dtype='float')
  3265. chns['location'] = range(1,127)
  3266. chns['native_channel_name'] = channel_names
  3267. chns['native_order'] = range(126)
  3268. chns['port_prefix'] = np.full(126,'A')
  3269. chns['ref_channel'] = np.zeros(shape=126,dtype='float')
  3270. # Create the channel information data
  3271. chdata = pd.DataFrame(chns)
  3273. MP = BSPM(recorded_data=recdata,channel_data=chdata,layout=EMG_126CH_lyt,metadata=meta)
  3274. # Save BSPM object
  3275. MP.save(savepath=r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\Scalability\\'+filename+'.bspm')
  3277. # Call DataSaver function
  3278. apply2all(read_directory=readpath,extension='.mat',func=func)
  3281. # %% [markdown]
  3282. # ### Preprocess data and obtain relevant EMG envelopes
  3284. # %% [markdown]
  3285. # Loading and data preprocessing
  3287. # %%
  3288. from biolab.bspm import BSPM
  3290. DATA = {
  3291. 'BICEPS': r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\Scalability\Ruben_electrode_biceps_ordered.bspm',
  3292. 'BTD': r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\Scalability\Ruben_electrode_btd_ordered.bspm',
  3293. 'FIST1': r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\Scalability\Ruben_electrode_fist_ordered.bspm',
  3294. 'FIST2': r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Data\Scalability\Ruben_electrode_fist_ordered2.bspm'
  3295. }
  3297. # Load BSPM data
  3298. MP = BSPM.load(loadpath=DATA['BTD'])
  3299. # Preprocess data
  3300. rMP = MP.resample(newfs=800)
  3301. rMP.preprocess(bw=(5,399))
  3303. # %% [markdown]
  3304. # EMG envelope
  3306. # %%
  3307. from scipy.signal import find_peaks
  3309. window_size = 500
  3310. envelopes = rMP.preprocessed_data.abs().rolling(window=window_size,min_periods=1).mean().loc[70:,:]
  3311. channel = 'A126'
  3313. peaks_dict = {}
  3314. for channel in envelopes.columns:
  3315. # Find peaks in the current channel data
  3316. peaks, _ = find_peaks(envelopes[channel],distance=9*800)
  3317. peaks_dict[channel] = peaks
  3319. plt.plot(rMP.preprocessed_data[channel][70:].abs())
  3320. plt.plot(envelopes[channel])
  3321. # Mark the peaks with vertical lines (vlines) on the signal
  3322. plt.vlines(envelopes.index[peaks_dict[channel]], ymin=envelopes[channel].min(), ymax=envelopes[channel].max(),
  3323. colors='r', linestyles='--', label='Peaks')
  3324. plt.show()
  3326. # %% [markdown]
  3327. # Compute average EMG values and generate mapping
  3329. # %%
  3330. average_peaks_dict = {}
  3332. for channel, indices in peaks_dict.items():
  3333. # Get the values at the specified peak indices for the current channel
  3334. peak_values = envelopes[channel].iloc[indices]
  3336. # Compute the average of those values
  3337. average_peaks_dict[channel] = peak_values.mean()
  3340. # %% [markdown]
  3341. # Generate mapping
  3343. # %%
  3344. from biolab.bspm import layout2map
  3346. new_keys = rMP.channel_data['location'].values
  3347. # Modify dictionary to relate to locations
  3348. updated_dict = {new_keys[i]: value for i, (_,value) in enumerate(average_peaks_dict.items())}
  3350. vmap = layout2map(layout=rMP.layout,mapping=updated_dict)
  3352. # %% [markdown]
  3353. # Display BSPM
  3355. # %%
  3356. plt.imshow(vmap[:,:,0])
  3357. plt.colorbar()
  3359. # %% [markdown]
  3360. # Interpolate BSPM and save as figure
  3362. # %%
  3363. from biolab.bspm import interpolate
  3365. intmap = interpolate(valmap=vmap,method='cubic')
  3368. # %%
  3369. import matplotlib.pyplot as plt
  3371. plt.imshow(intmap,cmap='viridis')
  3372. plt.colorbar()
  3373. plt.savefig(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Results\Figure_5\Fist1_viridis.pdf')
  3374. plt.show()
  3376. # %% [markdown]
  3377. # Study differences between biceps, triceps and deltoid activity
  3379. # %%
  3380. from biolab.bspm import layout2map
  3382. muscle_peaks_dict = {}
  3384. for channel, indices in peaks_dict.items():
  3385. # Get the values at the specified peak indices for the current channel
  3386. peak_values = envelopes[channel].iloc[indices]
  3387. muscle_peaks_dict[channel] = peak_values.values[0] # 0 - Biceps, 1 - Triceps, 2 - Deltoids
  3389. new_keys = rMP.channel_data['location'].values
  3390. # Modify dictionary to relate to locations
  3391. updated_dict = {new_keys[i]: value for i, (_,value) in enumerate(muscle_peaks_dict.items())}
  3393. vmap = layout2map(layout=rMP.layout,mapping=updated_dict)
  3395. # %%
  3396. plt.imshow(vmap[:,:,0])
  3397. plt.colorbar()
  3399. # %% [markdown]
  3400. # Interpolate BSPM and save as figure
  3402. # %%
  3403. from biolab.bspm import interpolate
  3405. intmap = interpolate(valmap=vmap,method='cubic')
  3407. # %%
  3408. import matplotlib.pyplot as plt
  3410. plt.imshow(intmap,vmin=0,vmax=200,cmap='viridis')
  3411. plt.colorbar()
  3412. plt.savefig(r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Results\Figure_5\Deltoid_BTD_viridis.pdf',dpi=1000)
  3413. plt.show()
  3415. # %%
  3416. import numpy as np
  3417. import plotly.graph_objects as go
  3418. import matplotlib
  3419. import matplotlib.pyplot as plt
  3421. # --- Example Data ---
  3422. testdata = intmap
  3424. # --- Downsampling the Array for Efficiency ---
  3425. downsample_factor = (10, 5) # (rows, cols)
  3426. data_ds = testdata[::downsample_factor[0], ::downsample_factor[1]]
  3428. # --- Circular Loop Parameters ---
  3429. rows, cols = data_ds.shape
  3431. # Limit the angle range to 300 degrees, which is approximately 5.24 radians
  3432. theta = np.linspace(0, np.deg2rad(30), cols) # 300 degrees in radians
  3433. base_radius = 10 # Base radius of the loop
  3435. # --- Create 2D Grid for Circular Loop ---
  3436. theta_grid, z_grid = np.meshgrid(theta, np.arange(rows))
  3438. # --- Apply Height Scaling ---
  3439. height_factor = 0.3 # Make the shape 30% as tall
  3440. z_scaled = z_grid * height_factor
  3442. # --- Apply Width Scaling ---
  3443. width_factor = 3 # Make the shape 50% wider
  3444. bulge_factor = 0.1 # How much wider it gets at the center
  3446. # Bulging radius with width scaling applied
  3447. radius_grid = base_radius * (1 + bulge_factor * np.sin(np.pi * z_grid / rows))
  3448. radius_scaled = radius_grid * width_factor # Scale the width
  3450. # Circular coordinates with varying radius and width scaling
  3451. x_grid = radius_scaled * np.cos(theta_grid)
  3452. y_grid = radius_scaled * np.sin(theta_grid)
  3454. # --- Preprocess Data ---
  3455. # Handle NaN values by replacing them with -1
  3456. data_with_nan_replacement = np.copy(data_ds)
  3457. data_with_nan_replacement[np.isnan(data_with_nan_replacement)] = -1 # Replace NaNs with -1
  3459. # --- Create a custom jet colormap ---
  3460. jet = matplotlib.colormaps['viridis'] # Use the new Matplotlib API
  3462. # --- Create a colorscale (with -1 mapped to transparent) ---
  3463. jet_colors = jet(np.linspace(0, 1, 256)) # Get the jet colormap's RGBA values
  3464. jet_colors[0] = [0, 0, 0, 0] # Set the first color to fully transparent (for NaNs)
  3466. # Convert the RGBA colors to Plotly's 'rgba' format
  3467. jet_colors_rgba = [f'rgba({int(c[0]*255)}, {int(c[1]*255)}, {int(c[2]*255)}, {c[3]})' for c in jet_colors]
  3469. # --- Plotting with Plotly (3D Surface) ---
  3470. fig = go.Figure(data=[go.Surface(
  3471. z=z_scaled, # Apply the scaled height
  3472. x=x_grid, # Circular x coordinates with bulge and width effect
  3473. y=y_grid, # Circular y coordinates with bulge and width effect
  3474. surfacecolor=data_with_nan_replacement, # Color data
  3475. colorscale=jet_colors_rgba, # Custom colorscale with NaN handling
  3476. colorbar=dict(title="Intensity"), # Colorbar for intensity values
  3477. )])
  3479. # --- Adjusting Aspect Ratio and Axis Range ---
  3480. fig.update_layout(
  3481. scene=dict(
  3482. xaxis=dict(title='X', range=[-base_radius * width_factor * (1 + bulge_factor), base_radius * width_factor * (1 + bulge_factor)]),
  3483. yaxis=dict(title='Y', range=[-base_radius * width_factor * (1 + bulge_factor), base_radius * width_factor * (1 + bulge_factor)]),
  3484. zaxis=dict(title='Layer', range=[0, rows * height_factor]), # Adjust z-axis range for the scaled height
  3485. camera=dict(eye=dict(x=1.5, y=1.5, z=1.5)), # Adjust camera angle for better visualization
  3486. ),
  3487. margin=dict(l=0, r=0, b=0, t=50), # Adjust margins for a clean view
  3488. )
  3490. # --- Display the plot ---
  3491. fig.show()
  3494. # %%
  3495. import numpy as np
  3496. import plotly.graph_objects as go
  3497. import matplotlib
  3498. import matplotlib.pyplot as plt
  3500. # --- Example Data ---
  3501. testdata = intmap.T
  3503. # --- Downsampling the Array for Efficiency ---
  3504. downsample_factor = (5, 5) # (rows, cols)
  3505. data_ds = testdata[::downsample_factor[0], ::downsample_factor[1]]
  3507. # --- Circular Loop Parameters ---
  3508. rows, cols = data_ds.shape
  3510. # Limit the angle range to 300 degrees, which is approximately 5.24 radians
  3511. theta = np.linspace(0, np.deg2rad(300), cols) # 300 degrees in radians
  3512. base_radius = 10 # Base radius of the loop
  3514. # --- Create 2D Grid for Circular Loop ---
  3515. theta_grid, z_grid = np.meshgrid(theta, np.arange(rows))
  3517. # --- Apply Height Scaling ---
  3518. height_factor = 0.15 # Make the shape 30% as tall
  3519. z_scaled = z_grid * height_factor
  3521. # --- Apply Width Scaling ---
  3522. width_factor = 3 # Make the shape 50% wider
  3523. bulge_factor = 0.15 # How much wider it gets at the center
  3525. # Bulging radius with width scaling applied
  3526. radius_grid = base_radius * (1 + bulge_factor * np.sin(np.pi * z_grid / rows))
  3527. radius_scaled = radius_grid * width_factor # Scale the width
  3529. # Circular coordinates with varying radius and width scaling
  3530. x_grid = radius_scaled * np.cos(theta_grid)
  3531. y_grid = radius_scaled * np.sin(theta_grid)
  3533. # --- Preprocess Data ---
  3534. # Handle NaN values by replacing them with -1
  3535. data_with_nan_replacement = np.copy(data_ds)
  3536. data_with_nan_replacement[np.isnan(data_with_nan_replacement)] = -1 # Replace NaNs with -1
  3538. # --- Create a custom jet colormap ---
  3539. jet = matplotlib.colormaps['viridis'] # Use the new Matplotlib API
  3541. # --- Create a colorscale (with -1 mapped to transparent) ---
  3542. jet_colors = jet(np.linspace(0, 1, 256)) # Get the jet colormap's RGBA values
  3543. jet_colors[0] = [0, 0, 0, 0] # Set the first color to fully transparent (for NaNs)
  3545. # Convert the RGBA colors to Plotly's 'rgba' format
  3546. jet_colors_rgba = [f'rgba({int(c[0]*255)}, {int(c[1]*255)}, {int(c[2]*255)}, {c[3]})' for c in jet_colors]
  3548. # --- Plotting with Plotly (3D Surface) ---
  3549. fig = go.Figure(data=[go.Surface(
  3550. z=z_scaled, # Apply the scaled height
  3551. x=x_grid, # Circular x coordinates with bulge and width effect
  3552. y=y_grid, # Circular y coordinates with bulge and width effect
  3553. surfacecolor=data_with_nan_replacement, # Color data
  3554. colorscale=jet_colors_rgba, # Custom colorscale with NaN handling
  3555. colorbar=dict(title="Intensity"), # Colorbar for intensity values
  3556. )])
  3558. # --- Adjusting Aspect Ratio and Axis Range ---
  3559. fig.update_layout(
  3560. scene=dict(
  3561. xaxis=dict(title='X', range=[-base_radius * width_factor * (1 + bulge_factor), base_radius * width_factor * (1 + bulge_factor)]),
  3562. yaxis=dict(title='Y', range=[-base_radius * width_factor * (1 + bulge_factor), base_radius * width_factor * (1 + bulge_factor)]),
  3563. zaxis=dict(title='Layer', range=[0, rows * height_factor]), # Adjust z-axis range for the scaled height
  3564. camera=dict(eye=dict(x=1.5, y=1.5, z=1.5)), # Adjust camera angle for better visualization
  3565. ),
  3566. margin=dict(l=0, r=0, b=0, t=50), # Adjust margins for a clean view
  3567. )
  3569. # --- Display the plot ---
  3570. fig.show()
  3573. # %%
  3574. from biolab.bspm import interpolate
  3576. intmaptest = interpolate(valmap=vmaptest,method='cubic')
  3577. intmappred = interpolate(valmap=vmappred,method='cubic')
  3579. # %% [markdown]
  3580. # ## Convert `jet` colorspace to `viridis`
  3582. # %%
  3583. import fitz # PyMuPDF
  3584. import numpy as np
  3585. import cv2
  3586. import matplotlib.pyplot as plt
  3587. import matplotlib.cm as cm
  3588. from matplotlib.colors import Normalize
  3589. from PIL import Image
  3591. def extract_image_from_pdf(pdf_path):
  3592. doc = fitz.open(pdf_path)
  3593. for page_index in range(len(doc)):
  3594. page = doc.load_page(page_index)
  3595. images = page.get_images(full=True)
  3596. for img_index, img in enumerate(images):
  3597. xref = images[img_index][0]
  3598. base_image = doc.extract_image(xref)
  3599. image_bytes = base_image["image"]
  3600. img = cv2.imdecode(np.frombuffer(image_bytes, np.uint8), cv2.IMREAD_COLOR)
  3601. return img
  3602. return None
  3604. def convert_jet_to_viridis(image):
  3605. image_rgb = cv2.cvtColor(image, cv2.COLOR_BGR2RGB)
  3607. # Convert to float32 and normalize
  3608. img_float = image_rgb.astype(np.float32) / 255.0
  3610. # Simulate inverse of `jet` colormap
  3611. jet = cm.get_cmap('jet', 256)
  3612. jet_colors = (jet(np.linspace(0, 1, 256))[:, :3]) # ignore alpha
  3614. # Reshape image and jet LUT
  3615. pixels = img_float.reshape(-1, 3)
  3616. distances = np.linalg.norm(pixels[:, None] - jet_colors[None, :], axis=2)
  3617. jet_indices = np.argmin(distances, axis=1)
  3619. # Normalize back to data values
  3620. normed_data = jet_indices / 255.0
  3622. # Apply viridis colormap
  3623. viridis = cm.get_cmap('viridis')
  3624. viridis_img = viridis(normed_data)[:, :3]
  3625. viridis_img = (viridis_img.reshape(image.shape[0], image.shape[1], 3) * 255).astype(np.uint8)
  3627. return cv2.cvtColor(viridis_img, cv2.COLOR_RGB2BGR)
  3629. def save_image_as_pdf(image, output_pdf_path):
  3630. img_pil = Image.fromarray(cv2.cvtColor(image, cv2.COLOR_BGR2RGB))
  3631. img_pil.save(output_pdf_path, "PDF")
  3633. # === USAGE ===
  3634. filepath = r'C:\Users\Ruben\OneDrive - University of Cambridge\PhD\BSPM_for_ECGi\Publications\Wireless High-Density E-textile Eutectogel Electrode Arrays for Spatio-Temporal Machine Learning in Cutaneous Electrophysiology\Experiments\INTAN\Results\Figure_4\EMG_BSPM_regression_test_UP'
  3635. input_pdf = filepath+'.pdf'
  3636. output_pdf = filepath+'_viridis.pdf'
  3638. image = extract_image_from_pdf(input_pdf)
  3639. if image is not None:
  3640. viridis_image = convert_jet_to_viridis(image)
  3641. save_image_as_pdf(viridis_image, output_pdf)
  3642. print(f"Saved converted PDF as {output_pdf}")
  3643. else:
  3644. print("No image found in PDF.")
  3647. # %%
  3648. import numpy as np
  3649. import matplotlib.pyplot as plt
  3650. import seaborn as sns
  3652. def plot_confusion_matrix(cm, labels=None, title='Confusion Matrix', pdf_path='confusion_matrix.pdf'):
  3653. fig, ax = plt.subplots(figsize=(6, 5), constrained_layout=True)
  3655. sns.heatmap(cm, annot=True, fmt='d', cmap='Blues',
  3656. xticklabels=labels, yticklabels=labels,
  3657. cbar=True, square=True, linewidths=0.5, vmin=0, vmax=20, ax=ax)
  3659. ax.set_title(title)
  3660. ax.set_xlabel('Predicted')
  3661. ax.set_ylabel('Actual')
  3663. # Save to PDF
  3664. fig.savefig(pdf_path, format='pdf')
  3665. print(f"Saved confusion matrix to: {pdf_path}")
  3667. # === Example Usage ===
  3668. conf_matrix = np.array([[13, 4, 3],
  3669. [10, 11, 0],
  3670. [8, 7, 5]])
  3672. class_labels = ['Class A', 'Class B', 'Class C']
  3674. plot_confusion_matrix(conf_matrix, labels=class_labels, pdf_path='my_confusion_matrix.pdf')

Wireless_BSPM.ipynb at commit c1ccfee, no license · at the source

Overview

Authors: Ruben Ruiz-Mateos Serrano1,2, Charlie Brunt1,2, Xudong Tao1,2, Maciej Zajaczkowski3, Antonio Dominguez-Alfaro2,4, Matias L. Picchio5,6, Daniele Mantione7,5, Emmanuel M. Drakakis3, David Mecerreyes7,5, George G. Malliaras1,2
ORCID iDs: Ruben Ruiz-Mateos Serrano, Xudong Tao, Maciej Zajaczkowski, Matias L. Picchio, Daniele Mantione, George G. Malliaras
  1. Institute for Biomedical Innovation, University of Cambridge,Cambridge, UK
  2. Electrical Engineering Division, Department of Engineering, University of Cambridge,Cambridge, UK
  3. Department of Bioengineering, Faculty of Engineering, Imperial College London,London, UK
  4. Present Address: Instituto de Microelectrónica, IMSE-CNM, (CSIC Universidad de Sevilla),Av. Américo Vespucio 28, 41092 Sevilla, Spain
  5. IKERBASQUE, Basque Foundation for Science,Bilbao, Spain
  6. POLYMAT, Department of Mining-Metallurgy Engineering and Materials Science, School of Engineering, University of the Basque Country (UPV/EHU),Bilbao, Spain
  7. POLYMAT, University of the Basque Country UPV/EHU,Av.Tolosa 72, 20018 Donostia-San Sebastian, Gipuzkoa Spain
Institutions: University of Cambridge (United Kingdom); Imperial College London (United Kingdom); Universidad de Sevilla (Spain); Ikerbasque (Spain); University of the Basque Country (Spain)
Journal: Nature communications, volume 17, issue 1, article 8441
Dates: received 14 October 2025; accepted 17 June 2026; published online 8 July 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41467-026-75134-1 · PMID 42420293 · PMCID PMC13478180 · OpenAlex W7167728318
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: other (modality), human (organism)
Methods: Spectral & time-frequency, Preprocessing, Connectivity, Smoothing, state filtering, decompositions, Machine learning, Evoked potentials, Physiology & signal measures
Keywords: Biomedical engineering, Electrophysiology, Preclinical research, Sensorimotor processing
MeSH: Body Surface Potential Mapping*, Muscle, Skeletal*, Textiles*, Algorithms, Electrodes, Humans, Machine Learning, Wearable Electronic Devices (* major topic)
Topic: Neuroscience and Neural Engineering (Cellular and Molecular Neuroscience, Neuroscience), according to OpenAlex
Funding: Engineering and Physical Sciences Research Council (EP/S022139/1); Ikerbasque, Basque Foundation for Science (RYC2024-048187-I); Eusko Jaurlaritza (Basque Government) (IT1766-22)
Citations: not cited yet (Europe PMC); 64 references in the paper

Abstract

Cutaneous electrophysiology is a fundamental non-invasive technique for assessing electrically active organs such as the brain, heart, and muscles. Standard approaches, however, are limited in spatial resolution, reducing sensitivity to certain pathological features. The development of body surface potential mapping using electrode arrays has helped overcome these limitations, enhancing the diagnostic power of cutaneous recordings, yet clinical adoption remains constrained by challenges in electrode performance, wiring complexity, wearability, data transmission, and interpretability. Here, we present a hybrid e-textile electrode array system that overcomes these barriers, enabling simultaneous mapping of electrical activity along the cortico-muscular axis. The system combines application-specific conducting polymer coatings to improve electrode performance, a flexible fabrication process for robust connectivity and wearability, and interpretable machine learning algorithms for data analysis. In controlled single-subject experiments, we demonstrate reliable muscle and brain recordings, enabling classification of grasped object shapes and somatosensory stimuli. Simultaneous multi-site recordings along the cortico-muscular axis provide spatial maps of reaction time distributions and allow prediction of muscle activation patterns from cortical activity. This platform establishes a framework for wearable, multi-modal electrophysiological mapping and non-invasive study of cortico-muscular dynamics, representing a step towards practical brain–body interfaces with applications in neurorehabilitation, prosthetics, and human–machine interaction.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 19 matches between paragraphs and lines of code.

Zenodo 20041637

License: CC-BY-4.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: the text, “Preprocessing of recorded electrophysiological s”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Matplotlib (2 files), NumPy (2 files), pandas (2 files), scikit-learn (2 files), SciPy (2 files), seaborn (2 files), h5py (1 file), OpenCV (1 file), Pillow (1 file), Plotly (1 file), specparam (formerly FOOOF) (1 file), WFDB Python (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
3 files
At the source:

rr1017/cortico-muscular-axis-mapping

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: c1ccfee5ca3314ef5c3c2c2fb8b7aa6a14008166, 28 April 2026
Languages: Jupyter (3)
Size: 3 files, 3 scripts
Software Heritage: not archived
Found in: the Zenodo archive record
Holds: 3 notebooks
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Matplotlib (2 files), NumPy (2 files), pandas (2 files), scikit-learn (2 files), SciPy (2 files), seaborn (2 files), h5py (1 file), OpenCV (1 file), Pillow (1 file), Plotly (1 file), specparam (formerly FOOOF) (1 file), WFDB Python (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
3 files
At the source: github.com/rr1017/cortico-muscular-axis-mapping/tree/zenodo

Code availability

The codes used for collecting data are available from the corresponding author upon reasonable request.

Reproduced under the paper's license (CC BY), from the paper cited above.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 6 scripts, each with its path and the digest of its content;
  • 19 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability

All data supporting the findings of this study are available within the article and its supplementary files. Any additional requests for information can be directed to, and will be fulfilled by, the corresponding author. Source data are provided with this paper. The code used to process and analyse the data in this study is publicly available at 10.5281/zenodo.20041637. The electrophysiological datasets generated in this study consist of human EEG and EMG recordings and are not publicly available due to ethical and data protection constraints, including participant consent limitations and compliance with applicable data privacy regulations. Access to a de-identified minimum dataset sufficient to reproduce the analyses may be granted upon request for non-commercial research purposes, subject to institutional data-sharing agreements and ethical approval where required. Requests should be directed to the corresponding author. Applicants must provide a brief research proposal outlining the intended use of the data. Requests are typically reviewed within 2–4 weeks. Approved users will be granted access to the data for a defined period and under conditions that prohibit re-identification or onward sharing. Source data are provided with this paper.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 10 authors, 4 keywords, 8 MeSH terms, 3 funders, 58 references.

Cite

This paper

Serrano, R. R.-M., Brunt, C., Tao, X., Zajaczkowski, M., Dominguez-Alfaro, A., Picchio, M. L., Mantione, D., Drakakis, E. M., Mecerreyes, D., & Malliaras, G. G. (2026). Body surface potential mapping of the cortico-muscular axis using smart textile electrode arrays. Nature communications, 17(1), 8441. https://doi.org/10.1038/s41467-026-75134-1

BibTeX

@article{serrano2026body,
author = {Serrano, Ruben Ruiz-Mateos and Brunt, Charlie and Tao, Xudong and Zajaczkowski, Maciej and Dominguez-Alfaro, Antonio and Picchio, Matias L. and Mantione, Daniele and Drakakis, Emmanuel M. and Mecerreyes, David and Malliaras, George G.},
title = {{Body surface potential mapping of the cortico-muscular axis using smart textile electrode arrays}},
journal = {Nature communications},
year = {2026},
month = jul,
volume = {17},
number = {1},
pages = {8441},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-75134-1},
url = {https://doi.org/10.1038/s41467-026-75134-1},
pmid = {42420293},
pmcid = {PMC13478180}
}

RIS

TY - JOUR
AU - Serrano, Ruben Ruiz-Mateos
AU - Brunt, Charlie
AU - Tao, Xudong
AU - Zajaczkowski, Maciej
AU - Dominguez-Alfaro, Antonio
AU - Picchio, Matias L.
AU - Mantione, Daniele
AU - Drakakis, Emmanuel M.
AU - Mecerreyes, David
AU - Malliaras, George G.
TI - Body surface potential mapping of the cortico-muscular axis using smart textile electrode arrays
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/07/08
VL - 17
IS - 1
SP - 8441
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-75134-1
UR - https://doi.org/10.1038/s41467-026-75134-1
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-75134-1",
"type": "article-journal",
"title": "Body surface potential mapping of the cortico-muscular axis using smart textile electrode arrays",
"container-title": "Nature communications",
"author": [
{
"family": "Serrano",
"given": "Ruben Ruiz-Mateos"
},
{
"family": "Brunt",
"given": "Charlie"
},
{
"family": "Tao",
"given": "Xudong"
},
{
"family": "Zajaczkowski",
"given": "Maciej"
},
{
"family": "Dominguez-Alfaro",
"given": "Antonio"
},
{
"family": "Picchio",
"given": "Matias L."
},
{
"family": "Mantione",
"given": "Daniele"
},
{
"family": "Drakakis",
"given": "Emmanuel M."
},
{
"family": "Mecerreyes",
"given": "David"
},
{
"family": "Malliaras",
"given": "George G."
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "8441",
"DOI": "10.1038/s41467-026-75134-1",
"PMID": "42420293",
"PMCID": "PMC13478180",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-75134-1",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
8
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1002/advs.202519893 [code]
NeuroSuite for Long-Term Functional and Structural Studies of Air-Liquid Interface Cerebral Organoids.
Journal: Advanced science (Weinheim, Baden-Wurttemberg, Germany)
In common: specparam (formerly FOOOF), Pillow, pandas, 3 other tools, 2 authors
[2] doi:10.7554/elife.100605 [code]
Age-related changes in ‘cortical’ 1/f dynamics are linked to cardiac activity
Journal: n/a
In common: WFDB Python, specparam (formerly FOOOF), h5py, 6 other tools, other
[3] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: WFDB Python, Plotly, OpenCV, 7 other tools
[4] doi:10.1016/j.xcrm.2026.102766 [code]
A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: Plotly, OpenCV, h5py, 7 other tools
[5] doi:10.1038/s41597-025-05174-7 [code]
A large-scale MEG and EEG dataset for object recognition in naturalistic scenes
Journal: n/a
In common: Plotly, OpenCV, h5py, 7 other tools
[6] doi:10.1016/j.isci.2026.116825 [code]
Social hierarchy shapes behavioral and transcriptional responses to chronic stress and ketamine in male mice.
Journal: iScience
In common: Plotly, OpenCV, h5py, 7 other tools
[7] doi:10.1162/imag.a.1269 [code]
From early to contemporary normative modeling: Mapping individual differences in neurophysiological signals.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: specparam (formerly FOOOF), Plotly, seaborn, 5 other tools, 1 reference
[8] doi:10.3389/fncom.2026.1786996 [code]
Schumann-anchored golden ratio organization of human neural oscillations.
Journal: Frontiers in computational neuroscience
In common: specparam (formerly FOOOF), h5py, Pillow, 6 other tools
[9] doi:10.1371/journal.pone.0348866 [code]
Using deep learning to identify inherited retinal diseases based on wide-field retinal imaging data.
Journal: PloS one
In common: Plotly, OpenCV, Pillow, 6 other tools, other
[10] doi:10.1038/s41467-026-75128-z
Surface circumferential spinal cord recording in freely moving rodents.
Journal: Nature communications
In common: 2 authors

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.