In [1]:

    

%matplotlib inline
from matplotlib import pylab as pl
import cPickle as pickle
import pandas as pd
import numpy as np
import os
import random


    

In [2]:

    

import sys
sys.path.append('..')


    

In [120]:

    

target = 'Dog_4'
data_type = 'test' # preictal interictal, test


    

In [121]:

    

channel = 0


    

In [122]:

    

import scipy.io

nchain = 0 # which chain you want

all_data = last_sequence = last_data_length_sec = last_Fs = last_channels = last_d_shape = None
for segment in range(10000):
    fname = '../seizure-data/%s/%s_%s_segment_%04d.mat'%(target,target,data_type,segment+1)
    try:
        data = scipy.io.loadmat(fname)
    except:
        break
    k = '%s_segment_%d'%(data_type,segment+1)
    data_length_sec = data[k]['data_length_sec'][0,0][0,0]
    try:
        sequence = data[k]['sequence'][0,0][0,0]
    except:
        sequence = 1
    Fs = float(data[k]['sampling_frequency'][0,0][0,0])
    channels = [t[0] for t in data[k]['channels'][0,0][0]]
    d = data[k]['data'][0,0]
    
    assert len(channels) == d.shape[0]
    assert int(Fs*data_length_sec + 0.5) == d.shape[1],int(Fs*data_length_sec + 0.5)
    assert last_data_length_sec is None or last_data_length_sec == data_length_sec
    last_data_length_sec = data_length_sec
    assert last_Fs is None or last_Fs == Fs
    last_Fs = Fs
    assert last_channels is None or all(c1==c2 for c1,c2 in zip(last_channels, channels))
    last_channels = channels
    assert last_d_shape is None or last_d_shape == d.shape
    last_d_shape = d.shape
    
    if last_sequence is None:
        all_data = d.astype(float)
    elif last_sequence < sequence:
        if Fs > 400:
            data_offset = d[:,:1].mean(axis=-1) - all_data[:,-1:].mean(axis=-1)
            d -= data_offset.reshape(-1,1)
        all_data = np.hstack((all_data,d.astype(float)))
    else:
        if nchain == 0:
            break
        nchain -= 1
        all_data = d.astype(float)
    last_sequence = sequence
all_data = all_data.copy()
data_length_sec, sequence, Fs, all_data.shape


    

    
Out[122]:





(600, 1, 399.609756097561, (16, 239766))

In [123]:

    

print 'next segment', segment, last_sequence


    

    




next segment 1 1

In [124]:

    

int(all_data.shape[1]*399.61/Fs + 0.5)


    

    
Out[124]:





239766

In [125]:

    

all_data.shape


    

    
Out[125]:





(16, 239766)

In [126]:

    

5000/399.61


    

    
Out[126]:





12.51219939440955

In [127]:

    

Nchannels = 4 # d.shape[0]
for i in range(Nchannels):
    pl.subplot(Nchannels,1,i+1)
    pl.plot(all_data[i,:int(12.5*40000)])
    pl.title(channels[i])
pl.gcf().set_size_inches(14,8);


    

In [128]:

    

Nfs = int(Fs*data_length_sec+0.5)
N = int(399.61*data_length_sec)
Nfs, N


    

    
Out[128]:





(239766, 239766)

In [129]:

    

data = all_data[:,:Nfs]
data.shape, data[:10]


    

    
Out[129]:





((16, 239766), array([[   0.,  -14.,  -30., ...,  201.,  167.,  137.],
        [  36.,   59.,   74., ...,  171.,  170.,  175.],
        [  15.,   45.,   74., ...,  314.,  297.,  279.],
        ..., 
        [ -30.,   12.,   71., ...,  223.,  190.,  157.],
        [  -7.,  -45.,  -75., ...,  -62.,  -92., -136.],
        [ -96., -100.,  -76., ..., -331., -329., -362.]]))

In [130]:

    

print data.mean(axis=-1)
data = data - data.mean(axis=-1).reshape(-1,1)


    

    




[-0.36146493  0.30739137 -0.01033091  0.34801431 -0.28038171 -0.08463669
  0.06846676  0.06031714 -0.26574243  0.31994528  0.34421478 -0.31781404
 -0.27251987  0.0457738   0.03407489 -0.15897166]

In [131]:

    

e = np.zeros(data.shape[0])
q = int(Fs/60.+0.5)
for i in range(q):
    s = np.mean(data[:,(np.arange(i,Nfs,Fs/60.)).astype(int)], axis=-1)
    e += s*s
e = np.mean(np.sqrt(e/q))
e


    

    
Out[131]:





0.072554564957357931

In [132]:

    

from scipy.signal import resample, hann
if Fs > 410:
    if e > 2:
        notchwidth = 8.
        notchwidths = notchwidth/2./Fs
        notchfreq = 60.
        notchfreqs = notchfreq/Fs
        import numpy.fft
        fftfreq = np.abs(numpy.fft.fftfreq(data.shape[-1]))
        notch = np.ones(Nfs)
        for h in range(1,int((Fs/2.)/notchfreq+0.5)):
            notch[np.abs(fftfreq - h*notchfreqs) < notchwidths] = 0.
    else:
        notch = None
    data = resample(data, 239766, axis=-1, window=notch)
    F = 399.61
else:
    F = Fs


    

In [133]:

    

if notch:
    pl.plot(fftfreq, notch)
    pl.gca().set_xlim(0,4*notchfreq/Fs)


    

In [134]:

    

data.shape, data[:10]


    

    
Out[134]:





((16, 239766),
 array([[  3.61464928e-01,  -1.36385351e+01,  -2.96385351e+01, ...,
           2.01361465e+02,   1.67361465e+02,   1.37361465e+02],
        [  3.56926086e+01,   5.86926086e+01,   7.36926086e+01, ...,
           1.70692609e+02,   1.69692609e+02,   1.74692609e+02],
        [  1.50103309e+01,   4.50103309e+01,   7.40103309e+01, ...,
           3.14010331e+02,   2.97010331e+02,   2.79010331e+02],
        ..., 
        [ -3.00603171e+01,   1.19396829e+01,   7.09396829e+01, ...,
           2.22939683e+02,   1.89939683e+02,   1.56939683e+02],
        [ -6.73425757e+00,  -4.47342576e+01,  -7.47342576e+01, ...,
          -6.17342576e+01,  -9.17342576e+01,  -1.35734258e+02],
        [ -9.63199453e+01,  -1.00319945e+02,  -7.63199453e+01, ...,
          -3.31319945e+02,  -3.29319945e+02,  -3.62319945e+02]]))

In [135]:

    

WINDOW=int(F*1.) # 10 sec window
# power of 2 FFT
FFT_SIZE=1024
while FFT_SIZE < WINDOW*2:
    FFT_SIZE*=2

Pxx, freqs, bins, im = pl.specgram(data[channel,:].astype(float), NFFT=WINDOW, Fs=F, noverlap=WINDOW/2, pad_to=FFT_SIZE)# , detrend='linear'
#pl.gca().set_ylim((0,200.))
pl.gcf().set_size_inches(18,8);
#pl.tight_layout();
pl.autoscale(tight=True);
data.shape, WINDOW


    

    
Out[135]:





((16, 239766), 399)






    

In [136]:

    

pl.plot(freqs[20:1500],Pxx.mean(axis=-1)[20:1500])


    

    
Out[136]:





[<matplotlib.lines.Line2D at 0x117006d90>]






    

In [197]:

    

Pxx.shape


    

    
Out[197]:





(513, 1197)

In [ ]: