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
Out[151]:
In [152]:
d.mean(axis=-1)
Out[152]:
In [153]:
print 'next segment', segment
In [154]:
N=d.shape[1]
N
Out[154]:
In [155]:
c=4
all_data[:,c*N:(c+1)*N].mean(axis=-1)
Out[155]:
In [156]:
Fs
Out[156]:
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]:
In [173]:
x2[-10:]
Out[173]:
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]:
In [325]:
c=3
pl.plot(all_data[0,c*N:(c+1)*N])
Out[325]:
In [315]:
all_data.shape[1]/N
Out[315]:
In [247]:
import scipy.signal
W=400
pl.plot(scipy.signal.correlate(all_data[3,-W:],all_data[2,-W:], mode='full'))
Out[247]:
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]:
In [237]:
N
Out[237]:
In [258]:
int(N*400./Fs)
Out[258]:
In [259]:
Fs
Out[259]:
In [260]:
N
Out[260]:
In [ ]: