ggw/chest_strap
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

updated analsysi

Browse Source
This commit is contained in:
ggw
2026-05-26 08:11:29 -07:00
parent 7d7297b12f
commit 293e7e85a7
3 changed files with 121 additions and 89 deletions
Download Patch File Download Diff File Expand all files Collapse all files
code/l452_code/packet_parser_ble.py
+4 -2
View File
@@ -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)
code/l452_code/packet_parser_helpers.py
+114 -84
View File
@@ -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
code/l452_code/packet_parser_log.py
+3 -3
View File
@@ -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()