The objective is to calculate power spectrum for each frequencies and output data to excel files.
First load packages:
import os
import sys
sys.path.append('Code/code/')
from load_data import load_MEG_dataset
import torch
import torchvision
import matplotlib.pyplot as plt
import numpy as np
from torch.autograd import Variable
from torch.cuda import amp
import torch.nn.functional as F
import math
from scipy.integrate import simps
from mne.time_frequency import psd_array_welch
from band_power import (
bandpower_multi_bands,
standard_scaling_sklearn,
)
import sklearn
import pandas as pd
Then calculate power spectrum for each subject and frequency:
# for sub in [sub for sub in range (1,29) if sub not in [6, 12, 14 ,23]]:
for sub in range(1,2):
exec('abs_{} = []'.format(sub))
exec('rel_{} = []'.format(sub))
Split = 0.90
X_train, y_train = load_MEG_dataset([str(i).zfill(3) for i in range(sub,sub+1)], mode = 'concatenate', output_format='numpy',shuffle = False, training=True, train_test_split=Split, batch_size=500)#, pca_n_components =30)
X_test, y_test = load_MEG_dataset([str(i).zfill(3) for i in range(sub,sub+1)], mode = 'concatenate', output_format='numpy',shuffle = False, training=False, train_test_split=Split, batch_size=500)#, pca_n_components =30)
X_train, X_test = (X_train-X_train.mean())/X_train.std(), (X_test-X_test.mean())/X_test.std()
X_train = X_train[:, None, ...]
X_test = X_test[:, None, ...]
y_train = (y_train / 2) - 1
y_test = (y_test / 2) - 1
X = np.swapaxes(X_train, 2, -1).squeeze()
data = X[X.shape[0]-1, 70, :]
psd_mne, freqs_mne = psd_array_welch(data, 100, 1., 70., n_per_seg=None,
n_overlap=0, n_jobs=1)
for low, high in [(0.5, 4), (4, 8), (8, 10), (10, 12), (12, 30),
(30, 70)]:
print("processing bands (low, high) : ({},{})".format(low, high))
# Find intersecting values in frequency vector
idx_delta = np.logical_and(freqs_mne >= low, freqs_mne <= high)
# Frequency resolution
freq_res = freqs_mne[1] - freqs_mne[0] # = 1 / 4 = 0.25
# Compute the absolute power by approximating the area under the curve
power = simps(psd_mne[idx_delta], dx=freq_res)
exec('abs_{}.append(power)'.format(sub))
# print('Absolute power: {:.4f} uV^2'.format(power))
total_power = simps(psd_mne, dx=freq_res)
rel_power = power / total_power
exec('rel_{}.append(rel_power)'.format(sub))
# print('Relative power: {:.4f}'.format(rel_power))
Last, output calculated power spectrum as xmlx files:
df1 = pd.DataFrame({'sub1':abs_1,
'sub2':abs_2,
'sub3':abs_3,
'sub4':abs_4,
'sub5':abs_5,
'sub7':abs_7,
'sub8':abs_8,
'sub9':abs_9,
'sub10':abs_10,
'sub11':abs_11,
'sub13':abs_13,
'sub15':abs_15,
'sub16':abs_16,
'sub17':abs_17,
'sub18':abs_18,
'sub19':abs_19,
'sub20':abs_20,
'sub21':abs_21,
'sub22':abs_22,
'sub24':abs_24,
'sub25':abs_25,
'sub26':abs_26,
'sub27':abs_27,
'sub28':abs_28},)
df1.to_excel('models/abs.xlsx', sheet_name='sheet1', index=False)
df2 = pd.DataFrame({'sub1':rel_1,
'sub2':rel_2,
'sub3':rel_3,
'sub4':rel_4,
'sub5':rel_5,
'sub7':rel_7,
'sub8':rel_8,
'sub9':rel_9,
'sub10':rel_10,
'sub11':rel_11,
'sub13':rel_13,
'sub15':rel_15,
'sub16':rel_16,
'sub17':rel_17,
'sub18':rel_18,
'sub19':rel_19,
'sub20':rel_20,
'sub21':rel_21,
'sub22':rel_22,
'sub24':rel_24,
'sub25':rel_25,
'sub26':rel_26,
'sub27':rel_27,
'sub28':rel_28},)
df2.to_excel('models/rel.xlsx', sheet_name='sheet1', index=False)
Hint: potential promotions are:
- to use wavelet instead of Fourier transformation;
- to use loop somehow pass data to DataFrame instead of listing them.