diff --git a/examples/example_03_ecg.py b/examples/example_03_ecg.py index b4c283b..a667e61 100644 --- a/examples/example_03_ecg.py +++ b/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. + diff --git a/physio/ecg.py b/physio/ecg.py index d99e769..e8e2393 100644 --- a/physio/ecg.py +++ b/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 diff --git a/physio/tools.py b/physio/tools.py index 2ab8478..59cdead 100644 --- a/physio/tools.py +++ b/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)