From 41df1f4abbb7a294fa3e1571897043c0d03fa064 Mon Sep 17 00:00:00 2001
From: CesarCaballeroGaudes <c.caballero@bcbl.eu>
Date: Mon, 11 May 2020 13:06:38 +0200
Subject: [PATCH 01/18] implementing retroicor

---
 phys2denoise/metrics/retroicor.py | 12 ++++++++++++
 1 file changed, 12 insertions(+)
 create mode 100644 phys2denoise/metrics/retroicor.py

diff --git a/phys2denoise/metrics/retroicor.py b/phys2denoise/metrics/retroicor.py
new file mode 100644
index 0000000..df9ee3c
--- /dev/null
+++ b/phys2denoise/metrics/retroicor.py
@@ -0,0 +1,12 @@
+"""implementation of retroicor"""
+
+import argparse as ap
+import numpy as np
+import scipy as sp
+import json
+import string
+import random
+import matplotlib as mpl; mpl.use('TkAgg')
+import matplotlib.pyplot as plt
+from scipy.signal import argrelmax
+

From 4ad18ba408801976ec0373ff813252499106b72d Mon Sep 17 00:00:00 2001
From: CesarCaballeroGaudes <c.caballero@bcbl.eu>
Date: Wed, 17 Jun 2020 13:31:14 +0200
Subject: [PATCH 02/18] draft implementation of retroicor

---
 phys2denoise/metrics/retroicor.py | 75 +++++++++++++++++++++++++++++++
 1 file changed, 75 insertions(+)

diff --git a/phys2denoise/metrics/retroicor.py b/phys2denoise/metrics/retroicor.py
index df9ee3c..e63e929 100644
--- a/phys2denoise/metrics/retroicor.py
+++ b/phys2denoise/metrics/retroicor.py
@@ -10,3 +10,78 @@
 import matplotlib.pyplot as plt
 from scipy.signal import argrelmax
 
+"""
+This function computes RETROICOR regressors (Glover et al. 2000) 
+"""
+
+
+def compute_card_phase(card_peaks_timings,slice_timings,nscans,TR):
+    """
+
+    This function creates cardiac phase from cardiac peaks.
+    Assumes that timing of cardiac events are given in same units
+    as slice timings, for example seconds.
+    
+    """
+
+    nscans = np.shape(slice_timings)
+    phase_card = np.zeros(nscans)
+    for ii in range(nscans):
+        # find previous cardiac peaks
+        previous_card_peaks = np.asarray(np.nonzero(card_peaks_timings < slice_timings[ii]))
+        if np.size(previous_card_peaks) == 0:
+            t1 = 0
+        else:
+            last_peak = previous_card_peaks[0][-1]
+            t1 = card_peaks_timings[last_peak]
+        
+        # find posterior cardiac peaks
+        next_card_peaks = np.asarray(np.nonzero(card_peaks_timings > slice_timings[ii]))
+        if np.size(next_card_peaks) == 0:
+            t2 = nscans * TR
+        else:
+            next_peak = next_card_peaks[0][0]
+            t2 = card_peaks_timings[next_peak]
+        
+        # compute cardiac phase
+        phase_card[ii] = (2*np.math.pi*(slice_timings[ii] - t1))/(t2-t1)
+
+    return phase_card
+
+def compute_resp_phase(resp,sampling_time):
+
+    """
+    This function creates respiration phase from resp trace. 
+    """
+
+
+def compute_retroicor_regressors(physio,TR,nscans,slice_timings,n_harmonics,card=FALSE,resp=FALSE):
+    nslices = np.shape(slice_timings) # number of slices
+
+    # if respiration, compute histogram and temporal derivative of respiration signal
+    if resp:
+        resp_hist, resp_hist_bins = plt.hist(physio,bins=100)
+        resp_diff = np.diff(physio,n=1)
+    
+    #initialize output variables
+    retroicor_regressors = []
+    phase = np.empty((nscans,nslices))
+
+    for jj in range(nslices):
+        # Initialize slice timings for current slice
+        crslice_timings = TR * np.arange(nscans)+slice_timings[jj]
+        # Compute physiological phases using the timings of physio events (e.g. peaks) slice sampling times
+        if card:
+            phase[,jj] = compute_phase_card(physio,crslice_timings)
+        if resp:
+            phase[,jj] = compute_phase_resp(resp_diff,resp_hist,resp_hist_bins,crslice_timings)
+        # Compute retroicor regressors
+        for nn in range(n_harmonics):
+            retroicor_regressors[jj][:,2*nn] = np.cos((nn+1)*phase[jj])
+            retricor_regressor[jj][:,2*nn+1] = np.sin((nn+1)*phase[jj])
+
+    return retroicor_regressors,phase
+
+
+        
+

From d877f83877f5fd53de3c53c4d75db32369534b3b Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Fri, 13 Nov 2020 09:41:08 -0500
Subject: [PATCH 03/18] Add duecredit and some references.

---
 phys2denoise/due.py        | 74 ++++++++++++++++++++++++++++++++++++++
 phys2denoise/references.py | 10 ++++++
 2 files changed, 84 insertions(+)
 create mode 100644 phys2denoise/due.py
 create mode 100644 phys2denoise/references.py

diff --git a/phys2denoise/due.py b/phys2denoise/due.py
new file mode 100644
index 0000000..9a1c4dd
--- /dev/null
+++ b/phys2denoise/due.py
@@ -0,0 +1,74 @@
+# emacs: at the end of the file
+# ex: set sts=4 ts=4 sw=4 et:
+# ## ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### #
+"""
+
+Stub file for a guaranteed safe import of duecredit constructs:  if duecredit
+is not available.
+
+To use it, place it into your project codebase to be imported, e.g. copy as
+
+    cp stub.py /path/tomodule/module/due.py
+
+Note that it might be better to avoid naming it duecredit.py to avoid shadowing
+installed duecredit.
+
+Then use in your code as
+
+    from .due import due, Doi, BibTeX, Text
+
+See  https://github.com/duecredit/duecredit/blob/master/README.md for examples.
+
+Origin:     Originally a part of the duecredit
+Copyright:  2015-2019  DueCredit developers
+License:    BSD-2
+"""
+
+__version__ = '0.0.8'
+
+
+class InactiveDueCreditCollector(object):
+    """Just a stub at the Collector which would not do anything"""
+    def _donothing(self, *args, **kwargs):
+        """Perform no good and no bad"""
+        pass
+
+    def dcite(self, *args, **kwargs):
+        """If I could cite I would"""
+        def nondecorating_decorator(func):
+            return func
+        return nondecorating_decorator
+
+    active = False
+    activate = add = cite = dump = load = _donothing
+
+    def __repr__(self):
+        return self.__class__.__name__ + '()'
+
+
+def _donothing_func(*args, **kwargs):
+    """Perform no good and no bad"""
+    pass
+
+
+try:
+    from duecredit import due, BibTeX, Doi, Url, Text
+    if 'due' in locals() and not hasattr(due, 'cite'):
+        raise RuntimeError(
+            "Imported due lacks .cite. DueCredit is now disabled")
+except Exception as e:
+    if not isinstance(e, ImportError):
+        import logging
+        logging.getLogger("duecredit").error(
+            "Failed to import duecredit due to %s" % str(e))
+    # Initiate due stub
+    due = InactiveDueCreditCollector()
+    BibTeX = Doi = Url = Text = _donothing_func
+
+# Emacs mode definitions
+# Local Variables:
+# mode: python
+# py-indent-offset: 4
+# tab-width: 4
+# indent-tabs-mode: nil
+# End:
diff --git a/phys2denoise/references.py b/phys2denoise/references.py
new file mode 100644
index 0000000..305d231
--- /dev/null
+++ b/phys2denoise/references.py
@@ -0,0 +1,10 @@
+"""
+References to be imported and injected throughout the package
+"""
+from .due import Doi
+
+CHANG_GLOVER_2009 = Doi('10.1016/j.neuroimage.2009.04.048')
+
+POWER_2018 = Doi('10.1073/pnas.1720985115')
+
+POWER_2020 = Doi('10.1016/j.neuroimage.2019.116234')

From 869ea49a8680fc78ebff76ff59c75ebe7d1d3809 Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Fri, 13 Nov 2020 09:41:15 -0500
Subject: [PATCH 04/18] Draft some metrics.

---
 phys2denoise/metrics/cardiac.py    |  62 ++++++++
 phys2denoise/metrics/chest_belt.py | 222 +++++++++++++++++++++++++++++
 phys2denoise/metrics/gas.py        |   2 +
 phys2denoise/metrics/utils.py      |  93 ++++++++++++
 4 files changed, 379 insertions(+)
 create mode 100644 phys2denoise/metrics/cardiac.py
 create mode 100644 phys2denoise/metrics/chest_belt.py
 create mode 100644 phys2denoise/metrics/gas.py
 create mode 100644 phys2denoise/metrics/utils.py

diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py
new file mode 100644
index 0000000..a0a0b2b
--- /dev/null
+++ b/phys2denoise/metrics/cardiac.py
@@ -0,0 +1,62 @@
+"""Denoising metrics for cardio recordings
+"""
+import numpy as np
+
+from ..due import due
+from .. import references
+
+
+def iht():
+    """Instantaneous heart rate
+    """
+    pass
+
+
+@due.dcite(references.CHANG_GLOVER_2009)
+def crf(samplerate, oversampling=50, time_length=32, onset=0., tr=2.):
+    """
+    Calculate the cardiac response function using the definition
+    supplied in Chang and Glover (2009).
+
+    Inputs
+    ------
+    samplerate : :obj:`float`
+        Sampling rate of data, in seconds.
+    oversampling : :obj:`int`, optional
+        Temporal oversampling factor, in seconds. Default is 50.
+    time_length : :obj:`int`, optional
+        RRF kernel length, in seconds. Default is 32.
+    onset : :obj:`float`, optional
+        Onset of the response, in seconds. Default is 0.
+
+    Outputs
+    -------
+    crf: array-like
+        Cardiac or "heart" response function
+
+    Notes
+    -----
+    This cardiac response function was defined in [1]_, Appendix A.
+
+    The core code for this function comes from metco2, while several of the
+    parameters, including oversampling, time_length, and onset, are modeled on
+    nistats' HRF functions.
+
+    References
+    ----------
+    .. [1] C. Chang & G. H. Glover, "Relationship between respiration,
+       end-tidal CO2, and BOLD signals in resting-state fMRI," Neuroimage,
+       issue 47, vol. 4, pp. 1381-1393, 2009.
+    """
+    def _crf(t):
+        rf = (0.6 * t ** 2.7 * np.exp(-t / 1.6) - 16 * (1 / np.sqrt(2 * np.pi * 9))
+              * np.exp(-0.5 * (((t - 12) ** 2) / 9)))
+        return rf
+
+    dt = tr / oversampling
+    time_stamps = np.linspace(0, time_length,
+                              np.rint(float(time_length) / dt).astype(np.int))
+    time_stamps -= onset
+    crf_arr = _crf(time_stamps)
+    crf_arr = crf_arr / max(abs(crf_arr))
+    return crf_arr
diff --git a/phys2denoise/metrics/chest_belt.py b/phys2denoise/metrics/chest_belt.py
new file mode 100644
index 0000000..f488a3d
--- /dev/null
+++ b/phys2denoise/metrics/chest_belt.py
@@ -0,0 +1,222 @@
+"""Denoising metrics for chest belt recordings
+"""
+import numpy as np
+import pandas as pd
+from scipy.ndimage.filters import convolve1d
+from scipy.signal import resample, detrend
+from scipy.stats import zscore
+
+from . import utils
+from ..due import due
+from .. import references
+
+
+@due.dcite(references.POWER_2018)
+def rpv(belt_ts, window):
+    """Respiratory pattern variability
+
+    Parameters
+    ----------
+    belt_ts
+    window
+
+    Returns
+    -------
+    rpv_arr
+
+    Notes
+    -----
+    This metric was first introduced in [1]_.
+
+    1. Z-score respiratory belt signal
+    2. Calculate upper envelope
+    3. Calculate standard deviation of envelope
+
+    References
+    ----------
+    .. [1] J. D. Power et al., "Ridding fMRI data of motion-related influences:
+       Removal of signals with distinct spatial and physical bases in multiecho
+       data," Proceedings of the National Academy of Sciences, issue 9, vol.
+       115, pp. 2105-2114, 2018.
+    """
+    # First, z-score respiratory traces
+    resp_z = zscore(belt_ts)
+
+    # Collect upper envelope
+    rpv_upper_env = utils.rms_envelope_1d(resp_z, window)
+
+    # Calculate standard deviation
+    rpv_arr = np.std(rpv_upper_env)
+    return rpv_arr
+
+
+@due.dcite(references.POWER_2020)
+def env(belt_ts, samplerate, out_samplerate, window=10, lags=(0,)):
+    """Respiratory pattern variability calculated across a sliding window
+
+    Parameters
+    ----------
+    belt_ts : (X,) :obj:`numpy.ndarray`
+        A 1D array with the respiratory belt time series.
+    samplerate : :obj:`float`
+        Sampling rate for belt_ts, in Hertz.
+    out_samplerate : :obj:`float`
+        Sampling rate for the output time series in seconds.
+        Corresponds to TR in fMRI data.
+    window : :obj:`int`, optional
+        Size of the sliding window, in the same units as out_samplerate.
+        Default is 6.
+    lags : (Y,) :obj:`tuple` of :obj:`int`, optional
+        List of lags to apply to the rv estimate. In the same units as
+        out_samplerate.
+
+    Returns
+    -------
+    env_arr
+
+    Notes
+    -----
+    This metric was first introduced in [1]_.
+
+    Across a sliding window, do the following:
+    1. Z-score respiratory belt signal
+    2. Calculate upper envelope
+    3. Calculate standard deviation of envelope
+
+    References
+    ----------
+    .. [1] J. D. Power et al., "Characteristics of respiratory measures in
+       young adults scanned at rest, including systematic changes and 'missed'
+       deep breaths," Neuroimage, vol. 204, 2020.
+    """
+    window = window * samplerate / out_samplerate
+    # Calculate RPV across a rolling window
+    env_arr = pd.Series(belt_ts).rolling(window=window, center=True).apply(
+        rpv, window=window)
+    env_arr[np.isnan(env_arr)] = 0.
+    return env_arr
+
+
+@due.dcite(references.CHANG_GLOVER_2009)
+def rv(belt_ts, samplerate, out_samplerate, window=6, lags=(0,)):
+    """Respiratory variance
+
+    Parameters
+    ----------
+    belt_ts : (X,) :obj:`numpy.ndarray`
+        A 1D array with the respiratory belt time series.
+    samplerate : :obj:`float`
+        Sampling rate for belt_ts, in Hertz.
+    out_samplerate : :obj:`float`
+        Sampling rate for the output time series in seconds.
+        Corresponds to TR in fMRI data.
+    window : :obj:`int`, optional
+        Size of the sliding window, in the same units as out_samplerate.
+        Default is 6.
+    lags : (Y,) :obj:`tuple` of :obj:`int`, optional
+        List of lags to apply to the rv estimate. In the same units as
+        out_samplerate.
+
+    Returns
+    -------
+    rv_out : (Z, 2Y) :obj:`numpy.ndarray`
+        Respiratory variance values, with lags applied, downsampled to
+        out_samplerate, convolved with an RRF, and detrended/normalized.
+        The first Y columns are not convolved with the RRF, while the second Y
+        columns are.
+
+    Notes
+    -----
+    Respiratory variance (RV) was introduced in [1]_, and consists of the
+    standard deviation of the respiratory trace within a 6-second window.
+
+    This metric is often lagged back and/or forward in time and convolved with
+    a respiratory response function before being included in a GLM.
+    Regressors also often have mean and linear trends removed and are
+    standardized prior to regressions.
+
+    References
+    ----------
+    .. [1] C. Chang & G. H. Glover, "Relationship between respiration,
+       end-tidal CO2, and BOLD signals in resting-state fMRI," Neuroimage,
+       issue 4, vol. 47, pp. 1381-1393, 2009.
+    """
+    # Raw respiratory variance
+    rv_arr = pd.Series(belt_ts).rolling(window=window, center=True).std()
+    rv_arr[np.isnan(rv_arr)] = 0.
+
+    # Apply lags
+    n_out_samples = int((belt_ts.shape[0] / samplerate) / out_samplerate)
+    # convert lags from out_samplerate to samplerate
+    delays = [abs(int(lag * samplerate)) for lag in lags]
+    rv_with_lags = utils.apply_lags(rv_arr, lags=delays)
+
+    # Downsample to out_samplerate
+    rv_with_lags = resample(rv_with_lags, num=n_out_samples, axis=0)
+
+    # Convolve with rrf
+    rrf_arr = rrf(out_samplerate, oversampling=1)
+    rv_convolved = convolve1d(rv_with_lags, rrf_arr, axis=0)
+
+    # Concatenate the raw and convolved versions
+    rv_combined = np.hstack((rv_with_lags, rv_convolved))
+
+    # Detrend and normalize
+    rv_combined = rv_combined - np.mean(rv_combined, axis=0)
+    rv_combined = detrend(rv_combined, axis=0)
+    rv_out = zscore(rv_combined, axis=0)
+    return rv_out
+
+
+def rvt(belt_ts, samplerate, out_samplerate, window=10, lags=(0,)):
+    """Respiratory volume-per-time
+    """
+    pass
+
+
+@due.dcite(references.CHANG_GLOVER_2009)
+def rrf(samplerate, oversampling=50, time_length=50, onset=0., tr=2.):
+    """
+    Calculate the respiratory response function using the definition
+    supplied in Chang and Glover (2009).
+
+    Inputs
+    ------
+    samplerate : :obj:`float`
+        Sampling rate of data, in seconds.
+    oversampling : :obj:`int`, optional
+        Temporal oversampling factor, in seconds. Default is 50.
+    time_length : :obj:`int`, optional
+        RRF kernel length, in seconds. Default is 50.
+    onset : :obj:`float`, optional
+        Onset of the response, in seconds. Default is 0.
+
+    Outputs
+    -------
+    rrf: array-like
+        respiratory response function
+
+    Notes
+    -----
+    This respiratory response function was defined in [1]_, Appendix A.
+
+    The core code for this function comes from metco2, while several of the
+    parameters, including oversampling, time_length, and onset, are modeled on
+    nistats' HRF functions.
+
+    References
+    ----------
+    .. [1] C. Chang & G. H. Glover, "Relationship between respiration,
+       end-tidal CO2, and BOLD signals in resting-state fMRI," Neuroimage,
+       issue 47, vol. 4, pp. 1381-1393, 2009.
+    """
+    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
+    dt = tr / oversampling
+    time_stamps = np.linspace(0, time_length,
+                              np.rint(float(time_length) / dt).astype(np.int))
+    time_stamps -= onset
+    rrf_arr = _rrf(time_stamps)
+    rrf_arr = rrf_arr / max(abs(rrf_arr))
+    return rrf_arr
diff --git a/phys2denoise/metrics/gas.py b/phys2denoise/metrics/gas.py
new file mode 100644
index 0000000..12ea419
--- /dev/null
+++ b/phys2denoise/metrics/gas.py
@@ -0,0 +1,2 @@
+"""Denoising metrics for gas recordings
+"""
diff --git a/phys2denoise/metrics/utils.py b/phys2denoise/metrics/utils.py
new file mode 100644
index 0000000..a55976b
--- /dev/null
+++ b/phys2denoise/metrics/utils.py
@@ -0,0 +1,93 @@
+import numpy as np
+
+
+def mirrorpad_1d(arr, buffer=250):
+    """
+    Pad both sides of array with flipped values from array of length 'buffer'.
+
+    Parameters
+    ----------
+    arr
+    buffer
+
+    Returns
+    -------
+    arr_out
+    """
+    mirror = np.flip(arr, axis=0)
+    idx = range(arr.shape[0] - buffer, arr.shape[0])
+    pre_mirror = np.take(mirror, idx, axis=0)
+    idx = range(0, buffer)
+    post_mirror = np.take(mirror, idx, axis=0)
+    arr_out = np.concatenate((pre_mirror, arr, post_mirror), axis=0)
+    return arr_out
+
+
+def rms_envelope_1d(arr, window=500):
+    """
+    Conceptual translation of MATLAB 2017b's envelope(X, x, 'rms') function.
+    https://www.mathworks.com/help/signal/ref/envelope.html
+
+    Parameters
+    ----------
+    arr
+    window
+
+    Returns
+    -------
+    rms_env : numpy.ndarray
+        The upper envelope.
+    """
+    assert arr.ndim == 1, 'Input data must be 1D'
+    assert window % 2 == 0, 'Window must be even'
+    n_t = arr.shape[0]
+    buf = int(window / 2)
+
+    # Pad array at both ends
+    arr = np.copy(arr).astype(float)
+    mean = np.mean(arr)
+    arr -= mean
+    arr = mirrorpad_1d(arr, buffer=buf)
+    rms_env = np.empty(n_t)
+    for i in range(n_t):
+        # to match matlab
+        start_idx = i + buf
+        stop_idx = i + buf + window
+
+        # but this is probably more appropriate
+        # start_idx = i + buf - 1
+        # stop_idx = i + buf + window
+        window_arr = arr[start_idx:stop_idx]
+        rms = np.sqrt(np.mean(window_arr ** 2))
+        rms_env[i] = rms
+    rms_env += mean
+    return rms_env
+
+
+def apply_lags(arr1d, lags):
+    """
+    Apply delays (lags) to an array.
+
+    Parameters
+    ----------
+    arr1d : (X,) :obj:`numpy.ndarray`
+        One-dimensional array to apply delays to.
+    lags : (Y,) :obj:`tuple` or :obj:`int`
+        Delays, in the same units as arr1d, to apply to arr1d. Can be negative,
+        zero, or positive integers.
+
+    Returns
+    -------
+    arr_with_lags : (X, Y) :obj:`numpy.ndarray`
+        arr1d shifted according to lags. Each column corresponds to a lag.
+    """
+    arr_with_lags = np.zeros((arr1d.shape[0], len(lags)))
+    for i_lag, lag in enumerate(lags):
+        if lag < 0:
+            arr_delayed = np.hstack((arr1d[lag:], np.zeros(lag)))
+        elif lag > 0:
+            arr_delayed = np.hstack((np.zeros(lag), arr1d[lag:]))
+        else:
+            arr_delayed = arr1d.copy()
+        arr_with_lags[:, i_lag] = arr_delayed
+    return arr_with_lags

From 23f851fb1f842b6a88bd5ac987c6a14a5d30002d Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Wed, 30 Dec 2020 12:16:10 -0500
Subject: [PATCH 05/18] Run black and isort.

---
 phys2denoise/due.py                | 19 +++++++++++-------
 phys2denoise/metrics/cardiac.py    | 18 +++++++++--------
 phys2denoise/metrics/chest_belt.py | 32 ++++++++++++++----------------
 phys2denoise/metrics/utils.py      |  4 ++--
 phys2denoise/references.py         |  6 +++---
 5 files changed, 42 insertions(+), 37 deletions(-)

diff --git a/phys2denoise/due.py b/phys2denoise/due.py
index 9a1c4dd..6013ed7 100644
--- a/phys2denoise/due.py
+++ b/phys2denoise/due.py
@@ -24,26 +24,29 @@
 License:    BSD-2
 """
 
-__version__ = '0.0.8'
+__version__ = "0.0.8"
 
 
 class InactiveDueCreditCollector(object):
     """Just a stub at the Collector which would not do anything"""
+
     def _donothing(self, *args, **kwargs):
         """Perform no good and no bad"""
         pass
 
     def dcite(self, *args, **kwargs):
         """If I could cite I would"""
+
         def nondecorating_decorator(func):
             return func
+
         return nondecorating_decorator
 
     active = False
     activate = add = cite = dump = load = _donothing
 
     def __repr__(self):
-        return self.__class__.__name__ + '()'
+        return self.__class__.__name__ + "()"
 
 
 def _donothing_func(*args, **kwargs):
@@ -52,15 +55,17 @@ def _donothing_func(*args, **kwargs):
 
 
 try:
-    from duecredit import due, BibTeX, Doi, Url, Text
-    if 'due' in locals() and not hasattr(due, 'cite'):
-        raise RuntimeError(
-            "Imported due lacks .cite. DueCredit is now disabled")
+    from duecredit import BibTeX, Doi, Text, Url, due
+
+    if "due" in locals() and not hasattr(due, "cite"):
+        raise RuntimeError("Imported due lacks .cite. DueCredit is now disabled")
 except Exception as e:
     if not isinstance(e, ImportError):
         import logging
+
         logging.getLogger("duecredit").error(
-            "Failed to import duecredit due to %s" % str(e))
+            "Failed to import duecredit due to %s" % str(e)
+        )
     # Initiate due stub
     due = InactiveDueCreditCollector()
     BibTeX = Doi = Url = Text = _donothing_func
diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py
index a0a0b2b..19ab81a 100644
--- a/phys2denoise/metrics/cardiac.py
+++ b/phys2denoise/metrics/cardiac.py
@@ -2,18 +2,17 @@
 """
 import numpy as np
 
-from ..due import due
 from .. import references
+from ..due import due
 
 
 def iht():
-    """Instantaneous heart rate
-    """
+    """Instantaneous heart rate"""
     pass
 
 
 @due.dcite(references.CHANG_GLOVER_2009)
-def crf(samplerate, oversampling=50, time_length=32, onset=0., tr=2.):
+def crf(samplerate, oversampling=50, time_length=32, onset=0.0, tr=2.0):
     """
     Calculate the cardiac response function using the definition
     supplied in Chang and Glover (2009).
@@ -48,14 +47,17 @@ def crf(samplerate, oversampling=50, time_length=32, onset=0., tr=2.):
        end-tidal CO2, and BOLD signals in resting-state fMRI," Neuroimage,
        issue 47, vol. 4, pp. 1381-1393, 2009.
     """
+
     def _crf(t):
-        rf = (0.6 * t ** 2.7 * np.exp(-t / 1.6) - 16 * (1 / np.sqrt(2 * np.pi * 9))
-              * np.exp(-0.5 * (((t - 12) ** 2) / 9)))
+        rf = 0.6 * t ** 2.7 * np.exp(-t / 1.6) - 16 * (
+            1 / np.sqrt(2 * np.pi * 9)
+        ) * np.exp(-0.5 * (((t - 12) ** 2) / 9))
         return rf
 
     dt = tr / oversampling
-    time_stamps = np.linspace(0, time_length,
-                              np.rint(float(time_length) / dt).astype(np.int))
+    time_stamps = np.linspace(
+        0, time_length, np.rint(float(time_length) / dt).astype(np.int)
+    )
     time_stamps -= onset
     crf_arr = _crf(time_stamps)
     crf_arr = crf_arr / max(abs(crf_arr))
diff --git a/phys2denoise/metrics/chest_belt.py b/phys2denoise/metrics/chest_belt.py
index f488a3d..3d0e905 100644
--- a/phys2denoise/metrics/chest_belt.py
+++ b/phys2denoise/metrics/chest_belt.py
@@ -3,12 +3,12 @@
 import numpy as np
 import pandas as pd
 from scipy.ndimage.filters import convolve1d
-from scipy.signal import resample, detrend
+from scipy.signal import detrend, resample
 from scipy.stats import zscore
 
-from . import utils
-from ..due import due
 from .. import references
+from ..due import due
+from . import utils
 
 
 @due.dcite(references.POWER_2018)
@@ -91,9 +91,10 @@ def env(belt_ts, samplerate, out_samplerate, window=10, lags=(0,)):
     """
     window = window * samplerate / out_samplerate
     # Calculate RPV across a rolling window
-    env_arr = pd.Series(belt_ts).rolling(window=window, center=True).apply(
-        rpv, window=window)
-    env_arr[np.isnan(env_arr)] = 0.
+    env_arr = (
+        pd.Series(belt_ts).rolling(window=window, center=True).apply(rpv, window=window)
+    )
+    env_arr[np.isnan(env_arr)] = 0.0
     return env_arr
 
 
@@ -143,7 +144,7 @@ def rv(belt_ts, samplerate, out_samplerate, window=6, lags=(0,)):
     """
     # Raw respiratory variance
     rv_arr = pd.Series(belt_ts).rolling(window=window, center=True).std()
-    rv_arr[np.isnan(rv_arr)] = 0.
+    rv_arr[np.isnan(rv_arr)] = 0.0
 
     # Apply lags
     n_out_samples = int((belt_ts.shape[0] / samplerate) / out_samplerate)
@@ -168,14 +169,8 @@ def rv(belt_ts, samplerate, out_samplerate, window=6, lags=(0,)):
     return rv_out
 
 
-def rvt(belt_ts, samplerate, out_samplerate, window=10, lags=(0,)):
-    """Respiratory volume-per-time
-    """
-    pass
-
-
 @due.dcite(references.CHANG_GLOVER_2009)
-def rrf(samplerate, oversampling=50, time_length=50, onset=0., tr=2.):
+def rrf(samplerate, oversampling=50, time_length=50, onset=0.0, tr=2.0):
     """
     Calculate the respiratory response function using the definition
     supplied in Chang and Glover (2009).
@@ -210,12 +205,15 @@ def rrf(samplerate, oversampling=50, time_length=50, onset=0., tr=2.):
        end-tidal CO2, and BOLD signals in resting-state fMRI," Neuroimage,
        issue 47, vol. 4, pp. 1381-1393, 2009.
     """
+
     def _rrf(t):
-        rf = (0.6 * t ** 2.1 * np.exp(-t / 1.6) - 0.0023 * t ** 3.54 * np.exp(-t / 4.25))
+        rf = 0.6 * t ** 2.1 * np.exp(-t / 1.6) - 0.0023 * t ** 3.54 * np.exp(-t / 4.25)
         return rf
+
     dt = tr / oversampling
-    time_stamps = np.linspace(0, time_length,
-                              np.rint(float(time_length) / dt).astype(np.int))
+    time_stamps = np.linspace(
+        0, time_length, np.rint(float(time_length) / dt).astype(np.int)
+    )
     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 a55976b..e44c2a5 100644
--- a/phys2denoise/metrics/utils.py
+++ b/phys2denoise/metrics/utils.py
@@ -38,8 +38,8 @@ def rms_envelope_1d(arr, window=500):
     rms_env : numpy.ndarray
         The upper envelope.
     """
-    assert arr.ndim == 1, 'Input data must be 1D'
-    assert window % 2 == 0, 'Window must be even'
+    assert arr.ndim == 1, "Input data must be 1D"
+    assert window % 2 == 0, "Window must be even"
     n_t = arr.shape[0]
     buf = int(window / 2)
 
diff --git a/phys2denoise/references.py b/phys2denoise/references.py
index 305d231..c4b557f 100644
--- a/phys2denoise/references.py
+++ b/phys2denoise/references.py
@@ -3,8 +3,8 @@
 """
 from .due import Doi
 
-CHANG_GLOVER_2009 = Doi('10.1016/j.neuroimage.2009.04.048')
+CHANG_GLOVER_2009 = Doi("10.1016/j.neuroimage.2009.04.048")
 
-POWER_2018 = Doi('10.1073/pnas.1720985115')
+POWER_2018 = Doi("10.1073/pnas.1720985115")
 
-POWER_2020 = Doi('10.1016/j.neuroimage.2019.116234')
+POWER_2020 = Doi("10.1016/j.neuroimage.2019.116234")

From a4454e296535bb1093e5328279239c5ade3c617a Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Wed, 30 Dec 2020 12:17:47 -0500
Subject: [PATCH 06/18] Minor changes.

---
 phys2denoise/due.py        | 13 +++++--------
 phys2denoise/references.py |  4 +---
 2 files changed, 6 insertions(+), 11 deletions(-)

diff --git a/phys2denoise/due.py b/phys2denoise/due.py
index 6013ed7..17071d4 100644
--- a/phys2denoise/due.py
+++ b/phys2denoise/due.py
@@ -1,10 +1,7 @@
 # emacs: at the end of the file
 # ex: set sts=4 ts=4 sw=4 et:
 # ## ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### #
-"""
-
-Stub file for a guaranteed safe import of duecredit constructs:  if duecredit
-is not available.
+"""Stub file for a guaranteed safe import of duecredit constructs if duecredit is not available.
 
 To use it, place it into your project codebase to be imported, e.g. copy as
 
@@ -28,14 +25,14 @@
 
 
 class InactiveDueCreditCollector(object):
-    """Just a stub at the Collector which would not do anything"""
+    """Just a stub at the Collector which would not do anything."""
 
     def _donothing(self, *args, **kwargs):
-        """Perform no good and no bad"""
+        """Perform no good and no bad."""
         pass
 
     def dcite(self, *args, **kwargs):
-        """If I could cite I would"""
+        """If I could cite I would."""
 
         def nondecorating_decorator(func):
             return func
@@ -50,7 +47,7 @@ def __repr__(self):
 
 
 def _donothing_func(*args, **kwargs):
-    """Perform no good and no bad"""
+    """Perform no good and no bad."""
     pass
 
 
diff --git a/phys2denoise/references.py b/phys2denoise/references.py
index c4b557f..297ba3b 100644
--- a/phys2denoise/references.py
+++ b/phys2denoise/references.py
@@ -1,6 +1,4 @@
-"""
-References to be imported and injected throughout the package
-"""
+"""References to be imported and injected throughout the package."""
 from .due import Doi
 
 CHANG_GLOVER_2009 = Doi("10.1016/j.neuroimage.2009.04.048")

From 4612d4fe4a557dc88f256053baa400b076c8519c Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Wed, 30 Dec 2020 12:25:47 -0500
Subject: [PATCH 07/18] More fixing.

---
 phys2denoise/due.py                |  1 +
 phys2denoise/metrics/cardiac.py    |  9 +++------
 phys2denoise/metrics/chest_belt.py | 13 +++++--------
 phys2denoise/metrics/gas.py        |  3 +--
 phys2denoise/metrics/utils.py      | 15 ++++++++-------
 5 files changed, 18 insertions(+), 23 deletions(-)

diff --git a/phys2denoise/due.py b/phys2denoise/due.py
index 17071d4..4ff2457 100644
--- a/phys2denoise/due.py
+++ b/phys2denoise/due.py
@@ -43,6 +43,7 @@ def nondecorating_decorator(func):
     activate = add = cite = dump = load = _donothing
 
     def __repr__(self):
+        """Perform magic."""
         return self.__class__.__name__ + "()"
 
 
diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py
index 19ab81a..98a5f23 100644
--- a/phys2denoise/metrics/cardiac.py
+++ b/phys2denoise/metrics/cardiac.py
@@ -1,5 +1,4 @@
-"""Denoising metrics for cardio recordings
-"""
+"""Denoising metrics for cardio recordings."""
 import numpy as np
 
 from .. import references
@@ -7,15 +6,13 @@
 
 
 def iht():
-    """Instantaneous heart rate"""
+    """Calculate instantaneous heart rate."""
     pass
 
 
 @due.dcite(references.CHANG_GLOVER_2009)
 def crf(samplerate, oversampling=50, time_length=32, onset=0.0, tr=2.0):
-    """
-    Calculate the cardiac response function using the definition
-    supplied in Chang and Glover (2009).
+    """Calculate the cardiac response function using Chang and Glover's definition.
 
     Inputs
     ------
diff --git a/phys2denoise/metrics/chest_belt.py b/phys2denoise/metrics/chest_belt.py
index 3d0e905..dc24794 100644
--- a/phys2denoise/metrics/chest_belt.py
+++ b/phys2denoise/metrics/chest_belt.py
@@ -1,5 +1,4 @@
-"""Denoising metrics for chest belt recordings
-"""
+"""Denoising metrics for chest belt recordings."""
 import numpy as np
 import pandas as pd
 from scipy.ndimage.filters import convolve1d
@@ -13,7 +12,7 @@
 
 @due.dcite(references.POWER_2018)
 def rpv(belt_ts, window):
-    """Respiratory pattern variability
+    """Calculate respiratory pattern variability.
 
     Parameters
     ----------
@@ -52,7 +51,7 @@ def rpv(belt_ts, window):
 
 @due.dcite(references.POWER_2020)
 def env(belt_ts, samplerate, out_samplerate, window=10, lags=(0,)):
-    """Respiratory pattern variability calculated across a sliding window
+    """Calculate respiratory pattern variability across a sliding window.
 
     Parameters
     ----------
@@ -100,7 +99,7 @@ def env(belt_ts, samplerate, out_samplerate, window=10, lags=(0,)):
 
 @due.dcite(references.CHANG_GLOVER_2009)
 def rv(belt_ts, samplerate, out_samplerate, window=6, lags=(0,)):
-    """Respiratory variance
+    """Calculate respiratory variance.
 
     Parameters
     ----------
@@ -171,9 +170,7 @@ def rv(belt_ts, samplerate, out_samplerate, window=6, lags=(0,)):
 
 @due.dcite(references.CHANG_GLOVER_2009)
 def rrf(samplerate, oversampling=50, time_length=50, onset=0.0, tr=2.0):
-    """
-    Calculate the respiratory response function using the definition
-    supplied in Chang and Glover (2009).
+    """Calculate the respiratory response function using Chang and Glover's definition.
 
     Inputs
     ------
diff --git a/phys2denoise/metrics/gas.py b/phys2denoise/metrics/gas.py
index 12ea419..a4e3a59 100644
--- a/phys2denoise/metrics/gas.py
+++ b/phys2denoise/metrics/gas.py
@@ -1,2 +1 @@
-"""Denoising metrics for gas recordings
-"""
+"""Denoising metrics for gas recordings."""
diff --git a/phys2denoise/metrics/utils.py b/phys2denoise/metrics/utils.py
index e44c2a5..5369081 100644
--- a/phys2denoise/metrics/utils.py
+++ b/phys2denoise/metrics/utils.py
@@ -1,9 +1,9 @@
+"""Miscellaneous utility functions for metric calculation."""
 import numpy as np
 
 
 def mirrorpad_1d(arr, buffer=250):
-    """
-    Pad both sides of array with flipped values from array of length 'buffer'.
+    """Pad both sides of array with flipped values from array of length 'buffer'.
 
     Parameters
     ----------
@@ -24,9 +24,7 @@ def mirrorpad_1d(arr, buffer=250):
 
 
 def rms_envelope_1d(arr, window=500):
-    """
-    Conceptual translation of MATLAB 2017b's envelope(X, x, 'rms') function.
-    https://www.mathworks.com/help/signal/ref/envelope.html
+    """Conceptual translation of MATLAB 2017b's envelope(X, x, 'rms') function.
 
     Parameters
     ----------
@@ -37,6 +35,10 @@ def rms_envelope_1d(arr, window=500):
     -------
     rms_env : numpy.ndarray
         The upper envelope.
+
+    Notes
+    -----
+    https://www.mathworks.com/help/signal/ref/envelope.html
     """
     assert arr.ndim == 1, "Input data must be 1D"
     assert window % 2 == 0, "Window must be even"
@@ -65,8 +67,7 @@ def rms_envelope_1d(arr, window=500):
 
 
 def apply_lags(arr1d, lags):
-    """
-    Apply delays (lags) to an array.
+    """Apply delays (lags) to an array.
 
     Parameters
     ----------

From 4202c2ec425d1c4eca34a35d1f956ef9af431d03 Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Wed, 30 Dec 2020 14:26:02 -0500
Subject: [PATCH 08/18] Fix crf/rrf docstrings.

---
 phys2denoise/metrics/cardiac.py    | 6 +++---
 phys2denoise/metrics/chest_belt.py | 8 ++++----
 2 files changed, 7 insertions(+), 7 deletions(-)

diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py
index 98a5f23..df13d4d 100644
--- a/phys2denoise/metrics/cardiac.py
+++ b/phys2denoise/metrics/cardiac.py
@@ -14,8 +14,8 @@ def iht():
 def crf(samplerate, oversampling=50, time_length=32, onset=0.0, tr=2.0):
     """Calculate the cardiac response function using Chang and Glover's definition.
 
-    Inputs
-    ------
+    Parameters
+    ----------
     samplerate : :obj:`float`
         Sampling rate of data, in seconds.
     oversampling : :obj:`int`, optional
@@ -25,7 +25,7 @@ def crf(samplerate, oversampling=50, time_length=32, onset=0.0, tr=2.0):
     onset : :obj:`float`, optional
         Onset of the response, in seconds. Default is 0.
 
-    Outputs
+    Returns
     -------
     crf: array-like
         Cardiac or "heart" response function
diff --git a/phys2denoise/metrics/chest_belt.py b/phys2denoise/metrics/chest_belt.py
index dc24794..26c3c3f 100644
--- a/phys2denoise/metrics/chest_belt.py
+++ b/phys2denoise/metrics/chest_belt.py
@@ -172,18 +172,18 @@ def rv(belt_ts, samplerate, out_samplerate, window=6, lags=(0,)):
 def rrf(samplerate, oversampling=50, time_length=50, onset=0.0, tr=2.0):
     """Calculate the respiratory response function using Chang and Glover's definition.
 
-    Inputs
-    ------
+    Parameters
+    ----------
     samplerate : :obj:`float`
         Sampling rate of data, in seconds.
     oversampling : :obj:`int`, optional
-        Temporal oversampling factor, in seconds. Default is 50.
+        Temporal oversampling factor. Default is 50.
     time_length : :obj:`int`, optional
         RRF kernel length, in seconds. Default is 50.
     onset : :obj:`float`, optional
         Onset of the response, in seconds. Default is 0.
 
-    Outputs
+    Returns
     -------
     rrf: array-like
         respiratory response function

From e2116dc0951a01b2ddad7f8e5f66d80d87dca545 Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Wed, 30 Dec 2020 14:26:58 -0500
Subject: [PATCH 09/18] Fix again.

---
 phys2denoise/metrics/cardiac.py    | 2 +-
 phys2denoise/metrics/chest_belt.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py
index df13d4d..f295848 100644
--- a/phys2denoise/metrics/cardiac.py
+++ b/phys2denoise/metrics/cardiac.py
@@ -27,7 +27,7 @@ def crf(samplerate, oversampling=50, time_length=32, onset=0.0, tr=2.0):
 
     Returns
     -------
-    crf: array-like
+    crf : array-like
         Cardiac or "heart" response function
 
     Notes
diff --git a/phys2denoise/metrics/chest_belt.py b/phys2denoise/metrics/chest_belt.py
index 26c3c3f..9ff8225 100644
--- a/phys2denoise/metrics/chest_belt.py
+++ b/phys2denoise/metrics/chest_belt.py
@@ -185,7 +185,7 @@ def rrf(samplerate, oversampling=50, time_length=50, onset=0.0, tr=2.0):
 
     Returns
     -------
-    rrf: array-like
+    rrf : array-like
         respiratory response function
 
     Notes

From 889d3e57b9f17eedf127703a57138741f0c04de6 Mon Sep 17 00:00:00 2001
From: Stefano Moia <s.moia@bcbl.eu>
Date: Thu, 25 Feb 2021 15:52:02 +0100
Subject: [PATCH 10/18] Move metric call function to metrics.utils

---
 phys2denoise/metrics/utils.py | 31 +++++++++++++++++++++++++++++++
 1 file changed, 31 insertions(+)

diff --git a/phys2denoise/metrics/utils.py b/phys2denoise/metrics/utils.py
index 5369081..7426a5c 100644
--- a/phys2denoise/metrics/utils.py
+++ b/phys2denoise/metrics/utils.py
@@ -1,6 +1,37 @@
 """Miscellaneous utility functions for metric calculation."""
+import logging
+
 import numpy as np
 
+LGR = logging.getLogger(__name__)
+LGR.setLevel(logging.INFO)
+
+
+def print_metric_call(metric, args):
+    """
+    Log a message to describe how a metric is being called.
+
+    Parameters
+    ----------
+    metric : function
+        Metric function that is being called
+    args : dict
+        Dictionary containing all arguments that are used to parametrise metric
+
+    Notes
+    -----
+    Outcome
+        An info-level message for the logger.
+    """
+    msg = f'The {metric} regressor will be computed using the following parameters:'
+
+    for arg in args:
+        msg = f'{msg}\n    {arg} = {args[arg]}'
+
+    msg = f'{msg}\n'
+
+    LGR.info(msg)
+
 
 def mirrorpad_1d(arr, buffer=250):
     """Pad both sides of array with flipped values from array of length 'buffer'.

From 6d8d0e90c1965c8b34f9d4eec50f4e42dbda9d5e Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Thu, 25 Feb 2021 11:08:59 -0500
Subject: [PATCH 11/18] Refactor and document RETROICOR code a bit.

---
 phys2denoise/metrics/retroicor.py | 147 ++++++++++++++++++++----------
 phys2denoise/references.py        |   2 +
 2 files changed, 101 insertions(+), 48 deletions(-)

diff --git a/phys2denoise/metrics/retroicor.py b/phys2denoise/metrics/retroicor.py
index e63e929..8b9fa6e 100644
--- a/phys2denoise/metrics/retroicor.py
+++ b/phys2denoise/metrics/retroicor.py
@@ -1,87 +1,138 @@
-"""implementation of retroicor"""
+"""These functions compute RETROICOR regressors (Glover et al. 2000)."""
 
-import argparse as ap
 import numpy as np
-import scipy as sp
-import json
-import string
-import random
-import matplotlib as mpl; mpl.use('TkAgg')
+import matplotlib as mpl
+
+mpl.use("TkAgg")
 import matplotlib.pyplot as plt
-from scipy.signal import argrelmax
 
-"""
-This function computes RETROICOR regressors (Glover et al. 2000) 
-"""
+from .. import references
+from ..due import due
 
 
-def compute_card_phase(card_peaks_timings,slice_timings,nscans,TR):
-    """
+def compute_phase_card(card_peaks_timings, slice_timings, n_scans, t_r):
+    """Calculate cardiac phase from cardiac peaks.
 
-    This function creates cardiac phase from cardiac peaks.
     Assumes that timing of cardiac events are given in same units
     as slice timings, for example seconds.
-    
-    """
 
-    nscans = np.shape(slice_timings)
-    phase_card = np.zeros(nscans)
-    for ii in range(nscans):
+    Parameters
+    ----------
+    peaks : 1D array_like
+    slice_timings : 1D array_like
+    n_scans : int
+    t_r : float
+
+    Returns
+    -------
+    phase_card : array_like
+        Cardiac phase signal.
+    """
+    n_scans = np.shape(slice_timings)
+    phase_card = np.zeros(n_scans)
+    for ii in range(n_scans):
         # find previous cardiac peaks
-        previous_card_peaks = np.asarray(np.nonzero(card_peaks_timings < slice_timings[ii]))
+        previous_card_peaks = np.asarray(
+            np.nonzero(card_peaks_timings < slice_timings[ii])
+        )
         if np.size(previous_card_peaks) == 0:
             t1 = 0
         else:
             last_peak = previous_card_peaks[0][-1]
             t1 = card_peaks_timings[last_peak]
-        
+
         # find posterior cardiac peaks
         next_card_peaks = np.asarray(np.nonzero(card_peaks_timings > slice_timings[ii]))
         if np.size(next_card_peaks) == 0:
-            t2 = nscans * TR
+            t2 = n_scans * t_r
         else:
             next_peak = next_card_peaks[0][0]
             t2 = card_peaks_timings[next_peak]
-        
+
         # compute cardiac phase
-        phase_card[ii] = (2*np.math.pi*(slice_timings[ii] - t1))/(t2-t1)
+        phase_card[ii] = (2 * np.math.pi * (slice_timings[ii] - t1)) / (t2 - t1)
 
     return phase_card
 
-def compute_resp_phase(resp,sampling_time):
 
-    """
-    This function creates respiration phase from resp trace. 
-    """
+def compute_phase_resp(resp, sampling_time):
+    """Calculate respiratory phase from respiratory signal.
 
+    Parameters
+    ----------
+    resp
+    sampling_time
 
-def compute_retroicor_regressors(physio,TR,nscans,slice_timings,n_harmonics,card=FALSE,resp=FALSE):
-    nslices = np.shape(slice_timings) # number of slices
+    Returns
+    -------
+    phase_resp : array_like
+        Respiratory phase signal.
+    """
+    pass
+
+
+@due.dcite(references.GLOVER_2000)
+def compute_retroicor_regressors(
+    physio, t_r, n_scans, slice_timings, n_harmonics, card=False, resp=False
+):
+    """Compute RETROICOR regressors.
+
+    Parameters
+    ----------
+    physio : array_like
+        1D array, whether cardiac or respiratory signal.
+    t_r : float
+    n_scans : int
+    slice_timings
+    n_harmonics : int
+    card : bool, optional
+    resp : bool, optional
+
+    Returns
+    -------
+    retroicor_regressors : list
+    phase : array_like
+        2D array of shape (n_scans, n_slices)
+
+    References
+    ----------
+    *   Glover, G. H., Li, T. Q., & Ress, D. (2000).
+        Image‐based method for retrospective correction of physiological
+        motion effects in fMRI: RETROICOR.
+        Magnetic Resonance in Medicine:
+        An Official Journal of the International Society for Magnetic Resonance in Medicine,
+        44(1), 162-167.
+    """
+    n_slices = np.shape(slice_timings)  # number of slices
 
     # if respiration, compute histogram and temporal derivative of respiration signal
     if resp:
-        resp_hist, resp_hist_bins = plt.hist(physio,bins=100)
-        resp_diff = np.diff(physio,n=1)
-    
-    #initialize output variables
+        # TODO: Replace with numpy.histogram
+        resp_hist, resp_hist_bins = plt.hist(physio, bins=100)
+        resp_diff = np.diff(physio, n=1)
+
+    # initialize output variables
     retroicor_regressors = []
-    phase = np.empty((nscans,nslices))
+    phase = np.empty((n_scans, n_slices))
 
-    for jj in range(nslices):
+    for i_slice in range(n_slices):
         # Initialize slice timings for current slice
-        crslice_timings = TR * np.arange(nscans)+slice_timings[jj]
+        crslice_timings = t_r * np.arange(n_scans) + slice_timings[i_slice]
         # Compute physiological phases using the timings of physio events (e.g. peaks) slice sampling times
         if card:
-            phase[,jj] = compute_phase_card(physio,crslice_timings)
+            phase[:, i_slice] = compute_phase_card(physio, crslice_timings)
         if resp:
-            phase[,jj] = compute_phase_resp(resp_diff,resp_hist,resp_hist_bins,crslice_timings)
-        # Compute retroicor regressors
-        for nn in range(n_harmonics):
-            retroicor_regressors[jj][:,2*nn] = np.cos((nn+1)*phase[jj])
-            retricor_regressor[jj][:,2*nn+1] = np.sin((nn+1)*phase[jj])
-
-    return retroicor_regressors,phase
-
-
-        
+            phase[:, i_slice] = compute_phase_resp(
+                resp_diff, resp_hist, resp_hist_bins, crslice_timings
+            )
 
+        # 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]
+            )
+
+    return retroicor_regressors, phase
diff --git a/phys2denoise/references.py b/phys2denoise/references.py
index 297ba3b..eeb914a 100644
--- a/phys2denoise/references.py
+++ b/phys2denoise/references.py
@@ -6,3 +6,5 @@
 POWER_2018 = Doi("10.1073/pnas.1720985115")
 
 POWER_2020 = Doi("10.1016/j.neuroimage.2019.116234")
+
+GLOVER_2000 = Doi("10.1002/1522-2594(200007)44:1<162::AID-MRM23>3.0.CO;2-E")

From b98e19c245c2ea10f20140a43ef973a9df32f816 Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Thu, 25 Feb 2021 11:37:10 -0500
Subject: [PATCH 12/18] Resolve missing arguments and some bugs.

---
 phys2denoise/metrics/retroicor.py | 105 +++++++++++++++++++++++-------
 1 file changed, 81 insertions(+), 24 deletions(-)

diff --git a/phys2denoise/metrics/retroicor.py b/phys2denoise/metrics/retroicor.py
index 8b9fa6e..deb57d9 100644
--- a/phys2denoise/metrics/retroicor.py
+++ b/phys2denoise/metrics/retroicor.py
@@ -10,7 +10,7 @@
 from ..due import due
 
 
-def compute_phase_card(card_peaks_timings, slice_timings, n_scans, t_r):
+def compute_phase_card(peaks, slice_timings, n_scans, t_r):
     """Calculate cardiac phase from cardiac peaks.
 
     Assumes that timing of cardiac events are given in same units
@@ -19,61 +19,108 @@ def compute_phase_card(card_peaks_timings, slice_timings, n_scans, t_r):
     Parameters
     ----------
     peaks : 1D array_like
+        Cardiac peak times, in seconds.
     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.
 
     Returns
     -------
     phase_card : array_like
         Cardiac phase signal.
     """
-    n_scans = np.shape(slice_timings)
+    n_slices = np.shape(slice_timings)
     phase_card = np.zeros(n_scans)
-    for ii in range(n_scans):
+
+    for i_slice in range(n_slices):
         # find previous cardiac peaks
-        previous_card_peaks = np.asarray(
-            np.nonzero(card_peaks_timings < slice_timings[ii])
-        )
+        previous_card_peaks = np.asarray(np.nonzero(peaks < slice_timings[i_slice]))
         if np.size(previous_card_peaks) == 0:
             t1 = 0
         else:
             last_peak = previous_card_peaks[0][-1]
-            t1 = card_peaks_timings[last_peak]
+            t1 = peaks[last_peak]
 
         # find posterior cardiac peaks
-        next_card_peaks = np.asarray(np.nonzero(card_peaks_timings > slice_timings[ii]))
+        next_card_peaks = np.asarray(np.nonzero(peaks > slice_timings[i_slice]))
         if np.size(next_card_peaks) == 0:
             t2 = n_scans * t_r
         else:
             next_peak = next_card_peaks[0][0]
-            t2 = card_peaks_timings[next_peak]
+            t2 = peaks[next_peak]
 
         # compute cardiac phase
-        phase_card[ii] = (2 * np.math.pi * (slice_timings[ii] - t1)) / (t2 - t1)
+        phase_card[i_slice] = (2 * np.math.pi * (slice_timings[i_slice] - t1)) / (
+            t2 - t1
+        )
 
     return phase_card
 
 
-def compute_phase_resp(resp, sampling_time):
+def compute_phase_resp(resp, sample_rate, n_scans, slice_timings, t_r):
     """Calculate respiratory phase from respiratory signal.
 
     Parameters
     ----------
-    resp
-    sampling_time
+    resp : 1D array_like
+        Respiratory signal.
+    sample_rate : float
+        Sample rate of physio, in Hertz.
+    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.
 
     Returns
     -------
     phase_resp : array_like
         Respiratory phase signal.
     """
-    pass
+    n_slices = np.shape(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 = plt.hist(resp, bins=100)
+
+    # first compute derivative of respiration signal
+    resp_diff = np.diff(resp, 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] * sample_rate)])
+            )  # 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))
+            numerator = np.sum(resp_hist[0:thisBin])
+            phase_resp_crSlice[j_scan] = (
+                np.math.pi * np.sign(resp_diff[iphys]) * (numerator / len(resp))
+            )
+
+        phase_resp[:, i_slice] = phase_resp_crSlice
+
+    return phase_resp
 
 
 @due.dcite(references.GLOVER_2000)
 def compute_retroicor_regressors(
-    physio, t_r, n_scans, slice_timings, n_harmonics, card=False, resp=False
+    physio,
+    sample_rate,
+    t_r,
+    n_scans,
+    slice_timings,
+    n_harmonics,
+    card=False,
+    resp=False,
 ):
     """Compute RETROICOR regressors.
 
@@ -81,12 +128,22 @@ def compute_retroicor_regressors(
     ----------
     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.
     t_r : float
+        Imaging sample rate, in seconds.
     n_scans : int
-    slice_timings
+        Number of volumes in the imaging run.
+    slice_timings : 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
     -------
@@ -105,12 +162,6 @@ def compute_retroicor_regressors(
     """
     n_slices = np.shape(slice_timings)  # number of slices
 
-    # if respiration, compute histogram and temporal derivative of respiration signal
-    if resp:
-        # TODO: Replace with numpy.histogram
-        resp_hist, resp_hist_bins = plt.hist(physio, bins=100)
-        resp_diff = np.diff(physio, n=1)
-
     # initialize output variables
     retroicor_regressors = []
     phase = np.empty((n_scans, n_slices))
@@ -118,12 +169,18 @@ def compute_retroicor_regressors(
     for i_slice in range(n_slices):
         # Initialize slice timings for current slice
         crslice_timings = t_r * np.arange(n_scans) + slice_timings[i_slice]
-        # Compute physiological phases using the timings of physio events (e.g. peaks) slice sampling times
+
+        # Compute physiological phases using the timings of physio events (e.g. peaks)
+        # slice sampling times
         if card:
             phase[:, i_slice] = compute_phase_card(physio, crslice_timings)
         if resp:
             phase[:, i_slice] = compute_phase_resp(
-                resp_diff, resp_hist, resp_hist_bins, crslice_timings
+                physio,
+                sample_rate,
+                n_scans,
+                slice_timings,
+                t_r,
             )
 
         # Compute retroicor regressors

From 58b908c9597c92e325ba60859e4a7e32dfc8fd32 Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Thu, 25 Feb 2021 11:37:15 -0500
Subject: [PATCH 13/18] Rename function.

---
 phys2denoise/metrics/retroicor.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/phys2denoise/metrics/retroicor.py b/phys2denoise/metrics/retroicor.py
index deb57d9..a99d3cc 100644
--- a/phys2denoise/metrics/retroicor.py
+++ b/phys2denoise/metrics/retroicor.py
@@ -112,7 +112,7 @@ def compute_phase_resp(resp, sample_rate, n_scans, slice_timings, t_r):
 
 
 @due.dcite(references.GLOVER_2000)
-def compute_retroicor_regressors(
+def retroicor(
     physio,
     sample_rate,
     t_r,

From d6e198469055839b00a6c6f76ba04b374b1a33dd Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Thu, 25 Feb 2021 11:51:07 -0500
Subject: [PATCH 14/18] Fix function call.

---
 phys2denoise/metrics/retroicor.py | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/phys2denoise/metrics/retroicor.py b/phys2denoise/metrics/retroicor.py
index a99d3cc..7072501 100644
--- a/phys2denoise/metrics/retroicor.py
+++ b/phys2denoise/metrics/retroicor.py
@@ -173,7 +173,13 @@ def retroicor(
         # Compute physiological phases using the timings of physio events (e.g. peaks)
         # slice sampling times
         if card:
-            phase[:, i_slice] = compute_phase_card(physio, crslice_timings)
+            phase[:, i_slice] = compute_phase_card(
+                physio,
+                crslice_timings,
+                n_scans,
+                t_r,
+            )
+
         if resp:
             phase[:, i_slice] = compute_phase_resp(
                 physio,

From 58f41ba7b58b1c0a6364bd80680e5cab930d65d6 Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Thu, 25 Feb 2021 12:33:16 -0500
Subject: [PATCH 15/18] Move phase functions to appropriate files.

---
 phys2denoise/metrics/cardiac.py    |  50 +++++++++++++
 phys2denoise/metrics/chest_belt.py |  55 ++++++++++++++
 phys2denoise/metrics/retroicor.py  | 111 ++---------------------------
 3 files changed, 109 insertions(+), 107 deletions(-)

diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py
index f295848..43c8079 100644
--- a/phys2denoise/metrics/cardiac.py
+++ b/phys2denoise/metrics/cardiac.py
@@ -59,3 +59,53 @@ def _crf(t):
     crf_arr = _crf(time_stamps)
     crf_arr = crf_arr / max(abs(crf_arr))
     return crf_arr
+
+
+def cardiac_phase(peaks, slice_timings, n_scans, t_r):
+    """Calculate cardiac phase from cardiac peaks.
+
+    Assumes that timing of cardiac events are given in same units
+    as slice timings, for example seconds.
+
+    Parameters
+    ----------
+    peaks : 1D array_like
+        Cardiac peak times, in seconds.
+    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.
+
+    Returns
+    -------
+    phase_card : array_like
+        Cardiac phase signal.
+    """
+    n_slices = np.shape(slice_timings)
+    phase_card = np.zeros(n_scans)
+
+    for i_slice in range(n_slices):
+        # find previous cardiac peaks
+        previous_card_peaks = np.asarray(np.nonzero(peaks < slice_timings[i_slice]))
+        if np.size(previous_card_peaks) == 0:
+            t1 = 0
+        else:
+            last_peak = previous_card_peaks[0][-1]
+            t1 = peaks[last_peak]
+
+        # find posterior cardiac peaks
+        next_card_peaks = np.asarray(np.nonzero(peaks > slice_timings[i_slice]))
+        if np.size(next_card_peaks) == 0:
+            t2 = n_scans * t_r
+        else:
+            next_peak = next_card_peaks[0][0]
+            t2 = peaks[next_peak]
+
+        # compute cardiac phase
+        phase_card[i_slice] = (2 * np.math.pi * (slice_timings[i_slice] - t1)) / (
+            t2 - t1
+        )
+
+    return phase_card
diff --git a/phys2denoise/metrics/chest_belt.py b/phys2denoise/metrics/chest_belt.py
index 9ff8225..7c5ae1c 100644
--- a/phys2denoise/metrics/chest_belt.py
+++ b/phys2denoise/metrics/chest_belt.py
@@ -1,10 +1,14 @@
 """Denoising metrics for chest belt recordings."""
+import matplotlib as mpl
 import numpy as np
 import pandas as pd
 from scipy.ndimage.filters import convolve1d
 from scipy.signal import detrend, resample
 from scipy.stats import zscore
 
+mpl.use("TkAgg")
+import matplotlib.pyplot as plt
+
 from .. import references
 from ..due import due
 from . import utils
@@ -215,3 +219,54 @@ def _rrf(t):
     rrf_arr = _rrf(time_stamps)
     rrf_arr = rrf_arr / max(abs(rrf_arr))
     return rrf_arr
+
+
+def respiratory_phase(resp, sample_rate, n_scans, slice_timings, t_r):
+    """Calculate respiratory phase from respiratory signal.
+
+    Parameters
+    ----------
+    resp : 1D array_like
+        Respiratory signal.
+    sample_rate : float
+        Sample rate of physio, in Hertz.
+    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.
+
+    Returns
+    -------
+    phase_resp : array_like
+        Respiratory phase signal.
+    """
+    n_slices = np.shape(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 = plt.hist(resp, bins=100)
+
+    # first compute derivative of respiration signal
+    resp_diff = np.diff(resp, 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] * sample_rate)])
+            )  # 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))
+            numerator = np.sum(resp_hist[0:thisBin])
+            phase_resp_crSlice[j_scan] = (
+                np.math.pi * np.sign(resp_diff[iphys]) * (numerator / len(resp))
+            )
+
+        phase_resp[:, i_slice] = phase_resp_crSlice
+
+    return phase_resp
diff --git a/phys2denoise/metrics/retroicor.py b/phys2denoise/metrics/retroicor.py
index 7072501..690b797 100644
--- a/phys2denoise/metrics/retroicor.py
+++ b/phys2denoise/metrics/retroicor.py
@@ -1,116 +1,13 @@
 """These functions compute RETROICOR regressors (Glover et al. 2000)."""
 
 import numpy as np
-import matplotlib as mpl
-
-mpl.use("TkAgg")
-import matplotlib.pyplot as plt
 
+from .cardiac import cardiac_phase
+from .chest_belt import respiratory_phase
 from .. import references
 from ..due import due
 
 
-def compute_phase_card(peaks, slice_timings, n_scans, t_r):
-    """Calculate cardiac phase from cardiac peaks.
-
-    Assumes that timing of cardiac events are given in same units
-    as slice timings, for example seconds.
-
-    Parameters
-    ----------
-    peaks : 1D array_like
-        Cardiac peak times, in seconds.
-    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.
-
-    Returns
-    -------
-    phase_card : array_like
-        Cardiac phase signal.
-    """
-    n_slices = np.shape(slice_timings)
-    phase_card = np.zeros(n_scans)
-
-    for i_slice in range(n_slices):
-        # find previous cardiac peaks
-        previous_card_peaks = np.asarray(np.nonzero(peaks < slice_timings[i_slice]))
-        if np.size(previous_card_peaks) == 0:
-            t1 = 0
-        else:
-            last_peak = previous_card_peaks[0][-1]
-            t1 = peaks[last_peak]
-
-        # find posterior cardiac peaks
-        next_card_peaks = np.asarray(np.nonzero(peaks > slice_timings[i_slice]))
-        if np.size(next_card_peaks) == 0:
-            t2 = n_scans * t_r
-        else:
-            next_peak = next_card_peaks[0][0]
-            t2 = peaks[next_peak]
-
-        # compute cardiac phase
-        phase_card[i_slice] = (2 * np.math.pi * (slice_timings[i_slice] - t1)) / (
-            t2 - t1
-        )
-
-    return phase_card
-
-
-def compute_phase_resp(resp, sample_rate, n_scans, slice_timings, t_r):
-    """Calculate respiratory phase from respiratory signal.
-
-    Parameters
-    ----------
-    resp : 1D array_like
-        Respiratory signal.
-    sample_rate : float
-        Sample rate of physio, in Hertz.
-    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.
-
-    Returns
-    -------
-    phase_resp : array_like
-        Respiratory phase signal.
-    """
-    n_slices = np.shape(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 = plt.hist(resp, bins=100)
-
-    # first compute derivative of respiration signal
-    resp_diff = np.diff(resp, 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] * sample_rate)])
-            )  # 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))
-            numerator = np.sum(resp_hist[0:thisBin])
-            phase_resp_crSlice[j_scan] = (
-                np.math.pi * np.sign(resp_diff[iphys]) * (numerator / len(resp))
-            )
-
-        phase_resp[:, i_slice] = phase_resp_crSlice
-
-    return phase_resp
-
-
 @due.dcite(references.GLOVER_2000)
 def retroicor(
     physio,
@@ -173,7 +70,7 @@ def retroicor(
         # Compute physiological phases using the timings of physio events (e.g. peaks)
         # slice sampling times
         if card:
-            phase[:, i_slice] = compute_phase_card(
+            phase[:, i_slice] = cardiac_phase(
                 physio,
                 crslice_timings,
                 n_scans,
@@ -181,7 +78,7 @@ def retroicor(
             )
 
         if resp:
-            phase[:, i_slice] = compute_phase_resp(
+            phase[:, i_slice] = respiratory_phase(
                 physio,
                 sample_rate,
                 n_scans,

From e0b4a1166228ba2d8c193f0a6e825fe5288bcd73 Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Thu, 25 Feb 2021 12:33:25 -0500
Subject: [PATCH 16/18] Fill out init file.

---
 phys2denoise/metrics/__init__.py | 15 +++++++++++++++
 1 file changed, 15 insertions(+)

diff --git a/phys2denoise/metrics/__init__.py b/phys2denoise/metrics/__init__.py
index e69de29..1337dec 100644
--- a/phys2denoise/metrics/__init__.py
+++ b/phys2denoise/metrics/__init__.py
@@ -0,0 +1,15 @@
+"""Metrics derived from physiological signals."""
+from .cardiac import crf, cardiac_phase
+from .chest_belt import rpv, env, rv, rrf, respiratory_phase
+from .retroicor import retroicor
+
+__all__ = [
+    "crf",
+    "cardiac_phase",
+    "rpv",
+    "env",
+    "rv",
+    "rrf",
+    "respiratory_phase",
+    "retroicor",
+]

From 6508a5efaee6ceaea3f773fba80d59a042230a74 Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Thu, 25 Feb 2021 12:46:32 -0500
Subject: [PATCH 17/18] Preallocate RETROICOR contents.

---
 phys2denoise/metrics/retroicor.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/phys2denoise/metrics/retroicor.py b/phys2denoise/metrics/retroicor.py
index 690b797..f4ae089 100644
--- a/phys2denoise/metrics/retroicor.py
+++ b/phys2denoise/metrics/retroicor.py
@@ -64,6 +64,8 @@ def retroicor(
     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]
 

From e729bc386bc329fc4d6d3bbe5ab25d347c2155ca Mon Sep 17 00:00:00 2001
From: Taylor Salo <tsalo006@fiu.edu>
Date: Thu, 25 Feb 2021 13:09:38 -0500
Subject: [PATCH 18/18] Reorg retroicor.

---
 phys2denoise/metrics/__init__.py                     | 2 +-
 phys2denoise/metrics/{retroicor.py => multimodal.py} | 5 +++++
 2 files changed, 6 insertions(+), 1 deletion(-)
 rename phys2denoise/metrics/{retroicor.py => multimodal.py} (94%)

diff --git a/phys2denoise/metrics/__init__.py b/phys2denoise/metrics/__init__.py
index 1337dec..8eb43c8 100644
--- a/phys2denoise/metrics/__init__.py
+++ b/phys2denoise/metrics/__init__.py
@@ -1,7 +1,7 @@
 """Metrics derived from physiological signals."""
 from .cardiac import crf, cardiac_phase
 from .chest_belt import rpv, env, rv, rrf, respiratory_phase
-from .retroicor import retroicor
+from .multimodal import retroicor
 
 __all__ = [
     "crf",
diff --git a/phys2denoise/metrics/retroicor.py b/phys2denoise/metrics/multimodal.py
similarity index 94%
rename from phys2denoise/metrics/retroicor.py
rename to phys2denoise/metrics/multimodal.py
index f4ae089..fb95eae 100644
--- a/phys2denoise/metrics/retroicor.py
+++ b/phys2denoise/metrics/multimodal.py
@@ -48,6 +48,11 @@ def retroicor(
     phase : array_like
         2D array of shape (n_scans, n_slices)
 
+    Notes
+    -----
+    RETROICOR regressors should be regressed from the imaging data *before*
+    any other preprocessing, including slice-timing correction and motion correction.
+
     References
     ----------
     *   Glover, G. H., Li, T. Q., & Ress, D. (2000).