Exercise 3: Filtering

Simulate a simple EEG signal

We start out with a simple signal:

from numpy import cos, sin, pi,  arange

sample_rate = 100.0
nsamples = 400 # fixed at 11.11.2020
t = arange(nsamples) / sample_rate
x = cos(2*pi*0.5*t) + 0.2*sin(2*pi*2.5*t+0.1) + \
        0.2*sin(2*pi*15.3*t) + 0.1*sin(2*pi*16.7*t + 0.1) + \
            0.1*sin(2*pi*23.45*t+.8)

Task: Plot the signal against time.

Transform to fourier space

Task
  • run the fourier analysis via fft.
  • plot the full log10-magnitude
Note

If you plot e.g. the magnitude you will notice it is mirrorsymmetric around the middle! This is because FFT is actually defined for complex signals. Because we only have real signals, practically half of the FFT is redundant. But this also means, all of our manual-filters in the next steps need to take the mirror-symmetry into account. For this reason, typically a rfft function exists, returning only the real part.

A simple filter

A very simple frequency filter is to set the unwanted amplitude of the fourier components to something close to zero. We have to be a tad careful not to remove the phase information though.

Your tasks are therefore
  • split the complex fourier result into angle and magnitude
  • set the respective magnitudes to zero (we start with a lowpass filter: magnitude[30:370] = 0)
  • combine angle and magnitude back to a complex fourier coefficient (\(m*e^{1j*ang}\))
  • apply the inverse FFT
  • plot the signal with what you started out

Highpass instead of lowpass

Task

Repeat the steps from above, but this time, remove the low frequency components

What happens to the frequency and time response if we add “artefacts”?

Let’s see how the FFT reacts to steps and spikes. A step-function could look like this:

Task

Add a DC-offset (a step-function) starting from x[200:] and investigate the fourier space. Filter it again (low or high pass) and transfer it back to the time domain and investigate the signal around the step.

Impulse Response Function

To get a bit deeper understanding of what is going on, have a look at the fourier transform of a impulse signal (e.g. 1:400 => 0. and 200 => 1.).

Question
  • What do you observe?
  • Why did we see ringing in the previous filtered step-function, if when filtering we put coefficients to 0 - after all, we are not adding something?

Filtering EEG data

Usually we would not built our own filters - but understanding the properties is still very important! We will load our data from last time again:

from mne_bids import (BIDSPath,read_raw_bids)
import mne_bids
import importlib
import ccs_eeg_utils

bids_root = "../local/bids"
subject_id = '030'


bids_path = BIDSPath(subject=subject_id,task="P3",session="P3",
                     datatype='eeg', suffix='eeg',
                     root=bids_root)
raw = read_raw_bids(bids_path)
ccs_eeg_utils.read_annotations_core(bids_path,raw)
raw.load_data()
Task
  • Choose the channel “Pz”, plot the channel
  • Plot the fourier space using raw.plot_psd()

Now we filter using raw.filter(), specify a highpass of 0.5Hz and a lowpass of 50Hz.

Task
  • Plot the fourier spectrum after filtering again.
  • Plot the channel again, did the filter work as indented?
Bonus-task

Bonus: comparing over-filtered ERPs

If you want, you can compare the ERP with and without filtering. You can also use “invalid” filter settings - HP up to 2-5Hz, lowpass until 10-20Hz. I say invalid here, because usually with such ranges, you would filter out results that you are actually interested in.

Bonus: Electrical Artefacts

Instead of notch filtering 50/60Hz artefacts, one can also try to regress it out in smarter ways. A good tool for this is Zap_Line with a python implementation here: https://github.com/nbara/python-meegkit/. There are also several robust detrending methods, which could potentially replace highpass filters in the future. But more work needs to be invested to see how results compare. These methods are not (yet) common.