Questions for calculate power

đź’¬ Support & Discussions
1

Hi, everyone

I want to calculate power (not psd) of EEG signal. Should I calculate psd first, then multiply psd with freq? Or if there are other methods? Thanks.

2

Instantaneous power: raw.get_data() ** 2
Average power (per channel): np.mean(raw.get_data() ** 2, axis=-1)

Is one of those what you’re after?

1 Like
3

Maybe I didn’t describe my question clearly. I want to calculate the power of the alpha band, not the power spectrum density (which is power/frequency). Can I just multiply the frequency to psd?

4

You could use raw.compute_psd(fmin=8, fmax=12, or the equivalent with epochs.compute_psd, to restrict your estimation to the alpha band. Does that help?

5

Thanks. But that doesn’t solve my problem. I can use compute_psd to calculate the power spectrum density.

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

But that psd is density, what I want to calculate the the power. So, should just multiply the freqs with psds to get the power?

6

Hmm… maybe I’m missing something but this doesn’t seem like it would work. Would it make sense to multiply 8 (Hz) x the psd (at 8Hz)?

What I would do is sum the psd values between 8-12Hz. But it’s not clear to me the difference between psd and power that you are referring to.

have you searched either online or in the literature for a conversion calculation to go from psd to power? If you can find one, it should be easy enough to do with the output of compute_psd

7

Alternatively, you could calculate the power of the signal as @drammock demonstrates (amplitude squared), and then bandpass filter your data between 8-12Hz?

8

Hello, you can either follow this tutorial here:

Or use the Yasa package, which does all the work for you:

Best wishes,
Richard

9

I think typically when people report “alpha band power” it is mean power in the alpha band. So expanding on @scott-huberty 's suggestion you could do something like raw.compute_psd(fmin=8, fmax=12).get_data().mean(axis=-1) to get the mean power between 8 and 12 Hz. See @richard’s answer below, it should be integrated, not a mean.

Do not do that. Computing instantaneous signal power is a non-linear operation and will change the frequency content of the signal, so band-passing afterward will not have the effect you expect/want.

import numpy as np
import matplotlib.pyplot as plt
import mne
t = np.linspace(0, 3, 3001)
freqs = [10, 16, 21, 25, 28, 30]
y = np.sum([np.sin(2 * np.pi * f * t) for f in freqs], axis=0)
alpha_band = mne.filter.filter_data(y, sfreq=1000, lfreq=8, hfreq=12)
inst_power = y ** 2
not_alpha = mne.filter.filter_data(inst_power, sfreq=1000, lfreq=8, hfreq=12)

fig, ax = plt.subplots()
ax.plot(t, alpha_band)  # only the 10Hz component, as expected
ax.plot(t, not_alpha)

Figure_1

2 Likes
10

I think you actually need to integrate, this is also what Yasa does:

Given that this seems to be all but trivial, but a common use case, I think we should maybe make it easier for users to get the power of a given frequency band…

1 Like
11

My mistake, thanks for clarifying!

1 Like
12

That solved my problem! Thank you all!