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
Next run the FFT algorithm. :::: callout-tip ## Task
plot the magnitude response (\(log10(abs(fft(x))))\)). ::::
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.
- 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
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 need to take the mirror-symmetry into account.
Highpass instead of lowpass
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”?
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 spike.
Impulse Response Function
To get a bit deeper understanding of what is going on, have a look at the fourier transform of a new impulse signal (e.g. 1:400 => 0. and 200 => 1.).
What do you observe?
Why would we see ringing if we put most of the coefficients to 0?
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()- Choose the channel “Pz”, plot the channel (previous HW it was “Cz”)
- 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.
- Plot the fourier spectrum after filtering again.
- Plot the channel again, did the filter work as indented?
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.