diff --git a/code/l452_code/packet_parser_helpers.py b/code/l452_code/packet_parser_helpers.py index 2d7339f..525d310 100644 --- a/code/l452_code/packet_parser_helpers.py +++ b/code/l452_code/packet_parser_helpers.py @@ -1,5 +1,6 @@ import numpy as np import matplotlib.pyplot as plt +import matplotlib import scipy as sp import time @@ -34,115 +35,189 @@ def get_type_list(lines): i += 4 return types -reds = [] -irs = [] -greens = [] +RIG_amp = np.zeros((0,3)) ppg_freq_Hz = 50 def process_ppg(d, t): - global greens, reds, irs + global RIG_amp 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'bytes[180]': - reds += [int.from_bytes(block[3 * i : 3 * i + 3], byteorder = 'big') for i in range(0,60,3)] - irs += [int.from_bytes(block[3 * i : 3 * i + 3], byteorder = 'big') for i in range(1,60,3)] - greens += [int.from_bytes(block[3 * i : 3 * i + 3], byteorder = 'big') for i in range(2,60,3)] - if len(reds) > 400: - reds = reds[-400:] - irs = irs[-400:] - greens = greens[-400:] - fig, axs = plt.subplots(3) - axs[0].set_title('red') - axs[1].set_title('ir') - axs[2].set_title('green') - axs[0].plot(reds) - axs[1].plot(irs) - axs[2].plot(greens) - plt.savefig("ppg.png") - plt.close() + RIG_amp = np.concatenate((RIG_amp, np.array([int.from_bytes(block[3 * i : 3 * i + 3], byteorder = 'big') for i in range(0,60)]).reshape(20,3))) -accs = [] -gyros = [] -imu_sparse = 0 + if len(ecgs) < ecg_freq_Hz * 100: # 500 + return + + # top is short again, bottom is long + fig, axs = plt.subplots(2) + + for i, ax in enumerate(axs): + + b, a = sp.signal.butter(2, 0.05 / (0.5 * ppg_freq_Hz), btype = 'highpass') + b1, a1 = sp.signal.butter(2, 0.05 / (0.5 * 0.1 * ecg_freq_Hz), btype = 'highpass') + + if i == 0: + lookback_s = 5 + else: + lookback_s = 60 + + breath_points = int(lookback_s * (0.1 * ecg_freq_Hz)) + if i == 1: + ax.plot((len(ecgs) - 10 * np.arange(breath_points)[::-1]) / (ecg_freq_Hz), sp.signal.filtfilt(b1, a1, t2s[-breath_points:]), color = 'black') + r_ax = ax.twinx() + i_ax = ax.twinx() + g_ax = ax.twinx() + tt = - ppg_freq_Hz * lookback_s + xs = len(ecgs) / ecg_freq_Hz - np.arange(len(RIG_amp))[::-1] / ppg_freq_Hz + r_ax.plot(xs[tt:], sp.signal.filtfilt(b, a, RIG_amp[tt:,0]), 'r', alpha = 0.5) + i_ax.plot(xs[tt:], sp.signal.filtfilt(b, a, RIG_amp[tt:,1]), 'm', alpha = 0.5) + g_ax.plot(xs[tt:], sp.signal.filtfilt(b, a, RIG_amp[tt:,2]), 'g', alpha = 0.5) + + r_ax.set_ylim(-750,750) + i_ax.set_ylim(-750,750) + g_ax.set_ylim(-750,750) + + if i == 0: + for peak in r_peaks: + peak_s = peak / ecg_freq_Hz + if peak_s < xs[-1] and peak_s > xs[tt]: + r_ax.axvline(x = peak_s, color = 'black', alpha = 0.5) + + plt.tight_layout() + plt.savefig("ppg.png") + plt.close() + +accs = np.zeros((0,3)) +gyros = np.zeros((0,3)) imu_freq_Hz = 240 +last_imu_graph = 0 +imu_ts = [] def process_imu(d, t): - global accs, gyros, imu_sparse + global accs, gyros, imu_sparse, last_imu_graph, imu_ts 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_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)]) + 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: - gyros.append(250 / (1<<16) * data) + gyros = np.concatenate((gyros, (250 / (1<<16) * data))) elif imu_reading_type == 2: - accs.append(4 / (1<<16) * data) + accs = np.concatenate((accs, (4 / (1<<16) * data))) else: pass #assert False + last_imu_graph += 1 + if len(ecgs) < ecg_freq_Hz * 100: + return + + if (last_imu_graph % 12 == 11):#True and time.time() - last_imu_graph > 0.5: + tt = int(5 * imu_freq_Hz) + if len(gyros) < 480: + return + fig, axs = plt.subplots(2) + b, a = sp.signal.butter(2, 15 / (0.5 * imu_freq_Hz), btype = 'lowpass') + g = sp.signal.filtfilt(b, a, np.abs(np.diff(np.array(gyros), axis = 0)), axis = 0) + a = sp.signal.filtfilt(b, a, np.abs(np.diff(np.array(accs), axis = 0)), axis = 0) + axs[0].set_ylabel("dps") + xs_g = len(ecgs) / ecg_freq_Hz - (np.arange(len(gyros) - 1) / imu_freq_Hz)[::-1] + xs_a = len(ecgs) / ecg_freq_Hz - (np.arange(len(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 r_peaks: + peak_s = peak / ecg_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() - # imu_sparse += 1 - # if imu_sparse % 5 == 4: - # tt = int(5 * imu_freq_Hz) - # if len(gyros) > 480: - # gyros = gyros[-1600:] - # accs = accs[-1600:] - # else: - # return - # fig, axs = plt.subplots(2) - # b, a = sp.signal.butter(2, 20 / (0.5 * imu_freq_Hz), btype = 'lowpass') - # g = sp.signal.filtfilt(b,a,np.abs(np.diff(np.array(gyros), axis = 0)), axis = 0) - # a = sp.signal.filtfilt(b,a,np.abs(np.diff(np.array(accs), axis = 0)), axis = 0) - # #np.save("gyros.npy", g) - # #np.save("accs.npy", a) - # a -= np.mean(a, axis = 0).reshape(1,3) - # axs[0].set_ylabel("dps") - # axs[0].plot(g[-tt:,0]) - # axs[0].plot(g[-tt:,1]) - # axs[0].plot(g[-tt:,2]) - # axs[1].set_ylabel("g") - # axs[1].plot(a[-tt:,0]) - # axs[1].plot(a[-tt:,1]) - # axs[1].plot(a[-tt:,2]) - # plt.savefig("acc_gyro.png") - # plt.close() + if len(r_peaks) > 30: + fig, axs = plt.subplots(2) + for i in range(1,30): + xs_g = len(ecgs) / ecg_freq_Hz - (np.arange(len(gyros)) / imu_freq_Hz)[::-1] + ys = np.interp(np.linspace(r_peaks[-i] / ecg_freq_Hz - 1.0, r_peaks[-i] / ecg_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(r_peaks[-i] / ecg_freq_Hz - 1.0, r_peaks[-i] / ecg_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(r_peaks[-i] / ecg_freq_Hz - 1.0, r_peaks[-i] / ecg_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(ecgs) / ecg_freq_Hz - (np.arange(len(accs)) / imu_freq_Hz)[::-1] + ys = np.interp(np.linspace(r_peaks[-i] / ecg_freq_Hz - 1.0, r_peaks[-i] / ecg_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(r_peaks[-i] / ecg_freq_Hz - 1.0, r_peaks[-i] / ecg_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(r_peaks[-i] / ecg_freq_Hz - 1.0, r_peaks[-i] / ecg_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) -ecgs = [] -t1s = [] -t2s = [] -strains = [] + plt.savefig("accs_ensemble.png") + plt.close() + + +ecgs = np.zeros(0) +t1s = np.zeros(0) +t2s = np.zeros(0) +strains = np.zeros(0) ts = [] adc_sparse = 0 # Should be running at 488Hz -freq_Hz = 488.28125 +ecg_freq_Hz = 488.28125 max_hr_Hz = 120 / 60 last_adc_graph = 0 r_peaks = [] +f_pwrs = [] def process_adc(d, t): - global ecgs, t1s, t2s, strains, adc_sparse, ts, last_adc_graph, r_peaks + global ecgs, t1s, t2s, strains, adc_sparse, ts, last_adc_graph, r_peaks, f_pwrs 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': - ts.append((1 / 2000) * int.from_bytes(block[:4], byteorder = 'little', signed = True)) + #if e['name'] == b't': + # ts.append((1 / 2000) * int.from_bytes(block[:4], byteorder = 'little', signed = True)) if e['name'] == b'ekg_readings_cnts[50]': - ecgs = ecgs + [(2.4 / (1<<24)) * int.from_bytes(block[4 * i : 4 * i + 4], byteorder = 'little', signed = True) for i in range(50)] + ecgs = np.concatenate((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]': - strains = strains + [(2.4 / (1<<24)) * int.from_bytes(block[4 * i : 4 * i + 4], byteorder = 'little', signed = True) for i in range(5)] + strains = np.concatenate((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]': - t1s = t1s + [(2.4 / (1<<24)) * int.from_bytes(block[4 * i : 4 * i + 4], byteorder = 'little', signed = True) for i in range(5)] + t1s = np.concatenate((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]': - t2s = t2s + [(2.4 / (1<<24)) * int.from_bytes(block[4 * i : 4 * i + 4], byteorder = 'little', signed = True) for i in range(5)] - - if True and time.time() - last_adc_graph > 0.5:# and len(ecgs) > 500: - last_adc_graph = time.time() - #ecgs = ecgs[- int(5 * 60 * freq_Hz):] - #b, a = sp.signal.bessel(2, [1 / (0.5 * freq_Hz), 120 / (0.5 * freq_Hz)], btype = 'bandpass') + t2s = np.concatenate((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)]))) + if len(ecgs) < ecg_freq_Hz * 100: + return + + last_adc_graph += 1 + # About every second + if last_adc_graph % 10 == 9: #True and time.time() - last_adc_graph > 0.5:# and len(ecgs) > 500: + #last_adc_graph = time.time() + #ecgs = ecgs[- int(5 * 60 * ecg_freq_Hz):] + #b, a = sp.signal.bessel(2, [1 / (0.5 * ecg_freq_Hz), 120 / (0.5 * ecg_freq_Hz)], btype = 'bandpass') #ecgs_ = sp.signal.filtfilt(b, a, ecgs) if len(r_peaks) > 0: @@ -154,45 +229,66 @@ def process_adc(d, t): last_ind = start for ind in inds: - if ind - last_ind < freq_Hz / max_hr_Hz: + if ind - last_ind < ecg_freq_Hz / max_hr_Hz: continue - region_start = ind - int(0.1 * freq_Hz) - region_end = ind + int(0.1 * freq_Hz) + region_start = ind - int(0.1 * ecg_freq_Hz) + region_end = ind + int(0.1 * ecg_freq_Hz) peak = region_start + np.argmax(ecgs[region_start : region_end]) r_peaks.append(peak) last_ind = ind - fig, axs = plt.subplots(3, 1, height_ratios = [1,1,1]) - axs[0].plot(np.arange(len(ecgs))[int(- 5 * freq_Hz):] / freq_Hz, ecgs[int(- 5 * freq_Hz):], 'k.', linestyle='--', alpha = 0.5) - + fig, axs = plt.subplots(2, 1, height_ratios = [1,2]) + axs[0].plot(np.arange(len(ecgs))[int(- 5 * ecg_freq_Hz):] / ecg_freq_Hz, ecgs[int(- 5 * ecg_freq_Hz):], 'k.', linestyle='--', alpha = 0.5) + axs[0].set_ylim(np.amin(ecgs[int(- 5 * ecg_freq_Hz):]) - 0.001, np.amax(ecgs[int(- 5 * ecg_freq_Hz):]) + 0.001) strain_ax = axs[1].twinx() - strain_ax.plot(np.arange(len(t2s)) * 10 / freq_Hz, t2s, 'g', alpha = 0.5) + strain_ax.plot(np.arange(len(t2s)) * 10 / ecg_freq_Hz, t2s, color = 'orange', alpha = 0.5) strain_ax.set_yticks([]) axs[1].plot([],'k.',label = 'R-R int') #alsho 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 + axs[1].plot([],'b.',label = '|delta(R-R int)|') axs[1].plot([],'g',label = 'chest') for peak in r_peaks: - axs[0].plot(peak / freq_Hz, ecgs[peak], 'ro') + axs[0].plot(peak / ecg_freq_Hz, ecgs[peak], 'ro') - #fn = sp.interpolate.CubicSpline(0.5 * (r_peaks[1:] + r_peaks[:-1]), np.diff(r_peaks) - # Xs = np.linspace(r_peaks[0], r_peaks[-1], 1000) - #np.log(np.abs(np.fft.rfft())) + # Poincare + #R = np.diff(r_peaks[-180:]) / freq_Hz + #axs[2].plot(R[:-1], R[1:], 'k.', alpha = 0.2) - axs[1].plot(np.array(r_peaks[:-1]) / freq_Hz, 1000 * np.diff(r_peaks) / freq_Hz, 'k.') + # 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))) - axs[0].set_xlim((len(ecgs) / freq_Hz) - 5, len(ecgs) / freq_Hz) + axs[1].plot(np.array(r_peaks[:-1]) / ecg_freq_Hz, 1000 * np.diff(r_peaks) / ecg_freq_Hz, 'k.', alpha = 0.25) + + axs[0].set_xlim((len(ecgs) / ecg_freq_Hz) - 5, len(ecgs) / ecg_freq_Hz) hrv_ax = axs[1].twinx() - hrv = np.abs(1000 * np.diff(np.diff(r_peaks)) / freq_Hz) - hrv_ax.plot(np.array(r_peaks[2:]) / freq_Hz, hrv, 'b.') + hrv = np.abs(1000 * np.diff(np.diff(r_peaks)) / ecg_freq_Hz) + hrv_ax.plot(np.array(r_peaks[2:]) / ecg_freq_Hz, hrv, 'b.', alpha = 0.25) if len(r_peaks) > 10: N = int(np.ceil(len(r_peaks) / 30)) for i in range(N): - block = 1000 * np.array(r_peaks[30 * i : min(len(r_peaks), 30 * i + 30)]) / freq_Hz + block = 1000 * np.array(r_peaks[30 * i : min(len(r_peaks), 30 * i + 30)]) / ecg_freq_Hz 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) @@ -216,7 +312,17 @@ def process_adc(d, t): #if (len(ecgs) / freq_Hz) > 6: # plt.show() # quit() - plt.savefig("hrv_biofeedback.png") + plt.savefig("hrv_biofeedback.png", dpi = 150) + plt.close() + + fig, ax = plt.subplots(1) + for i in range(1,30): + dat = ecgs[r_peaks[-i] - int(0.3 * ecg_freq_Hz) : + r_peaks[-i] + int(0.9 * ecg_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() def make_graphs(): diff --git a/code/l452_code/packet_parser_log.py b/code/l452_code/packet_parser_log.py index 1d90ed3..da2e4a7 100644 --- a/code/l452_code/packet_parser_log.py +++ b/code/l452_code/packet_parser_log.py @@ -3,7 +3,7 @@ import time import matplotlib.pyplot as plt -log = open('00140426.LOG','rb').read() +log = open('longish.LOG','rb').read() types = packet_parser_helpers.get_type_list(packet_parser_helpers.packet_definitions)