Looking at EEG power spectra¶

Here we examine a number of ways in which the frequency content of EEG data can be assessed.

We use the N170 dataset as an example, but the following can be used with minimal adjustments with any eeg or time series data.

Setup¶

In [2]:

import sys
eegnb_dir = "C:\\Users\\john_griffiths\\Code\\libraries_of_mine\\github\\eeg-notebooks"
sys.path.append(eegnb_dir)

In [1]:

# Generic imports
import os,sys,glob,numpy as np, pandas as pd

# we're currently in sandbox (i.e. test / dev area.)
# but we want the behaviour to be as if we're in one dir up; so move there now
os.chdir('../')

# eeg-notebooks utils
from utils import utils 

# mne functions
from mne import Epochs, find_events
from mne.time_frequency import psd_welch

# visualization stuff
%matplotlib inline 
from matplotlib import pyplot as plt
import seaborn as sns

# others
from scipy.signal import welch

Load the data¶

We will use session 1, subject 1 from the eeg notebooks example N170 dataset

In [2]:

subject = 1
session = 1
raw = utils.load_data('visual/N170', sfreq=256., 
                      subject_nb=subject, session_nb=session)

Creating RawArray with float64 data, n_channels=5, n_times=30732
    Range : 0 ... 30731 =      0.000 ...   120.043 secs
Ready.
Creating RawArray with float64 data, n_channels=5, n_times=30732
    Range : 0 ... 30731 =      0.000 ...   120.043 secs
Ready.
Creating RawArray with float64 data, n_channels=5, n_times=30732
    Range : 0 ... 30731 =      0.000 ...   120.043 secs
Ready.
Creating RawArray with float64 data, n_channels=5, n_times=30732
    Range : 0 ... 30731 =      0.000 ...   120.043 secs
Ready.
Creating RawArray with float64 data, n_channels=5, n_times=30720
    Range : 0 ... 30719 =      0.000 ...   119.996 secs
Ready.
Creating RawArray with float64 data, n_channels=5, n_times=30732
    Range : 0 ... 30731 =      0.000 ...   120.043 secs
Ready.

Calculate PSD for a raw (non-epoched) recording¶

This is the simplest option.

In [3]:

# compute the poweer spectral density (PSD) using 
# the MNE psd_welch function
# (this is simply a wrapper on scipy.signal.welch
#  that adds compatbility for MNE data types)

psd,freqs = psd_welch(raw)

# Place in a pandas dataframe for convenience

df_psd = pd.DataFrame(psd, columns=freqs).T
df_psd.columns = raw.ch_names[:4]
df_psd.index.names = ['freq']
df_psd.columns.names = ['chan']

Effective window size : 1.000 (s)

In [4]:

# (reminder: pandas dataframes are like a bit like excel spreadsheets. 
#  You can view their contents within notebooks with commands like the following: 

df_psd.iloc[:10,:]

Out[4]:

chan	TP9	AF7	AF8	TP10
freq
0.0	7.857495e-12	8.728770e-13	9.292256e-13	1.154006e-11
1.0	4.834027e-11	7.087248e-12	5.778365e-12	6.940097e-11
2.0	7.097305e-11	4.820972e-12	2.761128e-12	7.803431e-11
3.0	6.414010e-11	3.632046e-12	1.375676e-12	6.852216e-11
4.0	4.564189e-11	2.447174e-12	1.015633e-12	4.741013e-11
5.0	2.399177e-11	1.383455e-12	8.356244e-13	2.572205e-11
6.0	1.053377e-11	8.235806e-13	6.736721e-13	1.132447e-11
7.0	5.375513e-12	6.215222e-13	5.176218e-13	5.815875e-12
8.0	3.451319e-12	4.510599e-13	4.286293e-13	3.447529e-12
9.0	2.929786e-12	3.904274e-13	4.247347e-13	2.907828e-12

Plot the full power spectrum

In [6]:

fig, ax = plt.subplots(figsize=(8,4))
df_psd.plot(logy=True,ax=ax);

Double checking this agains the MNE plot_psd function, which should compute and plot the same numbers (with a few small differences in scaling and visualization):

In [7]:

raw.plot_psd();

Effective window size : 8.000 (s)

Remember, the 60Hz spike in the power spectrum is an artifact due to the pervasive mains electricity background noise.

Now, calculate the average power within a set of pre-defined canonical frequency bands

In [8]:

# These are the conventional EEG frequency band names and ranges

freqs = ['delta', 'theta', 'alpha', 'beta', 'lowgamma', 'midgamma']

freq_bands = dict(delta = [0.5,2],
                  theta = [4,8], 
                  alpha = [8,12], 
                  beta = [12,20],
                  lowgamma=[20,30],
                  midgamma=[30,50])

# Average the power within each of these bands

psd_fb = {}
for band_name, (rlow,rhigh) in freq_bands.items():
    psd_fb[band_name] = df_psd.loc[rlow:rhigh].mean(axis=0)

    
# Put in a pandas dataframe

df_psd_fb = pd.DataFrame(psd_fb).T.loc[freqs]
df_psd_fb.index.names = ['freq band']

Now plot the power in each frequency band for the four sensors:

In [9]:

fig, ax = plt.subplots(ncols=2, nrows=2,figsize=(12,8))

sens,thisax = 'AF7', ax[0][0]
df_psd_fb[sens].plot(kind='bar', ax=thisax,logy=True)
thisax.set_title(sens)

sens,thisax = 'AF8', ax[0][1]
df_psd_fb[sens].plot(kind='bar', ax=thisax,logy=True)
thisax.set_title(sens)

sens,thisax = 'TP9', ax[1][0]
df_psd_fb[sens].plot(kind='bar', ax=thisax,logy=True)
thisax.set_title(sens)

sens,thisax = 'TP10', ax[1][1]
df_psd_fb[sens].plot(kind='bar', ax=thisax,logy=True)
thisax.set_title(sens)


plt.tight_layout()

Complementary plot: power for each frequency band as a function of sensor:

In [10]:

fig, ax = plt.subplots(nrows=6, figsize=(8,20))


chans = ['AF7', 'AF8', 'TP9', 'TP10']

thisfreq,thisax = 'delta', ax[0]
df_psd_fb[chans].loc[thisfreq].plot(kind='bar',ax=thisax,logy=True)
thisax.set_title(thisfreq);

thisfreq,thisax = 'theta', ax[1]
df_psd_fb[chans].loc[thisfreq].plot(kind='bar',ax=thisax,logy=True)
thisax.set_title(thisfreq);

thisfreq,thisax = 'alpha', ax[2]
df_psd_fb[chans].loc[thisfreq].plot(kind='bar',ax=thisax,logy=True)
thisax.set_title(thisfreq);

thisfreq,thisax = 'beta', ax[3]
df_psd_fb[chans].loc[thisfreq].plot(kind='bar',ax=thisax,logy=True)
thisax.set_title(thisfreq);

thisfreq,thisax = 'lowgamma', ax[4]
df_psd_fb[chans].loc[thisfreq].plot(kind='bar',ax=thisax,logy=True)
thisax.set_title(thisfreq);

thisfreq,thisax = 'midgamma', ax[5]
df_psd_fb[chans].loc[thisfreq].plot(kind='bar',ax=thisax,logy=True)
thisax.set_title(thisfreq);


plt.tight_layout()