Obtaining average band power of EEG signals using mne Python library

Question

I'm looking for a way to obtain the gamma band's average frequency of a channel in an EEG signal from an edf file and I am unable to figure out how to do so. I checked various sources online and I found out in order to do that I need to obtain the PSD(Power Spectral Density) from the signal using the Welch's method but I was unable to find a way to do that using the mne library. All I have been able to accomplish so far has been attached below. I'd be thankful for any kind of help.

import mne
file = "H S1 EC.edf"
data = mne.io.read_raw_edf(file)
raw_data = data.get_data()
info = data.info
channels = data.ch_names

Answer 1

To get the power spectrum of your signal you can:

import mne

file = "H S1 EC.edf"
data = mne.io.read_raw_edf(file, preload=True)
psds, freqs = mne.time_frequency.psd_welch(data)

This gives you the spectral power for all frequencies, so to get average gamma power you would have to average a slice of the psds array that corresponds to your frequency band of interest (the frequencies corresponding to power values in psds can be found in freqs array).

However, you seem to be interested not in the average gamma power, but gamma frequency, and this is much trickier. Gamma power is generally weak in EEG and the signal in this frequency range can be dominated by muscle artifacts (and possible microsaccade artifacts) - beacause of this I doubt you would be able identify a reliable gamma peak (one that is not an artifact). I would suggest you to take a look at the "Dissociating neuronal gamma-band activity from cranial and ocular muscle activity in EEG" 2013 paper by Joerg Hipp and Markus Siegel, they give some good advice on preprocessing and analysing the data to get good estimates of neuronal gamma power.

Answer 2

To obtain the PSD using Welch's method in MNE, you can use the psd_welch() function. Here's an example code snippet to get you started:

import mne

file = "H S1 EC.edf"
data = mne.io.read_raw_edf(file)

# Define the frequency range of interest (gamma band)
freq_range = (30, 100)

# Obtain the PSD using Welch's method
psds, freqs = mne.time_frequency.psd_welch(data, fmin=freq_range[0], fmax=freq_range[1], n_fft=1024)

# Find the index of the frequency bin with the highest power
gamma_index = psds.mean(axis=0).argmax()

# Obtain the frequency corresponding to the highest power bin
gamma_freq = freqs[gamma_index]

# Print the average gamma frequency for each channel
for i, channel in enumerate(data.ch_names):
    print("Channel {}: {:.2f} Hz".format(channel, psds[i,gamma_index]))

In this code, we first define the frequency range of interest as the gamma band (30-100 Hz). We then use the psd_welch() function to obtain the PSD for this frequency range using the Welch's method with a window size of 1024 samples. The resulting PSDs are stored in the psds variable, and the corresponding frequencies are stored in the freqs variable.

We then find the index of the frequency bin with the highest power for each channel, and obtain the corresponding frequency. Finally, we print out the average gamma frequency for each channel.

Note that the psd_welch() function returns the PSD as a 2D array, where each row corresponds to a channel and each column corresponds to a frequency bin. To obtain the average gamma frequency for each channel, we first take the mean across the frequency axis (axis=1) to obtain the average power for each channel, and then find the index of the highest power bin for each channel.

Alternatively,To obtain the PSD of an EEG signal using the Welch's method and calculate the average frequency of the gamma band, you can use the YASA library. Here's some sample code that shows how to do this:

import mne
import yasa

# Load EEG data
file = "H S1 EC.edf"
data = mne.io.read_raw_edf(file)
raw_data = data.get_data()
sf = data.info['sfreq']  # Sampling frequency

# Compute the power spectral density using Welch's method
psd, freqs = yasa.psd_welch(raw_data, sf=sf, fmax=100)

# Calculate the average frequency of the gamma band (30-100 Hz)
gamma_freqs = freqs[(freqs >= 30) & (freqs <= 100)]
gamma_psd = psd[:, (freqs >= 30) & (freqs <= 100)]
avg_gamma_freq = (gamma_freqs * gamma_psd).sum() / gamma_psd.sum()

print("Average gamma frequency:", avg_gamma_freq)

This code loads the EEG data using MNE, computes the PSD using YASA's psd_welch() function, and then calculates the average frequency of the gamma band by taking the weighted average of the frequency bins in the gamma range. Note that the psd_welch() function takes care of windowing and overlapping the EEG signal, which are important steps in obtaining an accurate PSD estimate.

Additionally you can refer to my github code on how to Compute the average bandpower of an EEG signal