Band Power Continuous Data

MNE-Python version: 0.22.1
operating system: Windows 10

Hello,

I would like to calculate band power for resting-state data. I found a script from the sleep stage classification tutorial that is designed for epoch objects which generates the following error when applied to continuous data: “IndexError: too many indices for array: array is 2-dimensional, but 3 were indexed”

https://mne.tools/stable/auto_tutorials/clinical/60_sleep.html

# 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)

To adjust this so it accepts continuous data, do I just need to remove one of the colons in “psds_band=” as follows? It runs, but I want to make sure that this is correct/valid.

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))

Thank you for your time,
Dps

Hello,

like you noticed, Raw data has one less dimension than Epochs. So I would say that in principle, your adjusted indexing is correct. As for the rest of the code / algorithm, I’m not 100% certain, but it seems ok. Maybe @cbrnr or @Denis can provide a more solid judgment here?

Hi! It’s a bit difficult to see just by looking at the two screenshots instead of an entire script; but as far as I understand your situation you are right. If your psd variable is only two-dimensional, you will need to remove the middle comma/colon.

1 Like

I agree, this should be correct because psds has shape (n_channels, n_freqs) for continuous input data and (n_epochs, n_channels, n_freqs) for epochs.

1 Like

Hello,
Your responses are greatly appreciated. Here is the full code:

# Sampling frequency
sampling_freq = 500

# Import data
mat_data = read_mat('filepath')
data = np.array(mat_data['EEG']['data'])

#Info and raw objects
ch_names = mat_data['EEG']['chanlocs']['labels']
sampling_freq = 500
info = mne.create_info(ch_names, ch_types=['eeg']*64 , sfreq=sampling_freq)
raw = mne.io.RawArray(data,info)

# Define function
def eeg_power_band(epochs):
    """EEG relative power band feature extraction.

    This function takes an ``mne.Epochs`` object and creates EEG features based
    on relative power in specific frequency bands that are compatible with
    scikit-learn.

    Parameters
    ----------
    epochs : Epochs
        The data.

    Returns
    -------
    X : numpy array of shape [n_samples, 5]
        Transformed data.
    """
    # specific frequency bands
    FREQ_BANDS = {"delta": [0.5, 5.0],
                  "theta": [5.0, 8.0],
                  "alpha": [8.0, 13.0],
                  "sigma": [13.0, 16.0],
                  "beta": [16.0, 30.0]}

    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)

# Run calculation and average across electrodes 
power = eeg_power_band(raw)
avg_power = np.mean(power, axis=0).reshape(1,5)

A further question is that I see very different relationships between variables when running this script with or without psd normalization. I am quite new to this so I could be mistaken, but shouldn’t the scaling change but the relationships remain stable? Without modifying the code at all, I get a graph that looks like this:
image

After commenting out the psd normalization but keeping everything else the same, I see this:

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

image

Thank you

Yes I would say, it should retain the pattern. But it’s difficult to figure what’s going on without a reproducible example.