diff --git a/.all-contributorsrc b/.all-contributorsrc
index 46f9a58..1e8e6aa 100644
--- a/.all-contributorsrc
+++ b/.all-contributorsrc
@@ -110,6 +110,20 @@
         "review",
         "userTesting"
       ]
+    },
+    {
+      "login": "maestroque",
+      "name": "George Kikas",
+      "avatar_url": "https://avatars.githubusercontent.com/u/74024609?v=4",
+      "profile": "https://github.com/maestroque",
+      "contributions": [
+        "code",
+        "ideas",
+        "infra",
+        "docs",
+        "bug",
+        "review"
+      ]
     }
   ],
   "contributorsPerLine": 8,
diff --git a/.circleci/config.yml b/.circleci/config.yml
index fc36d30..dd45175 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -15,14 +15,14 @@ orbs:
 # Define a job to be invoked later in a workflow.
 # See: https://circleci.com/docs/2.0/configuration-reference/#jobs
 jobs:
-  test37: # This is the name of the job, feel free to change it to better match what you're trying to do!
+  test39: # This is the name of the job, feel free to change it to better match what you're trying to do!
     # These next lines defines a Docker executors: https://circleci.com/docs/2.0/executor-types/
     # You can specify an image from Dockerhub or use one of the convenience images from CircleCI's Developer Hub
     # A list of available CircleCI Docker convenience images are available here: https://circleci.com/developer/images/image/cimg/python
     # The executor is the environment in which the steps below will be executed - below will use a python 3.6.14 container
     # Change the version below to your required version of python
     docker:
-      - image: cimg/python:3.7
+      - image: cimg/python:3.9
     working_directory: /tmp/src/phys2denoise
     resource_class: medium
     # Checkout the code as the first step. This is a dedicated CircleCI step.
@@ -46,7 +46,7 @@ jobs:
           command: |
             pytest --cov=./phys2denoise
             mkdir /tmp/src/coverage
-            mv ./.coverage /tmp/src/coverage/.coverage.py37
+            mv ./.coverage /tmp/src/coverage/.coverage.py39
       - store_artifacts:
           path: /tmp/src/coverage
       # Persist the specified paths (workspace/echo-output) into the workspace for use in downstream job.
@@ -56,11 +56,11 @@ jobs:
           root: /tmp
           # Must be relative path from root
           paths:
-            - src/coverage/.coverage.py37
+            - src/coverage/.coverage.py39
 
-  test310:
+  test311:
     docker:
-      - image: cimg/python:3.10
+      - image: cimg/python:3.11
     working_directory: /tmp/src/phys2denoise
     resource_class: medium
     steps:
@@ -75,17 +75,17 @@ jobs:
           command: |
             pytest --cov=./phys2denoise
             mkdir /tmp/src/coverage
-            mv ./.coverage /tmp/src/coverage/.coverage.py310
+            mv ./.coverage /tmp/src/coverage/.coverage.py311
       - store_artifacts:
           path: /tmp/src/coverage
       - persist_to_workspace:
           root: /tmp
           paths:
-            - src/coverage/.coverage.py310
+            - src/coverage/.coverage.py311
 
   style_check:
     docker:
-      - image: cimg/python:3.7
+      - image: cimg/python:3.9
     working_directory: /tmp/src/phys2denoise
     resource_class: small
     steps:
@@ -105,7 +105,7 @@ jobs:
   merge_coverage:
     working_directory: /tmp/src/phys2denoise
     docker:
-      - image: cimg/python:3.10
+      - image: cimg/python:3.11
     resource_class: small
     steps:
       - attach_workspace:
@@ -133,13 +133,13 @@ workflows:
     # Inside the workflow, you define the jobs you want to run.
     jobs:
       - style_check
-      - test37:
+      - test39:
           requires:
             - style_check
-      - test310:
+      - test311:
           requires:
             - style_check
       - merge_coverage:
           requires:
-            - test37
-            - test310
+            - test39
+            - test311
diff --git a/docs/user_guide/metrics.rst b/docs/user_guide/metrics.rst
index e387600..779f31b 100644
--- a/docs/user_guide/metrics.rst
+++ b/docs/user_guide/metrics.rst
@@ -71,6 +71,10 @@ example shows how to compute the heart rate and the heart rate variability using
 
     from phys2denoise.metrics.chest_belt import respiratory_variance_time
 
-    # Given that the respiratory signal is stored in `data`, the peaks in `peaks`, the troughs in `troughs`
+    # Given that the respiratory signal is stored in `data` (which is not a physio.Physio instance), the peaks in `peaks`, the troughs in `troughs`
     # and the sample rate in `sample_rate`
-    _, rvt = respiratory_variance_time(data, peaks, troughs, sample_rate)
+    rvt = respiratory_variance_time(data, peaks, troughs, sample_rate)
+
+The computed respiratory variance time is stored in the variable ``rvt``. An internal check is performed to verify if the input data is a Physio object or not
+determining the appropriate output format. If the input is not a ``Physio`` object, the output will be a numpy array only containing the computed metric. Otherwise,
+the output will be a tuple with the updated Physio object and the computed metric.
diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py
index 40ebc47..5b7588d 100644
--- a/phys2denoise/metrics/cardiac.py
+++ b/phys2denoise/metrics/cardiac.py
@@ -1,14 +1,24 @@
 """Denoising metrics for cardio recordings."""
 import numpy as np
+from loguru import logger
+from physutils import io, physio
 
 from .. import references
 from ..due import due
 from .responses import crf
 from .utils import apply_function_in_sliding_window as afsw
-from .utils import convolve_and_rescale
-
-
-def _cardiac_metrics(card, peaks, samplerate, metric, window=6, central_measure="mean"):
+from .utils import convolve_and_rescale, return_physio_or_metric
+
+
+def _cardiac_metrics(
+    data,
+    metric,
+    peaks=None,
+    fs=None,
+    window=6,
+    central_measure="mean",
+    **kwargs,
+):
     """
     Compute cardiac metrics.
 
@@ -21,14 +31,14 @@ def _cardiac_metrics(card, peaks, samplerate, metric, window=6, central_measure=
 
     Parameters
     ----------
-    card : list or 1D numpy.ndarray
-        Timeseries of recorded cardiac signal
-    peaks : list or 1D numpy.ndarray
-        array of peak indexes for card.
-    samplerate : float
-        Sampling rate for card, in Hertz.
+    data : physutils.Physio, np.ndarray, or array-like object
+        Object containing the timeseries of the recorded cardiac signal
     metrics : "hbi", "hr", "hrv", string
         Cardiac metric(s) to calculate.
+    fs : float, optional if data is a physutils.Physio
+        Sampling rate of `data` in Hz
+    peaks : :obj:`numpy.ndarray`, optional if data is a physutils.Physio
+        Indices of peaks in `data`
     window : float, optional
         Size of the sliding window, in seconds.
         Default is 6.
@@ -71,8 +81,32 @@ def _cardiac_metrics(card, peaks, samplerate, metric, window=6, central_measure=
     Annual International Conference of the IEEE Engineering in Medicine and
     Biology Society (EMBC), doi: 10.1109/EMBC.2016.7591347.
     """
+    if isinstance(data, physio.Physio):
+        # Initialize physio object
+        data = physio.check_physio(data, ensure_fs=True, copy=True)
+    elif fs is not None and peaks is not None:
+        data = io.load_physio(data, fs=fs)
+        data._metadata["peaks"] = peaks
+    else:
+        raise ValueError(
+            """
+            To use this function you should either provide a Physio object
+            with existing peaks metadata (e.g. using the peakdet module), or
+            by providing the physiological data timeseries, the sampling frequency,
+            and the peak indices separately.
+            """
+        )
+    if data.peaks.size == 0:
+        raise ValueError(
+            """
+            Peaks must be a non-empty list.
+            Make sure to run peak detection on your physiological data first,
+            using the peakdet module, or other software of your choice.
+            """
+        )
+
     # Convert window to samples, but halves it.
-    halfwindow_samples = int(round(window * samplerate / 2))
+    halfwindow_samples = int(round(window * data.fs / 2))
 
     if central_measure in ["mean", "average", "avg"]:
         central_measure_operator = np.mean
@@ -85,14 +119,17 @@ def _cardiac_metrics(card, peaks, samplerate, metric, window=6, central_measure=
             f" {central_measure} is not a supported metric of centrality."
         )
 
-    idx_arr = np.arange(len(card))
+    idx_arr = np.arange(len(data))
     idx_min = afsw(idx_arr, np.min, halfwindow_samples)
     idx_max = afsw(idx_arr, np.max, halfwindow_samples)
 
-    card_met = np.empty_like(card)
+    card_met = np.empty_like(data)
     for n, i in enumerate(idx_min):
         diff = (
-            np.diff(peaks[np.logical_and(peaks >= i, peaks <= idx_max[n])]) / samplerate
+            np.diff(
+                data.peaks[np.logical_and(data.peaks >= i, data.peaks <= idx_max[n])]
+            )
+            / data.fs
         )
         if metric == "hbi":
             card_met[n] = central_measure_operator(diff) if diff.size > 0 else 0
@@ -109,13 +146,15 @@ def _cardiac_metrics(card, peaks, samplerate, metric, window=6, central_measure=
     card_met[np.isnan(card_met)] = 0.0
 
     # Convolve with crf and rescale
-    card_met = convolve_and_rescale(card_met, crf(samplerate), rescale="rescale")
+    card_met = convolve_and_rescale(card_met, crf(data.fs), rescale="rescale")
 
-    return card_met
+    return data, card_met
 
 
 @due.dcite(references.CHANG_CUNNINGHAM_GLOVER_2009)
-def heart_rate(card, peaks, samplerate, window=6, central_measure="mean"):
+@return_physio_or_metric()
+@physio.make_operation()
+def heart_rate(data, peaks=None, fs=None, window=6, central_measure="mean", **kwargs):
     """
     Compute average heart rate (HR) in a sliding window.
 
@@ -125,12 +164,12 @@ def heart_rate(card, peaks, samplerate, window=6, central_measure="mean"):
 
     Parameters
     ----------
-    card : list or 1D numpy.ndarray
-        Timeseries of recorded cardiac signal
-    peaks : list or 1D numpy.ndarray
-        array of peak indexes for card.
-    samplerate : float
-        Sampling rate for card, in Hertz.
+    data : physutils.Physio, np.ndarray, or array-like object
+        Object containing the timeseries of the recorded cardiac signal
+    fs : float, optional if data is a physutils.Physio
+        Sampling rate of `data` in Hz
+    peaks : :obj:`numpy.ndarray`, optional if data is a physutils.Physio
+        Indices of peaks in `data`
     window : float, optional
         Size of the sliding window, in seconds.
         Default is 6.
@@ -167,13 +206,25 @@ def heart_rate(card, peaks, samplerate, window=6, central_measure="mean"):
     Annual International Conference of the IEEE Engineering in Medicine and
     Biology Society (EMBC), doi: 10.1109/EMBC.2016.7591347.
     """
-    return _cardiac_metrics(
-        card, peaks, samplerate, metric="hrv", window=6, central_measure="mean"
+    data, hr = _cardiac_metrics(
+        data,
+        metric="hr",
+        peaks=peaks,
+        fs=fs,
+        window=window,
+        central_measure=central_measure,
+        **kwargs,
     )
 
+    return data, hr
+
 
 @due.dcite(references.PINHERO_ET_AL_2016)
-def heart_rate_variability(card, peaks, samplerate, window=6, central_measure="mean"):
+@return_physio_or_metric()
+@physio.make_operation()
+def heart_rate_variability(
+    data, peaks=None, fs=None, window=6, central_measure="mean", **kwargs
+):
     """
     Compute average heart rate variability (HRV) in a sliding window.
 
@@ -183,12 +234,12 @@ def heart_rate_variability(card, peaks, samplerate, window=6, central_measure="m
 
     Parameters
     ----------
-    card : list or 1D numpy.ndarray
-        Timeseries of recorded cardiac signal
-    peaks : list or 1D numpy.ndarray
-        array of peak indexes for card.
-    samplerate : float
-        Sampling rate for card, in Hertz.
+    data : physutils.Physio, np.ndarray, or array-like object
+        Object containing the timeseries of the recorded cardiac signal
+    fs : float, optional if data is a physutils.Physio
+        Sampling rate of `data` in Hz
+    peaks : :obj:`numpy.ndarray`, optional if data is a physutils.Physio
+        Indices of peaks in `data`
     window : float, optional
         Size of the sliding window, in seconds.
         Default is 6.
@@ -223,24 +274,36 @@ def heart_rate_variability(card, peaks, samplerate, window=6, central_measure="m
     Annual International Conference of the IEEE Engineering in Medicine and
     Biology Society (EMBC), doi: 10.1109/EMBC.2016.7591347.
     """
-    return _cardiac_metrics(
-        card, peaks, samplerate, metric="hrv", window=6, central_measure="std"
+    data, hrv = _cardiac_metrics(
+        data,
+        metric="hrv",
+        peaks=peaks,
+        fs=fs,
+        window=window,
+        central_measure=central_measure,
+        **kwargs,
     )
 
+    return data, hrv
+
 
 @due.dcite(references.CHEN_2020)
-def heart_beat_interval(card, peaks, samplerate, window=6, central_measure="mean"):
+@return_physio_or_metric()
+@physio.make_operation()
+def heart_beat_interval(
+    data, peaks=None, fs=None, window=6, central_measure="mean", **kwargs
+):
     """
     Compute average heart beat interval (HBI) in a sliding window.
 
     Parameters
     ----------
-    card : list or 1D numpy.ndarray
-        Timeseries of recorded cardiac signal
-    peaks : list or 1D numpy.ndarray
-        array of peak indexes for card.
-    samplerate : float
-        Sampling rate for card, in Hertz.
+    data : physutils.Physio, np.ndarray, or array-like object
+        Object containing the timeseries of the recorded cardiac signal
+    fs : float, optional if data is a physutils.Physio
+        Sampling rate of `data` in Hz
+    peaks : :obj:`numpy.ndarray`, optional if data is a physutils.Physio
+        Indices of peaks in `data`
     window : float, optional
         Size of the sliding window, in seconds.
         Default is 6.
@@ -272,12 +335,22 @@ def heart_beat_interval(card, peaks, samplerate, window=6, central_measure="mean
     .. [1] J. E. Chen et al., "Resting-state "physiological networks"", Neuroimage,
         vol. 213, pp. 116707, 2020.
     """
-    return _cardiac_metrics(
-        card, peaks, samplerate, metric="hbi", window=6, central_measure="mean"
+    data, hbi = _cardiac_metrics(
+        data,
+        metric="hbi",
+        peaks=peaks,
+        fs=fs,
+        window=window,
+        central_measure=central_measure,
+        **kwargs,
     )
 
+    return data, hbi
 
-def cardiac_phase(peaks, sample_rate, slice_timings, n_scans, t_r):
+
+@return_physio_or_metric()
+@physio.make_operation()
+def cardiac_phase(data, slice_timings, n_scans, t_r, peaks=None, fs=None, **kwargs):
     """Calculate cardiac phase from cardiac peaks.
 
     Assumes that timing of cardiac events are given in same units
@@ -285,27 +358,53 @@ def cardiac_phase(peaks, sample_rate, slice_timings, n_scans, t_r):
 
     Parameters
     ----------
-    peaks : 1D array_like
-        Cardiac peak times, in seconds.
-    sample_rate : float
-        Sample rate of physio, in Hertz.
+    data : physutils.Physio, np.ndarray, or array-like object
+        Object containing the timeseries of the recorded cardiac signal
     slice_timings : 1D array_like
         Slice times, in seconds.
     n_scans : int
         Number of volumes in the imaging run.
     t_r : float
         Sampling rate of the imaging run, in seconds.
+    fs : float, optional if data is a physutils.Physio
+        Sampling rate of `data` in Hz
+    peaks : :obj:`numpy.ndarray`, optional if data is a physutils.Physio
+        Indices of peaks in `data`
 
     Returns
     -------
     phase_card : array_like
         Cardiac phase signal, of shape (n_scans,)
     """
+    if isinstance(data, physio.Physio):
+        # Initialize physio object
+        data = physio.check_physio(data, ensure_fs=True, copy=True)
+    elif fs is not None and peaks is not None:
+        data = io.load_physio(data, fs=fs)
+        data._metadata["peaks"] = peaks
+    else:
+        raise ValueError(
+            """
+            To use this function you should either provide a Physio object
+            with existing peaks metadata (e.g. using the peakdet module), or
+            by providing the physiological data timeseries, the sampling frequency,
+            and the peak indices separately.
+            """
+        )
+    if data.peaks.size == 0:
+        raise ValueError(
+            """
+            Peaks must be a non-empty list.
+            Make sure to run peak detection on your physiological data first,
+            using the peakdet module, or other software of your choice.
+            """
+        )
+
     assert slice_timings.ndim == 1, "Slice times must be a 1D array"
     n_slices = np.size(slice_timings)
     phase_card = np.zeros((n_scans, n_slices))
 
-    card_peaks_sec = peaks / sample_rate
+    card_peaks_sec = data.peaks / data.fs
     for i_slice in range(n_slices):
         # generate slice acquisition timings across all scans
         times_crSlice = t_r * np.arange(n_scans) + slice_timings[i_slice]
@@ -332,4 +431,4 @@ def cardiac_phase(peaks, sample_rate, slice_timings, n_scans, t_r):
             ) / (t2 - t1)
         phase_card[:, i_slice] = phase_card_crSlice
 
-    return phase_card
+    return data, phase_card
diff --git a/phys2denoise/metrics/chest_belt.py b/phys2denoise/metrics/chest_belt.py
index eeeb08a..4e080fc 100644
--- a/phys2denoise/metrics/chest_belt.py
+++ b/phys2denoise/metrics/chest_belt.py
@@ -1,6 +1,7 @@
 """Denoising metrics for chest belt recordings."""
 import numpy as np
 import pandas as pd
+from physutils import io, physio
 from scipy.interpolate import interp1d
 from scipy.stats import zscore
 
@@ -8,11 +9,15 @@
 from ..due import due
 from .responses import rrf
 from .utils import apply_function_in_sliding_window as afsw
-from .utils import convolve_and_rescale, rms_envelope_1d
+from .utils import convolve_and_rescale, return_physio_or_metric, rms_envelope_1d
 
 
 @due.dcite(references.BIRN_2006)
-def respiratory_variance_time(resp, peaks, troughs, samplerate, lags=(0, 4, 8, 12)):
+@return_physio_or_metric()
+@physio.make_operation()
+def respiratory_variance_time(
+    data, peaks=None, troughs=None, fs=None, lags=(0, 4, 8, 12), **kwargs
+):
     """
     Implement the Respiratory Variance over Time (Birn et al. 2006).
 
@@ -20,14 +25,14 @@ def respiratory_variance_time(resp, peaks, troughs, samplerate, lags=(0, 4, 8, 1
 
     Parameters
     ----------
-    resp: 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
+    data : physutils.Physio, np.ndarray, or array-like object
+        Object containing the timeseries of the recorded respiratory signal
+    fs : float, optional if data is a physutils.Physio
+        Sampling rate of `data` in Hz
+    peaks : :obj:`numpy.ndarray`, optional if data is a physutils.Physio
+        Indices of peaks in `data`
+    troughs : :obj:`numpy.ndarray`, optional if data is a physutils.Physio
+        Indices of troughs in `data`
     lags: tuple
         lags in seconds of the RVT output. Default is 0, 4, 8, 12.
 
@@ -42,13 +47,38 @@ def respiratory_variance_time(resp, peaks, troughs, samplerate, lags=(0, 4, 8, 1
        respiratory-variation-related fluctuations from neuronal-activity-related
        fluctuations in fMRI”, NeuroImage, vol.31, pp. 1536-1548, 2006.
     """
-    timestep = 1 / samplerate
+    if isinstance(data, physio.Physio):
+        # Initialize physio object
+        data = physio.check_physio(data, ensure_fs=True, copy=True)
+    elif fs is not None and peaks is not None and troughs is not None:
+        data = io.load_physio(data, fs=fs)
+        data._metadata["peaks"] = peaks
+        data._metadata["troughs"] = troughs
+    else:
+        raise ValueError(
+            """
+            To use this function you should either provide a Physio object
+            with existing peaks and troughs metadata (e.g. using the peakdet module), or
+            by providing the physiological data timeseries, the sampling frequency,
+            and the peak and trough indices separately.
+            """
+        )
+    if data.peaks.size == 0 or data.troughs.size == 0:
+        raise ValueError(
+            """
+            Peaks and troughs must be non-empty lists.
+            Make sure to run peak/trough detection on your physiological data first,
+            using the peakdet module, or other software of your choice.
+            """
+        )
+
+    timestep = 1 / data.fs
     # respiration belt timing
-    time = np.arange(0, len(resp) * timestep, timestep)
-    peak_vals = resp[peaks]
-    trough_vals = resp[troughs]
-    peak_time = time[peaks]
-    trough_time = time[troughs]
+    time = np.arange(0, len(data) * timestep, timestep)
+    peak_vals = data[data.peaks]
+    trough_vals = data[data.troughs]
+    peak_time = time[data.peaks]
+    trough_time = time[data.troughs]
     mid_peak_time = (peak_time[:-1] + peak_time[1:]) / 2
     period = np.diff(peak_time)
     # interpolate peak values over all timepoints
@@ -77,17 +107,19 @@ def respiratory_variance_time(resp, peaks, troughs, samplerate, lags=(0, 4, 8, 1
         )
         rvt_lags[:, ind] = temp_rvt
 
-    return rvt_lags
+    return data, rvt_lags
 
 
 @due.dcite(references.POWER_2018)
-def respiratory_pattern_variability(resp, window):
+@return_physio_or_metric()
+@physio.make_operation()
+def respiratory_pattern_variability(data, window, **kwargs):
     """Calculate respiratory pattern variability.
 
     Parameters
     ----------
-    resp : str or 1D numpy.ndarray
-        Tiemseries representing respiration activity.
+    data : physutils.Physio, np.ndarray, or array-like object
+        Object containing the timeseries of the recorded respiratory signal
     window : int
         Window length in samples.
 
@@ -111,27 +143,33 @@ def respiratory_pattern_variability(resp, window):
        data," Proceedings of the National Academy of Sciences, issue 9, vol.
        115, pp. 2105-2114, 2018.
     """
+    # Initialize physio object
+    data = physio.check_physio(data, ensure_fs=False, copy=True)
+
     # First, z-score respiratory traces
-    resp_z = zscore(resp)
+    resp_z = zscore(data.data)
 
     # Collect upper envelope
     rpv_upper_env = rms_envelope_1d(resp_z, window)
 
     # Calculate standard deviation
     rpv_val = np.std(rpv_upper_env)
-    return rpv_val
+
+    return data, rpv_val
 
 
 @due.dcite(references.POWER_2020)
-def env(resp, samplerate, window=10):
+@return_physio_or_metric()
+@physio.make_operation()
+def env(data, fs=None, window=10, **kwargs):
     """Calculate respiratory pattern variability across a sliding window.
 
     Parameters
     ----------
-    resp : (X,) :obj:`numpy.ndarray`
-        A 1D array with the respiratory belt time series.
-    samplerate : :obj:`float`
-        Sampling rate for resp, in Hertz.
+    data : physutils.Physio, np.ndarray, or array-like object
+        Object containing the timeseries of the recorded respiratory signal
+    fs : float, optional if data is a physutils.Physio
+        Sampling rate of `data` in Hz
     window : :obj:`int`, optional
         Size of the sliding window, in seconds.
         Default is 10.
@@ -155,30 +193,67 @@ def env(resp, samplerate, window=10):
        young adults scanned at rest, including systematic changes and 'missed'
        deep breaths," Neuroimage, vol. 204, 2020.
     """
+
+    def _respiratory_pattern_variability(data, window):
+        """
+        Respiratory pattern variability function without utilizing
+        the physutils.Physio object history, only to be used within
+        the context of this function. This is done to only store the
+        chest_belt.env operation call to the history and not all the
+        subsequent sub-operations
+        """
+        # First, z-score respiratory traces
+        resp_z = zscore(data)
+
+        # Collect upper envelope
+        rpv_upper_env = rms_envelope_1d(resp_z, window)
+
+        # Calculate standard deviation
+        rpv_val = np.std(rpv_upper_env)
+        return rpv_val
+
+    if isinstance(data, physio.Physio):
+        # Initialize physio object
+        data = physio.check_physio(data, ensure_fs=True, copy=True)
+    elif fs is not None:
+        data = io.load_physio(data, fs=fs)
+    else:
+        raise ValueError(
+            """
+            To use this function you should either provide a Physio object
+            with the sampling frequency encapsulated, or
+            by providing the physiological data timeseries and the sampling
+            frequency separately.
+            """
+        )
+
     # Convert window to Hertz
-    window = int(window * samplerate)
+    window = int(window * data.fs)
 
     # Calculate RPV across a rolling window
 
     env_arr = (
-        pd.Series(resp)
+        pd.Series(data.data)
         .rolling(window=window, center=True)
-        .apply(respiratory_pattern_variability, args=(window,))
+        .apply(_respiratory_pattern_variability, args=(window,))
     )
     env_arr[np.isnan(env_arr)] = 0.0
-    return env_arr
+
+    return data, env_arr
 
 
 @due.dcite(references.CHANG_GLOVER_2009)
-def respiratory_variance(resp, samplerate, window=6):
+@return_physio_or_metric()
+@physio.make_operation()
+def respiratory_variance(data, fs=None, window=6, **kwargs):
     """Calculate respiratory variance.
 
     Parameters
     ----------
-    resp : (X,) :obj:`numpy.ndarray`
-        A 1D array with the respiratory belt time series.
-    samplerate : :obj:`float`
-        Sampling rate for resp, in Hertz.
+    data : physutils.Physio, np.ndarray, or array-like object
+        Object containing the timeseries of the recorded respiratory signal
+    fs : float, optional if data is a physutils.Physio
+        Sampling rate of `data` in Hz
     window : :obj:`int`, optional
         Size of the sliding window, in seconds.
         Default is 6.
@@ -206,49 +281,81 @@ def respiratory_variance(resp, samplerate, window=6):
        end-tidal CO2, and BOLD signals in resting-state fMRI," Neuroimage,
        issue 4, vol. 47, pp. 1381-1393, 2009.
     """
+    if isinstance(data, physio.Physio):
+        # Initialize physio object
+        data = physio.check_physio(data, ensure_fs=True, copy=True)
+    elif fs is not None:
+        data = io.load_physio(data, fs=fs)
+    else:
+        raise ValueError(
+            """
+            To use this function you should either provide a Physio object
+            with the sampling frequency encapsulated, or
+            by providing the physiological data timeseries and the sampling
+            frequency separately.
+            """
+        )
+
     # Convert window to Hertz
-    halfwindow_samples = int(round(window * samplerate / 2))
+    halfwindow_samples = int(round(window * data.fs / 2))
 
     # Raw respiratory variance
-    rv_arr = afsw(resp, np.std, halfwindow_samples)
+    rv_arr = afsw(data, np.std, halfwindow_samples)
 
     # Convolve with rrf
-    rv_out = convolve_and_rescale(rv_arr, rrf(samplerate), rescale="zscore")
+    rv_out = convolve_and_rescale(rv_arr, rrf(data.fs), rescale="zscore")
 
-    return rv_out
+    return data, rv_out
 
 
-def respiratory_phase(resp, sample_rate, n_scans, slice_timings, t_r):
+@return_physio_or_metric()
+@physio.make_operation()
+def respiratory_phase(data, n_scans, slice_timings, t_r, fs=None, **kwargs):
     """Calculate respiratory phase from respiratory signal.
 
     Parameters
     ----------
-    resp : 1D array_like
-        Respiratory signal.
-    sample_rate : float
-        Sample rate of physio, in Hertz.
+    data : physutils.Physio, np.ndarray, or array-like object
+        Object containing the timeseries of the recorded respiratory signal
     n_scans
         Number of volumes in the imaging run.
     slice_timings
         Slice times, in seconds.
     t_r
         Sample rate of the imaging run, in seconds.
+    fs : float, optional if data is a physutils.Physio
+        Sampling rate of `data` in Hz
 
     Returns
     -------
     phase_resp : array_like
         Respiratory phase signal, of shape (n_scans, n_slices).
     """
+    if isinstance(data, physio.Physio):
+        # Initialize physio object
+        data = physio.check_physio(data, ensure_fs=True, copy=True)
+    elif fs is not None:
+        data = io.load_physio(data, fs=fs)
+    else:
+        raise ValueError(
+            """
+            To use this function you should either provide a Physio object
+            with the sampling frequency encapsulated, or
+            by providing the physiological data timeseries and the sampling
+            frequency separately.
+            """
+        )
+
     assert slice_timings.ndim == 1, "Slice times must be a 1D array"
     n_slices = np.size(slice_timings)
     phase_resp = np.zeros((n_scans, n_slices))
 
     # generate histogram from respiratory signal
     # TODO: Replace with numpy.histogram
-    resp_hist, resp_hist_bins = np.histogram(resp, bins=100)
+    resp_hist, resp_hist_bins = np.histogram(data, bins=100)
 
     # first compute derivative of respiration signal
-    resp_diff = np.diff(resp, n=1)
+    resp_diff = np.diff(data, n=1)
 
     for i_slice in range(n_slices):
         # generate slice acquisition timings across all scans
@@ -256,15 +363,15 @@ def respiratory_phase(resp, sample_rate, n_scans, slice_timings, t_r):
         phase_resp_crSlice = np.zeros(n_scans)
         for j_scan in range(n_scans):
             iphys = int(
-                max([1, round(times_crSlice[j_scan] * sample_rate)])
+                max([1, round(times_crSlice[j_scan] * data.fs)])
             )  # closest idx in resp waveform
             iphys = min([iphys, len(resp_diff)])  # cannot be longer than resp_diff
-            thisBin = np.argmin(abs(resp[iphys] - resp_hist_bins))
+            thisBin = np.argmin(abs(data[iphys] - resp_hist_bins))
             numerator = np.sum(resp_hist[0:thisBin])
             phase_resp_crSlice[j_scan] = (
-                np.math.pi * np.sign(resp_diff[iphys]) * (numerator / len(resp))
+                np.math.pi * np.sign(resp_diff[iphys]) * (numerator / len(data))
             )
 
         phase_resp[:, i_slice] = phase_resp_crSlice
 
-    return phase_resp
+    return data, phase_resp
diff --git a/phys2denoise/metrics/multimodal.py b/phys2denoise/metrics/multimodal.py
index 7a298b2..fa28bae 100644
--- a/phys2denoise/metrics/multimodal.py
+++ b/phys2denoise/metrics/multimodal.py
@@ -1,34 +1,35 @@
 """These functions compute RETROICOR regressors (Glover et al. 2000)."""
 
 import numpy as np
+from physutils import io, physio
 
 from .. import references
 from ..due import due
 from .cardiac import cardiac_phase
 from .chest_belt import respiratory_phase
+from .utils import return_physio_or_metric
 
 
 @due.dcite(references.GLOVER_2000)
+@return_physio_or_metric()
+@physio.make_operation()
 def retroicor(
-    physio,
-    sample_rate,
+    data,
     t_r,
     n_scans,
     slice_timings,
     n_harmonics,
-    card=False,
-    resp=False,
+    physio_type=None,
+    fs=None,
+    cardiac_peaks=None,
+    **kwargs,
 ):
     """Compute RETROICOR regressors.
 
     Parameters
     ----------
-    physio : array_like
-        1D array, whether cardiac or respiratory signal.
-        If cardiac, the array is a set of peaks in seconds.
-        If respiratory, the array is the actual respiratory signal.
-    sample_rate : float
-        Physio sample rate, in Hertz.
+    data : physutils.Physio, np.ndarray, or array-like object
+        Object containing the timeseries of the recorded respiratory or cardiac signal
     t_r : float
         Imaging sample rate, in seconds.
     n_scans : int
@@ -59,6 +60,44 @@ def retroicor(
        correction of physiological motion effects in fMRI: RETROICOR“, Magn. Reson. Med.,
        issue 1, vol. 44, pp. 162-167, 2000.
     """
+    if isinstance(data, physio.Physio):
+        # Initialize physio object
+        data = physio.check_physio(data, ensure_fs=True, copy=True)
+        if data.physio_type is None and physio_type is not None:
+            data._physio_type = physio_type
+        elif data.physio_type is None and physio_type is None:
+            raise ValueError(
+                """
+                Since the provided Physio object does not specify a `physio_type`,
+                this function's `physio_type` parameter must be specified as a
+                value from {'cardiac', 'respiratory'}
+                """
+            )
+
+    elif fs is not None and physio_type is not None:
+        data = io.load_physio(data, fs=fs)
+        data._physio_type = physio_type
+        if data.physio_type == "cardiac":
+            data._metadata["peaks"] = cardiac_peaks
+    else:
+        raise ValueError(
+            """
+            To use this function you should either provide a Physio object
+            with existing peaks metadata if it describes a cardiac signal
+            (e.g. using the peakdet module), or
+            by providing the physiological data timeseries, the sampling frequency,
+            the physio_type and the peak indices separately.
+            """
+        )
+    if not data.peaks and data.physio_type == "cardiac":
+        raise ValueError(
+            """
+            Peaks must be a non-empty list for cardiac data.
+            Make sure to run peak detection on your cardiac data first,
+            using the peakdet module, or other software of your choice.
+            """
+        )
+
     n_slices = np.shape(slice_timings)  # number of slices
 
     # initialize output variables
@@ -73,18 +112,17 @@ def retroicor(
 
         # Compute physiological phases using the timings of physio events (e.g. peaks)
         # slice sampling times
-        if card:
+        if data.physio_type == "cardiac":
             phase[:, i_slice] = cardiac_phase(
-                physio,
+                data,
                 crslice_timings,
                 n_scans,
                 t_r,
             )
 
-        if resp:
+        if data.physio_type == "respiratory":
             phase[:, i_slice] = respiratory_phase(
-                physio,
-                sample_rate,
+                data,
                 n_scans,
                 slice_timings,
                 t_r,
@@ -99,4 +137,9 @@ def retroicor(
                 (j_harm + 1) * phase[i_slice]
             )
 
-    return retroicor_regressors, phase
+    data._computed_metrics["retroicor_regressors"] = dict(
+        metrics=retroicor_regressors, phase=phase
+    )
+    retroicor_regressors = dict(metrics=retroicor_regressors, phase=phase)
+
+    return data, retroicor_regressors
diff --git a/phys2denoise/metrics/responses.py b/phys2denoise/metrics/responses.py
index aabcec9..1a2aeee 100644
--- a/phys2denoise/metrics/responses.py
+++ b/phys2denoise/metrics/responses.py
@@ -2,6 +2,7 @@
 import logging
 
 import numpy as np
+from loguru import logger
 
 from .. import references
 from ..due import due
@@ -56,10 +57,12 @@ def _crf(t):
         ) * np.exp(-0.5 * (((t - 12) ** 2) / 9))
         return rf
 
-    time_stamps = np.arange(0, time_length, 1 / samplerate)
+    time_stamps = np.arange(0, time_length, samplerate)
     time_stamps -= onset
+    logger.debug(f"Time stamps: {time_stamps}")
     crf_arr = _crf(time_stamps)
     crf_arr = crf_arr / max(abs(crf_arr))
+    logger.debug(f"CRF: {crf_arr}")
 
     if inverse:
         return -crf_arr
@@ -135,7 +138,7 @@ def _rrf(t):
         rf = 0.6 * t**2.1 * np.exp(-t / 1.6) - 0.0023 * t**3.54 * np.exp(-t / 4.25)
         return rf
 
-    time_stamps = np.arange(0, time_length, 1 / samplerate)
+    time_stamps = np.arange(0, time_length, samplerate)
     time_stamps -= onset
     rrf_arr = _rrf(time_stamps)
     rrf_arr = rrf_arr / max(abs(rrf_arr))
diff --git a/phys2denoise/metrics/utils.py b/phys2denoise/metrics/utils.py
index 34fb21b..e21a085 100644
--- a/phys2denoise/metrics/utils.py
+++ b/phys2denoise/metrics/utils.py
@@ -1,8 +1,12 @@
 """Miscellaneous utility functions for metric calculation."""
+import functools
+import inspect
 import logging
 
 import numpy as np
+from loguru import logger
 from numpy.lib.stride_tricks import sliding_window_view as swv
+from physutils.physio import Physio
 from scipy.interpolate import interp1d
 from scipy.stats import zscore
 
@@ -33,7 +37,7 @@ def print_metric_call(metric, args):
 
     msg = f"{msg}\n"
 
-    LGR.info(msg)
+    logger.info(msg)
 
 
 def mirrorpad_1d(arr, buffer=250):
@@ -58,7 +62,11 @@ def mirrorpad_1d(arr, buffer=250):
         post_mirror = np.take(mirror, idx, axis=0)
     except IndexError:
         len(arr)
-        LGR.warning(
+        # LGR.warning(
+        #     f"Requested buffer size ({buffer}) is longer than input array length "
+        #     f"({len(arr)}). Fixing buffer size to array length."
+        # )
+        logger.warning(
             f"Requested buffer size ({buffer}) is longer than input array length "
             f"({len(arr)}). Fixing buffer size to array length."
         )
@@ -332,3 +340,117 @@ def export_metric(
                 )
 
     return fileprefix
+
+
+def return_physio_or_metric(*, return_physio=True):
+    """
+    Decorator to check if the input is a Physio object.
+
+    Parameters
+    ----------
+    func : function
+        Function to be decorated
+
+    Returns
+    -------
+    function
+        Decorated function
+    """
+
+    def determine_return_type(func):
+        convolved_metrics = [
+            "respiratory_variance",
+            "heart_rate",
+            "heart_rate_variability",
+            "heart_beat_interval",
+        ]
+
+        @functools.wraps(func)
+        def wrapper(*args, **kwargs):
+            physio, metric = func(*args, **kwargs)
+            default_args = get_default_args(func)
+            if isinstance(args[0], Physio):
+                if "lags" in kwargs:
+                    has_lags = True if len(kwargs["lags"]) > 1 else False
+                elif "lags" in default_args:
+                    has_lags = True if len(default_args["lags"]) > 1 else False
+                else:
+                    has_lags = False
+
+                is_convolved = True if func.__name__ in convolved_metrics else False
+
+                physio._computed_metrics[func.__name__] = Metric(
+                    func.__name__,
+                    metric,
+                    kwargs,
+                    has_lags=has_lags,
+                    is_convolved=is_convolved,
+                )
+
+                return_physio_value = kwargs.get("return_physio", return_physio)
+                if return_physio_value:
+                    return physio
+                else:
+                    return metric
+            else:
+                return metric
+
+        return wrapper
+
+    return determine_return_type
+
+
+def get_default_args(func):
+    # Get the signature of the function
+    sig = inspect.signature(func)
+
+    # Extract default values for each parameter
+    defaults = {
+        k: v.default
+        for k, v in sig.parameters.items()
+        if v.default is not inspect.Parameter.empty
+    }
+
+    return defaults
+
+
+class Metric:
+    def __init__(self, name, data, args, has_lags=False, is_convolved=False):
+        self.name = name
+        self._data = data
+        self._args = args
+        self._has_lags = has_lags
+        self._is_convolved = is_convolved
+
+    def __array__(self):
+        return self.data
+
+    def __getitem__(self, slicer):
+        return self.data[slicer]
+
+    def __len__(self):
+        return len(self.data)
+
+    @property
+    def ndim(self):
+        return self.data.ndim
+
+    @property
+    def shape(self):
+        return self.data.shape
+
+    @property
+    def data(self):
+        return self._data
+
+    @property
+    def args(self):
+        return self._args
+
+    @property
+    def has_lags(self):
+        return self._has_lags
+
+    @property
+    def is_convolved(self):
+        return self._is_convolved
diff --git a/phys2denoise/tests/conftest.py b/phys2denoise/tests/conftest.py
index fcf05f3..757a071 100644
--- a/phys2denoise/tests/conftest.py
+++ b/phys2denoise/tests/conftest.py
@@ -1,5 +1,20 @@
 import numpy as np
 import pytest
+from _pytest.logging import LogCaptureFixture
+from loguru import logger
+
+
+@pytest.fixture
+def caplog(caplog: LogCaptureFixture):
+    handler_id = logger.add(
+        caplog.handler,
+        format="{message}",
+        level=20,
+        filter=lambda record: record["level"].no >= caplog.handler.level,
+        enqueue=False,
+    )
+    yield caplog
+    logger.remove(handler_id)
 
 
 @pytest.fixture(scope="module")
diff --git a/phys2denoise/tests/test_metrics_cardiac.py b/phys2denoise/tests/test_metrics_cardiac.py
index bb10d76..2d776cf 100644
--- a/phys2denoise/tests/test_metrics_cardiac.py
+++ b/phys2denoise/tests/test_metrics_cardiac.py
@@ -1,5 +1,7 @@
 """Tests for phys2denoise.metrics.cardiac."""
 import numpy as np
+from loguru import logger
+from physutils import physio
 
 from phys2denoise.metrics import cardiac
 
@@ -7,20 +9,15 @@
 def test_crf_smoke():
     """Basic smoke test for CRF calculation."""
     samplerate = 0.01  # in seconds
-    oversampling = 20
     time_length = 20
     onset = 0
     tr = 0.72
     crf_arr = cardiac.crf(
-        samplerate,
-        oversampling=oversampling,
-        time_length=time_length,
-        onset=onset,
-        tr=tr,
+        samplerate, time_length=time_length, onset=onset, inverse=False
     )
-    pred_len = np.rint(time_length / (tr / oversampling))
+    pred_len = np.rint(time_length * 1 / samplerate)
     assert crf_arr.ndim == 1
-    assert crf_arr.shape == pred_len
+    assert len(crf_arr) == pred_len
 
 
 def test_cardiac_phase_smoke():
@@ -30,12 +27,55 @@ def test_cardiac_phase_smoke():
     sample_rate = 1 / 0.01
     slice_timings = np.linspace(0, t_r, 22)[1:-1]
     peaks = np.array([0.534, 0.577, 10.45, 20.66, 50.55, 90.22])
+    data = np.zeros(peaks.shape)
     card_phase = cardiac.cardiac_phase(
-        peaks,
-        sample_rate=sample_rate,
+        data,
+        peaks=peaks,
+        fs=sample_rate,
         slice_timings=slice_timings,
         n_scans=n_scans,
         t_r=t_r,
     )
     assert card_phase.ndim == 2
     assert card_phase.shape == (n_scans, slice_timings.size)
+
+
+def test_cardiac_phase_smoke_physio_obj():
+    """Basic smoke test for cardiac phase calculation."""
+    t_r = 1.0
+    n_scans = 200
+    sample_rate = 1 / 0.01
+    slice_timings = np.linspace(0, t_r, 22)[1:-1]
+    peaks = np.array([0.534, 0.577, 10.45, 20.66, 50.55, 90.22])
+    data = np.zeros(peaks.shape)
+    phys = physio.Physio(data, sample_rate, physio_type="cardiac")
+    phys._metadata["peaks"] = peaks
+
+    # Test where the physio object is returned
+    phys = cardiac.cardiac_phase(
+        phys,
+        slice_timings=slice_timings,
+        n_scans=n_scans,
+        t_r=t_r,
+    )
+
+    assert phys.history[0][0] == "phys2denoise.metrics.cardiac.cardiac_phase"
+    assert phys.computed_metrics["cardiac_phase"].ndim == 2
+    assert phys.computed_metrics["cardiac_phase"].shape == (
+        n_scans,
+        slice_timings.size,
+    )
+    assert phys.computed_metrics["cardiac_phase"].args["slice_timings"] is not None
+    assert phys.computed_metrics["cardiac_phase"].args["n_scans"] is not None
+    assert phys.computed_metrics["cardiac_phase"].args["t_r"] is not None
+
+    # Test where the metric is returned
+    card_phase = cardiac.cardiac_phase(
+        phys,
+        slice_timings=slice_timings,
+        n_scans=n_scans,
+        t_r=t_r,
+        return_physio=False,
+    )
+    assert card_phase.ndim == 2
+    assert card_phase.shape == (n_scans, slice_timings.size)
diff --git a/phys2denoise/tests/test_metrics_chest_belt.py b/phys2denoise/tests/test_metrics_chest_belt.py
index 1250d85..a4257ec 100644
--- a/phys2denoise/tests/test_metrics_chest_belt.py
+++ b/phys2denoise/tests/test_metrics_chest_belt.py
@@ -14,17 +14,14 @@ def test_rrf_smoke():
     tr = 0.72
     rrf_arr = chest_belt.rrf(
         samplerate,
-        oversampling=oversampling,
         time_length=time_length,
         onset=onset,
-        tr=tr,
     )
-    pred_len = int(np.rint(time_length / (tr / oversampling)))
+    pred_len = int(np.rint(time_length * (1 / samplerate)))
     assert rrf_arr.ndim == 1
     assert rrf_arr.size == pred_len
 
 
-@mark.xfail
 def test_respiratory_phase_smoke():
     """Basic smoke test for respiratory phase calculation."""
     t_r = 1.0
@@ -35,7 +32,7 @@ def test_respiratory_phase_smoke():
     resp = np.random.normal(size=n_samples)
     resp_phase = chest_belt.respiratory_phase(
         resp,
-        sample_rate=sample_rate,
+        fs=sample_rate,
         slice_timings=slice_timings,
         n_scans=n_scans,
         t_r=t_r,
@@ -59,7 +56,7 @@ def test_env_smoke():
     resp = np.random.normal(size=n_samples)
     samplerate = 1 / 0.01
     window = 6
-    env_arr = chest_belt.env(resp, samplerate=samplerate, window=window)
+    env_arr = chest_belt.env(resp, fs=samplerate, window=window)
     assert env_arr.ndim == 1
     assert env_arr.shape == (n_samples,)
 
@@ -70,6 +67,6 @@ def test_respiratory_variance_smoke():
     resp = np.random.normal(size=n_samples)
     samplerate = 1 / 0.01
     window = 6
-    rv_arr = chest_belt.respiratory_variance(resp, samplerate=samplerate, window=window)
+    rv_arr = chest_belt.respiratory_variance(resp, fs=samplerate, window=window)
     assert rv_arr.ndim == 2
     assert rv_arr.shape == (n_samples, 2)
diff --git a/phys2denoise/tests/test_metrics_utils.py b/phys2denoise/tests/test_metrics_utils.py
index 916c485..1995aaf 100644
--- a/phys2denoise/tests/test_metrics_utils.py
+++ b/phys2denoise/tests/test_metrics_utils.py
@@ -16,12 +16,12 @@ def test_mirrorpad(short_arr):
     assert all(arr_mirror == expected_arr_mirror)
 
 
-def test_mirrorpad_exception(short_arr):
+def test_mirrorpad_exception(short_arr, caplog):
     """When passing array that is too short to perform mirrorpadding, the
     function should give an error."""
     arr = np.array(short_arr)
-    with pytest.raises(Exception) as e_info:
-        arr_mirror = mirrorpad_1d(short_arr)
+    arr_mirror = mirrorpad_1d(arr)
+    assert caplog.text.count("Fixing buffer size to array length.") > 0
 
 
 def test_rms_envelope():
diff --git a/phys2denoise/tests/test_rvt.py b/phys2denoise/tests/test_rvt.py
index 13e27d8..0bc4066 100644
--- a/phys2denoise/tests/test_rvt.py
+++ b/phys2denoise/tests/test_rvt.py
@@ -1,4 +1,5 @@
 import peakdet
+from physutils import io, physio
 
 from phys2denoise.metrics.chest_belt import respiratory_variance_time
 
@@ -15,8 +16,49 @@ def test_respiratory_variance_time(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 = test_respiratory_variance_time(
-        phys.data, phys.peaks, phys.troughs, samplerate=phys.fs
+
+    # TODO: Change to a simpler call once physutils are
+    # integrated to peakdet/prep4phys
+    r = respiratory_variance_time(
+        phys.data, fs=phys.fs, peaks=phys.peaks, troughs=phys.troughs
+    )
+    assert r is not None
+    assert len(r) == 18750
+
+
+def test_respiratory_variance_time(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)
+
+    # TODO: Change to a simpler call once physutils are
+    # integrated to peakdet/prep4phys
+    r = respiratory_variance_time(
+        phys.data, fs=phys.fs, peaks=phys.peaks, troughs=phys.troughs
     )
     assert r is not None
     assert len(r) == 18750
+
+
+def test_respiratory_variance_time_physio_obj(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)
+
+    # TODO: Change to a simpler call once physutils are
+    # integrated to peakdet/prep4phys
+    physio_obj = physio.Physio(phys.data, phys.fs)
+    physio_obj._metadata["peaks"] = phys.peaks
+    physio_obj._metadata["troughs"] = phys.troughs
+    physio_obj = respiratory_variance_time(physio_obj)
+
+    assert (
+        physio_obj.history[0][0]
+        == "phys2denoise.metrics.chest_belt.respiratory_variance_time"
+    )
+    assert physio_obj.computed_metrics["respiratory_variance_time"].ndim == 2
+    assert physio_obj.computed_metrics["respiratory_variance_time"].shape == (18750, 4)
+    assert physio_obj.computed_metrics["respiratory_variance_time"].has_lags == True
+
+    # assert r is not None
+    # assert len(r) == 18750
diff --git a/setup.cfg b/setup.cfg
index a3f8ff9..81268a1 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -19,13 +19,15 @@ provides =
     phys2denoise
 
 [options]
-python_requires = >=3.6.1
+python_requires = >=3.9
 install_requires =
-    numpy >=1.9.3
-    matplotlib >=3.1.1
+    numpy >=1.9.3, <2
+    matplotlib
     pandas
     scipy
     duecredit
+    loguru
+    physutils >=0.2.1
 tests_require =
     pytest >=5.3
 test_suite = pytest
@@ -48,7 +50,7 @@ test =
     %(style)s
     pytest >=5.3
     pytest-cov
-    peakdet
+    peakdet>=0.5.0
     coverage
 devtools =
     pre-commit