# Issue with EEG relative power estimation

**URL:** <https://mne.discourse.group/t/issue-with-eeg-relative-power-estimation/3174>\
**Category:** Support & Discussions\
**Tags:** preprocessing, eeg, epochs\
**Created:** [May 28, 2021, 11:40am UTC](https://mne.discourse.group/t/issue-with-eeg-relative-power-estimation/3174 "2021-05-28T11:40:56Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![eduardo](https://yyz2.discourse-cdn.com/free1/user_avatar/mne.discourse.group/eduardo/32/1157_2.png) [@eduardo](https://mne.discourse.group/u/eduardo)\
**Post date:** [May 28, 2021, 11:40am UTC](https://mne.discourse.group/t/issue-with-eeg-relative-power-estimation/3174/1 "2021-05-28T11:40:56Z")

</div>

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.](https://mne.tools/stable/auto_tutorials/clinical/60_sleep.html)

For my second attempt (`X_psd_2` in the code), I used `compute_pow_freq_bands` from [MNE features library](https://mne.tools/mne-features/generated/mne_features.univariate.compute_pow_freq_bands.html#mne_features.univariate.compute_pow_freq_bands), 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](https://global.discourse-cdn.com/free1/uploads/mne/original/1X/25c70877eac82718de54b1995c040a7a30d92e4a.png)

_Minimal working example_

```python
# 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()

```

---

<div class="post-metadata">

**Author:** ![eduardo](https://yyz2.discourse-cdn.com/free1/user_avatar/mne.discourse.group/eduardo/32/1157_2.png) [@eduardo](https://mne.discourse.group/u/eduardo)\
**Post date:** [May 28, 2021, 3:22pm UTC](https://mne.discourse.group/t/issue-with-eeg-relative-power-estimation/3174/2 "2021-05-28T15:22:31Z")

</div>

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:

```python
# 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](https://global.discourse-cdn.com/free1/uploads/mne/original/1X/ef6912d845448b537847facfc1e300a1d344bbcf.jpeg)