diff --git a/phys2denoise/due.py b/phys2denoise/due.py
new file mode 100644
index 0000000..4ff2457
--- /dev/null
+++ b/phys2denoise/due.py
@@ -0,0 +1,77 @@
+# 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):
+        """Perform magic."""
+        return self.__class__.__name__ + "()"
+
+
+def _donothing_func(*args, **kwargs):
+    """Perform no good and no bad."""
+    pass
+
+
+try:
+    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)
+        )
+    # 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/metrics/__init__.py b/phys2denoise/metrics/__init__.py
index e69de29..8eb43c8 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 .multimodal import retroicor
+
+__all__ = [
+    "crf",
+    "cardiac_phase",
+    "rpv",
+    "env",
+    "rv",
+    "rrf",
+    "respiratory_phase",
+    "retroicor",
+]
diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py
new file mode 100644
index 0000000..43c8079
--- /dev/null
+++ b/phys2denoise/metrics/cardiac.py
@@ -0,0 +1,111 @@
+"""Denoising metrics for cardio recordings."""
+import numpy as np
+
+from .. import references
+from ..due import due
+
+
+def iht():
+    """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 Chang and Glover's definition.
+
+    Parameters
+    ----------
+    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.
+
+    Returns
+    -------
+    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
+
+
+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
new file mode 100644
index 0000000..7c5ae1c
--- /dev/null
+++ b/phys2denoise/metrics/chest_belt.py
@@ -0,0 +1,272 @@
+"""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
+
+
+@due.dcite(references.POWER_2018)
+def rpv(belt_ts, window):
+    """Calculate 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,)):
+    """Calculate respiratory pattern variability 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.0
+    return env_arr
+
+
+@due.dcite(references.CHANG_GLOVER_2009)
+def rv(belt_ts, samplerate, out_samplerate, window=6, lags=(0,)):
+    """Calculate 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.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
+
+
+@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 Chang and Glover's definition.
+
+    Parameters
+    ----------
+    samplerate : :obj:`float`
+        Sampling rate of data, in seconds.
+    oversampling : :obj:`int`, optional
+        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.
+
+    Returns
+    -------
+    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
+
+
+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/gas.py b/phys2denoise/metrics/gas.py
new file mode 100644
index 0000000..a4e3a59
--- /dev/null
+++ b/phys2denoise/metrics/gas.py
@@ -0,0 +1 @@
+"""Denoising metrics for gas recordings."""
diff --git a/phys2denoise/metrics/multimodal.py b/phys2denoise/metrics/multimodal.py
new file mode 100644
index 0000000..fb95eae
--- /dev/null
+++ b/phys2denoise/metrics/multimodal.py
@@ -0,0 +1,105 @@
+"""These functions compute RETROICOR regressors (Glover et al. 2000)."""
+
+import numpy as np
+
+from .cardiac import cardiac_phase
+from .chest_belt import respiratory_phase
+from .. import references
+from ..due import due
+
+
+@due.dcite(references.GLOVER_2000)
+def retroicor(
+    physio,
+    sample_rate,
+    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.
+        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
+        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
+    -------
+    retroicor_regressors : list
+    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).
+        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
+
+    # 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]
+
+        # Compute physiological phases using the timings of physio events (e.g. peaks)
+        # slice sampling times
+        if card:
+            phase[:, i_slice] = cardiac_phase(
+                physio,
+                crslice_timings,
+                n_scans,
+                t_r,
+            )
+
+        if resp:
+            phase[:, i_slice] = respiratory_phase(
+                physio,
+                sample_rate,
+                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]
+            )
+
+    return retroicor_regressors, phase
diff --git a/phys2denoise/metrics/utils.py b/phys2denoise/metrics/utils.py
new file mode 100644
index 0000000..7426a5c
--- /dev/null
+++ b/phys2denoise/metrics/utils.py
@@ -0,0 +1,125 @@
+"""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'.
+
+    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.
+
+    Parameters
+    ----------
+    arr
+    window
+
+    Returns
+    -------
+    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"
+    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
diff --git a/phys2denoise/references.py b/phys2denoise/references.py
new file mode 100644
index 0000000..eeb914a
--- /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")
+
+GLOVER_2000 = Doi("10.1002/1522-2594(200007)44:1<162::AID-MRM23>3.0.CO;2-E")