In [9]:
%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 [149]:
target = 'Patient_2'
data_type = 'preictal' # preictal interictal, test

according to 140905-feature-importance the most important cannel is


In [150]:
channel = 3

In [151]:
import scipy.io
segment = 1 # start segment
all_data = last_sequence = last_data_length_sec = last_Fs = last_channels = last_d_shape = None
while True:
    fname = '../seizure-data/%s/%s_%s_segment_%04d.mat'%(target,target,data_type,segment)
    try:
        data = scipy.io.loadmat(fname)
    except:
        break
    k = '%s_segment_%d'%(data_type,segment)
    for k1 in data.keys():
        if k1 != k and data[k1]:
            print data[k1],
    print
    data_length_sec = data[k]['data_length_sec'][0,0][0,0]
    sequence = data[k]['sequence'][0,0][0,0]
    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]
    print segment,data_length_sec,sequence,Fs,d.shape
    
    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
        last_sequence = sequence
    elif last_sequence+1 == sequence:
        all_data = np.hstack((all_data,d-d.mean(axis=-1).reshape((-1,1))))
        last_sequence = sequence
    else:
        break
        # get the last sequence
        all_data = last_sequence = last_data_length_sec = last_Fs = last_channels = last_d_shape = None
    segment += 1
data_length_sec, sequence, Fs, all_data.shape


1.0 MATLAB 5.0 MAT-file, Platform: MACI64, Created on: Thu Aug 21 01:00:00 2014
1 600 1 5000.0 (24, 3000000)
1.0 MATLAB 5.0 MAT-file, Platform: MACI64, Created on: Thu Aug 21 01:00:00 2014
2 600 2 5000.0 (24, 3000000)
1.0 MATLAB 5.0 MAT-file, Platform: MACI64, Created on: Thu Aug 21 01:00:00 2014
3 600 3 5000.0 (24, 3000000)
1.0 MATLAB 5.0 MAT-file, Platform: MACI64, Created on: Thu Aug 21 01:00:00 2014
4 600 4 5000.0 (24, 3000000)
1.0 MATLAB 5.0 MAT-file, Platform: MACI64, Created on: Thu Aug 21 01:00:00 2014
5 600 5 5000.0 (24, 3000000)
1.0 MATLAB 5.0 MAT-file, Platform: MACI64, Created on: Thu Aug 21 01:00:00 2014
6 600 6 5000.0 (24, 3000000)
1.0 MATLAB 5.0 MAT-file, Platform: MACI64, Created on: Thu Aug 21 01:00:00 2014
7 600 1 5000.0 (24, 3000000)
Out[151]:
(600, 1, 5000.0, (24, 18000000))

In [152]:
d.mean(axis=-1)


Out[152]:
array([ 0.39343833, -0.17854733, -0.10525867, -0.19631467, -0.48245667,
        0.35039833, -0.20808667,  0.01547833,  0.388373  ,  0.244311  ,
       -0.331179  , -0.06954333,  0.313692  ,  0.056641  ,  0.13956433,
        0.13479633, -0.48531633, -0.34746833, -0.148391  , -0.04518367,
        0.12719   , -0.233886  , -0.00878833,  0.41231867])

In [153]:
print 'next segment', segment


next segment 7

In [154]:
N=d.shape[1]
N


Out[154]:
3000000

In [155]:
c=4
all_data[:,c*N:(c+1)*N].mean(axis=-1)


Out[155]:
array([  3.33511699e-09,  -1.14704797e-08,  -4.99722992e-09,
        -7.33716381e-09,  -1.04100460e-08,  -3.72182381e-09,
        -5.10724050e-09,  -4.33008874e-09,   1.45703233e-08,
        -8.32535712e-09,   1.48840326e-08,  -7.44664464e-09,
        -1.76307053e-08,   1.12835706e-08,  -4.68310021e-09,
        -9.33213636e-09,  -7.18143701e-09,  -7.65221068e-09,
         1.25766163e-08,  -1.55340651e-09,   3.14322368e-09,
         9.77262398e-09,   2.30453662e-09,   1.79559112e-08])

In [156]:
Fs


Out[156]:
5000.0

In [166]:
cutoff = Fs/100. # desired cutoff frequency of the filter, Hz
import scipy.signal
nyq = 0.5 * Fs
normal_cutoff = cutoff / nyq
order = 6
b, a = scipy.signal.butter(order, normal_cutoff, btype='low', analog=False)
# scipy.signal.lfilter(b, a, data)
h = scipy.signal.firwin(100, normal_cutoff)

In [167]:
Nbase = 8
Nchannels = 4 # d.shape[0]
W=300
W1 = max(int(Fs),W)
c=1
for i in range(Nchannels):
    pl.subplot(Nchannels,1,i+1)
    pl.plot(all_data[i+Nbase,c*N-W:c*N+W])
    # use filtfilt to get zero phase http://wiki.scipy.org/Cookbook/FiltFilt
    x1 = scipy.signal.filtfilt(b, a, all_data[i+Nbase,c*N-W1:c*N])
    x2 = scipy.signal.filtfilt(b, a, all_data[i+Nbase,c*N+W1:c*N:-1])[::-1]
    pl.plot(np.hstack((x1[-W:],x2[:W])))
#     x3 = scipy.signal.lfilter(h, 1, all_data[i+Nbase,c*N-W1:c*N+W1])
#     pl.plot(x3[W1-W:W1+W])

    pl.title(channels[i+Nbase])
pl.gcf().set_size_inches(14,8);



In [171]:
x2[10::-1]


Out[171]:
array([ 254.19976812,  254.19745329,  254.1954979 ,  254.19385798,
        254.19249314,  254.19136653,  254.19044469,  254.1896975 ,
        254.18909801,  254.18862229,  254.18824928])

In [173]:
x2[-10:]


Out[173]:
array([ 566.18886645,  566.07702448,  565.96704208,  565.85899896,
        565.7529759 ,  565.64905407,  565.54731453,  565.44783757,
        565.35070218,  565.25598556])

In [359]:
pl.plot(all_data[0,N-100:N+100]);
np.mean(np.abs(all_data[:,N-1]-all_data[:,N])/np.std(all_data[:,:N],axis=-1))


Out[359]:
3.285419563691018

In [325]:
c=3
pl.plot(all_data[0,c*N:(c+1)*N])


Out[325]:
[<matplotlib.lines.Line2D at 0x11aed1fd0>]

In [315]:
all_data.shape[1]/N


Out[315]:
6

In [247]:
import scipy.signal
W=400
pl.plot(scipy.signal.correlate(all_data[3,-W:],all_data[2,-W:], mode='full'))


Out[247]:
[<matplotlib.lines.Line2D at 0x11af8c2d0>]

In [248]:
NFFT=1024
Pxx, freqs, bins, im = pl.specgram(all_data[channel,:], NFFT=NFFT, pad_to=NFFT*2, Fs=Fs, noverlap=NFFT/2)
pl.gca().set_ylim((0,Fs/2.))
pl.show()



In [62]:
NFFT=1024
Pxx, freqs, bins, im = pl.gca().specgram(all_data[channel,:], NFFT=NFFT, pad_to=NFFT*2, Fs=Fs, noverlap=NFFT/2)
pl.gca().set_ylim((0,Fs/2.))
pl.show()



In [66]:
import matplotlib.mlab as mlab
Fc = 0.
Pxx, freqs, bins = mlab.specgram(all_data[channel,:], NFFT, Fs)

Z = 10. * np.log10(Pxx)
Z = np.flipud(Z)

xextent = 0, np.amax(bins)
xmin, xmax = xextent
freqs += Fc
extent = xmin, xmax, freqs[0], freqs[-1]
im = pl.gca().imshow(Z, None, extent=extent)
pl.gca().axis('auto')
pl.gca().set_ylim((0,Fs/2.))
pl.show()



In [262]:
import matplotlib.cbook as cbook

#This is a helper function that implements the commonality between the
#psd, csd, and spectrogram.  It is *NOT* meant to be used outside of mlab
def _spectral_helper(x, y, NFFT=256, Fs=2, detrend=mlab.detrend_none,
        window=mlab.window_hanning, noverlap=0, pad_to=None, sides='default',
        scale_by_freq=None):
    #The checks for if y is x are so that we can use the same function to
    #implement the core of psd(), csd(), and spectrogram() without doing
    #extra calculations.  We return the unaveraged Pxy, freqs, and t.
    same_data = y is x

    #Make sure we're dealing with a numpy array. If y and x were the same
    #object to start with, keep them that way
    x = np.asarray(x)
    if not same_data:
        y = np.asarray(y)
    else:
        y = x

    # zero pad x and y up to NFFT if they are shorter than NFFT
    if len(x)<NFFT:
        n = len(x)
        x = np.resize(x, (NFFT,))
        x[n:] = 0

    if not same_data and len(y)<NFFT:
        n = len(y)
        y = np.resize(y, (NFFT,))
        y[n:] = 0

    if pad_to is None:
        pad_to = NFFT

    if scale_by_freq is None:
        scale_by_freq = True

    # For real x, ignore the negative frequencies unless told otherwise
    if (sides == 'default' and np.iscomplexobj(x)) or sides == 'twosided':
        numFreqs = pad_to
        scaling_factor = 1.
    elif sides in ('default', 'onesided'):
        numFreqs = pad_to//2 + 1
        scaling_factor = 2.
    else:
        raise ValueError("sides must be one of: 'default', 'onesided', or "
            "'twosided'")

    if cbook.iterable(window):
        assert(len(window) == NFFT)
        windowVals = window
    else:
        windowVals = window(np.ones((NFFT,), x.dtype))

    step = NFFT - noverlap
    ind = np.arange(0, len(x) - NFFT + 1, step)
    n = len(ind)
    Pxy = np.zeros((numFreqs, n), np.complex_)

    # do the ffts of the slices
    for i in range(n):
        thisX = x[ind[i]:ind[i]+NFFT]
        thisX = windowVals * detrend(thisX)
        fx = np.fft.fft(thisX, n=pad_to)

        if same_data:
            fy = fx
        else:
            thisY = y[ind[i]:ind[i]+NFFT]
            thisY = windowVals * detrend(thisY)
            fy = np.fft.fft(thisY, n=pad_to)
        Pxy[:,i] = np.conjugate(fx[:numFreqs]) * fy[:numFreqs]

    # Scale the spectrum by the norm of the window to compensate for
    # windowing loss; see Bendat & Piersol Sec 11.5.2.
    Pxy /= (np.abs(windowVals)**2).sum()

    # Also include scaling factors for one-sided densities and dividing by the
    # sampling frequency, if desired. Scale everything, except the DC component
    # and the NFFT/2 component:
    Pxy[1:-1] *= scaling_factor

    # MATLAB divides by the sampling frequency so that density function
    # has units of dB/Hz and can be integrated by the plotted frequency
    # values. Perform the same scaling here.
    if scale_by_freq:
        Pxy /= Fs

    t = 1./Fs * (ind + NFFT / 2.)
    freqs = float(Fs) / pad_to * np.arange(numFreqs)

    if (np.iscomplexobj(x) and sides == 'default') or sides == 'twosided':
        # center the frequency range at zero
        freqs = np.concatenate((freqs[numFreqs//2:] - Fs, freqs[:numFreqs//2]))
        Pxy = np.concatenate((Pxy[numFreqs//2:, :], Pxy[:numFreqs//2, :]), 0)

    return Pxy, freqs, t

def myspecgram(x, y, NFFT=256, Fs=2, detrend=mlab.detrend_none, window=mlab.window_hanning,
        noverlap=128, pad_to=None, sides='default', scale_by_freq=None):
    """
    Compute a spectrogram of data in *x*.  Data are split into *NFFT*
    length segments and the PSD of each section is computed.  The
    windowing function *window* is applied to each segment, and the
    amount of overlap of each segment is specified with *noverlap*.

    If *x* is real (i.e. non-complex) only the spectrum of the positive
    frequencie is returned.  If *x* is complex then the complete
    spectrum is returned.

    %(PSD)s

      *noverlap*: integer
          The number of points of overlap between blocks.  The default value
          is 128.

    Returns a tuple (*Pxx*, *freqs*, *t*):

         - *Pxx*: 2-D array, columns are the periodograms of
           successive segments

         - *freqs*: 1-D array of frequencies corresponding to the rows
           in Pxx

         - *t*: 1-D array of times corresponding to midpoints of
           segments.

    .. seealso::

        :func:`psd`
            :func:`psd` differs in the default overlap; in returning
            the mean of the segment periodograms; and in not returning
            times.
    """
    assert(NFFT > noverlap)

    Pxx, freqs, t = _spectral_helper(x, y, NFFT, Fs, detrend, window,
        noverlap, pad_to, sides, scale_by_freq)
    Pxx = Pxx.real #Needed since helper implements generically

    return Pxx, freqs, t

import matplotlib.mlab as mlab
Fc = 0.
Pxx, freqs, bins = myspecgram(x[1,:], x[4,:], NFFT, Fs)

Z = 10. * np.log10(Pxx)
Z = np.flipud(Z)

xextent = 0, np.amax(bins)
xmin, xmax = xextent
freqs += Fc
extent = xmin, xmax, freqs[0], freqs[-1]
im = pl.gca().imshow(Z, None, extent=extent)
pl.gca().axis('auto')
pl.gca().set_ylim((0,Fs/2.))
pl.show()



In [256]:
x = all_data[:,-N:]
if Fs > 400.:
    x = scipy.signal.resample(x, int(N*400./Fs), axis=-1)
Cxy, f = mlab.cohere(x[channel,:], x[10,:], NFFT=400, Fs=Fs, pad_to=512)

In [257]:
pl.plot(f,Cxy)


Out[257]:
[<matplotlib.lines.Line2D at 0x113dd4590>]

In [237]:
N


Out[237]:
3000000

In [258]:
int(N*400./Fs)


Out[258]:
240000

In [259]:
Fs


Out[259]:
5000.0

In [260]:
N


Out[260]:
3000000

In [ ]: