added dead time feature for ECG

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]
+        start = max(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:
-
     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)
+    # axs[1].plot(ts[:-1], np.diff(adc_vars.ecgs[-ecg_points:]), 'k.', linestyle='--', alpha = 0.25)
-    for peak in adc_vars.r_peaks:
+    # for peak in adc_vars.r_peaks:
-        if peak < ts[-1] * adc_vars.freq_Hz and peak > ts[0] * adc_vars.freq_Hz:
+    #     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].axvline(peak / adc_vars.freq_Hz, color = 'red', alpha = 0.1)
     
 
-    plt.tight_layout()
+    # plt.tight_layout()
-    plt.savefig("hrv_biofeedback.png", dpi = 150)
+    # plt.savefig("hrv_biofeedback.png", dpi = 150)
-    plt.close()
+    # 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:
+        if len(r_peaks_) > 10:
-        N = int(np.ceil(len(r_peaks) / 30))
+            N = int(np.ceil(len(r_peaks_) / 30))
-        for i in range(N):
+            for i in range(N):
-            block = 1000 * np.array(r_peaks[30 * i : min(len(r_peaks), 30 * i + 30)]) / adc_vars.freq_Hz
+                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)
+                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))
+                sdnn = np.std(np.diff(block))
-            pNN50 = 100 * np.mean(np.abs(np.diff(np.diff(block))) > 50)
+                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], [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], [sdnn, sdnn], 'm', alpha = 0.5)
-            hrv_ax.plot([block[0] / 1000, block[-1] / 1000], [pNN50, pNN50], color = 'pink', 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) :