Issue with EEG relative power estimation

💬 Support & Discussions
, ,
1

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.

Figure_1
Figure_1916×693 410 KB

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')
raw = mne.io.read_raw_fif(sample_data_raw_file, verbose=False).crop(tmax=60)
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()
2

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:

Figure_2
Figure_2825×700 147 KB

1 Like