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,
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
# 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 = 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')
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)
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?
Thank you in advance.
Kind regards,
Laura