Develop new HRV metrics

examples/example_03_ecg.py

@@ -167,6 +167,23 @@
 #      the square root of the mean of the squared differences between successive RR intervals.
 #      Conceptually, it is similar to a second derivative of the RR intervals (if RR intervals are considered as a
 #      first derivative). RMSSD is very sensitive to outliers, which can artificially increase its value.
+#    * `HRV_pNN50`: (units = `%`) Percentage of successive RR interval differences greater than 50 ms.
+#      It is sensitive to outliers and ectopic beats.
+#    * `HRV_pNN20`: (units = `%`) Percentage of successive RR interval differences greater than 20 ms.
+#      It is sensitive to outliers and ectopic beats.
+#    * `HRV_SD1`: (units = `ms`) Standard deviation of (RR[n] - RR[n+1]) / sqrt(2).
+#      This metric is computed from successive differences between RR intervals.
+#      It is mathematically related to RMSSD by a constant scaling factor.
+#    * `HRV_SD2`: (units = `ms`) Standard deviation of (RR[n] + RR[n+1]) / sqrt(2).
+#      This metric is computed from successive sums of RR intervals.
+#    * `HRV_SD1SD2`: (units = `AU`) Ratio of `HRV_SD1` to `HRV_SD2`.
+#      This dimensionless metric compares variability derived from successive differences
+#      to variability derived from successive sums of RR intervals.
+#    * `HRV_S`: (units = `ms^2`) Area of the ellipse defined by `HRV_SD1` and `HRV_SD2`
+#      (computed as pi * SD1 * SD2). This scalar combines both SD1 and SD2 into a single
+#      dispersion measure.
+#    * `HRV_ShannonEntropy`: (units = `AU`) Shannon entropy of the probability distribution
+#      of RR intervals estimated using histogram binning.
 #
 # Some of these metrics can be visualized on the RR interval distribution below, which provides
 # a simple way to identify potential outliers in the detection.
@@ -177,6 +194,10 @@
 #    3) By interpreting results using robust metrics such as `HRV_Median`, `HRV_Mad`, or `HRV_MCV`
 #
 # While these three steps can reduce the impact of outliers, careful ECG data recording is no substitute for quality optimization.
+#
+# WARNING : Many of the metrics (e.g. `HRV_Mean`, `HRV_SD`, `HRV_RMSSD`, `HRV_pNN50`, `HRV_pNN20`, `HRV_SD1`, `HRV_SD2`, `HRV_S`, `HRV_ShannonEntropy`) are highly sensitive to outliers and artifacts in RR interval sequences.
+# Furthermore, the physiological interpretability of these metrics remains limited and should be considered with caution.
+
 
 
 

physio/ecg.py

@@ -4,7 +4,7 @@
 
 import scipy.interpolate
 
-from .tools import detect_peak, compute_median_mad
+from .tools import detect_peak, compute_median_mad, sd1_sd2_hrv, Shannon_Entropy
 from .preprocess import preprocess
 from .parameters import get_ecg_parameters, recursive_update
 
@@ -137,6 +137,13 @@ def compute_ecg_metrics(ecg_peaks, min_interval_ms=500., max_interval_ms=2000.,
     * HRV_MCV : HRV_Mad / HRV_Median = MAD Coefficient of Variation = Robust version of Coefficient of Variation
     * HRV_Asymmetry = HRV_Median - HRV_Mean = Difference between Median and Mean that diverges from 0 in case of outliers / non normal distribution of RR intervals.
     * HRV_RMSSD = Root-Mean Square of Successive Differences ~ like a 2nd derivative of RR intervals. Sensitive to fine variations of RR intervals but also very sensitive to outliers.
+    * HRV_pNN50 : Percentage of successive RR interval differences greater than 50 ms. Very sensitive to outliers.
+    * HRV_pNN20 : Percentage of successive RR interval differences greater than 20 ms. Very sensitive to outliers.
+    * HRV_SD1 : Standard deviation of (RR[n] - RR[n+1]) / sqrt(2). Quantifies variability based on successive differences between RR intervals. Directly related to RMSSD by a constant scaling factor.
+    * HRV_SD2 : Standard deviation of (RR[n] + RR[n+1]) / sqrt(2). Quantifies variability along the sum of consecutive RR intervals.
+    * HRV_SD1SD2 : Ratio HRV_SD1 / HRV_SD2. Compares the dispersion of successive differences to the dispersion of successive sums of RR intervals. Dimensionless metric derived from SD1 and SD2
+    * HRV_S : Area of the ellipse defined by SD1 and SD2 (pi * SD1 * SD2). Derived scalar combining SD1 and SD2. Expresses dispersion of RR intervals in the SD1 to SD2 representation.
+    * HRV_ShannonEntropy : Shannon entropy computed from histogram-based probability distribution of RR intervals. Could measures uncertainty of the RR interval distribution.
 
     These metrics are a bit more robust than others toolboxes because are computed after a cleaning of RR intervals based on min and max intervals as set.
 
@@ -174,12 +181,7 @@ def compute_ecg_metrics(ecg_peaks, min_interval_ms=500., max_interval_ms=2000.,
 
 
     delta_ms = np.diff(peak_ms)
-
-    # keep = delta_ms < max_interval_ms
-
-    # delta_ms = delta_ms[keep]
-
-
+
     metrics = pd.Series(dtype=float)
 
     metrics['HRV_Mean'] = np.nanmean(delta_ms)
@@ -190,10 +192,22 @@ def compute_ecg_metrics(ecg_peaks, min_interval_ms=500., max_interval_ms=2000.,
     metrics['HRV_Asymmetry'] = metrics['HRV_Median'] - metrics['HRV_Mean']
 
 
-    # TODO
-    metrics['HRV_RMSSD'] = np.sqrt(np.nanmean(np.diff(delta_ms)**2))
+    # TO DO : more robust outlier mask while working on diff_delta_ms
+    diff_delta_ms = np.diff(delta_ms)
+    metrics['HRV_RMSSD'] = np.sqrt(np.nanmean(diff_delta_ms**2))
+    metrics["HRV_pNN50"] = np.sum(np.abs(diff_delta_ms) > 50) / (len(diff_delta_ms) + 1) * 100
+    metrics["HRV_pNN20"] = np.sum(np.abs(diff_delta_ms) > 20) / (len(diff_delta_ms) + 1) * 100
+
+    # poincare
+    sd1, sd2, sd1_sd2, s = sd1_sd2_hrv(delta_ms)
+    metrics['HRV_SD1'] = sd1
+    metrics['HRV_SD2'] = sd2
+    metrics['HRV_SD1SD2'] = sd1_sd2
+    metrics['HRV_S'] = s
+
+    # entropy 
+    metrics['HRV_ShannonEntropy'] = Shannon_Entropy(delta_ms)
 
-    # return pd.DataFrame(metrics).T
     return metrics
 
 

physio/tools.py

@@ -129,8 +129,31 @@ def crosscorrelogram(a, b, bins):
     count, bins = np.histogram(diff, bins)
     return count, bins
 
+def sd1_sd2_hrv(rr_intervals):
+    rri_n = rr_intervals[:-1]
+    rri_plus = rr_intervals[1:]
 
+    x1 = (rri_n - rri_plus) / np.sqrt(2) 
+    x2 = (rri_n + rri_plus) / np.sqrt(2)
+    sd1 = np.std(x1, ddof=1)
+    sd2 = np.std(x2, ddof=1)
+    sd1_sd2 = sd1 / sd2
+    s = np.pi * sd1 * sd2 # Area of ellipse described by SD1 and SD2
+    return sd1, sd2, sd1_sd2, s
 
+def Shannon_Entropy(signal):
+    signal = np.asarray(signal)
+
+    n_bins = int(np.sqrt(len(signal)))
+
+    hist, _ = np.histogram(signal, bins=n_bins, density=False)
+
+    prob = hist / np.sum(hist)
+    prob = prob[prob > 0]
+
+    entropy = -np.sum(prob * np.log(prob))
+
+    return entropy
 
 # convolution stuff to keep in mind
 # sweep = np.arange(-60, 60)