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']