# Taking average of Power Spectral Densities

**URL:** <https://mne.discourse.group/t/taking-average-of-power-spectral-densities/5894>\
**Category:** Support & Discussions\
**Created:** [November 9, 2022, 10:37am UTC](https://mne.discourse.group/t/taking-average-of-power-spectral-densities/5894 "2022-11-09T10:37:02Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![WhosePolter](https://avatars.discourse-cdn.com/v4/letter/w/b4bc9f/32.png) [@WhosePolter](https://mne.discourse.group/u/WhosePolter)\
**Post date:** [November 9, 2022, 10:37am UTC](https://mne.discourse.group/t/taking-average-of-power-spectral-densities/5894/1 "2022-11-09T10:37:02Z")

</div>

I have extracted 10 PSDs (using [mne.Epochs.compute\_PSD()](https://mne.tools/stable/generated/mne.Epochs.html#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?

---

<div class="post-metadata">

**Author:** ![giuliagennari](https://yyz2.discourse-cdn.com/free1/user_avatar/mne.discourse.group/giuliagennari/32/1914_2.png) [@giuliagennari](https://mne.discourse.group/u/giuliagennari)\
**Post date:** [November 11, 2022, 10:08pm UTC](https://mne.discourse.group/t/taking-average-of-power-spectral-densities/5894/2 "2022-11-11T22:08:52Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![WhosePolter](https://avatars.discourse-cdn.com/v4/letter/w/b4bc9f/32.png) [@WhosePolter](https://mne.discourse.group/u/WhosePolter)\
**Post date:** [November 12, 2022, 5:30am UTC](https://mne.discourse.group/t/taking-average-of-power-spectral-densities/5894/3 "2022-11-12T05:30:35Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![richard](https://yyz2.discourse-cdn.com/free1/user_avatar/mne.discourse.group/richard/32/15_2.png) [@richard](https://mne.discourse.group/u/richard)\
**Post date:** [November 12, 2022, 7:52am UTC](https://mne.discourse.group/t/taking-average-of-power-spectral-densities/5894/4 "2022-11-12T07:52:08Z")

</div>

I think you should consider using [`mne.grand_average()`](https://mne.tools/stable/generated/mne.grand_average.html)

Best wishes,  
Richard

---

<div class="post-metadata">

**Author:** ![WhosePolter](https://avatars.discourse-cdn.com/v4/letter/w/b4bc9f/32.png) [@WhosePolter](https://mne.discourse.group/u/WhosePolter)\
**Post date:** [November 12, 2022, 8:10am UTC](https://mne.discourse.group/t/taking-average-of-power-spectral-densities/5894/5 "2022-11-12T08:10:35Z")

</div>

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`.

---

<div class="post-metadata">

**Author:** ![mscheltienne](https://yyz2.discourse-cdn.com/free1/user_avatar/mne.discourse.group/mscheltienne/32/827_2.png) [@mscheltienne](https://mne.discourse.group/u/mscheltienne)\
**Post date:** [November 12, 2022, 9:57am UTC](https://mne.discourse.group/t/taking-average-of-power-spectral-densities/5894/6 "2022-11-12T09:57:56Z")

</div>

@WhosePotter As mentioned, you can just retrieve the underlying numpy array with the [`get_data`](https://mne.tools/dev/generated/mne.time_frequency.EpochsSpectrum.html#mne.time_frequency.EpochsSpectrum.get_data) method and then average as you want.

```python
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

---

<div class="post-metadata">

**Author:** ![WhosePolter](https://avatars.discourse-cdn.com/v4/letter/w/b4bc9f/32.png) [@WhosePolter](https://mne.discourse.group/u/WhosePolter)\
**Post date:** [November 12, 2022, 10:03am UTC](https://mne.discourse.group/t/taking-average-of-power-spectral-densities/5894/7 "2022-11-12T10:03:11Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![mscheltienne](https://yyz2.discourse-cdn.com/free1/user_avatar/mne.discourse.group/mscheltienne/32/827_2.png) [@mscheltienne](https://mne.discourse.group/u/mscheltienne)\
**Post date:** [November 12, 2022, 10:49am UTC](https://mne.discourse.group/t/taking-average-of-power-spectral-densities/5894/8 "2022-11-12T10:49:27Z")

</div>

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:

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

```

---

<div class="post-metadata">

**Author:** ![richard](https://yyz2.discourse-cdn.com/free1/user_avatar/mne.discourse.group/richard/32/15_2.png) [@richard](https://mne.discourse.group/u/richard)\
**Post date:** [November 12, 2022, 8:14pm UTC](https://mne.discourse.group/t/taking-average-of-power-spectral-densities/5894/10 "2022-11-12T20:14:43Z")

</div>

> [@WhosePolter](#):
>
> But `mne.grand_average()` works only on `mne.time_frequency.AverageTFR` and `mne.Evoked`.

You are right, I got confused there, sorry!

---

<div class="post-metadata">

**Author:** ![WhosePolter](https://avatars.discourse-cdn.com/v4/letter/w/b4bc9f/32.png) [@WhosePolter](https://mne.discourse.group/u/WhosePolter)\
**Post date:** [November 13, 2022, 3:26pm UTC](https://mne.discourse.group/t/taking-average-of-power-spectral-densities/5894/11 "2022-11-13T15:26:19Z")

</div>

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.