Calculate spectral density

Hello everyone,

I concucted an EEG study with about 30 participants about VR creativity and now I am struggling Iâ€™m struggling with computing the average power of each frequency band.

I followed all tutorials that I could find online but still my values seem to be incorrect. Is it just a scaling problem? Or am I doing something wrong with the preprocessing?

``````import os
import numpy as np
import pandas as pd
import mne
import matplotlib.pyplot as plt
from mne.preprocessing import (ICA, create_eog_epochs, create_ecg_epochs,
corrmap)
import mat73
from mne.time_frequency import psd_multitaper
from numpy import pi as PI
from scipy import signal
from scipy.integrate import simps
import seaborn as sns
from scipy.signal import welch
from scipy.integrate import simps
from mne.time_frequency import psd_array_multitaper
from scipy.signal import spectrogram, welch

file = "C:\\sub13-2D.mat"

# extract values
data = np.array(dat['y'])

# create info
sf = 250
n_channels = 29
info = mne.create_info(ch_names=ch_names, sfreq=sf ,ch_types='eeg')
ch_types = ['misc','misc','eeg','eeg','eeg','eeg','eeg','eeg','eeg','eeg','eeg','eeg','eeg','eeg','eeg','eeg','eeg','eeg','misc','misc','misc','misc','misc','misc','misc','misc','misc','misc','misc']
ch_map = dict(zip(ch_names, ch_types))

# create raw
raw = mne.io.RawArray(data,info)
raw.set_channel_types(ch_map)

# set channel locations for montage
biosemi_montage = mne.channels.make_standard_montage('biosemi16')
raw.set_montage(biosemi_montage)

filt_raw = raw.copy().filter(l_freq=1., h_freq=60.,picks='eeg')

# ICA
ica = ICA(n_components=n_channels, random_state=97)
ica.fit(filt_raw,picks='eeg')

# plot the ICA to check components
ica.plot_sources(filt_raw, show_scrollbars=True)
ica.plot_components(outlines='skirt')

ica.exclude = list(comps)
ica.apply(filt_raw)

# call the function on the preprocessed data
eeg = eeg_power_band(filt_raw)

# this is the first function I could find:
def bandpower(data, sf, band):
# frequencies of interest
fmin = 1.
fmax = 60.

band = np.asarray(band)
low, high = band
win = 4 * sf

psd, freqs = psd_array_multitaper(data, sf, adaptive=True, fmin=fmin, fmax=fmax,
normalization='full', verbose=0)

psd = psd.T

# psd = 10. * np.log10(psd)

# Frequency resolution
freq_res = freqs[1] - freqs[0]

# Find index of band in frequency vector
idx_band = np.logical_and(freqs >= low, freqs <= high)

# Integral approximation of the spectrum using parabola (Simpson's rule)
bp = simps(psd[idx_band], dx=freq_res)
return bp

def eeg_power_band(eeg):
# frequencies of interest
fmin = 1.
fmax = 60.

# 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],
"gamma": [30,60]}

psds, freqs =  mne.time_frequency.psd_multitaper(eeg,
picks='eeg',
fmin=fmin, fmax=fmax)
# Normalize the PSDs
psds /= np.sum(psds, axis=1, keepdims=True)

X = []
for fmin, fmax in FREQ_BANDS.values():
freq_res = freqs[1] - freqs[0]

idx_band = np.logical_and(freqs >= fmin, freqs <= fmax)

psds_band = psds[:, (freqs >= fmin) & (freqs < fmax)].mean(axis=-1)
X.append(psds_band.reshape(len(psds), -1))

return np.concatenate(X, axis=1)

data = filt_raw[2][1]
print("Absolute power of delta in Fp1: ", float(bandpower(data, sf, [0.5,4])))
eeg_power_band(filt_raw)

``````

This is the output of the two functions:

All values are really low, and especially for alpha and delta Iâ€™m expecting values between 5 and 20 Hz. Does someone know what the problem is? Iâ€™m already stuck with this problem for a week and I really donâ€™t know anymore. Can someone help me please?

Kind regards,
Laura

Hi Laura,

I guess the line you need to change is

`````` psds_band = psds[:, (freqs >= fmin) & (freqs < fmax)].mean(axis=-1)
``````

I think it should be `sum` instead of `mean`. As you normalized the values, they will sum to 1.
Then when you take the average values of freqs in a band, they will not sum to 1 anymore.

Hi @lrajos,

The values look ok, but to judge if something is wrong it is better to look at the whole spectrum. Remember that these values represent power, not frequency (the frequency coordinates can be found in `freqs` output of psd_welch), and they are in Volt (squared) units. Although the name of the variable is `psd` (power spectral density), I donâ€™t remember if the power values are normalized per Hz (then the units would be `V ** 2 / Hz`).

BTW I see that you normalize the spectrum dividing by its sum - this likely causes your values to be small. As @peara suggests in such case you could sum the spectral values in given range instead of averaging.

Hi Laura, I just discovered this package today - you might find it useful for your purpose yasa.bandpower â€” yasa 0.5.1 documentation

The previous answers are good and I just wanted to add that you might want to consider the unit you have. MNE uses Volts, where you might want microvolts or something else.