541 lines
22 KiB
Python
541 lines
22 KiB
Python
import numpy as np
|
|
import matplotlib
|
|
import matplotlib.pyplot as plt
|
|
import scipy as sp
|
|
import time
|
|
|
|
matplotlib.rcParams['legend.framealpha'] = 1.0
|
|
|
|
packet_definitions = b"""packet_rtc 1 28 uint32_t t 0 4 RTC_TimeTypeDef sTime 4 20 RTC_DateTypeDef sDate 24 4
|
|
packet_vbatt 2 8 uint32_t t 0 4 uint16_t vbatt_cnts 4 2
|
|
packet_imu 3 148 uint32_t t 0 4 uint8_t data[141] 4 1
|
|
packet_adc 4 268 uint32_t t 0 4 uint8_t index 4 1 int32_t ekg_readings_cnts[50] 8 4 int32_t str_readings_cnts[5] 208 4 int32_t oT_readings_cnts[5] 228 4 int32_t iT_readings_cnts[5] 248 4
|
|
packet_spo2 5 184 uint32_t t 0 4 uint8_t bytes[180] 4 1
|
|
packet_msg 6 36 uint32_t t 0 4 char buff[32] 4 1""".split(b'\n')
|
|
|
|
def arr_sizes(s):
|
|
if s[-1:] != b']' or b'[' not in s:
|
|
return 1
|
|
v = int(s[s.rindex(b'[') + 1:-1])
|
|
return v
|
|
|
|
def get_type_list(lines):
|
|
|
|
types = []
|
|
|
|
for line in lines:
|
|
if line != b'':
|
|
L = line.split(b" ")
|
|
types.append({'type_name' : L[0],
|
|
'type_code' : int(L[1]),
|
|
'size' : int(L[2]),
|
|
'elements' : []
|
|
})
|
|
i = 3
|
|
while i < len(L):
|
|
types[-1]['elements'].append({'type_name' : L[i], 'name' : L[i + 1], 'offset' : int(L[i + 2]), 'n_elements' : arr_sizes(L[i + 1]), 'size' : int(L[i + 3]) * arr_sizes(L[i + 1])})
|
|
i += 4
|
|
return types
|
|
|
|
class vbatt_vars():
|
|
|
|
charges = []
|
|
ts = []
|
|
|
|
def process_vbatt(d, t):
|
|
global vbatt_vars
|
|
for e in t['elements']:
|
|
block = d[e['offset']:e['offset'] + e['size']]
|
|
element_size = int(len(block) / e['n_elements'])
|
|
if e['name'] == b't':
|
|
vbatt_vars.ts.append((1 / 2000) * int.from_bytes(block[:4], byteorder = 'little', signed = True))
|
|
if e['name'] == b'vbatt_cnts':
|
|
vbatt_vars.charges.append(2 * 3.3 * int.from_bytes(block, byteorder = 'little', signed = False) / (1 << 12))
|
|
|
|
if len(vbatt_vars.charges) % 10 == 9:
|
|
plt.plot(vbatt_vars.ts, vbatt_vars.charges)
|
|
plt.savefig("vbatt.png")
|
|
plt.close()
|
|
|
|
class ppg_vars():
|
|
|
|
# Red, IR, Green values
|
|
RIGs = np.zeros((0,3))
|
|
ts = []
|
|
freq_Hz = 50
|
|
|
|
def process_ppg(d, t):
|
|
global ppg_vars
|
|
for e in t['elements']:
|
|
block = d[e['offset']:e['offset'] + e['size']]
|
|
element_size = int(len(block) / e['n_elements'])
|
|
if e['name'] == b't':
|
|
ppg_vars.ts.append((1 / 2000) * int.from_bytes(block[:4], byteorder = 'little', signed = True))
|
|
if e['name'] == b'bytes[180]':
|
|
ppg_vars.RIGs = np.concatenate((ppg_vars.RIGs, np.array([int.from_bytes(block[3 * i : 3 * i + 3], byteorder = 'big') for i in range(0,60)]).reshape(20,3)))
|
|
|
|
if len(ppg_vars.ts) % 3 == 2:
|
|
ppg_graph()
|
|
|
|
def ppg_graph():
|
|
|
|
#ratio of pulsatile / DC per wavelength gives the SPO2 measure
|
|
|
|
# Plots of short term PPG vs R peaks and long term PPG vs breathing
|
|
fig, axs = plt.subplots(2, 1, height_ratios = [1,2])
|
|
|
|
for i, ax in enumerate(axs):
|
|
|
|
if i == 0:
|
|
b_ppg, a_ppg = sp.signal.butter(2, [1 / (0.5 * ppg_vars.freq_Hz)], btype = 'highpass')
|
|
if i == 1:
|
|
b_ppg, a_ppg = sp.signal.butter(2, 0.05 / (0.5 * ppg_vars.freq_Hz), btype = 'highpass')
|
|
|
|
b_breath, a_breath = sp.signal.butter(2, 0.05 / (0.5 * 0.1 * adc_vars.freq_Hz), btype = 'highpass')
|
|
|
|
if i == 0:
|
|
lookback_s = 5
|
|
if i == 1:
|
|
lookback_s = 60
|
|
|
|
breath_points = int((0.1 * adc_vars.freq_Hz) * lookback_s)
|
|
ppg_points = int(ppg_vars.freq_Hz * lookback_s)
|
|
|
|
ts = len(adc_vars.ecgs) / adc_vars.freq_Hz - np.arange(len(ppg_vars.RIGs))[::-1] / ppg_vars.freq_Hz
|
|
if i == 0: # plot R wave peaks
|
|
for peak in adc_vars.r_peaks:
|
|
peak_s = peak / adc_vars.freq_Hz
|
|
if peak_s < ts[-1] and peak_s > ts[max(-ppg_points, -len(ts))]:
|
|
ax.axvline(x = peak_s, color = 'black', alpha = 0.5)
|
|
if i == 1: # plot chest strap strain
|
|
ax.plot((len(adc_vars.ecgs) - 10 * np.arange(min(breath_points, len(adc_vars.t2s)))[::-1]) / (adc_vars.freq_Hz), sp.signal.filtfilt(b_breath, a_breath, adc_vars.t2s[-breath_points:]), color = 'black')
|
|
|
|
ppg_ax = ax.twinx()
|
|
|
|
if False:
|
|
smoothed_ppg = sp.signal.filtfilt(b_ppg, a_ppg, ppg_vars.RIGs[-ppg_points:,:], axis = 0)
|
|
else:
|
|
smoothed_ppg = ppg_vars.RIGs[-ppg_points:,:] - np.mean(ppg_vars.RIGs[-ppg_points:,:],axis = 0)
|
|
ppg_ax.plot(ts[-ppg_points:], smoothed_ppg[:,0], 'r', alpha = 0.5)
|
|
ppg_ax.plot(ts[-ppg_points:], smoothed_ppg[:,1], 'm', alpha = 0.5)
|
|
ppg_ax.plot(ts[-ppg_points:], smoothed_ppg[:,2], 'g', alpha = 0.5)
|
|
|
|
#if i == 0:
|
|
# lim = 75
|
|
#else:
|
|
# lim = 250
|
|
#ppg_ax.set_ylim(-lim, lim)
|
|
|
|
ax.set_yticks([])
|
|
|
|
ppg_ax.set_yticks([])
|
|
|
|
ax.set_xlabel("t (s)")
|
|
|
|
if i == 1:
|
|
|
|
ax.plot([],[],'r',label='Red Val', alpha = 0.5)
|
|
ax.plot([],[],'m',label='IR Val', alpha = 0.5)
|
|
ax.plot([],[],'g',label='Green Val', alpha = 0.5)
|
|
ax.plot([],[],'k',label='Chest Diameter')
|
|
ax.legend(loc = 'lower center',
|
|
bbox_to_anchor = (0.5, 1.02),
|
|
ncol = 4)
|
|
|
|
plt.tight_layout()
|
|
plt.savefig("ppg.png")
|
|
plt.close()
|
|
|
|
class imu_vars:
|
|
|
|
accs = np.zeros((0,3))
|
|
gyros = np.zeros((0,3))
|
|
ts = []
|
|
freq_Hz = 240
|
|
|
|
def process_imu(d, t):
|
|
global imu_vars
|
|
for e in t['elements']:
|
|
block = d[e['offset']:e['offset'] + e['size']]
|
|
element_size = int(len(block) / e['n_elements'])
|
|
if e['name'] == b't':
|
|
imu_vars.ts.append((1 / 2000) * int.from_bytes(block[:4], byteorder = 'little', signed = True))
|
|
if e['name'] == b'data[141]':
|
|
for i in range(20):
|
|
reading = block[1 + 7 * i : 8 + 7 * i]
|
|
#print(reading)
|
|
imu_reading_type = reading[0] >> 3
|
|
imu_reading_tag_cnt = (reading[0] >> 1) & 3
|
|
data = np.array([int.from_bytes(reading[2 * i + 1 : 2 * i + 3], byteorder = 'little', signed = True) for i in range(3)]).reshape(1,3)
|
|
if imu_reading_type == 1:
|
|
imu_vars.gyros = np.concatenate((imu_vars.gyros, (250 / (1<<16) * data)))
|
|
elif imu_reading_type == 2:
|
|
imu_vars.accs = np.concatenate((imu_vars.accs, (4 / (1<<16) * data)))
|
|
else:
|
|
assert False
|
|
|
|
def imu_graph():
|
|
return
|
|
tt = int(5 * imu_vars.freq_Hz)
|
|
if len(gyros) < 480:
|
|
return
|
|
fig, axs = plt.subplots(2)
|
|
b, a = sp.signal.butter(2, 15 / (0.5 * imu_vars.freq_Hz), btype = 'lowpass')
|
|
g = sp.signal.filtfilt(b, a, np.abs(np.diff(np.array(imu_vars.gyros), axis = 0)), axis = 0)
|
|
a = sp.signal.filtfilt(b, a, np.abs(np.diff(np.array(imu_vars.accs), axis = 0)), axis = 0)
|
|
axs[0].set_ylabel("dps")
|
|
xs_g = len(adc_vars.ecgs) / adc_vars.freq_Hz - (np.arange(len(imu_vars.gyros) - 1) / imu_freq_Hz)[::-1]
|
|
xs_a = len(adc_vars.ecgs) / adc_vars.freq_Hz - (np.arange(len(imu_vars.accs) - 1) / imu_freq_Hz)[::-1]
|
|
axs[0].plot(xs_g[-tt:], np.sum(g[-tt:,:], axis = 1))
|
|
#axs[0].plot(xs_g[-tt:], g[-tt:,1])
|
|
#axs[0].plot(xs_g[-tt:], g[-tt:,2])
|
|
axs[1].set_ylabel("g")
|
|
axs[1].plot(xs_a[-tt:], np.sum(np.power(a[-tt:,:], 2.0), axis = 1))
|
|
#axs[1].plot(xs_a[-tt:], a[-tt:,1])
|
|
#axs[1].plot(xs_a[-tt:], a[-tt:,2])
|
|
for peak in adc_vars.r_peaks:
|
|
peak_s = peak / adc_vars.freq_Hz
|
|
if peak_s < xs_a[-1] and peak_s > xs_a[-tt]:
|
|
axs[1].axvline(x = peak_s, color = 'black', alpha = 0.5)
|
|
axs[0].axvline(x = peak_s, color = 'black', alpha = 0.5)
|
|
#axs[0].set_ylim(0,1.0)
|
|
#axs[1].set_ylim(0,0.05)
|
|
plt.tight_layout()
|
|
plt.savefig("acc_gyro.png")
|
|
plt.close()
|
|
|
|
if len(adc_vars.r_peaks) > 30:
|
|
fig, axs = plt.subplots(2)
|
|
for i in range(1,30):
|
|
xs_g = len(adc_vars.ecgs) / adc_vars.freq_Hz - (np.arange(len(gyros)) / imu_freq_Hz)[::-1]
|
|
ys = np.interp(np.linspace(adc_vars.r_peaks[-i] / adc_vars.freq_Hz - 1.0, adc_vars.r_peaks[-i] / adc_vars.freq_Hz + 1.0, 240),
|
|
xs_g,
|
|
gyros[:,0])
|
|
axs[0].plot(np.convolve(np.abs(np.diff(ys)),np.ones(10),mode='valid'), 'k.', alpha = 0.05)
|
|
ys = np.interp(np.linspace(adc_vars.r_peaks[-i] / adc_vars.freq_Hz - 1.0, adc_vars.r_peaks[-i] / adc_vars.freq_Hz + 1.0, 240),
|
|
xs_g,
|
|
gyros[:,1])
|
|
axs[0].plot(np.convolve(np.abs(np.diff(ys)),np.ones(10),mode='valid'), 'b.', alpha = 0.05)
|
|
ys = np.interp(np.linspace(adc_vars.r_peaks[-i] / adc_vars.freq_Hz - 1.0, adc_vars.r_peaks[-i] / adc_vars.freq_Hz + 1.0, 240),
|
|
xs_g,
|
|
gyros[:,2])
|
|
axs[0].plot(np.convolve(np.abs(np.diff(ys)),np.ones(10),mode='valid'), 'r.', alpha = 0.05)
|
|
for i in range(1,30):
|
|
xs_a = len(adc_vars.ecgs) / adc_vars.freq_Hz - (np.arange(len(accs)) / imu_freq_Hz)[::-1]
|
|
ys = np.interp(np.linspace(adc_vars.r_peaks[-i] / adc_vars.freq_Hz - 1.0, adc_vars.r_peaks[-i] / adc_vars.freq_Hz + 1.0, 240),
|
|
xs_a,
|
|
accs[:,0])
|
|
axs[1].plot(np.convolve(np.abs(np.diff(ys)),np.ones(10),mode='valid'), 'k.', alpha = 0.05)
|
|
ys = np.interp(np.linspace(adc_vars.r_peaks[-i] / adc_vars.freq_Hz - 1.0, adc_vars.r_peaks[-i] / adc_vars.freq_Hz + 1.0, 240),
|
|
xs_a,
|
|
accs[:,1])
|
|
axs[1].plot(np.convolve(np.abs(np.diff(ys)),np.ones(10),mode='valid'), 'b.', alpha = 0.05)
|
|
ys = np.interp(np.linspace(adc_vars.r_peaks[-i] / adc_vars.freq_Hz - 1.0, adc_vars.r_peaks[-i] / adc_vars.freq_Hz + 1.0, 240),
|
|
xs_a,
|
|
accs[:,2])
|
|
axs[1].plot(np.convolve(np.abs(np.diff(ys)),np.ones(10),mode='valid'), 'r.', alpha = 0.05)
|
|
|
|
plt.savefig("accs_ensemble.png")
|
|
plt.close()
|
|
|
|
|
|
class adc_vars():
|
|
|
|
ecgs = np.zeros(0)
|
|
t1s = np.zeros(0)
|
|
t2s = np.zeros(0)
|
|
strains = np.zeros(0)
|
|
ts = []
|
|
|
|
freq_Hz = 488.28125
|
|
max_hr_Hz = 120 / 60
|
|
|
|
r_peaks = []
|
|
|
|
on_off_times = [0]
|
|
|
|
def process_adc(d, t):
|
|
global adc_vars
|
|
|
|
for e in t['elements']:
|
|
block = d[e['offset']:e['offset'] + e['size']]
|
|
element_size = int(len(block) / e['n_elements'])
|
|
|
|
if e['name'] == b't':
|
|
adc_vars.ts.append((1 / 2000) * int.from_bytes(block[:4], byteorder = 'little', signed = True))
|
|
if e['name'] == b'ekg_readings_cnts[50]':
|
|
adc_vars.ecgs = np.concatenate((adc_vars.ecgs, np.array([(2.4 / (1<<24)) * int.from_bytes(block[4 * i : 4 * i + 4], byteorder = 'little', signed = True) for i in range(50)])))
|
|
if e['name'] == b'str_readings_cnts[5]':
|
|
adc_vars.strains = np.concatenate((adc_vars.strains, np.array([(2.4 / (1<<24)) * int.from_bytes(block[4 * i : 4 * i + 4], byteorder = 'little', signed = True) for i in range(5)])))
|
|
if e['name'] == b'oT_readings_cnts[5]':
|
|
adc_vars.t1s = np.concatenate((adc_vars.t1s, np.array([(2.4 / (1<<24)) * int.from_bytes(block[4 * i : 4 * i + 4], byteorder = 'little', signed = True) for i in range(5)])))
|
|
if e['name'] == b'iT_readings_cnts[5]':
|
|
adc_vars.t2s = np.concatenate((adc_vars.t2s, np.array([(2.4 / (1<<24)) * int.from_bytes(block[4 * i : 4 * i + 4], byteorder = 'little', signed = True) for i in range(5)])))
|
|
|
|
update_ecg_peaks()
|
|
|
|
if len(adc_vars.ts) % 10 == 9:
|
|
adc_graphs()
|
|
|
|
def dead_time_cond(arr):
|
|
return np.amax(arr) - np.amin(arr) > 0.005 or np.amax(np.abs(np.diff(arr))) > 0.001
|
|
|
|
def update_ecg_peaks():
|
|
global adc_vars
|
|
|
|
start = adc_vars.on_off_times[-1]
|
|
|
|
if len(adc_vars.r_peaks) > 0:
|
|
start = max(start, adc_vars.r_peaks[-1])
|
|
|
|
v = np.diff(adc_vars.ecgs[start:])
|
|
inds = [start + e for e in np.where(v < -0.00025)[0]]
|
|
|
|
last_ind = start
|
|
for ind in inds:
|
|
if ind - last_ind < adc_vars.freq_Hz / adc_vars.max_hr_Hz:
|
|
continue
|
|
if dead_time_cond(adc_vars.ecgs[ind - 200 : max(ind + 200, len(adc_vars.ecgs))]):
|
|
# If not ok, and was ok, turn off
|
|
if len(adc_vars.on_off_times) % 2 == 1:
|
|
assert(adc_vars.on_off_times[-1] < ind - 200)
|
|
adc_vars.on_off_times.append(ind - 200)
|
|
continue
|
|
else:
|
|
# If ok, and was not ok, turn on
|
|
if len(adc_vars.on_off_times) % 2 == 0:
|
|
assert(adc_vars.on_off_times[-1] < ind - 200)
|
|
adc_vars.on_off_times.append(ind - 200)
|
|
region_start = ind - int(0.1 * adc_vars.freq_Hz)
|
|
region_end = ind + int(0.1 * adc_vars.freq_Hz)
|
|
if region_end > len(adc_vars.ecgs):
|
|
break
|
|
peak = region_start + np.argmax(adc_vars.ecgs[region_start : region_end])
|
|
adc_vars.r_peaks.append(peak)
|
|
last_ind = ind
|
|
|
|
def adc_graphs():
|
|
|
|
fig, axs = plt.subplots(2, 1, height_ratios = [2,1])
|
|
|
|
ecg_lookback_s = 5
|
|
RR_lookback_s = 600
|
|
ecg_points = int(ecg_lookback_s * adc_vars.freq_Hz)
|
|
ts = np.arange(max(len(adc_vars.ecgs) - ecg_points, 0), len(adc_vars.ecgs)) / adc_vars.freq_Hz
|
|
|
|
axs[0].plot(ts,
|
|
adc_vars.ecgs[-ecg_points:], 'k.', linestyle='--', alpha = 0.5)
|
|
axs[0].set_ylim(np.amin(adc_vars.ecgs[-ecg_points:]) - 0.001, np.amax(adc_vars.ecgs[-ecg_points:]) + 0.001)
|
|
for peak in adc_vars.r_peaks:
|
|
if peak < ts[-1] * adc_vars.freq_Hz and peak > ts[0] * adc_vars.freq_Hz:
|
|
axs[0].plot(peak / adc_vars.freq_Hz, adc_vars.ecgs[peak], 'ro')
|
|
|
|
axs[0].set_ylabel("V")
|
|
axs[0].set_xlabel("t (s)")
|
|
|
|
# axs[1].plot(ts[:-1], np.diff(adc_vars.ecgs[-ecg_points:]), 'k.', linestyle='--', alpha = 0.25)
|
|
# for peak in adc_vars.r_peaks:
|
|
# if peak < ts[-1] * adc_vars.freq_Hz and peak > ts[0] * adc_vars.freq_Hz:
|
|
# axs[1].axvline(peak / adc_vars.freq_Hz, color = 'red', alpha = 0.1)
|
|
|
|
|
|
# plt.tight_layout()
|
|
# plt.savefig("hrv_biofeedback.png", dpi = 150)
|
|
# plt.close()
|
|
|
|
# return
|
|
|
|
strain_ax = axs[1].twinx()
|
|
|
|
|
|
#b_breath, a_breath = sp.signal.butter(2, 0.05 / (0.5 * 0.1 * adc_vars.freq_Hz), btype = 'highpass')
|
|
sos_breath = sp.signal.butter(1, 0.05 / (0.5 * 0.1 * adc_vars.freq_Hz), btype = 'highpass', output = 'sos')
|
|
|
|
filtered_strain = sp.signal.sosfiltfilt(sos_breath, adc_vars.t2s[-int(RR_lookback_s * adc_vars.freq_Hz * 0.1):])
|
|
strain_ax.plot(np.arange(max(len(adc_vars.t2s) - len(filtered_strain), 0), len(adc_vars.t2s)) * 10 / adc_vars.freq_Hz, filtered_strain, color = 'orange', alpha = 0.5)
|
|
strain_ax.set_yticks([])
|
|
|
|
first_valid_RR = len(adc_vars.ecgs) - RR_lookback_s * adc_vars.freq_Hz
|
|
r_peaks = [e for e in adc_vars.r_peaks if e > first_valid_RR]
|
|
|
|
hrv_ax = axs[1].twinx()
|
|
|
|
live_times = adc_vars.on_off_times[:]
|
|
if len(adc_vars.on_off_times) % 2 == 1:
|
|
live_times.append(len(adc_vars.ecgs))
|
|
live_times = np.array(live_times).reshape(-1,2)
|
|
for start, end in live_times:
|
|
r_peaks_ = [e for e in r_peaks if e > start and e < end]
|
|
axs[1].plot(np.array(r_peaks_[:-1]) / adc_vars.freq_Hz, 1000 * np.diff(r_peaks_) / adc_vars.freq_Hz, 'k.', alpha = 0.5)
|
|
hrv = np.abs(1000 * np.diff(np.diff(r_peaks_)) / adc_vars.freq_Hz)
|
|
hrv_ax.plot(np.array(r_peaks_[2:]) / adc_vars.freq_Hz, hrv, 'b.', alpha = 0.5)
|
|
|
|
if len(r_peaks_) > 10:
|
|
|
|
N = int(np.floor(len(r_peaks_) / 30))
|
|
for i in range(N):
|
|
block = 1000 * np.array(r_peaks_[30 * i : min(len(r_peaks_), 30 * i + 30)]) / adc_vars.freq_Hz
|
|
if len(block) > 2:
|
|
bpm = 60 * len(block) / ((block[-1] - block[0]) / 1000)
|
|
rmssd = np.power(np.sum(np.power(np.diff(np.diff(block)), 2.0)) / (len(block) - 2), 0.5)
|
|
sdnn = np.std(np.diff(block))
|
|
pNN50 = 100 * np.mean(np.abs(np.diff(np.diff(block))) > 50)
|
|
hrv_ax.plot([block[0] / 1000, block[-1] / 1000], [rmssd, rmssd], 'r', alpha = 0.5)
|
|
hrv_ax.plot([block[0] / 1000, block[-1] / 1000], [sdnn, sdnn], 'm', alpha = 0.5)
|
|
hrv_ax.plot([block[0] / 1000, block[-1] / 1000], [pNN50, pNN50], color = 'pink', alpha = 0.5)
|
|
|
|
hrv_ax.minorticks_on()
|
|
hrv_ax.grid(alpha = 0.5)
|
|
hrv_ax.tick_params(axis = 'y', colors = 'blue')
|
|
hrv_ax.yaxis.tick_right()
|
|
|
|
axs[1].set_ylim(0, 1200)
|
|
hrv_ax.set_ylim(0, 240)
|
|
axs[1].set_ylabel("R-R delta T (ms)")
|
|
hrv_ax.set_ylabel("delta R-R (ms)", color = 'blue')
|
|
axs[1].plot([],'k.',label = 'R-R int', alpha = 0.5)
|
|
axs[1].plot([],'b.',label = '|delta(R-R int)|', alpha = 0.5)
|
|
axs[1].plot([],color = 'orange', label = 'Breath (tension)', alpha = 0.5)
|
|
|
|
try:
|
|
axs[1].plot([],color = 'red', label = f'RMSSD {int(rmssd)}ms', alpha = 0.5)
|
|
except NameError:
|
|
axs[1].plot([],color = 'red', label = f'RMSSD', alpha = 0.5)
|
|
try:
|
|
axs[1].plot([],color = 'magenta', label = f'SSDN {int(sdnn)}ms', alpha = 0.5)
|
|
except NameError:
|
|
axs[1].plot([],color = 'magenta', label = 'SSDN', alpha = 0.5)
|
|
try:
|
|
axs[1].plot([],color = 'black', label = f'BPM {int(bpm)}', alpha = 0.5)
|
|
except:
|
|
pass
|
|
axs[1].plot([],color = 'pink', label = 'pNN50', alpha = 0.5)
|
|
axs[1].legend(loc = 'lower center', bbox_to_anchor = (0.5, 1.02), ncol = 4)
|
|
axs[1].set_xlabel("t (s)")
|
|
|
|
hrv_ax.set_axisbelow(True)
|
|
plt.tight_layout()
|
|
plt.savefig("hrv_biofeedback.png", dpi = 150)
|
|
plt.close()
|
|
|
|
fig, ax = plt.subplots(1)
|
|
for i in range(1, min(len(adc_vars.r_peaks), 30)):
|
|
dat = adc_vars.ecgs[adc_vars.r_peaks[-i] - int(0.3 * adc_vars.freq_Hz) :
|
|
adc_vars.r_peaks[-i] + int(0.9 * adc_vars.freq_Hz)]
|
|
dat = dat - np.mean(np.sort(dat[40:-40]))
|
|
ax.plot(dat, '.', alpha = 0.02, color = matplotlib.colormaps['viridis'](i / 30.))
|
|
ax.set_ylim(-0.001, 0.002)
|
|
plt.savefig("pqrs_ensemble.png")
|
|
plt.close()
|
|
|
|
return
|
|
|
|
#also can do poincare plot
|
|
# for fft, resample RR interpolation to a 4Hz grid wich cubic spline
|
|
# welches method or FFT->PSD
|
|
# there's also Lomb-Scargle periodogram
|
|
|
|
# Poincare
|
|
#R = np.diff(r_peaks[-180:]) / freq_Hz
|
|
#axs[2].plot(R[:-1], R[1:], 'k.', alpha = 0.2)
|
|
|
|
# lombscargle
|
|
#Ys = sp.signal.lombscargle(np.diff(r_peaks[-240:])[1:] / freq_Hz,
|
|
# np.diff(np.diff(r_peaks[-240:])) / freq_Hz,
|
|
# np.linspace(0.01,2.0,100),
|
|
# floating_mean = False)
|
|
#axs[2].plot(np.linspace(0.01,2.0,100), Ys)
|
|
|
|
# Power spectrum
|
|
# if len(r_peaks) > 120:
|
|
# R = np.array(r_peaks[-120:]) / freq_Hz
|
|
# fn = sp.interpolate.CubicSpline(0.5 * (R[1:] + R[:-1]), np.diff(R))
|
|
# end = 0.5 * (R[-2] + R[-1])
|
|
# Xs = np.linspace(end - 60, end, 1000)
|
|
# #np.log(np.abs(np.fft.fft(fn(Xs))))[1:])
|
|
# X, Y = sp.signal.welch(fn(Xs), fs = 1000 / 60)
|
|
# f_pwrs.append(Y[1:10])
|
|
# if len(f_pwrs) > 60:
|
|
# f_pwrs = f_pwrs[-60:]
|
|
# axs[2].imshow(f_pwrs, aspect = 'auto', extent = (X[1], X[9], 0, len(f_pwrs)))
|
|
|
|
|
|
|
|
|
|
# def make_graphs():
|
|
|
|
# fig, axs = plt.subplots(2,2)
|
|
|
|
# b, a = sp.signal.butter(2, [1 / (0.5 * freq_Hz), 120 / (0.5 * freq_Hz)], btype = 'bandpass')
|
|
# ecgs_ = sp.signal.filtfilt(b, a, ecgs)
|
|
|
|
# r_peaks = sp.signal.find_peaks(ecgs, height = None, threshold = None, distance = freq_Hz / max_hr_Hz, width = 5, prominence = 0.001)[0]
|
|
# axs[0][0].plot(np.arange(len(ecgs)) / freq_Hz, ecgs_, 'k.', linestyle='--')
|
|
# for peak in r_peaks:
|
|
# axs[0][0].plot(peak / freq_Hz, ecgs_[peak], 'ro')
|
|
# axs[1][0].plot(r_peaks[:-1] / freq_Hz, 1 / (np.diff(r_peaks) /freq_Hz), 'k.')
|
|
|
|
# b, a = sp.signal.butter(2, [0.05 / (0.5 * ppg_freq_Hz), 4 / (0.5 * ppg_freq_Hz)], btype = 'bandpass')
|
|
# reds_ = sp.signal.filtfilt(b, a, reds)
|
|
# #irs_ = sp.signal.filtfilt(b, a, irs)
|
|
# #greens_ = sp.signal.filtfilt(b, a, greens)
|
|
|
|
# P_r = np.log(np.abs(np.fft.rfft(reds_)))
|
|
# P_i = np.log(np.abs(np.fft.rfft(irs)))
|
|
# P_g = np.log(np.abs(np.fft.rfft(greens)))
|
|
|
|
# #axs[0][1].plot(np.arange(P_r.shape[0]) * 1 / 200, P_r)
|
|
# #axs[0][1].plot(P_i)
|
|
# #axs[0][1].plot(P_g)
|
|
|
|
# axs[0][1].plot(np.arange(len(reds)) / 50, reds_, color = 'red')
|
|
# #ax2 = axs[0][1].twinx()
|
|
# #ax2.plot(np.arange(len(reds)) / 50, irs_, color = 'magenta')
|
|
# #ax3 = ax2.twinx()
|
|
# #ax3.plot(np.arange(len(reds)) / 50, greens_, color = 'green')
|
|
|
|
# #axs[1][1].plot(np.array(gyros)[:,0])
|
|
# #axs[1][1].plot(np.array(gyros)[:,1])
|
|
# #axs[1][1].plot(np.array(gyros)[:,2])
|
|
# gs = np.array(gyros)
|
|
# acs = np.array(accs)
|
|
# f_gs = np.log(np.abs(np.fft.rfft(acs, axis = 0)))
|
|
# axs[1][1].plot(f_gs[:,0], alpha = 0.25)
|
|
# axs[1][1].plot(f_gs[:,1], alpha = 0.25)
|
|
# axs[1][1].plot(f_gs[:,2], alpha = 0.25)
|
|
|
|
|
|
# #axs[1][1].plot(strains)
|
|
# #axs[1][1].plot(t1s)
|
|
# #axs[1][1].plot(t2s)
|
|
|
|
# #plt.show()
|
|
|
|
def read_and_process(types, cons, size):
|
|
index = 0
|
|
while (index < size):
|
|
packet_type = cons[index]
|
|
print(packet_type)
|
|
try:
|
|
t = [t for t in types if t['type_code'] == packet_type][0]
|
|
except:
|
|
print("HERE")
|
|
print(cons[index-5:index+5])
|
|
quit()
|
|
return
|
|
print(index, packet_type, t['type_name'])
|
|
d = cons[index + 1 : index + 1 + t['size']]
|
|
if t['type_name'] == b'packet_imu':
|
|
process_imu(d, t)
|
|
if t['type_name'] == b'packet_msg':
|
|
print(d)
|
|
if t['type_name'] == b'packet_adc':
|
|
process_adc(d, t)
|
|
if t['type_name'] == b'packet_spo2':
|
|
process_ppg(d, t)
|
|
if t['type_name'] == b'packet_vbatt':
|
|
process_vbatt(d, t)
|
|
index += 1 + t['size']
|
|
|
|
|