In [1]:
#sample frequency = ?
import pandas as pd
import numpy as np
df = pd.read_csv('/Users/John/Dropbox/LLU/ROP/aEEG/CrossEEGTest.txt',
parse_dates={'datetime': ['Date','Time']},
index_col='datetime',
usecols=['Date', 'Time', 'EEG'],
na_values=['0'],
infer_datetime_format='True'
)
In [2]:
df[:3]
Out[2]:
In [6]:
#asymmetrical filtering
#need a flat, band-pass filter with asymmetrical
#gain with a slope of 12 dB per decade and in the frequency range of 2-15
#so therefore use Butterworth filter, because its AKA as a maximally flat band pass filter
#documentation: scipy.signal.butter(N, Wn, btype='low', analog=False, output='ba')
# modified following code from BASS.py
from scipy.signal import butter, lfilter
# Sample rate and desired cutoff frequencies (in Hz).
fs = 100.0
lowcut = 2.0
highcut = 20.0
def butter_bandpass(lowcut, highcut, fs, order=5):
"""
Wrapper that put sorts parameters into the butterworth band pass function.
Parameters
----------
lowcut : float
lowerbond, Hz.
highcut : float
upperbound, Hz.
Returns
-------
a: type
varible for lfilter
b: type
varible for lfilter
References
----------
.. [1] http://wiki.scipy.org/Cookbook/ButterworthBandpass
.. [2] http://stackoverflow.com/questions/12093594/how-to-implement-band-pass-butterworth-filter-with-scipy-signal-butter
"""
nyq = 0.5 * fs
low = lowcut / nyq
high = highcut / nyq
b, a = butter(order, [low, high], btype='band')
return b, a
def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
"""
Description
Parameters
----------
param1 : type, shape (N,) optional
description.
param2 : type, shape (N,) optional
description.
Returns
-------
value : type, shape (N) optional
description.
Notes
-----
more note about usage and philosophy.
Examples
--------
?
References
----------
.. [1] http://wiki.scipy.org/Cookbook/ButterworthBandpass
.. [2] http://stackoverflow.com/questions/12093594/how-to-implement-band-pass-butterworth-filter-with-scipy-signal-butter
"""
b, a = butter_bandpass(lowcut, highcut, fs, order=order)
y = lfilter(b, a, data)
return y
#how do I know which order to use?
In [3]:
import numpy as np
from scipy.interpolate import UnivariateSpline
from scipy.signal import wiener, filtfilt, butter, gaussian, freqz
from scipy.ndimage import filters
import scipy.optimize as op
import matplotlib.pyplot as plt
def ssqe(sm, s, npts):
return np.sqrt(np.sum(np.power(s-sm,2)))/npts
#documentation: scipy.signal.butter(N, Wn, btype='low', analog=False, output='ba')
def testButterworth(nyf, x, y, s, npts):
b, a = butter(4, 1.5/nyf)
fl = filtfilt(b, a, y)
plt.plot(x,fl)
print "flerr", ssqe(fl, s, npts)
return fl
if __name__ == '__main__':
npts = 1024
end = 8
dt = end/float(npts)
nyf = 0.5/dt
sigma = 0.5
x = np.linspace(0,end,npts)
r = np.random.normal(scale = sigma, size=(npts))
s = np.sin(2*np.pi*x)#+np.sin(4*2*np.pi*x)
y = s + r
plt.plot(x,s)
plt.plot(x,y,ls='none',marker='.')
ga = testGauss(x, y, s, npts)
fl = testButterworth(nyf, x, y, s, npts)
wi = testWiener(x, y, s, npts)
sp = testSpline(x, y, s, npts)
plt.legend(['true','meas','gauss','iir','wie','spl'], loc='upper center')
plt.savefig("signalvsnoise.png")
plt.clf()
w = np.fft.fftfreq(npts, d=dt)
w = np.abs(w[:npts/2+1]) #only freqs for real fft
plotPowerSpectrum(s, w)
plotPowerSpectrum(y, w)
In [ ]:
In [ ]: