Taking average of Power Spectral Densities

đź’¬ Support & Discussions
1

I have extracted 10 PSDs (using mne.Epochs.compute_PSD()) for 10 test subjects corresponding to one “Task 1” epoch. Now, how to produce a single average PSD that will represent all subjects?

2

You could call .get_data() for each subject and then use NumPy to average together the resulting arrays.

1 Like
3

Thanks, @giuliagennari. I figured that out but shall I average along axis=0?

4

I think you should consider using mne.grand_average()

Best wishes,
Richard

5

How will that work? mne.Epochs.compute_PSD() returns mne.time_frequency.EpochsSpectrum. But mne.grand_average() works only on mne.time_frequency.AverageTFR and mne.Evoked.

6

@WhosePotter As mentioned, you can just retrieve the underlying numpy array with the get_data method and then average as you want.

epochs = ...
spectrum = epochs.compute_psd()
data = spectrum.get_data()

@richard I had no idea this grand_average function existed and could be applied to AverageTFR objects. Let’s say you have 30 epochs, is it equivalent to:

  • compute the TFR 3 times on 10 epochs and then use grand_average
  • compute the TFR once on all 30 epochs
1 Like
7

@mscheltienne To repeat my concern, shall I take mean of the NumPy arrays along axis=0?

8

I missed this message. It depends on what you want to do and on what is the information of interest for you.
A couple of examples:

import numpy as np

from mne import create_info, make_fixed_length_epochs
from mne.io import RawArray
from scipy.integrate import simpson


info = create_info(["AF7", "Fp1", "Fpz"], 1024, "eeg")
data = np.random.randn(3, 10240)
raw = RawArray(data, info)
epochs = make_fixed_length_epochs(raw, preload=True)  # (10, 3, 1024)
spectrum = epochs.compute_psd()
data = spectrum.get_data()  # (10, 3, 513) - (n_epochs, n_channels, n_freqs)

# average across epochs
np.average(data, axis=0)  # (3, 513)

# average across channels
np.average(data, axis=1)  # (10, 513)

# bandpower by integrating on the frequency dimension, e.g. alpha band
data = spectrum.get_data(fmin=8, fmax=13)
fq_res = spectrum.freqs[1] - spectrum.freqs[0]
bp = simpson(data, dx=fq_res, axis=-1)  # (10, 3)
2 Likes
10

You are right, I got confused there, sorry!

11

Nah @richard - you saved me!
I computed PSD on the grand-average of epochs['desired_event'].average()! Solved my question without going through Scheltienne’s approach which made plotting difficult.