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.externals.pymatreader import read_mat
from mne.preprocessing import (ICA, create_eog_epochs, create_ecg_epochs,
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

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

dat = mat73.loadmat(file)    # load the matlab file 

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

# create info 
sf = 250
n_channels = 29
ch_names = ['time','event','Fp1','Fp2','F3','Fz','F4','T7','C3','Cz','C4','T8','P3','Pz','P4','O1','O2','Oz','bad1','bad2','bad3','bad4','bad5','bad6','bad7','bad8','bad9','bad10','bad11']
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 =,info)

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

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

raw.copy().pick_types(eeg=True,exclude='bads').plot(n_channels = 16,scalings={"eeg": 950e-1}, remove_dc=True)

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

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

ica.exclude = list(comps) 

# 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, 
                                                      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])))

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?

Thank you in advance.

Kind regards,

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.