How to correctly compute the EEG Frequency Bands with Python?

Question

so I am trying to compute the EEG(25 channels, 512 sampling rate, 248832/channel) bands (alpha, beta, gamma, etc.) with Python. I managed to do so by: firstly filtering the signal with a butterworth filter that looks like this:

def butter_bandpass_filter(data, lowcut, highcut, fs, order=2):
    nyq = 0.5 * fs
    low = lowcut /nyq
    high = highcut/nyq
    b, a = butter(order, [low, high], btype='band')
    #print(b,a)
    y = lfilter(b, a, data)
    return y

I used the filter with the low set on 0.1 and high on 80. Than I compute the fft of the signal and store it in fft1, on which I use again the butterwort filter to extract the frequencies of each band and it looks like this:

for i in np.arange(n):
    alpha1 = butter_bandpass_filter(fft1[i, :], 8.1, 12.0, 256)
    beta1 = butter_bandpass_filter(fft1[i, :], 16.0, 36.0, 256)
    gamma1 = butter_bandpass_filter(fft1[i, :], 36.1, 80, 256)
    delta1 = butter_bandpass_filter(fft1[i, :], 0.0, 4.0, 256)
    sigma1 = butter_bandpass_filter(fft1[i, :], 12.1, 16.0, 256)
    theta1 = butter_bandpass_filter(fft1[i, :], 4.1, 8.0, 256)
    sumalpha1 = sum(abs(alpha1))
    sumbeta1 = sum(abs(beta1))
    sumgamma1 = sum(abs(gamma1))
    sumdelta1 = sum(abs(delta1))
    sumsigma1 = sum(abs(sigma1))
    sumtheta1 = sum(abs(theta1))
    objects = [sumalpha1, sumbeta1, sumgamma1, sumdelta1, sumsigma1, sumtheta1]
    N = len(objects)
    ra = range(N)
    plt.title(signal_labels[i])
    plt.autoscale
    somestuffneeded = np.arange(6)
    ticks = ['alpha','beta','gamma','delta','sigma','theta']
    plt.xticks(somestuffneeded, ticks)
    plt.bar(ra, objects)
    plt.show()

The problem is that I get a high value for gamma( which represents a high cognitive activity, which the person clearly didn't have). The result for the first channel is:

Can anyone point out a better solution/ has any ideas about what is wrong/ can tell me what is wrong? Am I using the filter incorect? Thanks!

Answer 1

I believe there is a much simpler way to do this with numpy.fft.rfft and numpy.fft.rfftfreq.

In the below example, I have two seconds of random data between 0.0 and 100.0 sampled at 512 Hz. I'm plotting the band amplitude below, but if you wanted power, it would just be the square of the values.

I think the reason you saw overly large values for Gamma, was that your gamma range is much larger than the others, and you're taking the sum of all the values in that range. This is fine if you are going to compare one EEG set's gamma value to another (relative amp/power), but not if you're just looking at one set and want to compare gamma to say alpha.

If you just want to know what EEG state the person is in, you could use np.max instead of np.mean. Then whichever band has the max will tell you about the person's state.

import numpy as np

fs = 512                                # Sampling rate (512 Hz)
data = np.random.uniform(0, 100, 1024)  # 2 sec of data b/w 0.0-100.0

# Get real amplitudes of FFT (only in postive frequencies)
fft_vals = np.absolute(np.fft.rfft(data))

# Get frequencies for amplitudes in Hz
fft_freq = np.fft.rfftfreq(len(data), 1.0/fs)

# Define EEG bands
eeg_bands = {'Delta': (0, 4),
             'Theta': (4, 8),
             'Alpha': (8, 12),
             'Beta': (12, 30),
             'Gamma': (30, 45)}

# Take the mean of the fft amplitude for each EEG band
eeg_band_fft = dict()
for band in eeg_bands:  
    freq_ix = np.where((fft_freq >= eeg_bands[band][0]) & 
                       (fft_freq <= eeg_bands[band][1]))[0]
    eeg_band_fft[band] = np.mean(fft_vals[freq_ix])

# Plot the data (using pandas here cause it's easy)
import pandas as pd
df = pd.DataFrame(columns=['band', 'val'])
df['band'] = eeg_bands.keys()
df['val'] = [eeg_band_fft[band] for band in eeg_bands]
ax = df.plot.bar(x='band', y='val', legend=False)
ax.set_xlabel("EEG band")
ax.set_ylabel("Mean band Amplitude")

Answer 2

Ok, so for those interested, I've computed the frequency bands of an eeg by using the butterworth filter described in the problem description. During the eeg analysis class I came to the conclusion that the frequency bands were computed from the fft of the eeg which was not enough because the fft should have been multiplied with its conjugate! so here is the code in python which computes the total power, the relative and the absolute frequency bands.

alpha1 = np.zeros((n, g.getNSamples()[0]))
beta1 = np.zeros((n, g.getNSamples()[0]))
gamma1 = np.zeros((n, g.getNSamples()[0]))
delta1 = np.zeros((n, g.getNSamples()[0]))
sigma1 = np.zeros((n, g.getNSamples()[0]))
theta1 = np.zeros((n, g.getNSamples()[0]))
#I only used the channels with the following indexes
for i in [1,5,10,15,18]:
    #note that fft1 is precomputed in another part of the code i used
    ccFFT1 = fft1[i,:]*np.conjugate(fft1[i,:]) 
    #the butterworth filter is used to compute each frequency band according to its featuring freq.
    alpha1 = butter_bandpass_filter(ccFFT1[:], 8.1, 12.0, 512)
    beta1 = butter_bandpass_filter(ccFFT1[:], 16.0, 36.0, 512)
    gamma1 = butter_bandpass_filter(ccFFT1[:], 36.1, 80, 512)
    delta1 = butter_bandpass_filter(ccFFT1[:], 0.0, 4.0, 512)
    sigma1 = butter_bandpass_filter(ccFFT1[:], 12.1, 16.0, 512)
    theta1 = butter_bandpass_filter(ccFFT1[:], 4.1, 8.0, 512)
    #this sums frequencies for each band 
    sumalpha1 = sum(abs(alpha1))
    sumbeta1 = sum(abs(beta1))
    sumgamma1 = sum(abs(gamma1))
    sumdelta1 = sum(abs(delta1))
    sumsigma1 = sum(abs(sigma1))
    sumtheta1 = sum(abs(theta1))
    #compute the total power
    totalPower = sumalpha1+sumbeta1+sumgamma1+sumdelta1+sumsigma1+sumtheta1
    #storing them in objects to plot
    #this computes the relativ power!
    objects = [sumalpha1/totalPower, sumbeta1/totalPower, sumgamma1/totalPower, sumdelta1/totalPower, sumsigma1/totalPower, sumtheta1/totalPower]
    #to plot the absolute power, don't divide each sum through the total power
    N = len(objects)
    ra = range(N)
    plt.title(str(signal_labels[i])+" relativ")
    plt.autoscale
    somestuffneeded = np.arange(6)
    ticks = ['alpha','beta','gamma','delta','sigma','theta']
    plt.xticks(somestuffneeded, ticks)
    plt.bar(ra, objects)
    plt.show()

Hope this helps! I am glad to answer questions if there are any regarding this task of computing the frequency bands with python!

Answer 3

I would use scipy.signal.welch for that:

def calc_bands_power(x, dt, bands):
    from scipy.signal import welch
    f, psd = welch(x, fs=1. / dt)
    power = {band: np.mean(psd[np.where((f >= lf) & (f <= hf))]) for band, (lf, hf) in bands.items()}
    return power

Also, mne-python is a great package for EEG/MEG analysis, it's worth taking a look!