-
Notifications
You must be signed in to change notification settings - Fork 0
/
h5_exploitation.py
115 lines (94 loc) · 3.09 KB
/
h5_exploitation.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
# -*- coding: utf-8 -*-
"""Formating data from first experiment.
Conditions:
* 8 minutes lying down in bed with eyes closed
* 4 two-minutes sessions: silence, monaural 1f 247Hz,
monaural 2f 240Hz + 254 Hz, binaural 240 Hz + 254 Hz.
Hypothesis of ear symetry in binaural beats.
"""
import h5py
import numpy as np
from matplotlib import pyplot as plt
from math import floor, sqrt
import array_functions as array
f = h5py.File("acquisition_1.h5", "r")
fs = 250
f_acc = 50
f2 = f[u'channel4'][u'filtered']
sec = 20
x_lims = [51000+8*250*60-sec*250, 51000+8*250*60+sec*250]
# Ok let's use channel 4
# Now look at the accelerometer
acc_x = np.array(f[u'accelerometer'][u'x'])
"""
plt.figure(3)
plt.plot(acc_x)
plt.title("acc_x")
plt.ylim([-5, 5])
plt.xlim([x_lims[0]/5, x_lims[1]/5])
plt.show()
"""
def SO_fft_window_rel(EEG_data, fs=250):
"""Return alpha power after fft on the EEG data.
Args:
EEG_data (np array)
fs (int): sampling frequency in Hz.
Returns:
alpha (np array): alpha power of the EEG_data
"""
n = len(EEG_data)
fft = abs(np.fft.fft(EEG_data, n=None, axis=-1)[:n / 2 + 1])
pow_fft = fft**2
# relative band power
tot_pow = np.sum(pow_fft)
alpha = sqrt(
np.sum(pow_fft[int(8 * n * 2 / fs):int(12 * n * 2 / fs)])/tot_pow)
theta = sqrt(
np.sum(pow_fft[int(4 * n * 2 / fs):int(8 * n * 2 / fs)])/tot_pow)
return alpha, theta
def alpha_power(EEG_data, win_len, exp_smoothing_coeff, fs=250):
"""Return alpha relative power on a sliding window at 1 Hz.
Args:
EEG_data (np array)
win_len (int): floating window length in sec
exp_smoothing_coeff (float)
fs (int): EEG sampling frequency in Hz
Returns:
alpha_smoothed (np array): relative alpha power at 1 Hz.
"""
npts_win = fs * win_len
alpha = []
# For computational reasons we calculate the powers at each second
for i in range(int(floor(len(EEG_data)/fs))-win_len):
sig_win = EEG_data[i*fs: i*fs+npts_win-1]
fft_w = SO_fft_window_rel(sig_win, fs)[0]
alpha.append(fft_w)
# Now append zeros at the beginning of all the channels
zeros = np.zeros((win_len, ))
alpha = np.concatenate((zeros, alpha), axis=0)
# Smooth it a bit
alpha_smoothed = array.mean_smoothing(
array.exp_smoothing(alpha, exp_smoothing_coeff),
win_len)
return alpha_smoothed
# Lionel veut que je lui calcule une fft
# Ok data is between 51000 and 51000+8*60*250
EEG_data = f2[51000:51000 + 8*250*60]
plt.figure(2)
plt.plot(EEG_data)
plt.ylim([-50, 50])
plt.title("EEG signal")
# plt.xlim(x_lims)
plt.show()
silence = SO_fft_window_rel(EEG_data[0:2*60*250])
monaural_1f = SO_fft_window_rel(EEG_data[2*60*250:4*60*250])
monaural_2f = SO_fft_window_rel(EEG_data[4*60*250:6*60*250])
binaural = SO_fft_window_rel(EEG_data[6*60*250:8*60*250])
alpha_slide = alpha_power(EEG_data, 10, 0.1)
plt.figure(4)
plt.plot(alpha_slide)
plt.axvline(x=2*60, linewidth=4, color="r")
plt.axvline(x=4*60, linewidth=4, color="r")
plt.axvline(x=6*60, linewidth=4, color="r")
plt.axvline(x=8*60, linewidth=4, color="r")
plt.show()