 # Issue with EEG relative power estimation

Dear forum, I am trying to estimate the relative power in the five main EEG bands from an mne.Epochs object (minimal working example at the bottom of this post).

For my first attempt (`X_psd_1` in the code), I used the function `eeg_power_band` that can be found in this MNE tutorial on sleep stage classification.

For my second attempt (`X_psd_2` in the code), I used `compute_pow_freq_bands` from MNE features library, as it seemed a more straight-forward approach.

After this, I expected X_psd_1 and X_psd_2 to be very similar as both functions estimate the relative power per band. Nonetheless, when I represented the results, they looked very different (see figure below these lines). I would appreciate if someone could shed some light on the matter.

Minimal working example

``````# Modules
import mne
import os
import numpy as np
from matplotlib import pyplot as plt
from mne.time_frequency import psd_welch
from mne_features.feature_extraction import extract_features

# Load an mne.Epochs object sample
sample_data_folder = mne.datasets.sample.data_path()
sample_data_raw_file = os.path.join(sample_data_folder, 'MEG', 'sample','sample_audvis_raw.fif')
events = mne.find_events(raw, stim_channel='STI 014')
epochs = mne.Epochs(raw, events, tmin=-0.3, tmax=0.7, picks='eeg')

# Function from sleep stage classification tutorial
def eeg_power_band(epochs):

# Specific frequency bands
FREQ_BANDS = {"delta": [0.5, 4.5],
"theta": [4.5, 8.5],
"alpha": [8.5, 11.5],
"sigma": [11.5, 15.5],
"beta": [15.5, 30]}

psds, freqs = psd_welch(epochs, picks='eeg', fmin=0.5, fmax=30.)

# Normalize the PSDs
psds /= np.sum(psds, axis=-1, keepdims=True)

X = []
for fmin, fmax in FREQ_BANDS.values():
psds_band = psds[:, :, (freqs >= fmin) & (freqs < fmax)].mean(axis=-1)
X.append(psds_band.reshape(len(psds), -1))

return np.concatenate(X, axis=1)

# Estimate relative power using eeg_power_band
X_psd_1 = eeg_power_band(epochs)

# Estimate relative power using compute_energy_freq_bands through extract_features
bands   = {"delta": [0.5, 4.5], "theta": [4.5, 8.5], "alpha": [8.5, 11.5],
"sigma": [11.5, 15.5], "beta":  [15.5, 30]}
X_psd_2 = extract_features(X = epochs.get_data(),
sfreq = epochs.info['sfreq'],
selected_funcs = ['pow_freq_bands'],
funcs_params = {'pow_freq_bands__freq_bands': bands,
'pow_freq_bands__normalize' : True})
# Graph both PSDs
plt.figure()
plt.subplot(211), plt.imshow(X_psd_1),
plt.xlabel('Channels x bands'),
plt.ylabel('Epoch #'),
plt.title('PSD with eeg_power_band')
plt.subplot(212), plt.imshow(X_psd_2),
plt.xlabel('Channels x bands'),
plt.ylabel('Epoch #'),
plt.title('PSD with compute_pow_freq_bands')
plt.tight_layout()
``````

I think I found the solution by myself, so I leave it here in case it can be helpful for anyone.

• Firstly, `extract_features` returns the relative power across epochs (rows) in the following manner in the columns:
[ch1_band1, ch1_band2, …, ch1_band5, …, …, ch60_band1, …, ch60_band5].

• On the other hand, `eeg_power_band` outputs the relative power across epochs (rows) in the following fashion in the columns:
[ch1_band1, ch2_band1, …, ch60_band1, …, ch1_band5, ch2_band_5, …, ch60_band5]

So the issue was caused by the order of the columns returned by both functions. I believe the approach followed by `eeg_power_band` is more intuitive, thus, to solve the issue you can do the following:

``````# Reorder the index of X_psd_2 to show all channels in a band together
new_order = [list(range(idx, len(bands) * len(epochs.info['chs']), len(bands))) for idx in range(len(bands))]

# Flatten list
new_order = [item for sublist in new_order for item in sublist]

# Graph both PSDs
plt.figure()
plt.subplot(211), plt.imshow(X_psd_1),
plt.xlabel('Channels x bands'),
plt.ylabel('Epoch #'),
plt.title('PSD with eeg_power_band')
plt.subplot(212), plt.imshow(X_psd_2[:, new_order]),
plt.xlabel('Channels x bands'),
plt.ylabel('Epoch #'),
plt.title('PSD with compute_pow_freq_bands REORDERED')
plt.tight_layout()
``````

And the new plot looks like this:

1 Like