diff --git a/phys2denoise/bids.py b/phys2denoise/bids.py
new file mode 100644
index 0000000..73546ae
--- /dev/null
+++ b/phys2denoise/bids.py
@@ -0,0 +1,2 @@
+"""Functions for indexing and manipulating BIDS datasets.
+"""
diff --git a/phys2denoise/cli/app.py b/phys2denoise/cli/app.py
new file mode 100644
index 0000000..46a0c1b
--- /dev/null
+++ b/phys2denoise/cli/app.py
@@ -0,0 +1,2 @@
+"""Command-line interface for the phys2denoise BIDS App
+"""
diff --git a/phys2denoise/cli/phys2denoise.py b/phys2denoise/cli/phys2denoise.py
new file mode 100644
index 0000000..9fd3fdf
--- /dev/null
+++ b/phys2denoise/cli/phys2denoise.py
@@ -0,0 +1,2 @@
+"""General-purpose phys2denoise command-line interface
+"""
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/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
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')
diff --git a/phys2denoise/utils.py b/phys2denoise/utils.py
new file mode 100644
index 0000000..924d3c9
--- /dev/null
+++ b/phys2denoise/utils.py
@@ -0,0 +1,2 @@
+"""General utility functions for phys2denoise.
+"""
diff --git a/requirements.txt b/requirements.txt
deleted file mode 100644
index 3d6f17a..0000000
--- a/requirements.txt
+++ /dev/null
@@ -1,2 +0,0 @@
-numpy>=1.9.3
-matplotlib>=3.1.1
\ No newline at end of file
diff --git a/setup.cfg b/setup.cfg
index 324284b..30d6607 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -23,7 +23,9 @@ provides =
 python_requires = >=3.6.1
 install_requires =
     numpy >=1.9.3
-    matplotlib >=3.1.1
+    pandas
+    duecredit
+    scipy
 tests_require =
     pytest >=3.6
 test_suite = pytest
@@ -53,7 +55,6 @@ all =
 [options.package_data]
 abagen =
     phys2denoise/data/*
-    phys2denoise/heuristics/*
     phys2denoise/tests/data/*
 
 [options.entry_points]