Import necessary libraries and functions
In [1]:
import numpy as np, cmath, scipy as sp
import scipy.io
from matplotlib import pyplot as plt
#import basic functions from numpy that we'll need
from numpy import pi, sin, cos, exp, sqrt, log, log10, random, angle, real, imag
from numpy import zeros, ceil, floor, absolute, linspace
from numpy.fft import fft, ifft
from scipy import signal as sig
from scipy.signal import hilbert
from matplotlib.pyplot import *
%matplotlib inline
Import optional library for prettier plots
In [2]:
import seaborn as sns
sns.set_palette('muted')
sns.set_style('darkgrid')
Import EEGdata from mat file using scipy.io
In [4]:
data = scipy.io.loadmat('sampleEEGdata')
EEGdata = data["EEG"][0,0]["data"]
EEGpnts = data["EEG"][0,0]["pnts"][0,0] #number of points in EEG data
EEGtimes = data["EEG"][0,0]["times"][0]
EEGsrate = float(data["EEG"][0,0]["srate"][0]) #make float for division purposes later
EEGtrials = data["EEG"][0,0]["trials"][0,0]
EEGnbchan = data["EEG"][0,0]["nbchan"][0,0]
EEGchanlocslabels=data["EEG"][0,0]["chanlocs"][0]["labels"]
In [13]:
a=[0,1]
a[:1]
Out[13]:
In [131]:
timewin = 500 # in ms
# convert ms to idx
timewinidx = int(np.round(timewin/(1000/EEGsrate)));
# create hann taper function
hann_win = .5*(1-cos(2*pi*np.arange(timewinidx)/(timewinidx-1)));
#detrend data (useful to attentuate super-low frequency artifacts in FFT
# from sampled data)
d = sig.detrend(EEGdata[19,:,15]);
subplot(311)
plot(EEGtimes,d)
title('One trial of data')
stime = np.argmin(absolute(EEGtimes-(-50)))
subplot(323)
plot(EEGtimes[stime:stime+timewinidx],d[stime:stime+timewinidx])
plot(EEGtimes[stime:stime+timewinidx],d[stime:stime+timewinidx]*hann_win,'r')
setp(gca(),'xlim',[-50, -50+timewin])
title('One short-time window of data, windowed')
tight_layout()
dfft = fft(d[stime:stime+timewinidx]*hann_win)
f = np.linspace(0,EEGsrate/2.,np.floor(len(hann_win)/2)+1) # frequencies of FFT
subplot(313)
plot(f[1:],absolute(dfft[1:int(floor(len(hann_win)/2.))+1])**2,'.-');
title('power spectrum from that time window')
_=setp(gca(),'xlim',[1 ,128],'ylim',[-1000, 25000],'xticks',np.arange(0,EEGsrate/2.,10))
In [64]:
# create TF matrix and input column of data at selected time point
tf = np.zeros([np.floor(len(hann_win)/2),EEGpnts])
tf[:,int(stime+timewinidx/2-11):int(stime+timewinidx/2+10)] = np.tile(absolute(dfft[1:int(floor(len(hann_win)/2.)+1)])*2,[21,1]).T
imshow(log10(tf+1),
aspect='auto',
extent=[0,700,80,0],
cmap=cm.jet)
grid(False)
In [171]:
timewin = 400 # in ms, for stFFT
times2save = np.arange(-300,1000 + 50,50) # in ms
channel2plot = 'P7'
frequency2plot = 15 # in Hz
timepoint2plot = 200 # ms
# convert from ms to index
times2saveidx = np.zeros(times2save.shape);
for i in xrange(len(times2save)):
times2saveidx[i]=np.argmin(absolute(EEGtimes-times2save[i]))
timewinidx = np.round(timewin/(1000/EEGsrate));
chan2useidx = EEGchanlocslabels == channel2plot
#create hann taper
hann_win = .5*(1-cos(2*pi*(np.arange(timewinidx))/(timewinidx-1)));
hann_win = np.reshape(hann_win,[len(hann_win),1])
#define frequencies
frex = linspace(0,EEGsrate/2,floor(timewinidx/2)+1);
#initialize power output matrix
tf = np.zeros([len(frex),len(times2save)])
#loop over time points and perform FFT
for timepointi in xrange(len(times2save)):
#extract time series data for this center time point
# note: the 'mod' function here corrects for even or odd number of points
selectedData = EEGdata[chan2useidx, int(times2saveidx[timepointi]-floor(timewinidx/2)-1):int(times2saveidx[timepointi]+floor(timewinidx/2)-((timewinidx+1)%2)),:]
tempdat = np.squeeze(selectedData)
#multiply the two arrays together by broadcasting
taperdat = tempdat * hann_win
fdat = fft(taperdat,axis=0)/timewinidx #
tf[:,timepointi] = np.mean(absolute(fdat[:int(floor(timewinidx/2)+1),:])**2,axis=1)
subplot(121)
freq2plotidx=np.argmin(absolute(frex-frequency2plot));
plot(times2save,np.mean(log10(tf[freq2plotidx-3:freq2plotidx+3,:]),axis=0))
title( 'Sensor ' + channel2plot + ', ' + str(frequency2plot) +' Hz' )
_=setp(gca(),'xlim',[times2save[0] ,times2save[-1]])
subplot(122)
time2plotidx=np.argmin(absolute(times2save-timepoint2plot));
plot(frex,log10(tf[:,time2plotidx]))
title( 'Sensor ' + channel2plot + ', ' + str(timepoint2plot) +' ms' )
_=setp(gca(),'xlim',[frex[0], 40])
In [174]:
contourf(times2save,frex,log10(tf),40,cmap=cm.jet)
title([ 'Sensor ' + channel2plot + ', power plot (no baseline correction)' ])
print 'Overlap of ' + str(100*(1-np.mean(np.diff(times2save))/timewin)) + '%'
In [120]:
#create hamming
hamming_win = .54 - .46*cos(2*pi*(np.arange(timewinidx))/(timewinidx-1))
hann_win = .5*(1-cos(2*pi*(np.arange(timewinidx))/(timewinidx-1)))
#create gaussian
gaus_win = exp(-.5*(2.5*(np.arange(-timewinidx/2,timewinidx/2))/(timewinidx/2))**2)
#plot together
plot(hann_win)
plot(hamming_win,'r')
plot(gaus_win,'k')
legend(['Hann','Hamming','Gaussian'])
_=setp(gca(),'xlim',[-5, timewinidx+5],'ylim',[-.1, 1.1],'yticks',np.arange(0,1.2,.2))
In [175]:
chan2use = 'O1';
chan2useidx = EEGchanlocslabels == chan2use
frequency2plot = 10 # in Hz
freq2plotidx=np.argmin(absolute(frex-frequency2plot))
# initialize ITPC output matrix
itpc = np.zeros([len(frex),len(times2save)])
for timepointi in xrange(len(times2save)):
# extract time series data for this center time point
selectedData = EEGdata[chan2useidx, int(times2saveidx[timepointi]-floor(timewinidx/2)-1):int(times2saveidx[timepointi]+floor(timewinidx/2)-((timewinidx+1)%2)),:]
tempdat = np.squeeze(selectedData)
taperdat = tempdat * hann_win
fdat = fft(taperdat,axis=0)/timewinidx;
itpc[:,timepointi] = absolute(np.mean(exp(1j*angle(fdat[:int(floor(timewinidx/2)+1),:])),axis=1)) # average over trials
contourf(times2save,frex,itpc,40,cmap=cm.jet)
_=setp(gca(),'xlim',[-200, 1000])
In [167]:
plot(times2save,np.mean(itpc[freq2plotidx-2:freq2plotidx+2,:],axis=0))
_=title([ 'ITPC at sensor ' + chan2use +', ' +str(frequency2plot) + ' Hz' ])