Add RVT algorithm (Birn et al. 2006)

phys2denoise/metrics/chest_belt.py

@@ -5,6 +5,7 @@
 from scipy.ndimage.filters import convolve1d
 from scipy.signal import detrend
 from scipy.stats import zscore
+from scipy.interpolate import interp1d
 
 mpl.use("TkAgg")
 import matplotlib.pyplot as plt
@@ -14,6 +15,68 @@
 from . import utils
 
 
+def rvt(belt_ts, peaks, troughs, samplerate, lags=(0, 4, 8, 12)):
+    """
+    Implement the Respiratory Variance over Time (Birn et al. 2006).
+
+    Procedural choices influenced by RetroTS
+
+    Parameters
+    ----------
+    belt_ts: array_like
+        respiratory belt data - samples x 1
+    peaks: array_like
+        peaks found by peakdet algorithm
+    troughs: array_like
+        troughs found by peakdet algorithm
+    samplerate: float
+        sample rate in hz of respiratory belt
+    lags: tuple
+        lags in seconds of the RVT output. Default is 0, 4, 8, 12.
+
+    Outputs
+    -------
+    rvt: array_like
+        calculated RVT and associated lags.
+    """
+    timestep = 1 / samplerate
+    # respiration belt timing
+    time = np.arange(0, len(belt_ts) * timestep, timestep)
+    peak_vals = belt_ts[peaks]
+    trough_vals = belt_ts[troughs]
+    peak_time = time[peaks]
+    trough_time = time[troughs]
+    mid_peak_time = (peak_time[:-1] + peak_time[1:]) / 2
+    period = peak_time[1:] - peak_time[:-1]
+    # interpolate peak values over all timepoints
+    peak_interp = interp1d(
+        peak_time, peak_vals, bounds_error=False, fill_value="extrapolate"
+    )(time)
+    # interpolate trough values over all timepoints
+    trough_interp = interp1d(
+        trough_time, trough_vals, bounds_error=False, fill_value="extrapolate"
+    )(time)
+    # interpolate period over  all timepoints
+    period_interp = interp1d(
+        mid_peak_time, period, bounds_error=False, fill_value="extrapolate"
+    )(time)
+    # full_rvt is (peak-trough)/period
+    full_rvt = (peak_interp - trough_interp) / period_interp
+    # calculate lags for RVT
+    rvt_lags = np.zeros((len(full_rvt), len(lags)))
+    for ind, lag in enumerate(lags):
+        start_index = np.argmin(np.abs(time - lag))
+        temp_rvt = np.concatenate(
+            (
+                np.full((start_index), full_rvt[0]),
+                full_rvt[: len(full_rvt) - start_index],
+            )
+        )
+        rvt_lags[:, ind] = temp_rvt
+
+    return rvt_lags
+
+
 @due.dcite(references.POWER_2018)
 def rpv(resp, window=6):
     """Calculate respiratory pattern variability.

phys2denoise/tests/conftest.py

@@ -0,0 +1,14 @@
+import numpy as np
+import pytest
+
+
+@pytest.fixture(scope="module")
+def fake_phys():
+    f = 0.3
+    fs = 62.5  # sampling rate
+    t = 300
+    samples = np.arange(t * fs) / fs
+    noise = np.random.normal(0, 0.5, len(samples))
+    fake_phys = 10 * np.sin(2 * np.pi * f * samples) + noise
+    return fake_phys
+

phys2denoise/tests/test_rvt.py

@@ -0,0 +1,19 @@
+import peakdet
+from phys2denoise.metrics.chest_belt import rvt
+
+
+def test_peakdet(fake_phys):
+    phys = peakdet.Physio(fake_phys, fs=62.5)
+    phys = peakdet.operations.filter_physio(phys, cutoffs=3, method="lowpass")
+    phys = peakdet.operations.peakfind_physio(phys)
+    assert phys.troughs is not None
+    assert phys.peaks is not None
+
+
+def test_rvt(fake_phys):
+    phys = peakdet.Physio(fake_phys, fs=62.5)
+    phys = peakdet.operations.filter_physio(phys, cutoffs=3, method="lowpass")
+    phys = peakdet.operations.peakfind_physio(phys)
+    r = rvt(phys.data, phys.peaks, phys.troughs, samplerate=phys.fs)
+    assert r is not None
+    assert len(r) == 18750

setup.cfg

@@ -49,6 +49,7 @@ test =
     %(style)s
     pytest >=5.3
     pytest-cov
+    peakdet
     coverage
 devtools =
     pre-commit