How to correctly compute the EEG Frequency Bands with Python?

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!

• Hello, Can I ask a question? Wouldn't the FFT1 function be called after you filtered each channel with the set of bandpass filters, rather than filtering the FFT1 result with the Butterworth filter? Convolution in the time-domain (i.e. filtering with the Butterworth filter and the filter function) versus multiplication in the frequency-domain. Thanks. May 2 '18 at 5:42

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]) &
(fft_freq <= eeg_bands[band]))
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") • great solution! good to know there is a simpler way to do it! Thanks a lot :D!!! Dec 11 '17 at 11:09
• how did you get the different bands to show up in different colors? Aug 20 '19 at 20:29
• also, what's the unit for the y axis? isn't uV anymore, isn't Hz anymore, what is it? Aug 21 '19 at 5:27
• fwiw i worked this into github.com/sixtyfive/tgam1/blob/master/script/python/… but it looks weird as there is always a lot of delta, even in a person who's engaged in programming. any comment welcome! Aug 21 '19 at 6:19
• @RainerVerteidiger - this is random data, not real acquired uV data. The y axis is the magnitude of the real component of the fft, meaned across the frequency band. The different colors show automatically using pandas bar plot. Aug 21 '19 at 21:33

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()))
beta1 = np.zeros((n, g.getNSamples()))
gamma1 = np.zeros((n, g.getNSamples()))
delta1 = np.zeros((n, g.getNSamples()))
sigma1 = np.zeros((n, g.getNSamples()))
theta1 = np.zeros((n, g.getNSamples()))
#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!

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!