diff --git a/phys2denoise/metrics/chest_belt.py b/phys2denoise/metrics/chest_belt.py index e8dafea..85be77c 100644 --- a/phys2denoise/metrics/chest_belt.py +++ b/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. diff --git a/phys2denoise/tests/conftest.py b/phys2denoise/tests/conftest.py new file mode 100644 index 0000000..1022b3e --- /dev/null +++ b/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 + diff --git a/phys2denoise/tests/test_rvt.py b/phys2denoise/tests/test_rvt.py new file mode 100644 index 0000000..7923ac6 --- /dev/null +++ b/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 diff --git a/setup.cfg b/setup.cfg index 031f8eb..a2f20fe 100644 --- a/setup.cfg +++ b/setup.cfg @@ -49,6 +49,7 @@ test = %(style)s pytest >=5.3 pytest-cov + peakdet coverage devtools = pre-commit