EEG Signal Processing in Python and Scipy.Signal (1): Spectrum Estimation, FIR Filter Design, Convolution and Windowing

by Forrest Sheng Bao http://fsbao.net

(Update 2012-12-22: People have been asking about the input format of the EEG data. Well, there is no input - because I didn't use any function. My code loads data from a text file (a TSV one), where rows correspond to sampling instant sequentially and columns correspond to channels. Therefore, you just need to ensure that the variable eeg is a 2-D list or numpy array such that any row is sublist/subarray representing the value from all channels at a sampling time. For example, if you have only two sampling instants for a 4-channel EEG, you would expect eeg to be like [[1,2,3,4],[5,6,7,8]].)

I am doing a take-home midterm test of a class I am taking. So, I decided to use Python to to it. Below is a code for one problem. It includes several frequency used functions in classical signal spectral analysis and FIR filter design. The problem itself is to design bandpass filters over alpha to theta bands and apply them onto a EEG series, and plot the time domain and frequency domain signal, as well as the frequency response of filters. I used window methods to design FIR bandpass filters. The filters coefficients are smoothed by a Kaiser window. Filtering is implemented by convolving original signal with coefficients of filters.

The result of running this code is given below. Black lines in the bottom plot is the amplitude property of different filters. As you can see, signals outside filter bands are kicked out respectively.

# This software is a free software that comes with NO WARRANTY
# You have the right to obtain, modify and redistribute this code 
# under the terms of  GNU General Public License ver. 3 or later
# Copyleft 2009 Forrest Sheng Bao http://fsbao.net

from numpy import *
from scipy.signal import *
from numpy.fft import * 
from matplotlib import *
from scipy import *
from pylab import *

f=open('./EEG1_1c31.txt','r')
chno = 16   # total number of channels
eeg = [] 

while True:
 testline = f.readline()
 if len(testline) == 0:
  break #EOF
 testline = testline.split()
 eeg.append([])
 for i in xrange(0,chno):
  eeg[-1].append(float(testline[i])) 

ch = 1 # particualr channel to study 
eeg = array(eeg)
y = eeg[:,ch]         # the signal, study channel 'ch'
L = len(y)            # signal length
fs = 500.0              # sampling rate
T = 1/fs                # sample time
t= linspace(1,L,L)*T   # time vector

f = fs*linspace(0,L/10,L/10)/L  # single side frequency vector, real frequency up to fs/2
Y = fft(y)

figure()
filtered = []
b= [] # store filter coefficient
cutoff = [0.5,4.0,7.0,12.0,30.0]

for band in xrange(0, len(cutoff)-1):
 wl = 2*cutoff[band]/fs*pi
 wh = 2*cutoff[band+1]/fs*pi
 M = 512      # Set number of weights as 128
 bn = zeros(M)
 
 for i in xrange(0,M):     # Generate bandpass weighting function
  n = i-  M/2       # Make symmetrical
  if n == 0:
   bn[i] = wh/pi - wl/pi;
  else:
   bn[i] = (sin(wh*n))/(pi*n) - (sin(wl*n))/(pi*n)   # Filter impulse response
 
 bn = bn*kaiser(M,5.2)  # apply Kaiser window, alpha= 5.2
 b.append(bn)
 
 [w,h]=freqz(bn,1)
 filtered.append(convolve(bn, y)) # filter the signal by convolving the signal with filter coefficients

figure(figsize=[16, 10])
subplot(2, 1, 1)
plot(y)
for i in xrange(0, len(filtered)):
  y_p = filtered[i]
  plot(y_p[ M/2:L+M/2])
axis('tight')
title('Time domain')
xlabel('Time (seconds)')

subplot(2, 1, 2)
plot(f,2*abs(Y[0:L/10]))
for i in xrange(0, len(filtered)):
  Y = filtered[i]
  Y = fft(Y [ M/2:L+M/2])
  plot(f,abs(Y[0:L/10]))
axis('tight')
legend(['original','delta band, 0-4 Hz','theta band, 4-7 Hz','alpha band, 7-12 Hz','beta band, 12-30 Hz'])

for i in xrange(0, len(filtered)):   # plot filter's frequency response
  H = abs(fft(b[i], L))
  H = H*1.2*(max(Y)/max(H))
  plot(f, 3*H[0:L/10], 'k')    
axis('tight')
title('Frequency domain')
xlabel('Frequency (Hz)')
subplots_adjust(left=0.04, bottom=0.04, right=0.99, top=0.97)
savefig('filtered.png')

Reference: Signal Processing Toolbox of Scipy, http://docs.scipy.org/doc/scipy/reference/signal.html