In [1]:
import scipy.io as sio
from matplotlib import pyplot as plt
%matplotlib inline
import numpy as np
import os
import mne
import numpy as np
import scipy.io as sio
In [2]:
# coding: utf-8
def _loadftfile(path):
filecontents = sio.whosmat(path)
strucname = filecontents[0][0]
mat = sio.loadmat(path, squeeze_me=True, struct_as_record=False)
matstruct = mat[strucname]
return matstruct
def _matstruc2mne(matstruct,ch_names=None):
if ch_names is None:
ch_names=list(matstruct.label)
myinfo = mne.create_info(ch_names=ch_names,sfreq=1/(matstruct.time[1] - matstruct.time[0]),ch_types='eeg')
ev_arr = mne.EvokedArray(matstruct.individual.mean(axis=0),myinfo,tmin=-0.2) ### Specific to this dataset, 200ms baseline
ev_arr.set_montage(mne.channels.read_montage("EGI_256"))
return ev_arr
def _matstruc2latency(matstruct,peak_tmin,peak_tmax,ch_names=None):
if ch_names is None:
ch_names=list(matstruct.label)
myinfo = mne.create_info(ch_names=ch_names,sfreq=1/(matstruct.time[1] - matstruct.time[0]),ch_types='eeg')
all_chpeaks = []
all_lat = []
all_amp = []
for mat in matstruct.individual:
ev_arr = mne.EvokedArray(mat,myinfo,tmin=-0.2) ### Specific to this dataset, 200ms baseline
ev_arr.set_montage(mne.channels.read_montage("EGI_256"))
chpeak,lat,amp = ev_arr.get_peak(tmin=peak_tmin,tmax=peak_tmax,return_amplitude=True,mode='neg')
all_lat.append(lat)
all_chpeaks.append(chpeak)
all_amp.append(amp)
return all_chpeaks,all_lat,all_amp
def _matstruc2mne_epochs(matstruct,ch_names=None):
if ch_names is None:
ch_names=list(matstruct.label)
myinfo = mne.create_info(ch_names=ch_names,sfreq=1/(matstruct.time[1] - matstruct.time[0]),ch_types='eeg')
ev_arr = mne.EpochsArray(matstruct.individual,myinfo,tmin=-0.2) ### Specific to this dataset, 200ms baseline
ev_arr.set_montage(mne.channels.read_montage("EGI_256"))
return ev_arr
In [3]:
tcrop = 0.7
In [5]:
matfile = '/home/nfarrugi/datasets/eeg4sounds/result-eeg4sounds/oddball/grav/grav_bin_dev.mat'
mat_bin_dev = _loadftfile(matfile)
matfile = '/home/nfarrugi/datasets/eeg4sounds/result-eeg4sounds/oddball/grav/grav_bin_std.mat'
mat_bin_std = _loadftfile(matfile)
matfile = '/home/nfarrugi/datasets/eeg4sounds/result-eeg4sounds/oddball/grav/grav_ste_dev.mat'
mat_ste_dev = _loadftfile(matfile)
matfile = '/home/nfarrugi/datasets/eeg4sounds/result-eeg4sounds/oddball/grav/grav_ste_std.mat'
mat_ste_std = _loadftfile(matfile)
In [6]:
ev_bin_dev = _matstruc2mne(mat_bin_dev).crop(tmax=tcrop)
ev_bin_std = _matstruc2mne(mat_bin_std).crop(tmax=tcrop)
ev_ste_dev = _matstruc2mne(mat_ste_dev).crop(tmax=tcrop)
ev_ste_std = _matstruc2mne(mat_ste_std).crop(tmax=tcrop)
mne.equalize_channels([ev_bin_dev,ev_ste_dev,ev_bin_std,ev_ste_std])
In [6]:
epochs_bin_dev = _matstruc2mne_epochs(mat_bin_dev).crop(tmax=tcrop)
epochs_bin_std = _matstruc2mne_epochs(mat_bin_std).crop(tmax=tcrop)
epochs_ste_dev = _matstruc2mne_epochs(mat_ste_dev).crop(tmax=tcrop)
epochs_ste_std = _matstruc2mne_epochs(mat_ste_std).crop(tmax=tcrop)
mne.equalize_channels([epochs_bin_dev,epochs_bin_std,epochs_ste_dev,epochs_ste_std])
In [7]:
peak_tmin = 0.30
peak_tmax = 0.35
ch_lat_bin_dev,lat_bin_dev,amp_bin_dev = _matstruc2latency(mat_bin_dev,peak_tmin=peak_tmin,peak_tmax=peak_tmax)
#ev_bin_std = _matstruc2latency(mat_bin_std)
ch_lat_ste_dev,lat_ste_dev,amp_ste_dev = _matstruc2latency(mat_ste_dev,peak_tmin=peak_tmin,peak_tmax=peak_tmax)
#ev_ste_std = _matstruc2latency(mat_ste_std)
In [9]:
plt.hist(lat_bin_dev,alpha=0.5)
plt.hist(lat_ste_dev,alpha=0.5)
plt.legend(['binaural','stereo'])
Out[9]:
In [10]:
plt.hist(amp_bin_dev,alpha=0.5)
plt.hist(amp_ste_dev,alpha=0.5)
plt.legend(['binaural','stereo'])
Out[10]:
In [11]:
plt.boxplot([lat_bin_dev,lat_ste_dev])
Out[11]:
In [7]:
T_obs_bin,clusters_bin,clusters_pb_bin,H0_bin = mne.stats.spatio_temporal_cluster_test(X_bin,threshold=None,n_permutations=1000,out_type='mask')
In [91]:
T_obs_ste,clusters_ste,clusters_pb_ste,H0_ste = mne.stats.spatio_temporal_cluster_test(X_ste,threshold=None,n_permutations=1000,out_type='mask')
In [92]:
T_obs_dev,clusters_dev,clusters_pb_dev,H0_dev = mne.stats.spatio_temporal_cluster_test(X_dev,threshold=12,n_permutations=1000,out_type='mask')
In [93]:
T_obs_stan,clusters_stan,clusters_pb_stan,H0_stan = mne.stats.spatio_temporal_cluster_test(X_stan,threshold=None,n_permutations=1000,out_type='mask')
In [109]:
evoked = mne.combine_evoked([ev_bin_dev, -ev_ste_dev],
weights='equal') # calculate difference wave
time_unit = dict(time_unit="ms")
evoked.plot_joint(title="Deviants Binaural vs. Stereo")#, ts_args=time_unit,topomap_args=time_unit) # show difference wave
plt.show()
for clust,pval in zip(clusters_dev,clusters_pb_dev):
if pval < 0.1:
evoked.plot_image(mask=clust.T,mask_style='contour')
plt.show()
In [106]:
evoked = mne.combine_evoked([ev_bin_dev, -ev_bin_std],
weights='equal') # calculate difference wave
time_unit = dict(time_unit="ms")
evoked.plot_joint(title="Binaural Deviant vs. Std")#, ts_args=time_unit,topomap_args=time_unit) # show difference wave
plt.show()
for clust,pval in zip(clusters_bin,clusters_pb_bin):
if pval < 0.005:
evoked.plot_image(mask=clust.T,mask_style='contour')
plt.show()
clustersize = np.sum(clust)
print(clustersize)