diff --git a/code/l452_code/packet_parser_helpers.py b/code/l452_code/packet_parser_helpers.py index e83a45d..f2484c7 100644 --- a/code/l452_code/packet_parser_helpers.py +++ b/code/l452_code/packet_parser_helpers.py @@ -4,6 +4,8 @@ import matplotlib.pyplot as plt import scipy as sp import time +# Try using filtfilt instead of filtfilt + 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 @@ -233,7 +235,7 @@ class adc_vars(): r_peaks = [] - #dead_times = [] + on_off_times = [0] def process_adc(d, t): global adc_vars @@ -258,16 +260,17 @@ def process_adc(d, t): 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 = adc_vars.r_peaks[-1] - else: - start = 0 - - # TODO: add dead time recognition - #if np.amax(adc_vars.ecgs[start:]) - np.amin(adc_vars.ecgs[start:start + 60]) > 0.005: - + 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]] @@ -275,8 +278,21 @@ def update_ecg_peaks(): 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 @@ -300,45 +316,53 @@ def adc_graphs(): 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) + # 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() + # plt.tight_layout() + # plt.savefig("hrv_biofeedback.png", dpi = 150) + # plt.close() - return + # 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') - filtered_strain = sp.signal.filtfilt(b_breath, a_breath, adc_vars.t2s[-int(RR_lookback_s * adc_vars.freq_Hz * 0.1):]) + + #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] - 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_ax = axs[1].twinx() - 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) + + 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.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)]) / adc_vars.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 = 'pink', alpha = 0.5) + 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)]) / adc_vars.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 = 'pink', alpha = 0.5) hrv_ax.minorticks_on() hrv_ax.grid(alpha = 0.5) @@ -362,7 +386,7 @@ def adc_graphs(): 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) :