Filtering EEG signals using SciPy¶

I am part of a team at the Illinois Geometry Lab on a project called Visual Cliffs, Virtual Reality and Movement Disorders. The goal of this project is to better understand human walking patterns in stressful situations.

Our setup is as follows:

the subject walks on a self-pacing treadmill,
while wearing an HTC Vive,
hooked up to an EEG.

Essentially, we let people walk through virtual environments that elicit anxiety (e.g. over a narrow bridge next to a steep cliff), and we measure their response. This semester, we are working on creating a feedback loop, where we change the virtual world depending on how anxious the user is.

I lead the development of the virtual reality software, and I'm not usually involved with EEG signal processing. Today, however, I wanted to give a very quick example of how you can filter an EEG signal to only get the relevant frequencies. Really, this is just an example of how to use the function scipy.signal.firwin. I originally wrote this to clarify the filtering for team members who were struggling with it, but maybe someone else finds it useful.

We'll start by importing some of the modules needed.

In [1]:

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
from math import pi

We'll be working with a signal that's sampled at 1000 Hz.

In [2]:

sampling_freq = 1000

Our sample is going to be 1.5 seconds long.

In [3]:

duration = 1.5

We set up the list of times.

In [4]:

t = np.arange(0.0, duration, 1/sampling_freq)

Here is the completely clean (noiseless) signal. It has a component at 10 Hz, and a component at 25 Hz.

In [5]:

y_clean = np.sin(2*pi*10*t) + 0.5 * np.sin(2*pi*25*t)

Now we add some noise to it. We'll add two sources of noise:

a sinusoidal signal at 60Hz (maybe there's some noise coming off power lines)
Gaussian white noise

In [6]:

y = y_clean + 0.1 * np.sin(2*pi*60*t) + 0.2 * np.random.normal(size=len(t))
plt.figure(figsize=(18,4))
plt.plot(t, y)
#plt.plot(t, y_clean)
plt.xlabel('time (s)')
plt.ylabel('signal')
plt.show()

We design a FIR filter using the window method. We'll pass in 400 as the first argument (we're using 400 taps). Higher numbers mean more precise filtering (the cutoff at our chosen frequencies will be steeper), but more processing and, in a real-time setting, increased latency. The latency when filtering using $n$ stops is about $n/2$ samples, so in this case we would get a latency of about 200 milliseconds. For our purposes, the precise choice of this first parameter is not very critical:

as far as processing power goes, a 400-tap filter is a piece of cake for any recent computer
a 200 millisecond delay is acceptable in our experiment
I suspect that a 400-tap filter is more than precise enough for our filtering

In [7]:

filter = signal.firwin(400, [0.01, 0.06], pass_zero=False)
plt.plot(filter)
plt.show()

The picture above shows the filter, but you can ignore the picture if you don't know how to interpret it.

Let's filter!

In [8]:

y2 = signal.convolve(y, filter, mode='same')

We compare three signals:

the noisy signal (semi-transparent blue)
the original clean signal (dashed green)
the filtered signal (solid orange)

In [9]:

plt.figure(figsize=(18,4))
plt.plot(t, y, alpha=0.4)
plt.plot(t, y_clean, '--', color='green')
plt.plot(t, y2)
plt.show()

Looks pretty good! Let's zoom in for a closer look.

In [10]:

plt.figure(figsize=(18,4))
plt.plot(t[:300], y[:300], alpha=0.4)
plt.plot(t[:300], y_clean[:300], '--', color='green')
plt.plot(t[:300], y2[:300])
plt.show()

Pretty close to the original signal!

And that's it. Very basic signal filtering using SciPy.