In this file, we're going to lowpass, bandpass, and highpass filter some EEG signals. The signals we'll be using are samples provided with the MNE package. It is recommended, but not required, that you first install MME. You can do so using the instructions on the MNE homepage. The key step is to go to a command window, and type
conda env update --file environment.yml in order to create the
mne environment. If you don't want to install MNE, you can do the whole MP without it, by simply avoiding the blocks of code marked Requires MNE below.
The first thing you should do is to download and install the template, from https://courses.engr.illinois.edu/ece401/fa2020/ece401_20fall_mp3.zip. When you unzip it, you will find the following files, among others:
import os import numpy as np import matplotlib.figure import matplotlib.pyplot %matplotlib inline
Here's how the debugging works. Every time one of your functions from mp3.py is called, we will first use the code
importlib.reload(mp3) in order to reload the most recent version of mp3. That way, if a block of code fails, you can revise the file mp3.py, and try again.
import mp3 import importlib importlib.reload(mp3)
<module 'mp3' from '/Users/jhasegaw/Dropbox/mark/teaching/ece401/ece401labs/20fall/mp3/src/mp3.py'>
Here's how the data works. If you want to use MNE to visualize the data, you first need to install the
mne environment as described above. Then, every time you want to run the MNE code on this page, you'll need to do the following steps:
conda activate mne, and then
jupyter notebook. That should open up your Jupyter browser, just as if you had opened it from the GUI, except it will make the following steps possible (they wouldn't be possible if you used the GUI).
Kernelmenu. Under that menu, choose
Python [conda env:mne].
Now you should be able to run the following block of code without errors:
## REQUIRES MNE import mne,os sample_data_folder = mne.datasets.sample.data_path() sample_data_raw_file = os.path.join(sample_data_folder,'MEG','sample','sample_audvis_raw.fif') raw = mne.io.read_raw_fif(sample_data_raw_file, preload=True, verbose=False) data, times = raw.get_data(picks='eeg',return_times=True) channel_indices_by_type = mne.io.pick.channel_indices_by_type(raw.info) len(channel_indices_by_type['eeg']) info = mne.pick_info(raw.info, channel_indices_by_type['eeg']) print('According to the MNE Info, the Sampling Freq is %g'%(info['sfreq']))
According to the MNE Info, the Sampling Freq is 600.615
If you don't want to use MNE, you can run the following block of code to load the repackaged data. You should run the following block whether or not you're using MNE, because some of the blocks that follow will use this version of the data.
import h5py with h5py.File('data.hdf5','r') as f: onesecond = f['onesecond'][:] sfreq = f['sfreq'] solutions = h5py.File('solutions.hdf5','r') print('According to data.hdf5, sfreq is %g, nchannels is %d, nsamps is %d'% (sfreq, onesecond.shape, onesecond.shape))
According to data.hdf5, sfreq is 600.615, nchannels is 60, nsamps is 600
EEG data consists of 60 different voltage signals, recorded in parallel from different electrodes on the scalp. Here is a function that you can use to choose three of the signals, and plot both the signals and their power spectra:
def plot_waves_and_spectra(signal, sfreq, nrows=3): nchannels, nsamps = signal.shape t = np.linspace(0,nsamps/sfreq,nsamps,endpoint=False) f = np.linspace(0,sfreq,nsamps,endpoint=False) fig = matplotlib.figure.Figure(figsize=(14,6)) axs = fig.subplots(nrows,2) channels = np.arange(0,nchannels,int(nchannels/nrows)) for plotrow in range(nrows): channel = channels[plotrow] axs[plotrow,0].plot(t,signal[channel,:]) axs[plotrow,0].set_title('Signal, Channel %d'%(channel)) axs[plotrow,0].set_ylabel('Volts') axs[plotrow,0].set_xlabel('Seconds') axs[plotrow,1].plot(f,20*np.log10(1e-3+np.abs(np.fft.fft(signal[channel,:])))) axs[plotrow,1].set_title('Spectrum, Channel %d'%(channel)) axs[plotrow,1].set_xlabel('Hertz') axs[plotrow,1].set_ylabel('dB') return(fig)
If you have installed MNE, you can use the
plot_topomap function to see how the voltages are distributed over the scalp at each time step. Here's a function that plots the voltages as a function of location, for 25 different time steps:
def plot_topomaps(data, info): nchannels, nsamps = data.shape fig = matplotlib.figure.Figure(figsize=(10,8)) axs = fig.subplots(5,5) for row in range(5): for col in range(5): samp = int(nsamps*(row*5+col)/25) t = samp/info['sfreq'] mne.viz.plot_topomap(data[:,samp],info,axes=axs[row,col],show=False) axs[row,col].set_title('%dms'%(int(np.round(1000*t)))) return(fig)
Running it requires MNE. It requires that you have already run the MNE sample-loading code, up above, to load the
info object that contains the spatial locations of the electrodes.
# REQUIRES MNE plot_topomaps(onesecond,info)
<Figure size 432x288 with 0 Axes>
Now you're ready to run the MP. Each block of code, below, will first show you the distributed solution (distributed to you in the file
solutions.hdf5), and will then run the code in your own file
mp3.py, to test whether or not your code works.
def plot_digital_filter(h, sfreq): fig = matplotlib.figure.Figure(figsize=(14,6)) axs = fig.subplots(3,1) axs.plot(h) axs.set_title('Impulse response $h[n]$') H = np.abs(np.fft.fft(h)) omega = np.linspace(0,2*np.pi,len(h)) axs.plot(omega,H) axs.set_title('Magnitude response $|H(\omega)|$') f = omega * sfreq / (2*np.pi) axs.loglog(f[1:],1e-3+H[1:]) axs.set_title('Log-log plot of magnitude vs. frequency in Hertz') return(fig)
This function tests whether or not you can create a lowpass filter with a specified cutoff frequency, and with even length. The distributed solution contains a 200-sample lowpass filter with a cutoff of pi/4. Notice that the loglog plot doesn't quite know what to do with the aliased part of the filter, up near $\omega=2\pi$. It isn't really designed to handle the frequency responses of digital filters like this one.