Repository: PennLINC/xcp_d

Issue:
Adding to the Methods boilerplate that cosine regressors are added in aCompCor
When the aCompCor denoising strategy is chosen, xcp_d adds cosine regressors to the confounds, but this is not mentioned in the Methods boilerplate in the final report. This is equivalent to applying a high-pass filter to the data.

I think it would be good to add this, because users might not be aware that these extra regressors are added. 

See this relevant thread: https://neurostars.org/t/high-pass-filtering-using-xcp-d/24528/6



Discussion:
codecov[bot]: # [Codecov](https://codecov.io/gh/PennLINC/xcp_d/pull/788?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennLINC) Report
Base: **70.31**% // Head: **70.31**% // No change to project coverage :thumbsup:
> Coverage data is based on head [(`c39e1f5`)](https://codecov.io/gh/PennLINC/xcp_d/pull/788?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennLINC) compared to base [(`39f949e`)](https://codecov.io/gh/PennLINC/xcp_d/commit/39f949ee09164e09527e680427121607f6d43d5d?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennLINC).
> Patch has no changes to coverable lines.

<details><summary>Additional details and impacted files</summary>


```diff
@@           Coverage Diff           @@
##             main     #788   +/-   ##
=======================================
  Coverage   70.31%   70.31%           
=======================================
  Files          71       71           
  Lines        6465     6465           
=======================================
  Hits         4546     4546           
  Misses       1919     1919           
```


| [Impacted Files](https://codecov.io/gh/PennLINC/xcp_d/pull/788?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennLINC) | Coverage Δ | |
|---|---|---|
| [xcp\_d/utils/confounds.py](https://codecov.io/gh/PennLINC/xcp_d/pull/788?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennLINC#diff-eGNwX2QvdXRpbHMvY29uZm91bmRzLnB5) | `78.53% <ø> (ø)` | |

Help us with your feedback. Take ten seconds to tell us [how you rate us](https://about.codecov.io/nps?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennLINC). Have a feature suggestion? [Share it here.](https://app.codecov.io/gh/feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennLINC)

</details>

[:umbrella: View full report at Codecov](https://codecov.io/gh/PennLINC/xcp_d/pull/788?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennLINC).   
:loudspeaker: Do you have feedback about the report comment? [Let us know in this issue](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennLINC).

tsalo: @leotozzi88 WDYT of the proposed change?
leotozzi88: I think it would be beneficial to also explain when the data is filtered (at which stage of the denoising pipeline) and what filter is used (i.e. how much is the bandpass?). Also, it might be useful to say explicitly that adding the cosine regressors is equivalent to doing this bandpass filtering.  
tsalo: >I think it would be beneficial to also explain when the data is filtered (at which stage of the denoising pipeline) and what filter is used (i.e. how much is the bandpass?). 

That is mentioned when the bandpass filter is applied. As an example, [here](https://output.circle-artifacts.com/output/job/5841e567-6a78-458e-b733-a73beec8f070/artifacts/0/src/xcp_d/.circleci/out/test_ds001419_nifti/xcp_d/sub-01.html#boilerplate) is the boilerplate for one of the tests. Filtering is mentioned after interpolation.

>Also, it might be useful to say explicitly that adding the cosine regressors is equivalent to doing this bandpass filtering.

AFAIK, the cosine regressors act as a high-pass filter, rather than a bandpass filter, so they're not equivalent. I'd love to put the cutoff value in the boilerplate, but I don't know that it's something that can be extracted from the fMRIPrep derivatives. The only place I've found reference to the exact cutoff is in the fMRIPrep code, where it appears to be hardcoded. I'm not really comfortable hardcoding a value into the boilerplate, when it could be changed in the fMRIPrep source code without me noticing.

Patch:
diff --git a/xcp_d/interfaces/ants.py b/xcp_d/interfaces/ants.py
index abacdb2..953e480 100644
--- a/xcp_d/interfaces/ants.py
+++ b/xcp_d/interfaces/ants.py
@@ -12,6 +12,12 @@ from nipype.interfaces.base import (
     TraitedSpec,
     traits,
 )
+from niworkflows.interfaces.fixes import (
+    FixHeaderApplyTransforms,
+    _FixTraitApplyTransformsInputSpec,
+)
+
+from xcp_d.utils.filemanip import fname_presuffix
 
 LOGGER = logging.getLogger("nipype.interface")
 
@@ -147,3 +153,45 @@ class CompositeInvTransformUtil(ANTSCommand):
         if self.inputs.process == "assemble":
             outputs["out_file"] = os.path.abspath(self.inputs.out_file)
         return outputs
+
+
+class _ApplyTransformsInputSpec(_FixTraitApplyTransformsInputSpec):
+    # Nipype's version doesn't have GenericLabel
+    interpolation = traits.Enum(
+        "Linear",
+        "NearestNeighbor",
+        "CosineWindowedSinc",
+        "WelchWindowedSinc",
+        "HammingWindowedSinc",
+        "LanczosWindowedSinc",
+        "MultiLabel",
+        "Gaussian",
+        "BSpline",
+        "GenericLabel",
+        argstr="%s",
+        usedefault=True,
+    )
+
+
+class ApplyTransforms(FixHeaderApplyTransforms):
+    """A modified version of FixHeaderApplyTransforms from niworkflows.
+
+    The niworkflows version of ApplyTransforms "fixes the resampled image header
+    to match the xform of the reference image".
+    This modification overrides the allowed interpolation values,
+    since FixHeaderApplyTransforms doesn't support GenericLabel,
+    which is preferred over MultiLabel.
+    """
+
+    input_spec = _ApplyTransformsInputSpec
+
+    def _run_interface(self, runtime):
+        # Run normally
+        self.inputs.output_image = fname_presuffix(
+            self.inputs.input_image,
+            suffix="_trans.nii.gz",
+            newpath=runtime.cwd,
+            use_ext=False,
+        )
+        runtime = super(ApplyTransforms, self)._run_interface(runtime)
+        return runtime
diff --git a/xcp_d/interfaces/connectivity.py b/xcp_d/interfaces/connectivity.py
index 0d42c05..9de2f80 100644
--- a/xcp_d/interfaces/connectivity.py
+++ b/xcp_d/interfaces/connectivity.py
@@ -7,7 +7,6 @@ import numpy as np
 from nilearn.input_data import NiftiLabelsMasker
 from nilearn.plotting import plot_matrix
 from nipype import logging
-from nipype.interfaces.ants.resampling import ApplyTransforms, ApplyTransformsInputSpec
 from nipype.interfaces.base import (
     BaseInterfaceInputSpec,
     File,
@@ -149,33 +148,6 @@ class NiftiConnect(SimpleInterface):
         return runtime
 
 
-class _ApplyTransformsInputSpec(ApplyTransformsInputSpec):
-    transforms = InputMultiObject(
-        traits.Either(File(exists=True), "identity"),
-        argstr="%s",
-        mandatory=True,
-        desc="transform files",
-    )
-
-
-class ApplyTransformsx(ApplyTransforms):
-    """ApplyTransforms from nipype as workflow.
-
-    This is a modification of the ApplyTransforms interface,
-    with an updated set of inputs and a different default output image name.
-    """
-
-    input_spec = _ApplyTransformsInputSpec
-
-    def _run_interface(self, runtime):
-        # Run normally
-        self.inputs.output_image = fname_presuffix(
-            self.inputs.input_image, suffix="_trans.nii.gz", newpath=runtime.cwd, use_ext=False
-        )
-        runtime = super(ApplyTransformsx, self)._run_interface(runtime)
-        return runtime
-
-
 class _ConnectPlotInputSpec(BaseInterfaceInputSpec):
     in_file = File(exists=True, mandatory=True, desc="bold file")
     atlas_names = InputMultiObject(
diff --git a/xcp_d/interfaces/filtering.py b/xcp_d/interfaces/filtering.py
index 1658cc2..bd18478 100644
--- a/xcp_d/interfaces/filtering.py
+++ b/xcp_d/interfaces/filtering.py
@@ -43,12 +43,12 @@ class FilteringData(SimpleInterface):
         # filter the data
         if self.inputs.bandpass_filter:
             filt_data = butter_bandpass(
-                data=data_matrix,
+                data=data_matrix.T,
                 fs=1 / self.inputs.TR,
                 lowpass=self.inputs.lowpass,
                 highpass=self.inputs.highpass,
                 order=self.inputs.filter_order,
-            )
+            ).T
         else:
             filt_data = data_matrix  # no filtering!
 
diff --git a/xcp_d/utils/confounds.py b/xcp_d/utils/confounds.py
index aadc207..71b351c 100644
--- a/xcp_d/utils/confounds.py
+++ b/xcp_d/utils/confounds.py
@@ -8,7 +8,7 @@ import numpy as np
 import pandas as pd
 from nilearn.interfaces.fmriprep import load_confounds
 from nipype import logging
-from scipy.signal import filtfilt, firwin, iirnotch
+from scipy.signal import butter, filtfilt, iirnotch
 
 from xcp_d.utils.doc import fill_doc
 
@@ -128,7 +128,8 @@ def load_motion(
     ----------
     .. footbibliography::
     """
-    assert motion_filter_type in ("lp", "notch", None), motion_filter_type
+    if motion_filter_type not in ("lp", "notch", None):
+        raise ValueError(f"Motion filter type '{motion_filter_type}' not supported.")
 
     # Select the motion columns from the overall confounds DataFrame
     motion_confounds_df = confounds_df[
@@ -137,18 +138,14 @@ def load_motion(
 
     # Apply LP or notch filter
     if motion_filter_type in ("lp", "notch"):
-        # TODO: Eliminate need for transpose. We control the filter function,
-        # so we can make it work on RxT data instead of TxR.
-        motion_confounds = motion_confounds_df.to_numpy().T
         motion_confounds = motion_regression_filter(
-            data=motion_confounds,
+            data=motion_confounds_df.to_numpy(),
             TR=TR,
             motion_filter_type=motion_filter_type,
             band_stop_min=band_stop_min,
             band_stop_max=band_stop_max,
             motion_filter_order=motion_filter_order,
         )
-        motion_confounds = motion_confounds.T  # Transpose motion confounds back to RxT
         motion_confounds_df = pd.DataFrame(
             data=motion_confounds,
             columns=motion_confounds_df.columns,
@@ -284,6 +281,9 @@ def describe_regression(params, custom_confounds_file):
             "were selected as nuisance regressors [@behzadi2007component], "
             "along with the six motion parameters and their temporal derivatives "
             "[@benchmarkp;@satterthwaite_2013]. "
+            "As the aCompCor regressors were generated on high-pass filtered data, "
+            "the associated cosine basis regressors were included. "
+            "This has the effect of high-pass filtering the data as well."
         ),
         "acompcor_gsr": (
             "The top 5 aCompCor principal components from the WM and CSF compartments "
@@ -291,6 +291,9 @@ def describe_regression(params, custom_confounds_file):
             "along with the six motion parameters and their temporal derivatives, "
             "mean white matter signal, mean CSF signal, and mean global signal "
             "[@benchmarkp;@satterthwaite_2013]. "
+            "As the aCompCor regressors were generated on high-pass filtered data, "
+            "the associated cosine basis regressors were included. "
+            "This has the effect of high-pass filtering the data as well."
         ),
         "aroma": (
             "AROMA motion-labeled components [@pruim2015ica], mean white matter signal, "
@@ -348,6 +351,73 @@ def describe_regression(params, custom_confounds_file):
     return desc
 
 
+@fill_doc
+def describe_censoring(
+    motion_filter_type,
+    motion_filter_order,
+    band_stop_min,
+    band_stop_max,
+    head_radius,
+    fd_thresh,
+):
+    """Build a text description of the motion parameter filtering and FD censoring process.
+
+    Parameters
+    ----------
+    %(motion_filter_type)s
+    %(motion_filter_order)s
+    %(band_stop_min)s
+    %(band_stop_max)s
+    %(head_radius)s
+    %(fd_thresh)s
+
+    Returns
+    -------
+    desc : str
+        A text description of the censoring procedure.
+    """
+    from num2words import num2words
+
+    filter_str, filter_post_str = "", ""
+    if motion_filter_type:
+        if motion_filter_type == "notch":
+            filter_sub_str = (
+                f"band-stop filtered to remove signals between {band_stop_min} and "
+                f"{band_stop_max} breaths-per-minute using a(n) "
+                f"{num2words(motion_filter_order, ordinal=True)}-order notch filter, "
+                "based on @fair2020correction"
+            )
+        else:  # lp
+            filter_sub_str = (
+                f"low-pass filtered below {band_stop_min} breaths-per-minute using a(n) "
+                f"{num2words(motion_filter_order, ordinal=True)}-order Butterworth filter, "
+                "based on @gratton2020removal"
+            )
+
+        filter_str = (
+            f"the six translation and rotation head motion traces were {filter_sub_str}. Next, "
+        )
+        filter_post_str = (
+            "The filtered versions of the motion traces and framewise displacement were not used "
+            "for denoising."
+        )
+
+    if isinstance(head_radius, float):
+        fd_substr = f"{head_radius} mm"
+    else:
+        fd_substr = "estimated from the preprocessed brain mask"
+
+    desc = (
+        f"In order to identify high-motion outlier volumes, {filter_str}"
+        "framewise displacement was calculated using the formula from @power_fd_dvars, "
+        f"with a head radius {fd_substr}. "
+        f"Volumes with {'filtered ' if motion_filter_type else ''}framewise displacement "
+        f"greater than {fd_thresh} mm were flagged as outliers and excluded from nuisance "
+        f"regression [@power_fd_dvars]. {filter_post_str}"
+    )
+    return desc
+
+
 @fill_doc
 def load_confound_matrix(params, img_file, custom_confounds=None):
     """Load a subset of the confounds associated with a given file.
@@ -509,40 +579,48 @@ def motion_regression_filter(
 
     Parameters
     ----------
-    data : numpy.ndarray
-        Data to filter.
+    data : (T, R) numpy.ndarray
+        Data to filter. T = time, R = motion regressors
+        The filter will be applied independently to each variable, across time.
     TR : float
         Repetition time of the data.
-    motion_filter_type : {"lp", "notch"}
-        The type of motion filter to apply.
-        If "notch", the frequencies between ``band_stop_min`` and ``band_stop_max`` will be
-        removed.
-        If "lp", the frequencies above ``band_stop_min`` will be removed.
+    %(motion_filter_type)s
         If not "notch" or "lp", an exception will be raised.
     %(band_stop_min)s
     %(band_stop_max)s
-    motion_filter_order : int, optional
-        Default is 4.
+    %(motion_filter_order)s
 
     Returns
     -------
-    data : numpy.ndarray
-        Filtered data.
+    data : (T, R) numpy.ndarray
+        Filtered data. Same shape as the original data.
+
+    Notes
+    -----
+    Low-pass filtering (``motion_filter_type = "lp"``) is performed with a Butterworth filter,
+    as in :footcite:t:`gratton2020removal`.
+    The order of the Butterworth filter is determined by ``motion_filter_order``,
+    although the original paper used a first-order filter.
+    The original paper also used zero-padding with a padding size of 100.
+    We use constant-padding, with the default padding size determined by
+    :func:`scipy.signal.filtfilt`.
+
+    Band-stop filtering (``motion_filter_type = "notch"``) is performed with a notch filter,
+    as in :footcite:t:`fair2020correction`.
+    This filter uses the mean of the stopband frequencies as the target frequency,
+    and the range between the two frequencies as the bandwidth.
+    The filter is applied with constant-padding, using the default padding size determined by
+    :func:`scipy.signal.filtfilt`.
 
     References
     ----------
     .. footbibliography::
     """
-    assert motion_filter_type in ("lp", "notch")
-
-    # casting all variables
-    TR = float(TR)
-    order = float(motion_filter_order)
-
-    filtered_data = data.copy()
+    if motion_filter_type not in ("lp", "notch"):
+        raise ValueError(f"Motion filter type '{motion_filter_type}' not supported.")
 
-    sampling_frequency = 1.0 / TR
-    nyquist_frequency = sampling_frequency / 2.0
+    sampling_frequency = 1 / TR
+    nyquist_frequency = sampling_frequency / 2
 
     if motion_filter_type == "lp":  # low-pass filter
         # Remove any frequencies above band_stop_min.
@@ -551,19 +629,29 @@ def motion_regression_filter(
         if band_stop_max:
             warnings.warn("The parameter 'band_stop_max' will be ignored.")
 
-        low_pass_freq_hertz = band_stop_min / 60  # change BPM to right time unit
+        lowpass_hz = band_stop_min / 60  # change BPM to right time unit
 
-        # cutting frequency
-        cutting_frequency = np.abs(
-            low_pass_freq_hertz
-            - (np.floor((low_pass_freq_hertz + nyquist_frequency) / sampling_frequency))
+        # Adjust frequency in case Nyquist is below cutoff.
+        # This won't have an effect if the data have a fast enough sampling rate.
+        lowpass_hz_adjusted = np.abs(
+            lowpass_hz
+            - (np.floor((lowpass_hz + nyquist_frequency) / sampling_frequency))
             * sampling_frequency
         )
-        # cutting frequency normalized between 0 and nyquist
-        Wn = np.amin(cutting_frequency) / nyquist_frequency  # cutoffs
-        filt_num = firwin(int(order) + 1, Wn, pass_zero="lowpass")  # create b_filt
-        filt_denom = 1.0
-        num_f_apply = 1  # num of times to apply
+        if lowpass_hz_adjusted != lowpass_hz:
+            warnings.warn(
+                f"Low-pass filter frequency is above Nyquist frequency ({nyquist_frequency} Hz), "
+                f"so it has been changed ({lowpass_hz} --> {lowpass_hz_adjusted} Hz)."
+            )
+
+        b, a = butter(
+            motion_filter_order / 2,
+            lowpass_hz_adjusted,
+            btype="lowpass",
+            output="ba",
+            fs=sampling_frequency,
+        )
+        filtered_data = filtfilt(b, a, data, axis=0, padtype="constant", padlen=data.shape[0] - 1)
 
     elif motion_filter_type == "notch":  # notch filter
         # Retain any frequencies *outside* the band_stop_min-band_stop_max range.
@@ -573,27 +661,43 @@ def motion_regression_filter(
         assert band_stop_min > 0
         assert band_stop_min < band_stop_max
 
-        # bandwidth as an array
-        bandstop_band = np.array([band_stop_min, band_stop_max])
-        bandstop_band_hz = bandstop_band / 60  # change BPM to Hertz
-        cutting_frequencies = np.abs(
-            bandstop_band_hz
-            - (np.floor((bandstop_band_hz + nyquist_frequency) / sampling_frequency))
+        stopband = np.array([band_stop_min, band_stop_max])
+        stopband_hz = stopband / 60  # change BPM to Hertz
+
+        # Adjust frequencies in case Nyquist is within/below band.
+        # This won't have an effect if the data have a fast enough sampling rate.
+        stopband_hz_adjusted = np.abs(
+            stopband_hz
+            - (np.floor((stopband_hz + nyquist_frequency) / sampling_frequency))
             * sampling_frequency
         )
-
-        # normalize cutting frequency
-        W_notch = cutting_frequencies / nyquist_frequency
-        Wn = np.mean(W_notch)
-        Wd = np.diff(W_notch)
-        bandwidth = np.abs(Wd)  # bandwidth
-        # create filter coefficients
-        filt_num, filt_denom = iirnotch(Wn, Wn / bandwidth)
-        num_f_apply = np.int(np.floor(order / 2))  # how many times to apply filter
-
-    for i_iter in range(num_f_apply):
-        for j_row in range(data.shape[0]):  # apply filters across columns
-            filtered_data[j_row, :] = filtfilt(filt_num, filt_denom, filtered_data[j_row, :])
+        if not np.array_equal(stopband_hz, stopband_hz_adjusted):
+            warnings.warn(
+                "One or both filter frequencies are above Nyquist frequency "
+                f"({nyquist_frequency} Hz), "
+                "so they have been changed "
+                f"({stopband_hz[0]} --> {stopband_hz_adjusted[0]}, "
+                f"{stopband_hz[1]} --> {stopband_hz_adjusted[1]} Hz)."
+            )
+
+        # Convert stopband to a single notch frequency.
+        freq_to_remove = np.mean(stopband_hz_adjusted)
+        bandwidth = np.abs(np.diff(stopband_hz_adjusted))
+
+        # Create filter coefficients.
+        b, a = iirnotch(freq_to_remove, freq_to_remove / bandwidth, fs=sampling_frequency)
+        n_filter_applications = int(np.floor(motion_filter_order / 2))
+
+        filtered_data = data.copy()
+        for _ in range(n_filter_applications):
+            filtered_data = filtfilt(
+                b,
+                a,
+                filtered_data,
+                axis=0,
+                padtype="constant",
+                padlen=data.shape[0] - 1,
+            )
 
     return filtered_data
 
diff --git a/xcp_d/utils/doc.py b/xcp_d/utils/doc.py
index 5926445..7489776 100644
--- a/xcp_d/utils/doc.py
+++ b/xcp_d/utils/doc.py
@@ -165,19 +165,23 @@ docdict[
     "motion_filter_type"
 ] = """
 motion_filter_type : {None, "lp", "notch"}
-    Type of band-stop filter to use for removing respiratory artifact from motion regressors.
+    Type of filter to use for removing respiratory artifact from motion regressors.
+
     If None, no filter will be applied.
 
-    If the filter type is set to "notch", then both ``band_stop_min`` and ``band_stop_max``
-    must be defined.
-    If the filter type is set to "lp", then only ``band_stop_max`` must be defined.
+    If the filter type is set to "notch", frequencies between
+    ``band_stop_min`` and ``band_stop_max`` will be removed with a notch filter.
+    In this case, both ``band_stop_min`` and ``band_stop_max`` must be defined.
+
+    If "lp", frequencies above ``band_stop_min`` will be removed with a Butterworth filter.
+    In this case, only ``band_stop_min`` must be defined.
 """
 
 docdict[
     "motion_filter_order"
 ] = """
 motion_filter_order : :obj:`int`
-    Number of filter coefficients for the band-stop filter, for filtering motion regressors.
+    Number of filter coefficients for the motion parameter filter.
     Motion filtering is only performed if ``motion_filter_type`` is not None.
     This parameter is used in conjunction with ``band_stop_max`` and ``band_stop_min``.
 """
@@ -186,7 +190,7 @@ docdict[
     "band_stop_min"
 ] = """
 band_stop_min : :obj:`float` or None
-    Lower frequency for the band-stop motion filter, in breaths-per-minute (bpm).
+    Lower frequency for the motion parameter filter, in breaths-per-minute (bpm).
     Motion filtering is only performed if ``motion_filter_type`` is not None.
     If used with the "lp" ``motion_filter_type``, this parameter essentially corresponds to a
     low-pass filter (the maximum allowed frequency in the filtered data).
@@ -215,7 +219,7 @@ docdict[
     "band_stop_max"
 ] = """
 band_stop_max : :obj:`float` or None
-    Upper frequency for the band-stop motion filter, in breaths-per-minute (bpm).
+    Upper frequency for the motion parameter filter, in breaths-per-minute (bpm).
     Motion filtering is only performed if ``motion_filter_type`` is not None.
     This parameter is only used if ``motion-filter-type`` is set to "notch".
     This parameter is used in conjunction with ``motion_filter_order`` and ``band_stop_min``.
diff --git a/xcp_d/utils/utils.py b/xcp_d/utils/utils.py
index 6c4b6f8..0d7e1c4 100644
--- a/xcp_d/utils/utils.py
+++ b/xcp_d/utils/utils.py
@@ -7,7 +7,6 @@ import tempfile
 from pathlib import Path
 
 import numpy as np
-from nipype.interfaces.ants import ApplyTransforms
 from scipy.signal import butter, filtfilt
 from templateflow.api import get as get_template
 
@@ -64,6 +63,8 @@ def get_segfile(bold_file):
     -----
     Only used in concatenation code and should be dropped in favor of BIDSLayout methods ASAP.
     """
+    from xcp_d.interfaces.ants import ApplyTransforms
+
     # get transform files
     dd = Path(os.path.dirname(bold_file))
     anatdir = str(dd.parent) + "/anat"
@@ -376,36 +377,40 @@ def butter_bandpass(data, fs, lowpass, highpass, order=2):
 
     Parameters
     ----------
-    data : numpy.ndarray
-        Voxels/vertices by timepoints dimension.
+    data : (T, S) numpy.ndarray
+        Time by voxels/vertices array of data.
     fs : float
         Sampling frequency. 1/TR(s).
     lowpass : float
-        frequency
+        frequency, in Hertz
     highpass : float
-        frequency
+        frequency, in Hertz
     order : int
         The order of the filter. This will be divided by 2 when calling scipy.signal.butter.
 
     Returns
     -------
-    filtered_data : numpy.ndarray
+    filtered_data : (T, S) numpy.ndarray
         The filtered data.
     """
-    nyq = 0.5 * fs  # nyquist frequency
-
-    # normalize the cutoffs
-    lowcut = np.float(highpass) / nyq
-    highcut = np.float(lowpass) / nyq
-
-    b, a = butter(order / 2, [lowcut, highcut], btype="band")  # get filter coeff
+    b, a = butter(
+        order / 2,
+        [highpass, lowpass],
+        btype="bandpass",
+        output="ba",
+        fs=fs,  # eliminates need to normalize cutoff frequencies
+    )
 
-    filtered_data = np.zeros(data.shape)  # create something to populate filtered values with
+    filtered_data = np.zeros_like(data)  # create something to populate filtered values with
 
     # apply the filter, loop through columns of regressors
-    for ii in range(filtered_data.shape[0]):
-        filtered_data[ii, :] = filtfilt(
-            b, a, data[ii, :], padtype="odd", padlen=3 * (max(len(b), len(a)) - 1)
+    for i_voxel in range(filtered_data.shape[1]):
+        filtered_data[:, i_voxel] = filtfilt(
+            b,
+            a,
+            data[:, i_voxel],
+            padtype="constant",
+            padlen=data.shape[0] - 1,
         )
 
     return filtered_data
diff --git a/xcp_d/workflows/anatomical.py b/xcp_d/workflows/anatomical.py
index 7d64f39..4a07e9c 100644
--- a/xcp_d/workflows/anatomical.py
+++ b/xcp_d/workflows/anatomical.py
@@ -10,10 +10,13 @@ from niworkflows.engine.workflows import LiterateWorkflow as Workflow
 from pkg_resources import resource_filename as pkgrf
 from templateflow.api import get as get_template
 
-from xcp_d.interfaces.ants import CompositeInvTransformUtil, ConvertTransformFile
+from xcp_d.interfaces.ants import (
+    ApplyTransforms,
+    CompositeInvTransformUtil,
+    ConvertTransformFile,
+)
 from xcp_d.interfaces.bids import DerivativesDataSink
 from xcp_d.interfaces.c3 import C3d  # TM
-from xcp_d.interfaces.connectivity import ApplyTransformsx
 from xcp_d.interfaces.nilearn import BinaryMath, Merge
 from xcp_d.interfaces.utils import FilterUndefined
 from xcp_d.interfaces.workbench import (  # MB,TM
@@ -134,7 +137,7 @@ def init_warp_anats_to_template_wf(
     else:
         # Warp the native T1w-space T1w and T1w segmentation files to the selected standard space.
         warp_t1w_to_template = pe.Node(
-            ApplyTransformsx(
+            ApplyTransforms(
                 num_threads=2,
                 reference_image=template_file,
                 interpolation="LanczosWindowedSinc",
@@ -156,10 +159,10 @@ def init_warp_anats_to_template_wf(
         # fmt:on
 
         warp_t1seg_to_template = pe.Node(
-            ApplyTransformsx(
+            ApplyTransforms(
                 num_threads=2,
                 reference_image=template_file,
-                interpolation="MultiLabel",
+                interpolation="GenericLabel",
                 input_image_type=3,
                 dimension=3,
             ),
diff --git a/xcp_d/workflows/bold.py b/xcp_d/workflows/bold.py
index 50fe747..ee59ec9 100644
--- a/xcp_d/workflows/bold.py
+++ b/xcp_d/workflows/bold.py
@@ -24,6 +24,7 @@ from xcp_d.interfaces.resting_state import DespikePatch
 from xcp_d.utils.bids import collect_run_data
 from xcp_d.utils.confounds import (
     consolidate_confounds,
+    describe_censoring,
     describe_regression,
     get_customfile,
 )
@@ -207,41 +208,17 @@ def init_boldpostprocess_wf(
         run_data["confounds"],
     )
     regression_description = describe_regression(params, custom_confounds_file)
+    censoring_description = describe_censoring(
+        motion_filter_type=motion_filter_type,
+        motion_filter_order=motion_filter_order,
+        band_stop_min=band_stop_min,
+        band_stop_max=band_stop_max,
+        head_radius=head_radius,
+        fd_thresh=fd_thresh,
+    )
 
     workflow = Workflow(name=name)
 
-    filter_str, filter_post_str = "", ""
-    if motion_filter_type:
-        if motion_filter_type == "notch":
-            filter_sub_str = (
-                f"band-stop filtered to remove signals between {band_stop_min} and "
-                f"{band_stop_max} breaths-per-minute using a notch filter, based on "
-                "@fair2020correction"
-            )
-        else:  # lp
-            filter_sub_str = (
-                f"low-pass filtered below {band_stop_min} breaths-per-minute, "
-                "based on @fair2020correction and @gratton2020removal"
-            )
-
-        filter_str = (
-            f"the six translation and rotation head motion traces were {filter_sub_str}. Next, "
-        )
-        filter_post_str = (
-            "The filtered versions of the motion traces and framewise displacement were not used "
-            "for denoising."
-        )
-
-    if isinstance(head_radius, float):
-        fd_substr = f"with a head radius of {head_radius} mm"
-    else:
-        fd_substr = "with a head radius estimated from the preprocessed brain mask"
-
-    fd_str = (
-        f"{filter_str}framewise displacement was calculated using the formula from "
-        f"@power_fd_dvars, {fd_substr}"
-    )
-
     if dummy_scans == 0 and dummytime != 0:
         dummy_scans = int(np.ceil(dummytime / TR))
 
@@ -267,18 +244,16 @@ def init_boldpostprocess_wf(
     bandpass_str = ""
     if bandpass_filter:
         bandpass_str = (
-            "The interpolated timeseries were then band-pass filtered to retain signals within "
-            f"the {lower_bpf}-{upper_bpf} Hz frequency band."
+            "The interpolated timeseries were then band-pass filtered using a(n) "
+            f"{num2words(bpf_order, ordinal=True)}-order Butterworth filter, "
+            f"in order to retain signals within the {lower_bpf}-{upper_bpf} Hz frequency band."
         )
 
     workflow.__desc__ = f"""\
-For each of the {num2words(n_runs)} BOLD series found per subject (across all tasks and sessions),
+For each of the {num2words(n_runs)} BOLD runs found per subject (across all tasks and sessions),
 the following post-processing was performed.
 First, {dummy_scans_str}outlier detection was performed.
-In order to identify high-motion outlier volumes, {fd_str}.
-Volumes with {'filtered ' if motion_filter_type else ''}framewise displacement greater than
-{fd_thresh} mm were flagged as outliers and excluded from nuisance regression [@power_fd_dvars].
-{filter_post_str}
+{censoring_description}
 {despike_str}
 Next, the BOLD data and confounds were mean-centered and linearly detrended.
 {regression_description}
diff --git a/xcp_d/workflows/cifti.py b/xcp_d/workflows/cifti.py
index 1c471bc..0ecb767 100644
--- a/xcp_d/workflows/cifti.py
+++ b/xcp_d/workflows/cifti.py
@@ -25,6 +25,7 @@ from xcp_d.interfaces.workbench import CiftiConvert
 from xcp_d.utils.bids import collect_run_data
 from xcp_d.utils.confounds import (
     consolidate_confounds,
+    describe_censoring,
     describe_regression,
     get_customfile,
 )
@@ -185,41 +186,17 @@ def init_ciftipostprocess_wf(
         run_data["confounds"],
     )
     regression_description = describe_regression(params, custom_confounds_file)
+    censoring_description = describe_censoring(
+        motion_filter_type=motion_filter_type,
+        motion_filter_order=motion_filter_order,
+        band_stop_min=band_stop_min,
+        band_stop_max=band_stop_max,
+        head_radius=head_radius,
+        fd_thresh=fd_thresh,
+    )
 
     workflow = Workflow(name=name)
 
-    filter_str, filter_post_str = "", ""
-    if motion_filter_type:
-        if motion_filter_type == "notch":
-            filter_sub_str = (
-                f"band-stop filtered to remove signals between {band_stop_min} and "
-                f"{band_stop_max} breaths-per-minute using a notch filter, based on "
-                "@fair2020correction"
-            )
-        else:  # lp
-            filter_sub_str = (
-                f"low-pass filtered below {band_stop_min} breaths-per-minute, "
-                "based on @fair2020correction and @gratton2020removal"
-            )
-
-        filter_str = (
-            f"the six translation and rotation head motion traces were {filter_sub_str}. Next, "
-        )
-        filter_post_str = (
-            "The filtered versions of the motion traces and framewise displacement were not used "
-            "for denoising."
-        )
-
-    if isinstance(head_radius, float):
-        fd_substr = f"with a head radius of {head_radius} mm"
-    else:
-        fd_substr = "with a head radius estimated from the preprocessed brain mask"
-
-    fd_str = (
-        f"{filter_str}framewise displacement was calculated using the formula from "
-        f"@power_fd_dvars, {fd_substr}"
-    )
-
     if dummy_scans == 0 and dummytime != 0:
         dummy_scans = int(np.ceil(dummytime / TR))
 
@@ -246,18 +223,16 @@ def init_ciftipostprocess_wf(
     bandpass_str = ""
     if bandpass_filter:
         bandpass_str = (
-            "The interpolated timeseries were then band-pass filtered to retain signals within "
-            f"the {lower_bpf}-{upper_bpf} Hz frequency band."
+            "The interpolated timeseries were then band-pass filtered using a(n) "
+            f"{num2words(bpf_order, ordinal=True)}-order Butterworth filter, "
+            f"in order to retain signals within the {lower_bpf}-{upper_bpf} Hz frequency band."
         )
 
     workflow.__desc__ = f"""\
-For each of the {num2words(n_runs)} BOLD series found per subject (across all tasks and sessions),
+For each of the {num2words(n_runs)} BOLD runs found per subject (across all tasks and sessions),
 the following post-processing was performed.
 First, {dummy_scans_str}outlier detection was performed.
-In order to identify high-motion outlier volumes, {fd_str}.
-Volumes with {'filtered ' if motion_filter_type else ''}framewise displacement greater than
-{fd_thresh} mm were flagged as outliers and excluded from nuisance regression [@power_fd_dvars].
-{filter_post_str}
+{censoring_description}
 {despike_str}
 Next, the BOLD data and confounds were mean-centered and linearly detrended.
 {regression_description}
diff --git a/xcp_d/workflows/connectivity.py b/xcp_d/workflows/connectivity.py
index 610a4d5..6128372 100644
--- a/xcp_d/workflows/connectivity.py
+++ b/xcp_d/workflows/connectivity.py
@@ -8,8 +8,9 @@ from nipype.interfaces import utility as niu
 from nipype.pipeline import engine as pe
 from niworkflows.engine.workflows import LiterateWorkflow as Workflow
 
+from xcp_d.interfaces.ants import ApplyTransforms
 from xcp_d.interfaces.bids import DerivativesDataSink
-from xcp_d.interfaces.connectivity import ApplyTransformsx, ConnectPlot, NiftiConnect
+from xcp_d.interfaces.connectivity import ConnectPlot, NiftiConnect
 from xcp_d.interfaces.prepostcleaning import CiftiPrepareForParcellation
 from xcp_d.interfaces.workbench import (
     CiftiCorrelation,
@@ -150,8 +151,8 @@ or were set to zero,  when the parcel had <{min_coverage * 100}% coverage.
 
     # Using the generated transforms, apply them to get everything in the correct MNI form
     warp_atlases_to_bold_space = pe.MapNode(
-        ApplyTransformsx(
-            interpolation="MultiLabel",
+        ApplyTransforms(
+            interpolation="GenericLabel",
             input_image_type=3,
             dimension=3,
         ),
diff --git a/xcp_d/workflows/plotting.py b/xcp_d/workflows/plotting.py
index 2167d25..ab3f0bc 100644
--- a/xcp_d/workflows/plotting.py
+++ b/xcp_d/workflows/plotting.py
@@ -3,9 +3,9 @@ from nipype import Function
 from nipype.interfaces import utility as niu
 from nipype.pipeline import engine as pe
 from niworkflows.engine.workflows import LiterateWorkflow as Workflow
-from niworkflows.interfaces.fixes import FixHeaderApplyTransforms as ApplyTransforms
 from templateflow.api import get as get_template
 
+from xcp_d.interfaces.ants import ApplyTransforms
 from xcp_d.interfaces.bids import DerivativesDataSink
 from xcp_d.interfaces.qc_plot import CensoringPlot, QCPlot
 from xcp_d.interfaces.report import FunctionalSummary
@@ -302,7 +302,7 @@ def init_qc_report_wf(
             ApplyTransforms(
                 dimension=3,
                 input_image=dseg_file,
-                interpolation="MultiLabel",
+                interpolation="GenericLabel",
             ),
             name="warp_dseg_to_bold",
             n_procs=omp_nthreads,


Test patch:
diff --git a/xcp_d/tests/data/band_passed_bold_data.txt b/xcp_d/tests/data/band_passed_bold_data.txt
deleted file mode 100644
index 533f4a7..0000000
--- a/xcp_d/tests/data/band_passed_bold_data.txt
+++ /dev/null
@@ -1,256 +0,0 @@
--0.6334
--0.38641
--0.15557
-0.066915
-0.27613
-0.36001
-0.30095
-0.24997
-0.29013
-0.39268
-0.4349
-0.42727
-0.45328
-0.48037
-0.55162
-0.66022
-0.7477
-0.78652
-0.65389
-0.32299
--0.0081086
--0.19708
--0.28015
--0.28672
--0.19374
--0.0102
-0.14556
-0.20757
-0.28976
-0.41794
-0.49182
-0.49186
-0.44776
-0.34889
-0.22028
-0.11174
-0.024844
--0.0069693
-0.02643
-0.064267
-0.10522
-0.17323
-0.22767
-0.23069
-0.20557
-0.18
-0.15969
-0.18245
-0.26063
-0.30304
-0.28405
-0.22798
-0.099038
--0.055244
--0.15352
--0.16357
--0.067023
-0.10787
-0.25284
-0.28939
-0.2392
-0.09392
--0.076285
--0.12554
--0.10495
--0.13935
--0.20715
--0.2023
--0.13782
--0.10878
--0.13515
--0.18394
--0.15873
--0.042173
-0.023365
--0.080776
--0.30295
--0.47377
--0.51238
--0.46839
--0.41014
--0.40647
--0.47426
--0.48879
--0.40409
--0.3659
--0.43458
--0.49934
--0.50471
--0.45953
--0.35118
--0.25617
--0.25328
--0.25517
--0.14278
-0.067989
-0.26639
-0.35358
-0.36182
-0.34325
-0.31764
-0.37293
-0.48897
-0.53334
-0.47211
-0.35721
-0.25656
-0.19676
-0.17213
-0.1946
-0.23514
-0.24165
-0.21981
-0.20836
-0.1822
-0.14941
-0.17977
-0.2136
-0.17154
-0.11484
-0.090216
-0.082234
-0.089067
-0.1067
-0.14717
-0.20942
-0.23246
-0.17179
-0.011544
--0.19072
--0.30135
--0.33732
--0.36584
--0.32632
--0.21326
--0.12831
--0.11591
--0.15661
--0.20812
--0.20638
--0.14682
--0.042837
-0.10744
-0.20131
-0.17274
-0.10418
-0.065378
-0.047884
-0.024518
-0.028811
-0.087365
-0.18379
-0.29592
-0.43195
-0.52467
-0.50354
-0.48639
-0.46012
-0.34686
-0.24316
-0.20706
-0.21535
-0.20419
-0.13658
-0.037589
--0.069842
--0.24287
--0.47389
--0.59751
--0.61302
--0.64637
--0.71063
--0.76042
--0.75218
--0.70155
--0.64627
--0.57114
--0.46298
--0.31823
--0.17871
--0.087656
--0.049633
--0.09647
--0.14223
--0.071972
-0.050308
-0.10115
-0.05645
--0.035015
--0.17092
--0.32035
--0.39689
--0.35322
--0.1756
-0.056878
-0.23567
-0.2955
-0.21726
-0.16284
-0.22667
-0.24946
-0.19545
-0.17268
-0.16991
-0.13669
-0.068193
--0.021715
--0.099321
--0.1532
--0.16906
--0.16703
--0.21748
--0.3261
--0.44842
--0.50672
--0.47157
--0.39124
--0.32701
--0.34282
--0.44496
--0.55013
--0.54642
--0.4356
--0.28042
--0.12043
--0.030544
-0.0070595
-0.077645
-0.1403
-0.15628
-0.14108
-0.063542
--0.058478
--0.18029
--0.30387
--0.39413
--0.41883
--0.39786
--0.31673
--0.19816
--0.10516
--0.077077
--0.12238
--0.15151
--0.15386
--0.26185
--0.51181
--0.7626
--0.93277
--1.0334
--1.0153
--0.88919
--0.75013
--0.67746
--0.64635
--0.57929
diff --git a/xcp_d/tests/data/low_passed_MR_data.txt b/xcp_d/tests/data/low_passed_MR_data.txt
deleted file mode 100644
index 360d1e9..0000000
--- a/xcp_d/tests/data/low_passed_MR_data.txt
+++ /dev/null
@@ -1,256 +0,0 @@
--0.83047
--0.38746
--0.224
--0.30009
-0.86581
-1.9384
--0.76131
--0.17471
--0.98129
-1.0126
-1.0697
--0.95268
-1.243
--0.62586
-0.18668
-0.90501
-0.48036
-2.1073
-2.4064
-0.20352
--1.5124
--0.57612
--0.79365
--0.807
--0.95154
-0.15834
-1.5803
--0.36625
--0.16126
-1.0112
-1.1261
-0.62167
-0.94908
-0.75325
--0.16108
-0.2651
--0.28285
--0.66882
-0.40909
-0.097503
--0.10353
-0.31397
-0.846
-0.43358
-0.22645
-0.24371
--0.0080523
--0.51385
-0.6815
-1.1352
--0.06822
-1.4415
-0.032633
--0.29929
--0.81599
--0.81416
--0.79059
-0.34383
-1.344
-0.67553
-0.86185
-0.88162
--1.8331
--0.34832
-0.32087
-0.4589
--1.0493
--0.84032
-0.50015
-0.042021
-0.31824
--0.86386
--1.2958
-0.55323
-1.1745
-1.2423
--1.1716
--1.3566
--1.3072
--0.25093
--0.36003
-0.80798
--1.1871
--1.4879
-0.09569
-0.86246
--0.63691
--1.0934
--0.42291
--1.2985
--0.070543
-0.66167
-0.056041
--1.5286
--0.64473
--0.10699
-1.484
-0.90432
-0.19021
-0.91364
--0.71409
--0.22481
-1.227
-1.4023
-0.78375
-0.40776
--0.33091
-0.44014
--0.51266
-0.40006
-0.39557
-0.88257
--0.14476
-0.53802
-0.59362
--0.59746
-0.19868
-1.2763
-0.22489
--0.17793
-0.29907
--0.12331
-0.30622
--0.060659
-0.026008
-0.50278
-1.1693
-0.66323
-1.1311
--1.4629
--0.53659
-0.039327
--0.76289
--1.0203
-0.41266
-0.572
-0.33129
-0.24969
--0.69268
--0.49333
--0.20933
--0.50928
-0.43334
-1.8388
--0.06599
-0.3623
--0.60872
-0.84416
--0.60405
-0.0043525
--0.41315
-0.33865
-0.065954
-0.27425
-2.242
--0.18171
-0.18395
-1.9789
--0.56516
-0.21449
--0.68158
-0.73417
-0.4278
-0.52091
--0.11864
-0.31122
-0.93595
--1.7984
--1.472
-0.077823
--0.51196
--0.45879
--1.3029
--0.91162
--0.4955
--0.70647
--0.80401
--0.51755
--0.6496
-0.85545
--0.44971
-1.2828
--0.261
--1.2324
--0.53675
-0.7739
-0.94114
-0.27142
-0.30677
-0.37245
--1.1441
--0.79897
--1.5445
--0.44367
-0.77398
-0.89962
-1.882
-0.15289
--1.4846
-1.232
-1.1143
--0.29104
-0.20096
-0.67501
-0.45616
-0.57704
--0.018059
--0.033342
--0.012847
--0.46151
-1.0367
-0.21948
-0.45089
--0.71127
--0.94038
--0.20882
--0.084579
-0.97537
-0.65049
-0.33703
--1.1748
--0.88791
-0.020509
--0.11565
-1.114
-1.1276
--0.21456
-0.77712
-1.2228
-0.57374
-1.4231
-1.0292
-0.37677
-0.73246
-0.20268
--0.26439
-0.24094
-0.032516
--0.095158
-1.3834
-0.64126
-2.1413
-0.11815
-0.59598
-1.5849
-2.3362
-0.44889
--0.49912
-0.4097
--0.98149
--0.31613
-0.34054
-1.3454
-1.2838
-0.23675
-0.92393
diff --git a/xcp_d/tests/data/notched_MR_data.txt b/xcp_d/tests/data/notched_MR_data.txt
deleted file mode 100644
index 6c61b1d..0000000
--- a/xcp_d/tests/data/notched_MR_data.txt
+++ /dev/null
@@ -1,256 +0,0 @@
--0.82855
--0.51965
-0.12628
--0.019105
-0.35003
-1.6596
--1.1379
-1.0555
--1.1408
-0.59591
-1.0227
--1.2791
-2.0665
--1.2337
-0.20457
-1.4573
-0.52997
-1.6836
-1.814
-0.44652
--0.99018
--0.40694
--1.4011
--0.56341
--0.52107
--0.21425
-1.3578
--0.47522
-0.37204
-1.057
-0.81325
-0.45961
-1.1154
-0.87722
--0.42011
-0.37557
--0.1245
--0.73838
-0.41589
--0.058781
-0.043053
-0.39158
-0.81263
-0.28095
-0.1685
-0.25367
-0.24508
--0.54508
-0.60465
-1.2542
--0.42612
-1.8231
--0.33319
--0.21941
--0.80018
--0.59351
--0.74391
-0.32745
-1.4552
-0.33806
-0.36204
-1.326
--1.5585
-0.21188
--0.54531
-0.049503
--0.66228
--0.14845
-0.61177
--0.99802
-0.062624
--0.29438
--0.57948
-0.7244
-0.089125
-0.89917
--1.0729
--0.45596
--1.1129
--0.32477
--1.35
-0.84309
--0.78555
--0.82172
--0.11544
-0.28544
--0.93226
--0.86725
-0.23258
--1.2896
--0.21111
-0.013479
--0.0088532
--1.1739
-0.10916
--0.50266
-1.2226
-0.74243
--0.024794
-1.2967
--0.61433
-0.057278
-1.0635
-0.95255
-0.69041
-0.8195
--0.44292
-0.65973
--0.78617
-0.61725
-0.16401
-1.1279
--0.35746
-0.42649
-0.75632
--0.41171
-0.22247
-1.0361
--0.030996
-0.082278
-0.66712
--0.46702
-0.15701
--0.054265
-0.39173
-0.61979
-1.0078
--0.097806
-1.67
--1.3412
--0.34106
--0.34146
--0.8083
--0.73315
-0.65152
-0.20914
-0.06542
-0.4906
--0.71564
--0.52491
-0.021399
--0.31603
-0.133
-1.7293
--0.50294
-1.1628
--0.86698
-1.0232
--1.1235
-0.33796
--0.51055
-0.61764
-0.14708
--0.33306
-2.5903
--0.16539
-0.18192
-2.0789
--1.3093
-1.091
--0.75885
-0.6551
--0.042871
-0.98866
-0.11873
--0.29366
-0.84897
--1.6351
--0.85676
-0.047293
--1.3511
--0.29551
--0.96427
--0.69486
--0.70747
--0.96524
--0.79838
--0.19115
--0.78514
-1.3738
--1.2187
-1.2448
--0.36102
--0.62744
--0.28313
-0.35159
-0.62226
-0.42372
-0.49233
-0.41516
--1.6844
--0.49192
--1.3496
--0.030688
-0.79739
--0.055196
-1.664
-0.86013
--1.1448
-1.2171
-0.40193
--0.34587
-0.80097
-0.7789
-0.065959
-0.53156
--0.064863
-0.04987
-0.2906
--0.61417
-1.2413
--0.33668
-0.53057
--0.65318
--0.74705
-0.09265
--0.2551
-0.89163
-0.16509
-0.47409
--0.99835
--0.59121
-0.22994
--0.54473
-0.87528
-1.2591
--0.14985
-1.0201
-1.2046
-0.23595
-1.5459
-1.0545
-0.29846
-0.86314
-0.20589
--0.42189
-0.24596
-0.20491
--0.016665
-1.7294
--0.32411
-2.3656
-0.30242
-1.2475
-1.12
-1.9273
-0.43634
--0.36407
-0.75065
--1.211
-0.079034
-0.14667
-0.93305
-1.407
-0.55762
-0.9315
diff --git a/xcp_d/tests/data/raw_data.txt b/xcp_d/tests/data/raw_data.txt
deleted file mode 100644
index 9e1b35b..0000000
--- a/xcp_d/tests/data/raw_data.txt
+++ /dev/null
@@ -1,256 +0,0 @@
--0.830468498
--0.352252186
--0.17477504
--0.480653419
-0.836836711
-2.538349341
--1.323334251
-0.128340263
--1.442379153
-1.30250821
-1.4099115
--1.66254297
-1.943684478
--1.084698452
-0.226818975
-1.098929202
-0.147188752
-2.295665813
-2.752557878
-0.138317725
--1.907066208
--0.364992993
--0.848110012
--0.764753365
--1.127694815
-0.078188895
-2.106629656
--0.715847393
--0.280515662
-1.166475005
-1.212821419
-0.485540987
-1.026016486
-0.870726021
--0.381757879
-0.428893027
--0.29913052
--0.89986852
-0.634745461
-0.067453591
--0.187120545
-0.291727455
-0.987694692
-0.39293457
-0.194551374
-0.279784956
-0.051220312
--0.774466244
-0.78678171
-1.408906951
--0.534098584
-1.92775843
--0.176247547
--0.243750362
--0.897600659
--0.792336865
--0.952974691
-0.353905454
-1.597026324
-0.527470251
-0.854202301
-1.341846522
--2.499533442
--0.167559322
-0.353015315
-0.717253725
--1.304851639
--1.005868952
-0.79068347
--0.116571329
-0.553089888
--0.960644805
--1.633802396
-0.761200282
-1.193307152
-1.632057216
--1.532189603
--1.336852405
--1.473846498
--0.041663074
--0.615507335
-1.314154809
--1.455066582
--1.742349221
-0.205304683
-1.192930399
--0.802823098
--1.265636422
--0.149331353
--1.636446678
-0.017344346
-0.828387266
-0.217738348
--1.909244901
--0.536821821
--0.302032298
-1.813582131
-0.914851753
--0.057080716
-1.309362078
--1.044735848
--0.34826681
-1.412561161
-1.502382929
-0.730375995
-0.490752272
--0.58612614
-0.74489965
--0.82815497
-0.574520732
-0.281841383
-1.139306269
--0.425867846
-0.636139891
-0.793178081
--0.898377114
-0.156244844
-1.597253911
-0.112439711
--0.308624935
-0.456659609
--0.275100831
-0.443143609
--0.134765127
--0.01832823
-0.460789425
-1.362315464
-0.451874567
-1.648383657
--2.028361961
--0.449256811
-0.235993202
--0.835172964
--1.275955163
-0.617035081
-0.612701504
-0.289381155
-0.395316122
--0.870562842
--0.497688378
--0.106671817
--0.687829135
-0.331880877
-2.365224743
--0.482230718
-0.647448207
--1.034424717
-1.339554598
--0.969140356
-0.208715601
--0.618593355
-0.512015644
-0.011354171
--0.043988626
-2.949092541
--0.630046194
--0.0468794
-2.683025511
--1.146690686
-0.552998745
--1.076458409
-1.030639584
-0.327529812
-0.652124811
--0.278861116
-0.245191594
-1.472513479
--2.275101542
--1.633290725
-0.415468937
--0.654768852
--0.296348315
--1.496918876
--0.904834447
--0.404181545
--0.725798166
--0.866485025
--0.42184678
--0.942666324
-1.341883717
--0.988434712
-1.817942688
--0.374436607
--1.45174072
--0.618681808
-0.934500962
-1.055929256
-0.160227277
-0.28740032
-0.632905577
--1.459041869
--0.581709725
--1.830149322
--0.44910308
-0.949274704
-0.717441409
-2.287828606
-0.166727663
--2.156490752
-1.689399307
-1.282280747
--0.582630964
-0.222614481
-0.779451473
-0.384707158
-0.696366923
--0.112715559
--0.038823701
-0.088088629
--0.789655973
-1.422960619
-0.006331684
-0.686480848
--0.854933634
--1.075235317
--0.090967251
--0.252771833
-1.194823819
-0.606063709
-0.540514429
--1.444939419
--0.967693808
-0.202051338
--0.347877973
-1.29008799
-1.341153819
--0.580797689
-0.875135602
-1.395449887
-0.320985222
-1.623382458
-1.062433309
-0.214105238
-0.876803202
-0.194407011
--0.41489216
-0.358458734
-0.036222668
--0.364631168
-1.771019567
-0.221273068
-2.730377585
--0.296165383
-0.564295552
-1.582620663
-2.729230123
-0.303564281
--0.790258479
-0.803380233
--1.319902702
--0.27384576
-0.271866922
-1.489553508
-1.437141525
--0.027561397
-0.923931315
diff --git a/xcp_d/tests/test_filtering.py b/xcp_d/tests/test_filtering.py
index 02c3728..ebb7c6a 100644
--- a/xcp_d/tests/test_filtering.py
+++ b/xcp_d/tests/test_filtering.py
@@ -1,95 +1,166 @@
 """Tests for filtering methods."""
-import os.path as op
+import re
 
 import numpy as np
 import pytest
+from scipy import signal
 
 from xcp_d.interfaces.filtering import butter_bandpass
 from xcp_d.utils.confounds import motion_regression_filter
 
 
-@pytest.fixture(scope="module")
-def data_files():
-    """Grab files for testing filters.
+def test_motion_filtering_lp():
+    """Run lowpass filter on toy data, compare to simplified results."""
+    raw_data = np.random.random(500)
 
-    These files were kindly provided by the DCAN lab to directly compare their
-    MATLAB filtering results to our Python filtering results.
+    band_stop_min = 6
+    TR = 0.8
 
-    We have converted the single-column CSV files provided to TXT files for easier I/O.
-    """
-    data_path = op.abspath(op.join(op.dirname(__file__), "data"))
-    data_files_ = {
-        "raw_data": op.join(data_path, "raw_data.txt"),
-        "lowpass_filtered": op.join(data_path, "low_passed_MR_data.txt"),
-        "notch_filtered": op.join(data_path, "notched_MR_data.txt"),
-        "butterworth_filtered": op.join(data_path, "band_passed_bold_data.txt"),
-    }
-    return data_files_
-
-
-def test_motion_filtering_lp(data_files):
-    """Run lowpass filter on toy data, compare to results that have been verified."""
-    raw_data = np.loadtxt(data_files["raw_data"])
-    lowpass_data_true = np.loadtxt(data_files["lowpass_filtered"])
-    raw_data = raw_data[None, :]  # add singleton row dimension
-
-    band_stop_min = 12
-    band_stop_max = 20
+    lowpass = band_stop_min / 60
+    b, a = signal.butter(
+        1,
+        lowpass,
+        btype="lowpass",
+        output="ba",
+        fs=1 / TR,
+    )
+    lowpass_data_true = signal.filtfilt(
+        b,
+        a,
+        raw_data,
+        padtype="constant",
+        padlen=raw_data.size - 1,
+    )
 
     # Confirm the LP filter runs with reasonable parameters
+    raw_data = raw_data[:, None]  # add singleton row dimension
     with pytest.warns(match="The parameter 'band_stop_max' will be ignored."):
         lowpass_data_test = motion_regression_filter(
             raw_data,
-            TR=0.8,
+            TR=TR,
             motion_filter_type="lp",
             band_stop_min=band_stop_min,
-            band_stop_max=band_stop_max,
+            band_stop_max=18,
             motion_filter_order=2,
         )
 
     # What's the difference from the verified data?
-    assert np.allclose(lowpass_data_test, lowpass_data_true, atol=1e-4)
+    assert np.allclose(np.squeeze(lowpass_data_test), lowpass_data_true)
+
+    # Use freq above Nyquist, and function will automatically modify the filter.
+    with pytest.warns(
+        UserWarning,
+        match=re.escape(
+            "Low-pass filter frequency is above Nyquist frequency (0.625 Hz), "
+            "so it has been changed (0.7 --> 0.55 Hz)."
+        ),
+    ):
+        motion_regression_filter(
+            raw_data,
+            TR=TR,  # 1.25 Hz
+            motion_filter_type="lp",
+            band_stop_min=42,  # 0.7 Hz > (1.25 / 2)
+            band_stop_max=None,
+            motion_filter_order=2,
+        )
+
 
+def test_motion_filtering_notch():
+    """Run notch filter on toy data, compare to simplified results."""
+    raw_data = np.random.random(500)
 
-def test_motion_filtering_notch(data_files):
-    """Run notch filter on toy data, compare to results that have been verified.
+    band_stop_min, band_stop_max = 12, 20
+    TR = 0.8
 
-    Notes
-    -----
-    This test requires a much more liberal tolerance value, because the notch filter
-    doesn't replicate well across languages.
-    """
-    raw_data = np.loadtxt(data_files["raw_data"])
-    notch_data_true = np.loadtxt(data_files["notch_filtered"])
-    raw_data = raw_data[None, :]  # add singleton row dimension
+    lowcut, highcut = band_stop_min / 60, band_stop_max / 60
+    stopband_hz_adjusted = [lowcut, highcut]
+    freq_to_remove = np.mean(stopband_hz_adjusted)
+    bandwidth = np.abs(np.diff(stopband_hz_adjusted))
 
-    band_stop_min = 12
-    band_stop_max = 20
+    # Create filter coefficients.
+    b, a = signal.iirnotch(freq_to_remove, freq_to_remove / bandwidth, fs=1 / TR)
+    notch_data_true = signal.filtfilt(
+        b,
+        a,
+        raw_data,
+        padtype="constant",
+        padlen=raw_data.size - 1,
+    )
 
     # Repeat for notch filter
+    raw_data = raw_data[:, None]  # add singleton row dimension
     notch_data_test = motion_regression_filter(
         raw_data,
-        TR=0.8,
+        TR=TR,  # 1.25 Hz
         motion_filter_type="notch",
         band_stop_min=band_stop_min,
         band_stop_max=band_stop_max,
         motion_filter_order=2,
     )
-    assert np.allclose(notch_data_test, notch_data_true, atol=1e-1)
+    notch_data_test = np.squeeze(notch_data_test)
+    assert np.allclose(notch_data_test, notch_data_true)
+
+    # Use band above Nyquist, and function will automatically modify the filter.
+    # NOTE: In this case, the min and max end up flipped for some reason.
+    with pytest.warns(
+        UserWarning,
+        match=re.escape(
+            "One or both filter frequencies are above Nyquist frequency (0.625 Hz), "
+            "so they have been changed (0.7 --> 0.55, 0.75 --> 0.5 Hz)."
+        ),
+    ):
+        motion_regression_filter(
+            raw_data,
+            TR=TR,
+            motion_filter_type="notch",
+            band_stop_min=42,
+            band_stop_max=45,  # 0.7 Hz > (1.25 / 2)
+            motion_filter_order=2,
+        )
+
+    # Using a filter type other than notch or lp should raise an exception.
+    with pytest.raises(ValueError, match="Motion filter type 'fail' not supported."):
+        motion_regression_filter(
+            raw_data,
+            TR=TR,
+            motion_filter_type="fail",
+            band_stop_min=band_stop_min,
+            band_stop_max=band_stop_max,
+            motion_filter_order=2,
+        )
 
 
-def test_bandpass_filtering(data_files):
+def test_bandpass_filtering():
     """Run Butterworth on toy data, compare to results that have been verified."""
-    raw_data = np.loadtxt(data_files["raw_data"])
-    butterworth_data_true = np.loadtxt(data_files["butterworth_filtered"])
-    raw_data = raw_data[None, :]  # add singleton row dimension
+    raw_data = np.random.random(500)
+
+    highpass, lowpass = 0.009, 0.08
+    TR = 0.8
+    fs = 1 / TR
+
+    b, a = signal.butter(
+        1,
+        [highpass, lowpass],
+        btype="bandpass",
+        output="ba",
+        fs=fs,
+    )
+    butterworth_data_true = signal.filtfilt(
+        b,
+        a,
+        raw_data,
+        padtype="constant",
+        padlen=raw_data.size - 1,
+    )
 
     # Confirm the butterworth filter runs with reasonable parameters
+    raw_data = raw_data[:, None]  # add singleton row dimension
     butterworth_data_test = butter_bandpass(
         raw_data,
-        fs=1 / 0.8,
-        highpass=0.009,
-        lowpass=0.080,
+        fs=fs,
+        highpass=highpass,
+        lowpass=lowpass,
         order=2,
     )
-    assert np.allclose(butterworth_data_test, butterworth_data_true, atol=1e-4)
+    butterworth_data_test = np.squeeze(butterworth_data_test)
+    assert np.allclose(butterworth_data_test, butterworth_data_true)
diff --git a/xcp_d/tests/test_workflows_connectivity.py b/xcp_d/tests/test_workflows_connectivity.py
index f6fce7a..5df2261 100644
--- a/xcp_d/tests/test_workflows_connectivity.py
+++ b/xcp_d/tests/test_workflows_connectivity.py
@@ -99,8 +99,8 @@ def test_nifti_conn(fmriprep_with_freesurfer_data, tmp_path_factory):
     assert ground_truth.shape == (400, 400)
 
     # Parcels with <50% coverage should have NaNs
-    # We know that 10 of the parcels in the 400-parcel Schaefer are flagged
-    assert np.sum(np.isnan(np.diag(xcp_array))) == 10
+    # We know that 14 of the parcels in the 400-parcel Schaefer are flagged
+    assert np.sum(np.isnan(np.diag(xcp_array))) == 14
 
     # If we replace the bad parcels' results in the "ground truth" matrix with NaNs,
     # the resulting matrix should match the workflow-generated one.