From 293e7e85a7561477adf32777354b47c54a5eb3af Mon Sep 17 00:00:00 2001 From: ggw Date: Tue, 26 May 2026 08:11:29 -0700 Subject: [PATCH] updated analsysi --- code/l452_code/packet_parser_ble.py | 6 +- code/l452_code/packet_parser_helpers.py | 198 ++++++++++++++---------- code/l452_code/packet_parser_log.py | 6 +- 3 files changed, 121 insertions(+), 89 deletions(-) diff --git a/code/l452_code/packet_parser_ble.py b/code/l452_code/packet_parser_ble.py index 6cded59..565ed6d 100644 --- a/code/l452_code/packet_parser_ble.py +++ b/code/l452_code/packet_parser_ble.py @@ -10,13 +10,15 @@ from packet_parser_helpers import * async def scan(): return await BleakScanner.find_device_by_name("XX-STM32") +# 2 imu 4 adc 6 ppg async def connect(device): async with BleakClient(device) as client: + #await client.write_gatt_char(TX_UUID, b'2', response=True) await client.write_gatt_char(TX_UUID, b'4', response=True) - await client.write_gatt_char(TX_UUID, b'6', response=True) + #await client.write_gatt_char(TX_UUID, b'6', response=True) await client.stop_notify(RX_UUID) await client.start_notify(RX_UUID, cb) - await asyncio.sleep(60) + await asyncio.sleep(600) #await client.write_gatt_char(TX_UUID, b'2', response=True) #await client.write_gatt_char(TX_UUID, b'S', response=True) await client.stop_notify(RX_UUID) diff --git a/code/l452_code/packet_parser_helpers.py b/code/l452_code/packet_parser_helpers.py index e533fdc..2d7339f 100644 --- a/code/l452_code/packet_parser_helpers.py +++ b/code/l452_code/packet_parser_helpers.py @@ -1,6 +1,7 @@ import numpy as np import matplotlib.pyplot as plt import scipy as sp +import time 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 @@ -46,23 +47,24 @@ def process_ppg(d, t): 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() + 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() accs = [] gyros = [] imu_sparse = 0 +imu_freq_Hz = 240 def process_imu(d, t): global accs, gyros, imu_sparse for e in t['elements']: @@ -85,23 +87,27 @@ def process_imu(d, t): # imu_sparse += 1 # if imu_sparse % 5 == 4: - # if len(gyros) > 1600: + # tt = int(5 * imu_freq_Hz) + # if len(gyros) > 480: # gyros = gyros[-1600:] # accs = accs[-1600:] + # else: + # return # fig, axs = plt.subplots(2) - # g = np.array(gyros) - # a = np.array(accs) + # 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[:,0]) - # axs[0].plot(g[:,1]) - # axs[0].plot(g[:,2]) + # 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[:,0]) - # axs[1].plot(a[:,1]) - # axs[1].plot(a[:,2]) + # 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() @@ -113,9 +119,11 @@ ts = [] adc_sparse = 0 # Should be running at 488Hz freq_Hz = 488.28125 -max_hr_Hz = 240 / 60 +max_hr_Hz = 120 / 60 +last_adc_graph = 0 +r_peaks = [] def process_adc(d, t): - global ecgs, t1s, t2s, strains, adc_sparse, ts + global ecgs, t1s, t2s, strains, adc_sparse, ts, last_adc_graph, r_peaks for e in t['elements']: block = d[e['offset']:e['offset'] + e['size']] element_size = int(len(block) / e['n_elements']) @@ -124,72 +132,95 @@ def process_adc(d, t): 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)] - # adc_sparse += 1 - # if adc_sparse % 10 == 9: - # if len(ecgs) > 200: - # ecgs = ecgs[-4 * 4096:] - # t1s = t1s[-4096:] - # t2s = t2s[-4096:] - # strains = strains[-4096:] - # fig, axs = plt.subplots(2) - # axs[0].set_title("ECG") - # #axs[2].plot(np.array(ts[-200:-1]), 1000 * np.diff(ts[-200:]),'k.',linestyle='--') - # #if len(ts) > 200: - # # axs[2].set_title(str((ts[-1] - ts[-200]) / 200)) - # #axs[2].set_xlabel('S') - # #axs[2].set_ylabel('mS') - # #axs[1].set_title("Strain") - # #axs[2].set_title("oT") - # #axs[3].set_title("iT") - # # 0.25 * 244 Hz - # #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) - # #beats = np.where(np.diff(ecgs_) < -0.0003)[0] - # #last_beat = -1000000 - # #beats_ = [] - # #for beat in beats: - # # if (beat - last_beat) / 488 > 1 / (max_hr_Hz): - # # last_beat = beat - # # beats_.append(beat) - # axs[0].plot(ecgs, 'k') - # r_peaks = sp.signal.find_peaks(ecgs, height = None, threshold = None, distance = freq_Hz / max_hr_Hz, width = 5, prominence = 0.001)[0] - # all_peaks = sp.signal.find_peaks(ecgs, height = None, threshold = None, width = 10, prominence = 0.0001)[0] - # # for peak in all_peaks: - # # ls = [e for e in r_peaks if peak - e < 0 and peak - e > -0.2 * freq_Hz] - # # if len(ls) >= 1: - # # axs.plot(peak, ecgs[peak], 'bo') # P peaks - # # for peak in all_peaks: - # # ls = [e for e in r_peaks if peak - e > 0 and peak - e < 0.4 * freq_Hz] - # # if len(ls) >= 1: - # # axs.plot(peak, ecgs[peak], 'yo') # T peaks - # for peak in r_peaks: - # axs[0].plot(peak, ecgs[peak], 'ro') - # if len(r_peaks) > 5: - # hr_bpm = np.mean(60 * freq_Hz / np.diff(r_peaks)) - # #hrv_ms = 1000 * np.std(np.diff(r_peaks)) / freq_Hz - # hrv_rmssd_ms = 1000. * np.power(np.mean(np.power(np.diff(np.diff(r_peaks / freq_Hz)), 2)), 0.5) - # axs[0].set_title(f"ECG HR:{np.round(hr_bpm)}BPM HRV (RMSSD):{np.round(hrv_rmssd_ms)}ms") - # hist = np.histogram(np.diff(r_peaks), bins = 8) - # axs[1].plot(0.5 * (hist[1][1:] + hist[1][:-1]), hist[0]) - - # #axs[1].plot(strains) - # #axs[2].plot(t1s) - # #axs[3].plot(t2s) - # fig.tight_layout() - # plt.savefig("adcs.png") - # plt.close() - - # fig, axs = plt.subplots() - # axs.plot(ts[1:], np.diff(ts), 'k.') - # plt.savefig("adc_ts.png") 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)] 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)] 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') + #ecgs_ = sp.signal.filtfilt(b, a, ecgs) + if len(r_peaks) > 0: + start = r_peaks[-1] + else: + start = 0 + v = np.convolve(np.abs(np.diff(ecgs[start:])), np.ones(5) / 5, mode = 'valid') + inds = [start + e for e in np.where(v > 0.0002)[0]] + + last_ind = start + for ind in inds: + if ind - last_ind < freq_Hz / max_hr_Hz: + continue + region_start = ind - int(0.1 * freq_Hz) + region_end = ind + int(0.1 * 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) + + strain_ax = axs[1].twinx() + strain_ax.plot(np.arange(len(t2s)) * 10 / freq_Hz, t2s, 'g', 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') + + #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())) + + axs[1].plot(np.array(r_peaks[:-1]) / freq_Hz, 1000 * np.diff(r_peaks) / freq_Hz, 'k.') + + axs[0].set_xlim((len(ecgs) / freq_Hz) - 5, len(ecgs) / 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.') + + 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 + 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 = 'orange', 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[0].set_xlabel("t (s)") + axs[1].set_xlabel("t (s)") + axs[0].set_ylabel("V (mV)") + axs[1].set_ylabel("beat delta T (ms)") + hrv_ax.set_ylabel("BPS delta T (ms)", color = 'blue') + axs[1].legend(loc='upper left', framealpha = 1) + hrv_ax.set_axisbelow(True) + plt.tight_layout() + #if (len(ecgs) / freq_Hz) > 6: + # plt.show() + # quit() + plt.savefig("hrv_biofeedback.png") + plt.close() + 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') @@ -213,7 +244,6 @@ def make_graphs(): #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() @@ -236,7 +266,7 @@ def make_graphs(): #axs[1][1].plot(t1s) #axs[1][1].plot(t2s) - plt.show() + #plt.show() def read_and_process(types, cons, size): index = 0 diff --git a/code/l452_code/packet_parser_log.py b/code/l452_code/packet_parser_log.py index de3b165..1d90ed3 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('5-23-2026-first-10-min.LOG','rb').read() +log = open('00140426.LOG','rb').read() types = packet_parser_helpers.get_type_list(packet_parser_helpers.packet_definitions) @@ -17,7 +17,7 @@ while index < len(log): #time.sleep(0.1) except IndexError: break - index += 4 + size + index += 4096 #4 + size -packet_parser_helpers.make_graphs() +#packet_parser_helpers.make_graphs()