Band Power Continuous Data

Support & Discussions
1

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

2

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?

3

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
4

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
5

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

6

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

7

Hi,

just for clarification, while estimating absolute vs relative power.

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

  1. I imagine that the normalization above is the key when you want to estimate the relative power vs absolute power. Are you agree that without such normalization the script would later compute the absolute power?

  2. As shown in the image below, do you think is correct to report the unit measure of relative power as “PSD (%)” while for the absolute power as “PSD (dB/Hz)”?

  3. Is generally relative power a more common measure, compared to absolute power, while representing the power in a frequency band?

From this article

image
image597×533 43.6 KB

Thank you

8

First of all, relative PSD should not change the shape of the absolute PSD, because you are merely dividing by a constant. This is also evident from the cartoon you posted.

I don’t know what the three bars correspond to in your original figure, but I suspect it’s probably different conditions in the data. Now when you normalize all three conditions with their individual total power, in general you will get three different scaling factors. That’s why the second bar chart looks different from the first one, but if you only consider one color, you can see that the power changes are the same as in the absolute case.

Regarding which measure is more common, relative band power changes get rid of individual differences between participants, which is probably what you want in most cases.

9

thanks @cbrnr

regarding the unit measure, looking at the figure could I state the following?

  • Power represents here the “absolute power” and the unit measure of the figure is “V ** 2 / Hz (dB)” and not “uV ** 2 / Hz (dB)”.

This doubt because if I plot the PSD ,using for instance the line below, the resulting figure will show
“uV ** 2 / Hz (dB)”:

raw_filtered_new.plot_psd(area_mode='range', tmax=10.0, average=True)

but once I get data from the spectrum, using for instance

psds_noNorm, freqs = spectrum.get_data(return_freqs=True)

data shows really small numbers, so I believe that it is “V ** 2 / Hz (dB)” and not “uV ** 2 / Hz (dB)”.

image

thanks

10

@cbrnr

psds_noNorm, freqs = spectrum.get_data(return_freqs=True)
or maybe absolute power from MNE is simply (V ** 2) instead of (uV ** 2)

11

Indeed it looks like these numbers are in \frac{\text{V}^2}{\text{Hz}}. Please note that this is the unit and not dB, which you wrote in parentheses. It’s either or. Sometimes, “absolute” power is measured in dB, but this is also relative to some reference power P_0 according to \text{PSD}=10 \cdot \log \frac{P}{P_0}.

MNE always uses V for voltages, not µV.

2 Likes
12

@cbrnr just for curiosity:

  1. It is frequent to find the unit measure of power as µV** 2 instead of µV** 2/Hz. Is this happening only because it is obvious that “/Hz” is included and for quick writing it is not specified? Here an example.

  2. While getting data from the spectrum, for example using,
    psds_noNorm, freqs = spectrum.get_data(return_freqs=True)
    is it possible to get data in uV instead of V? Sometimes I’ve seen that MNE allows specifying " units=‘uV’ ", but this parameter doesn’t work with “spectrum.get_data”

Thanks

13

  1. The correct unit includes “/Hz”, because this is a power spectral density. The area under the curve corresponds to power (e.g. when you integrate between two frequencies, you get µV²).

  2. It looks like that’s not possible, but you can always multiply by 10^{12} to convert to µV².

1 Like