diff --git a/phys2denoise/metrics/__init__.py b/phys2denoise/metrics/__init__.py
index 57eeef1..3f91005 100644
--- a/phys2denoise/metrics/__init__.py
+++ b/phys2denoise/metrics/__init__.py
@@ -1,9 +1,8 @@
 """Metrics derived from physiological signals."""
-from .cardiac import cardiac_phase, heart_beat_interval, heart_rate_variability
+from .cardiac import heart_beat_interval, heart_rate_variability
 from .chest_belt import (
     env,
     respiratory_pattern_variability,
-    respiratory_phase,
     respiratory_variance,
     respiratory_variance_time,
 )
@@ -12,8 +11,6 @@
 
 __all__ = [
     "retroicor",
-    "cardiac_phase",
-    "respiratory_phase",
     "heart_beat_interval",
     "heart_rate_variability",
     "env",
diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py
index 5b7588d..5f957df 100644
--- a/phys2denoise/metrics/cardiac.py
+++ b/phys2denoise/metrics/cardiac.py
@@ -346,89 +346,3 @@ def heart_beat_interval(
     )
 
     return data, hbi
-
-
-@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
-    as slice timings, for example seconds.
-
-    Parameters
-    ----------
-    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 = 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]
-        phase_card_crSlice = np.zeros(n_scans)
-        for j_scan in range(n_scans):
-            previous_card_peaks = np.asarray(
-                np.nonzero(card_peaks_sec < times_crSlice[j_scan])
-            )
-            if np.size(previous_card_peaks) == 0:
-                t1 = 0
-            else:
-                last_peak = previous_card_peaks[0][-1]
-                t1 = card_peaks_sec[last_peak]
-            next_card_peaks = np.asarray(
-                np.nonzero(card_peaks_sec > times_crSlice[j_scan])
-            )
-            if np.size(next_card_peaks) == 0:
-                t2 = n_scans * t_r
-            else:
-                next_peak = next_card_peaks[0][0]
-                t2 = card_peaks_sec[next_peak]
-            phase_card_crSlice[j_scan] = (
-                2 * np.math.pi * (times_crSlice[j_scan] - t1)
-            ) / (t2 - t1)
-        phase_card[:, i_slice] = phase_card_crSlice
-
-    return data, phase_card
diff --git a/phys2denoise/metrics/chest_belt.py b/phys2denoise/metrics/chest_belt.py
index 4e080fc..00ba671 100644
--- a/phys2denoise/metrics/chest_belt.py
+++ b/phys2denoise/metrics/chest_belt.py
@@ -306,72 +306,3 @@ def respiratory_variance(data, fs=None, window=6, **kwargs):
     rv_out = convolve_and_rescale(rv_arr, rrf(data.fs), rescale="zscore")
 
     return data, rv_out
-
-
-@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
-    ----------
-    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(data, bins=100)
-
-    # first compute derivative of respiration signal
-    resp_diff = np.diff(data, n=1)
-
-    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]
-        phase_resp_crSlice = np.zeros(n_scans)
-        for j_scan in range(n_scans):
-            iphys = int(
-                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(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(data))
-            )
-
-        phase_resp[:, i_slice] = phase_resp_crSlice
-
-    return data, phase_resp
diff --git a/phys2denoise/metrics/multimodal.py b/phys2denoise/metrics/multimodal.py
index fa28bae..2975a71 100644
--- a/phys2denoise/metrics/multimodal.py
+++ b/phys2denoise/metrics/multimodal.py
@@ -1,145 +1,235 @@
 """These functions compute RETROICOR regressors (Glover et al. 2000)."""
 
+from copy import deepcopy as dc
+
 import numpy as np
+from peakdet import Physio as pkphysio
 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
+from .phase import cardiac, respiratory
 
 
 @due.dcite(references.GLOVER_2000)
-@return_physio_or_metric()
-@physio.make_operation()
 def retroicor(
-    data,
-    t_r,
-    n_scans,
     slice_timings,
-    n_harmonics,
-    physio_type=None,
-    fs=None,
-    cardiac_peaks=None,
-    **kwargs,
+    resp_data=None,
+    card_data=None,
+    resp_fs=None,
+    card_fs=None,
+    n_harmonics=2,
+    base_offset=0,
+    nbins="p2d",
+    compute_interaction=True,
 ):
     """Compute RETROICOR regressors.
 
     Parameters
     ----------
-    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
-        Number of volumes in the imaging run.
-    slice_timings : array_like
+    slice_timings : float or array_like
         Slice times, in seconds.
-    n_harmonics : int
-        ???
-    card : bool, optional
-        Whether the physio data correspond to cardiac or repiratory signal.
-    resp : bool, optional
-        Whether the physio data correspond to cardiac or repiratory signal.
-
-    Returns
-    -------
-    retroicor_regressors : list
-    phase : array_like
-        2D array of shape (n_scans, n_slices)
+    resp_data : None, physutils.Physio, or array-like object, optional
+        Respiratory data, either a physio object or an array-like object describing the
+        recorded respiratory changes *timeseries*
+    card_data : None, physutils.Physio, or array-like object, optional
+        Cardiac data, either a physio object or an array-like object describing the
+        recorded cardiac *peaks' indices* (!!!)
+    resp_fs : None or int, optional
+        Respiratory data sampling rate, optional if resp_data is a physio object,
+        necessary otherwise.
+    card_fs : None or int, optional
+        Cardiac data sampling rate, optional if card_data is a physio object,
+        necessary otherwise.
+    n_harmonics : int > 0, optional
+        Number of harmonics to compute (for all modalities). Must be > 0. Default is 2.
+    base_offset : int or float >=0, optional
+        Offset time for the whole neuroimaging data wrt physiological data, default is 0.
+        Must be >=0. Negative offsets are currently not supported.
+        Note that neuroimaging data offset is normally found in BIDS as a metadata of the
+        physiological recording, so the signal must be inverted compared to that field.
+    nbins : int or string, optional
+        Number of bins to consider in making the histogram or method to estimate such
+        number when computing respiratory phase. The default option is "p2d" and
+        corresponds to either a quarter of the data amount or 100, whatever is higher.
+        Any option supported by `np.histogram` is also supported here.
+    compute_interaction : bool, optional
+        If respiratory and cardiac data are both given, compute their interactive terms,
+        see [2]
+        Default is True.
 
     Notes
     -----
     RETROICOR regressors should be regressed from the imaging data *before*
-    any other preprocessing, including slice-timing correction and motion correction.
+    any other preprocessing, including slice-timing correction and motion correction,
+    and it should be voxel or slice specific.
+    Despite this, many people regress one single set of regressors in the GLM - in that
+    case, don't use any offset.
+
+    This function does not work like the others due to the issue of having multiple
+    possible physio objects in entrance. It will return all physio objects with all
+    computed retroicor regressors, so pay attention to repetitions!
+
+    Does not support voxelwise correction yet, thus coefficients are not given.
 
     References
     ----------
     .. [1] G. H. Glover & T. Q. L. Ress, “Image_based method for retrospective
        correction of physiological motion effects in fMRI: RETROICOR“, Magn. Reson. Med.,
        issue 1, vol. 44, pp. 162-167, 2000.
+    .. [2] A.K. Harvey & all, "Brainstem functional magnetic resonance imaging:
+       Disentangling signal from physiological noise", JMRI, issue 6, vol. 28, 2008.
+
+    Returns
+    -------
+    retroicor_regressors : dict
+        Retroicor regressors are organised in a dictionary of two dictionaries, one for
+        retroicor regressors and one for phases. These dictionaries have one entry per
+        slice.
+        All entries are lists, always organised: cardiac regressors first, then
+        respiratory, then interactions.
+
+        If no physutils.Physio object is given, they will be returned as they are.
+    resp_data : physutils.Physio object
+        If resp_data is given as physutils.Physio object, returns it with the computed
+        retroicor regressors in its metrics metadata.
+    card_data : physutils.Physio object
+        If card_data is given as physutils.Physio object, returns it with the computed
+        retroicor regressors in its metrics metadata.
     """
-    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:
+    if n_harmonics <= 0:
+        raise ValueError(
+            f"The number of harmonics must be > 0 but {n_harmonics} was given."
+        )
+
+    # Parse input data
+    return_physio = False
+
+    if resp_data is not None:
+        if isinstance(resp_data, physio.Physio):
+            # Initialize physio object
+            resp_data = physio.check_physio(resp_data, ensure_fs=True, copy=True)
+            return_physio = True
+        elif isinstance(resp_data, pkphysio):
+            # Retrocompatibility with peakdet
+            resp_data = io.load_physio(resp_data.data, fs=resp_data.fs)
+        elif resp_fs is not None:
+            resp_data = physio.Physio(resp_data, fs=resp_fs)
+        else:
             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'}
+                resp_data is not a Physio object but resp_fs was not specified.
+                To use this function with respiratory data you should either provide a
+                Physio object describing the physiological data timeseries, or an
+                array-like object *and* the sampling frequency.
+                """
+            )
+    elif card_data is None:
+        return None
+
+    if card_data is not None:
+        if isinstance(card_data, physio.Physio):
+            # Initialize physio object
+            card_data = physio.check_physio(card_data, ensure_fs=True, copy=True)
+            return_physio = True
+        elif isinstance(card_data, pkphysio):
+            # Retrocompatibility with peakdet
+            peaks = dc(card_data.peaks)
+            card_data = io.load_physio(card_data.data, fs=card_data.fs)
+            card_data._metadata["peaks"] = peaks
+        elif card_fs is not None:
+            peaks = dc(card_data)
+            card_data = physio.Physio(np.zeros((card_data[-1] + 1)), fs=card_fs)
+            card_data._metadata["peaks"] = peaks
+        else:
+            raise ValueError(
+                """
+                card_data is not a Physio object but card_fs was not specified.
+                To use this function with cardiac data you should either provide a
+                Physio object with cardiac data timeseries and peaks, or an
+                array-like object describing peaks *and* the sampling frequency.
                 """
             )
 
-    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":
+        if not (hasattr(card_data, "peaks") and np.any(card_data.peaks)):
+            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.
+                """
+            )
+
+    slice_timings = np.asarray(slice_timings)
+    if slice_timings.squeeze().ndim > 1:
         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.
-            """
+            "The provided slice timings are in a multidimensional array. Please "
+            "provide a single 1D array-like data."
         )
 
-    n_slices = np.shape(slice_timings)  # number of slices
-
-    # initialize output variables
-    retroicor_regressors = []
-    phase = np.empty((n_scans, n_slices))
-
-    for i_slice in range(n_slices):
-        retroicor_regressors.append(np.empty((n_scans, 2 * n_harmonics)))
-
-        # Initialize slice timings for current slice
-        crslice_timings = t_r * np.arange(n_scans) + slice_timings[i_slice]
+    # Initialize output variable
+    retroicor_regressors = {}
+    phases = {}
 
+    # Compute slice-dependent retroicor regressors.
+    for n, slice_time in enumerate(slice_timings):
+        phases[n] = {}
+        offset = base_offset + slice_time
         # Compute physiological phases using the timings of physio events (e.g. peaks)
         # slice sampling times
-        if data.physio_type == "cardiac":
-            phase[:, i_slice] = cardiac_phase(
-                data,
-                crslice_timings,
-                n_scans,
-                t_r,
-            )
-
-        if data.physio_type == "respiratory":
-            phase[:, i_slice] = respiratory_phase(
-                data,
-                n_scans,
-                slice_timings,
-                t_r,
-            )
-
-        # Compute retroicor regressors
-        for j_harm in range(n_harmonics):
-            retroicor_regressors[i_slice][:, 2 * j_harm] = np.cos(
-                (j_harm + 1) * phase[i_slice]
-            )
-            retroicor_regressors[i_slice][:, 2 * j_harm + 1] = np.sin(
-                (j_harm + 1) * phase[i_slice]
-            )
-
-    data._computed_metrics["retroicor_regressors"] = dict(
-        metrics=retroicor_regressors, phase=phase
-    )
-    retroicor_regressors = dict(metrics=retroicor_regressors, phase=phase)
-
-    return data, retroicor_regressors
+        if card_data is not None:
+            phases[n]["card"] = cardiac(card_data.peaks, card_data.fs, offset)
+
+        if resp_data is not None:
+            phases[n]["resp"] = respiratory(resp_data.data, resp_data.fs, offset, nbins)
+
+            if card_data is not None:
+                # Cut to same length. It'll be better later.
+                length = np.min((phases[n]["card"].size, phases[n]["resp"].size))
+                phases[n]["card"] = phases[n]["card"][:length]
+                phases[n]["resp"] = phases[n]["resp"][:length]
+
+        retroicor_regressors[n] = []
+
+        # It's not computationally efficient to loop multiple times per if statement,
+        # but organising regressors this way is better for users.
+        if card_data is not None:
+            for m in range(1, n_harmonics + 1):
+                retroicor_regressors[n] = retroicor_regressors[n] + [
+                    np.cos(m * phases[n]["card"]),
+                    np.sin(m * phases[n]["card"]),
+                ]
+
+        if resp_data is not None:
+            for m in range(1, n_harmonics + 1):
+                retroicor_regressors[n] = retroicor_regressors[n] + [
+                    np.cos(m * phases[n]["resp"]),
+                    np.sin(m * phases[n]["resp"]),
+                ]
+
+            if compute_interaction and card_data is not None:
+                for m in range(1, n_harmonics + 1):
+                    retroicor_regressors[n] = retroicor_regressors[n] + [
+                        np.cos(m * phases[n]["card"]) * np.cos(m * phases[n]["resp"]),
+                        np.cos(m * phases[n]["card"]) * np.sin(m * phases[n]["resp"]),
+                        np.sin(m * phases[n]["card"]) * np.cos(m * phases[n]["resp"]),
+                        np.sin(m * phases[n]["card"]) * np.sin(m * phases[n]["resp"]),
+                    ]
+
+    retroicor_regressors = {"retroicor": retroicor_regressors, "phases": phases}
+
+    if return_physio:
+        if card_data is not None:
+            card_data._computed_metrics["retroicor_regressors"] = retroicor_regressors
+
+        if resp_data is not None:
+            resp_data._computed_metrics["retroicor_regressors"] = retroicor_regressors
+
+            if card_data is not None:
+                return resp_data, card_data
+            else:
+                return resp_data
+        else:
+            return card_data
+    else:
+        return retroicor_regressors
diff --git a/phys2denoise/metrics/phase.py b/phys2denoise/metrics/phase.py
new file mode 100644
index 0000000..a0df5a3
--- /dev/null
+++ b/phys2denoise/metrics/phase.py
@@ -0,0 +1,179 @@
+"""Denoising metrics for cardio recordings."""
+
+import numpy as np
+from loguru import logger
+from scipy.interpolate import interp1d
+
+
+def cardiac(peaks, fs, offset=0):
+    """Calculate cardiac phase from cardiac peaks.
+
+    Assumes that timing of cardiac events are given in same units
+    as slice timings, for example seconds.
+
+    This function will always return just the metric, never a physio object with the
+    metric inside.
+
+    Parameters
+    ----------
+    peaks : 1D array_like
+        Indices of peaks in a cardiac array
+    fs : float
+        Sampling rate of physiological data in Hz
+    offset : int or float >=0, optional
+        Offset time for the neuroimaging data. This includes slice timing (normally
+        positive) and neuroimaging data offset wrt physiological data (normally positive).
+        Must be >=0. Negative offsets are currently not supported.
+        Note that neuroimaging data offset is normally found in BIDS as a metadata of the
+        physiological recording, so the signal must be inverted compared to that field.
+
+    Returns
+    -------
+    phase_card : array_like
+        Cardiac phase signal, sampled at the physiological signal sampling rate and of
+        length = peaks[-1]
+
+    Note
+    ----
+    This function will always return just the metric, never a physio object with the
+    metric inside, unlike every other metric function in this module.
+
+    This function won't fail, but you will not be able to export regressors, if the
+    physiological data collected ends before the neuroimaging data.
+
+    `time1` and `time2` refer to the original formula from [1]:
+        phi(t) = 2*pi*(t-t1)/(t2-t1)
+    They represent the beat right before the reference timepoint and the one after that respectively.
+
+    Raises
+    ------
+    ValueError
+        If the offset is negative.
+
+    References
+    ----------
+    .. [1] G. H. Glover & T. Q. L. Ress, “Image_based method for retrospective
+       correction of physiological motion effects in fMRI: RETROICOR“, Magn. Reson. Med.,
+       issue 1, vol. 44, pp. 162-167, 2000.
+    """
+    if offset < 0:
+        raise ValueError("Negative offsets are not supported yet.")
+
+    # Add a beat before and after the current beats
+    avg_peak_dist = int(np.diff(peaks).mean().round())
+    peaks = np.append(
+        np.append(peaks[0] - avg_peak_dist, peaks),
+        [peaks[-1] + avg_peak_dist, peaks[-1] + 2 * avg_peak_dist],
+    )
+
+    # Pre-append additional peaks until we get to <0
+    while peaks[0] >= 0:
+        peaks = np.append((peaks[0] - avg_peak_dist), peaks)
+
+    # Transform in seconds and offset
+    peaks_sec = peaks / fs - offset
+
+    # Create timeline reference for neuroimaging on the peaks time, adding offset,
+    # accounting for the two extra peaks added above.
+    time = np.arange(peaks[-3] + 1) / fs
+
+    time1 = interp1d(peaks_sec, peaks_sec, kind="previous", assume_sorted=True)(time)
+    time2 = interp1d(peaks_sec, peaks_sec, kind="next", assume_sorted=True)(time)
+
+    return 2 * np.pi * ((time[1:] - time1[:-1]) / (time2[1:] - time1[:-1]))
+
+
+def respiratory(data, fs, offset=0, nbins="p2d"):
+    """Calculate respiratory phase from respiratory signal.
+
+    Parameters
+    ----------
+    data : array-like object
+        Recorded respiratory signal's timeseries
+    fs : float
+        Sampling rate of physiological data in Hz
+    offset : int or float >=0, optional
+        Offset time for the neuroimaging data, default is 0. This includes slice timing
+        (normally positive) and neuroimaging data offset wrt physiological data
+        (normally positive).
+        Must be >=0. Negative offsets are currently not supported.
+        Note that neuroimaging data offset is normally found in BIDS as a metadata of the
+        physiological recording, so the signal must be inverted compared to that field.
+    nbins : int or string, optional
+        Number of bins to consider in making the histogram or method to estimate such
+        number. The default option is "p2d" and corresponds to either a quarter of the
+        data amount or 100, whatever is higher.
+        Any option supported by `np.histogram` is also supported here.
+
+    Returns
+    -------
+    phase_resp : array_like
+        Respiratory phase signal, sampled at the physiological signal sampling rate and
+        of length = data.size
+
+    Note
+    ----
+    This function will always return just the metric, never a physio object with the
+    metric inside, unlike every other metric function in this module.
+
+    This function won't fail, but you will not be able to export regressors, if the
+    physiological data collected ends before the neuroimaging data.
+
+    This is an adapted and vectorized version of the original formula from [1].
+    The original formula has basically three elements for each timepoint:
+    - pi
+    - the area under the curve (AUC) left of any bin that timepoint is in,
+      divided by the total AUC
+    - the sign of the (discrete) derivative of the signal at that timepoint
+
+    Here it is implemented as:
+    - pi (duh)
+    - the cumulative sum of the counts up to the bin of the timepoint, in a histogram
+      normalised by the amount of data entries (so that total AUC, i.e. cumsum, = 1)
+    - the sign of the one-hop difference of the signal.
+
+    It is, however, exactly the same thing (beside decimal point error).
+
+    Also, while the histogram is created on the original data, the search happens with
+    the offsetted data (if offset is specified). This is so that all slice time shifts
+    refer to the same histogram.
+
+    Raises
+    ------
+    ValueError
+        If offset is < 0
+
+    References
+    ----------
+    .. [1] G. H. Glover & T. Q. L. Ress, “Image_based method for retrospective
+       correction of physiological motion effects in fMRI: RETROICOR“, Magn. Reson. Med.,
+       issue 1, vol. 44, pp. 162-167, 2000.
+    """
+    if offset < 0:
+        raise ValueError("Negative offsets are not supported yet.")
+    elif offset > 0:
+        # Interpolate data in neuroimaging's offsetted time.
+        time = np.arange(data.size)
+        data_offsetted = interp1d(
+            time, data, kind="linear", fill_value="extrapolate", assume_sorted=True
+        )(time + offset)
+    else:
+        # Skip interpolation
+        data_offsetted = data
+
+    nbins = max(100, int(data.size / 4)) if nbins == "p2d" else nbins
+    counts, binedge = np.histogram(data, bins=nbins)
+    # Normalize counts by total to avoid denominator computation later
+    counts = counts / counts.sum()
+    # For each element, find its bin
+    bincenter = (binedge[:-1] + binedge[1:]) / 2
+    idx = np.searchsorted(bincenter, data_offsetted, side="left")
+    # Cumulative sum computes AUC left of any bin
+    cumsum = np.cumsum(counts)
+    # use the bin index to find the AUC left of it.
+    raw_phase = np.where(idx > 0, cumsum[idx - 1], 0)
+
+    sign = np.sign(np.diff(data))
+    sign = np.append(sign, sign[-1])
+
+    return np.pi * raw_phase * sign
diff --git a/phys2denoise/tests/test_metrics_cardiac.py b/phys2denoise/tests/test_metrics_cardiac.py
index 2d776cf..60bd211 100644
--- a/phys2denoise/tests/test_metrics_cardiac.py
+++ b/phys2denoise/tests/test_metrics_cardiac.py
@@ -20,62 +20,62 @@ def test_crf_smoke():
     assert len(crf_arr) == pred_len
 
 
-def test_cardiac_phase_smoke():
-    """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)
-    card_phase = cardiac.cardiac_phase(
-        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():
+#     """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)
+#     card_phase = cardiac.cardiac_phase(
+#         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
+# 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,
-    )
+#     # 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
+#     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)
+#     # 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 a4257ec..f3c8e77 100644
--- a/phys2denoise/tests/test_metrics_chest_belt.py
+++ b/phys2denoise/tests/test_metrics_chest_belt.py
@@ -22,23 +22,23 @@ def test_rrf_smoke():
     assert rrf_arr.size == pred_len
 
 
-def test_respiratory_phase_smoke():
-    """Basic smoke test for respiratory phase calculation."""
-    t_r = 1.0
-    n_scans = 200
-    sample_rate = 1 / 0.01
-    slice_timings = np.linspace(0, t_r, 22)[1:-1]
-    n_samples = int(np.rint((n_scans * t_r) * sample_rate))
-    resp = np.random.normal(size=n_samples)
-    resp_phase = chest_belt.respiratory_phase(
-        resp,
-        fs=sample_rate,
-        slice_timings=slice_timings,
-        n_scans=n_scans,
-        t_r=t_r,
-    )
-    assert resp_phase.ndim == 2
-    assert resp_phase.shape == (n_scans, slice_timings.size)
+# def test_respiratory_phase_smoke():
+#     """Basic smoke test for respiratory phase calculation."""
+#     t_r = 1.0
+#     n_scans = 200
+#     sample_rate = 1 / 0.01
+#     slice_timings = np.linspace(0, t_r, 22)[1:-1]
+#     n_samples = int(np.rint((n_scans * t_r) * sample_rate))
+#     resp = np.random.normal(size=n_samples)
+#     resp_phase = chest_belt.respiratory_phase(
+#         resp,
+#         fs=sample_rate,
+#         slice_timings=slice_timings,
+#         n_scans=n_scans,
+#         t_r=t_r,
+#     )
+#     assert resp_phase.ndim == 2
+#     assert resp_phase.shape == (n_scans, slice_timings.size)
 
 
 def test_respiratory_pattern_variability_smoke():