Repository: PennLINC/xcp_d

Issue:
Modified motion filter parameters are not reflected in the boilerplate or metadata
## Summary
<!--What is the nature of the bug?-->
The function `motion_regression_filter` will automatically adjust the motion filter parameters `band_stop_min` and `band_stop_max` if either is above the Nyquist frequency for the data (see below). However, the adjustments are not reflected anywhere except in a warning message. We should make any adjustments before the boilerplate is generated, as well as update the metadata for any generated files.

https://github.com/PennLINC/xcp_d/blob/4121f083b2debb722587d5e69e154323ab995727/xcp_d/utils/confounds.py#L599-L610

https://github.com/PennLINC/xcp_d/blob/4121f083b2debb722587d5e69e154323ab995727/xcp_d/utils/confounds.py#L632-L646


Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/PennLINC/xcp_d/pull/1114?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennLINC) Report
Attention: Patch coverage is `91.76471%` with `7 lines` in your changes are missing coverage. Please review.
> Project coverage is 83.02%. Comparing base [(`919ccc1`)](https://app.codecov.io/gh/PennLINC/xcp_d/commit/919ccc13b67bdff437265345dfa22dbc41b986b2?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennLINC) to head [(`4ab8527`)](https://app.codecov.io/gh/PennLINC/xcp_d/pull/1114?dropdown=coverage&src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennLINC).

| [Files](https://app.codecov.io/gh/PennLINC/xcp_d/pull/1114?dropdown=coverage&src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennLINC) | Patch % | Lines |
|---|---|---|
| [xcp\_d/utils/confounds.py](https://app.codecov.io/gh/PennLINC/xcp_d/pull/1114?src=pr&el=tree&filepath=xcp_d%2Futils%2Fconfounds.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennLINC#diff-eGNwX2QvdXRpbHMvY29uZm91bmRzLnB5) | 85.71% | [3 Missing and 4 partials :warning: ](https://app.codecov.io/gh/PennLINC/xcp_d/pull/1114?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennLINC) |

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


```diff
@@            Coverage Diff             @@
##             main    #1114      +/-   ##
==========================================
+ Coverage   82.98%   83.02%   +0.04%     
==========================================
  Files          52       52              
  Lines        6610     6656      +46     
  Branches      905      914       +9     
==========================================
+ Hits         5485     5526      +41     
- Misses        912      914       +2     
- Partials      213      216       +3     
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/PennLINC/xcp_d/pull/1114?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennLINC).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennLINC).


Patch:
diff --git a/xcp_d/interfaces/censoring.py b/xcp_d/interfaces/censoring.py
index 4532f78..05c4fb0 100644
--- a/xcp_d/interfaces/censoring.py
+++ b/xcp_d/interfaces/censoring.py
@@ -14,7 +14,12 @@ from nipype.interfaces.base import (
     traits,
 )
 
-from xcp_d.utils.confounds import _infer_dummy_scans, load_confound_matrix, load_motion
+from xcp_d.utils.confounds import (
+    _infer_dummy_scans,
+    _modify_motion_filter,
+    load_confound_matrix,
+    load_motion,
+)
 from xcp_d.utils.filemanip import fname_presuffix
 from xcp_d.utils.modified_data import _drop_dummy_scans, compute_fd
 
@@ -452,13 +457,20 @@ class GenerateConfounds(SimpleInterface):
 
     def _run_interface(self, runtime):
         fmriprep_confounds_df = pd.read_table(self.inputs.fmriprep_confounds_file)
+        band_stop_min_adjusted, band_stop_max_adjusted, _ = _modify_motion_filter(
+            motion_filter_type=self.inputs.motion_filter_type,
+            band_stop_min=self.inputs.band_stop_min,
+            band_stop_max=self.inputs.band_stop_max,
+            TR=self.inputs.TR,
+        )
+
         motion_df = load_motion(
             fmriprep_confounds_df.copy(),
             TR=self.inputs.TR,
             motion_filter_type=self.inputs.motion_filter_type,
             motion_filter_order=self.inputs.motion_filter_order,
-            band_stop_min=self.inputs.band_stop_min,
-            band_stop_max=self.inputs.band_stop_max,
+            band_stop_min=band_stop_min_adjusted,
+            band_stop_max=band_stop_max_adjusted,
         )
 
         # Add in framewise displacement
@@ -504,7 +516,52 @@ class GenerateConfounds(SimpleInterface):
             custom_confounds=self.inputs.custom_confounds_file,
         )
 
+        # Compile motion metadata from confounds metadata, adding in filtering info
+        motion_metadata = {}
+        for col in motion_df.columns.tolist():
+            col_metadata = confounds_metadata.get(col, {})
+            if col == "framewise_displacement":
+                col_metadata["Description"] = (
+                    "Framewise displacement calculated according to Power et al. (2012)."
+                )
+                col_metadata["Units"] = "mm"
+                col_metadata["HeadRadius"] = self.inputs.head_radius
+
+            if self.inputs.motion_filter_type == "lp":
+                filters = col_metadata.get("SoftwareFilters", {})
+                filters["Butterworth low-pass filter"] = {
+                    "cutoff": band_stop_min_adjusted / 60,
+                    "order": self.inputs.motion_filter_order,
+                    "cutoff units": "Hz",
+                    "function": "scipy.signal.filtfilt",
+                }
+                col_metadata["SoftwareFilters"] = filters
+
+            elif self.inputs.motion_filter_type == "notch":
+                filters = col_metadata.get("SoftwareFilters", {})
+                filters["IIR notch digital filter"] = {
+                    "cutoff": [
+                        band_stop_max_adjusted / 60,
+                        band_stop_min_adjusted / 60,
+                    ],
+                    "order": self.inputs.motion_filter_order,
+                    "cutoff units": "Hz",
+                    "function": "scipy.signal.filtfilt",
+                }
+                col_metadata["SoftwareFilters"] = filters
+
+            motion_metadata[col] = col_metadata
+
+        self._results["motion_metadata"] = motion_metadata
+
         if confounds_df is not None:
+            # Update confounds metadata with modified motion metadata
+            for col in motion_df.columns.tolist():
+                if col in confounds_df.columns:
+                    base_metadata = confounds_metadata.get(col, {})
+                    base_metadata.update(motion_metadata[col])
+                    confounds_metadata[col] = base_metadata
+
             # Orthogonalize full nuisance regressors w.r.t. any signal regressors
             signal_columns = [c for c in confounds_df.columns if c.startswith("signal__")]
             if signal_columns:
@@ -590,28 +647,6 @@ class GenerateConfounds(SimpleInterface):
             sep="\t",
         )
 
-        # Compile metadata to pass along to outputs.
-        motion_metadata = {
-            "framewise_displacement": {
-                "Description": (
-                    "Framewise displacement calculated according to Power et al. (2012)."
-                ),
-                "Units": "mm",
-                "HeadRadius": self.inputs.head_radius,
-            }
-        }
-        if self.inputs.motion_filter_type == "lp":
-            motion_metadata["LowpassFilter"] = self.inputs.band_stop_max
-            motion_metadata["LowpassFilterOrder"] = self.inputs.motion_filter_order
-        elif self.inputs.motion_filter_type == "notch":
-            motion_metadata["BandstopFilter"] = [
-                self.inputs.band_stop_min,
-                self.inputs.band_stop_max,
-            ]
-            motion_metadata["BandstopFilterOrder"] = self.inputs.motion_filter_order
-
-        self._results["motion_metadata"] = motion_metadata
-
         outliers_metadata = {
             "framewise_displacement": {
                 "Description": "Outlier time series based on framewise displacement.",
diff --git a/xcp_d/utils/confounds.py b/xcp_d/utils/confounds.py
index 1b217a3..2a5e335 100644
--- a/xcp_d/utils/confounds.py
+++ b/xcp_d/utils/confounds.py
@@ -253,16 +253,16 @@ def describe_regression(params, custom_confounds_file, motion_filter_type):
 
 
 @fill_doc
-def describe_censoring(
+def describe_motion_parameters(
+    *,
     motion_filter_type,
     motion_filter_order,
     band_stop_min,
     band_stop_max,
     head_radius,
-    fd_thresh,
-    exact_scans,
+    TR,
 ):
-    """Build a text description of the motion parameter filtering and FD censoring process.
+    """Build a text description of the motion parameter filtering and FD calculation process.
 
     Parameters
     ----------
@@ -271,60 +271,107 @@ def describe_censoring(
     %(band_stop_min)s
     %(band_stop_max)s
     %(head_radius)s
-    %(fd_thresh)s
-    %(exact_scans)s
+    %(TR)s
 
     Returns
     -------
     desc : :obj:`str`
-        A text description of the censoring procedure.
+        A text description of the motion parameters.
     """
     from num2words import num2words
 
-    filter_str = ""
+    desc = ""
     if motion_filter_type:
+        band_stop_min_adjusted, band_stop_max_adjusted, is_modified = _modify_motion_filter(
+            motion_filter_type=motion_filter_type,
+            band_stop_min=band_stop_min,
+            band_stop_max=band_stop_max,
+            TR=TR,
+        )
         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"
-            )
+            if is_modified:
+                desc = (
+                    "The six translation and rotation head motion traces were "
+                    f"band-stop filtered to remove signals between {band_stop_min_adjusted} and "
+                    f"{band_stop_max_adjusted} breaths-per-minute "
+                    f"(automatically modified from {band_stop_min} and {band_stop_max} BPM due "
+                    "to Nyquist frequency constraints) using a(n) "
+                    f"{num2words(motion_filter_order, ordinal=True)}-order notch filter, "
+                    "based on @fair2020correction. "
+                )
+            else:
+                desc = (
+                    "The six translation and rotation head motion traces were "
+                    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"
-            )
+            if is_modified:
+                desc = (
+                    "The six translation and rotation head motion traces were "
+                    f"low-pass filtered below {band_stop_min_adjusted} breaths-per-minute "
+                    f"(automatically modified from {band_stop_min} BPM due to Nyquist frequency "
+                    "constraints) using a(n) "
+                    f"{num2words(motion_filter_order, ordinal=True)}-order Butterworth filter, "
+                    "based on @gratton2020removal. "
+                )
+            else:
+                desc = (
+                    "The six translation and rotation head motion traces were "
+                    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. "
+                )
+
+        desc += "The Volterra expansion of these filtered motion parameters was then calculated. "
+
+    desc += (
+        "Framewise displacement was calculated from the "
+        f"{'filtered ' if motion_filter_type else ''}motion parameters using the formula from "
+        f"@power_fd_dvars, with a head radius of {head_radius} mm. "
+    )
 
-        filter_str = (
-            f"the six translation and rotation head motion traces were {filter_sub_str}. Next, "
-        )
+    return desc
+
+
+@fill_doc
+def describe_censoring(motion_filter_type, fd_thresh, exact_scans):
+    """Build a text description of the FD censoring process.
+
+    Parameters
+    ----------
+    %(motion_filter_type)s
+    %(fd_thresh)s
+    %(exact_scans)s
 
-    outlier_str = ""
+    Returns
+    -------
+    desc : :obj:`str`
+        A text description of the censoring procedure.
+    """
+    desc = ""
     if fd_thresh > 0:
-        outlier_str = (
-            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 of {head_radius} mm. "
+        desc += (
             f"Volumes with {'filtered ' if motion_filter_type else ''}framewise displacement "
             f"greater than {fd_thresh} mm were flagged as high-motion outliers for the sake of "
             "later censoring [@power_fd_dvars]."
         )
 
-    exact_str = ""
     if exact_scans and (fd_thresh > 0):
-        exact_str = (
+        desc += (
             " Additional sets of censoring volumes were randomly selected to produce additional "
             f"correlation matrices limited to {list_to_str(exact_scans)} volumes."
         )
     elif exact_scans:
-        exact_str = (
+        desc += (
             "Volumes were randomly selected for censoring, to produce additional correlation "
             f"matrices limited to {list_to_str(exact_scans)} volumes."
         )
 
-    return outlier_str + exact_str
+    return desc
 
 
 def _get_acompcor_confounds(confounds_file):
@@ -584,9 +631,73 @@ def motion_regression_filter(
     if motion_filter_type not in ("lp", "notch"):
         raise ValueError(f"Motion filter type '{motion_filter_type}' not supported.")
 
+    lowpass_hz = band_stop_min / 60
+
+    sampling_frequency = 1 / TR
+
+    if motion_filter_type == "lp":  # low-pass filter
+        b, a = butter(
+            motion_filter_order / 2,
+            lowpass_hz,
+            btype="lowpass",
+            output="ba",
+            fs=sampling_frequency,
+        )
+        filtered_data = filtfilt(b, a, data, axis=0, padtype="constant", padlen=data.shape[0] - 1)
+
+    else:  # notch filter
+        highpass_hz = band_stop_max / 60
+        stopband_hz = np.array([lowpass_hz, highpass_hz])
+        # Convert stopband to a single notch frequency.
+        freq_to_remove = np.mean(stopband_hz)
+        bandwidth = np.abs(np.diff(stopband_hz))
+
+        # 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
+
+
+def _modify_motion_filter(motion_filter_type, band_stop_min, band_stop_max, TR):
+    """Modify the motion filter parameters based on the TR.
+
+    Parameters
+    ----------
+    motion_filter_type : str
+        The type of motion filter to apply.
+    band_stop_min : float
+        The minimum frequency to stop in the filter, in breaths-per-minute.
+    band_stop_max : float
+        The maximum frequency to stop in the filter, in breaths-per-minute.
+    TR : float
+        The repetition time of the data.
+
+    Returns
+    -------
+    band_stop_min_adjusted : float
+        The adjusted low-pass filter frequency, in breaths-per-minute.
+    band_stop_max_adjusted : float
+        The adjusted high-pass filter frequency, in breaths-per-minute.
+    is_modified : bool
+        Whether the filter parameters were modified.
+    """
     sampling_frequency = 1 / TR
     nyquist_frequency = sampling_frequency / 2
+    nyquist_bpm = nyquist_frequency * 60
 
+    is_modified = False
     if motion_filter_type == "lp":  # low-pass filter
         # Remove any frequencies above band_stop_min.
         assert band_stop_min is not None
@@ -603,20 +714,15 @@ def motion_regression_filter(
             - (np.floor((lowpass_hz + nyquist_frequency) / sampling_frequency))
             * sampling_frequency
         )
-        if lowpass_hz_adjusted != lowpass_hz:
+        band_stop_min_adjusted = lowpass_hz_adjusted * 60  # change Hertz back to BPM
+        if band_stop_min_adjusted != band_stop_min:
             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)."
+                f"Low-pass filter frequency is above Nyquist frequency ({nyquist_bpm} BPM), "
+                f"so it has been changed ({band_stop_min} --> {band_stop_min_adjusted} BPM)."
             )
+            is_modified = True
 
-        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)
+        band_stop_max_adjusted = None
 
     elif motion_filter_type == "notch":  # notch filter
         # Retain any frequencies *outside* the band_stop_min-band_stop_max range.
@@ -624,7 +730,7 @@ def motion_regression_filter(
         assert band_stop_min is not None
         assert band_stop_max > 0
         assert band_stop_min > 0
-        assert band_stop_min < band_stop_max
+        assert band_stop_min < band_stop_max, f"{band_stop_min} >= {band_stop_max}"
 
         stopband = np.array([band_stop_min, band_stop_max])
         stopband_hz = stopband / 60  # change BPM to Hertz
@@ -636,35 +742,21 @@ def motion_regression_filter(
             - (np.floor((stopband_hz + nyquist_frequency) / sampling_frequency))
             * sampling_frequency
         )
-        if not np.array_equal(stopband_hz, stopband_hz_adjusted):
+        stopband_adjusted = stopband_hz_adjusted * 60  # change Hertz back to BPM
+        if not np.array_equal(stopband_adjusted, stopband):
             warnings.warn(
-                "One or both filter frequencies are above Nyquist frequency "
-                f"({nyquist_frequency} Hz), "
+                f"One or both filter frequencies are above Nyquist frequency ({nyquist_bpm} BPM), "
                 "so they have been changed "
-                f"({stopband_hz[0]} --> {stopband_hz_adjusted[0]}, "
-                f"{stopband_hz[1]} --> {stopband_hz_adjusted[1]} Hz)."
+                f"({stopband[0]} --> {stopband_adjusted[0]}, "
+                f"{stopband[1]} --> {stopband_adjusted[1]} BPM)."
             )
+            is_modified = True
 
-        # 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,
-            )
+        band_stop_min_adjusted, band_stop_max_adjusted = stopband_adjusted
+    else:
+        band_stop_min_adjusted, band_stop_max_adjusted, is_modified = None, None, False
 
-    return filtered_data
+    return band_stop_min_adjusted, band_stop_max_adjusted, is_modified
 
 
 def _infer_dummy_scans(dummy_scans, confounds_file=None):
diff --git a/xcp_d/utils/modified_data.py b/xcp_d/utils/modified_data.py
index aa78b70..9efb2c9 100644
--- a/xcp_d/utils/modified_data.py
+++ b/xcp_d/utils/modified_data.py
@@ -8,7 +8,7 @@ import numpy as np
 import pandas as pd
 from nipype import logging
 
-from xcp_d.utils.confounds import _infer_dummy_scans, load_motion
+from xcp_d.utils.confounds import _infer_dummy_scans, _modify_motion_filter, load_motion
 from xcp_d.utils.doc import fill_doc
 from xcp_d.utils.filemanip import fname_presuffix
 
@@ -183,13 +183,19 @@ def flag_bad_run(
     fmriprep_confounds_df = fmriprep_confounds_df.drop(np.arange(dummy_scans))
 
     # Calculate filtered FD
+    band_stop_min_adjusted, band_stop_max_adjusted, _ = _modify_motion_filter(
+        motion_filter_type=motion_filter_type,
+        band_stop_min=band_stop_min,
+        band_stop_max=band_stop_max,
+        TR=TR,
+    )
     motion_df = load_motion(
         fmriprep_confounds_df,
         TR=TR,
         motion_filter_type=motion_filter_type,
         motion_filter_order=motion_filter_order,
-        band_stop_min=band_stop_min,
-        band_stop_max=band_stop_max,
+        band_stop_min=band_stop_min_adjusted,
+        band_stop_max=band_stop_max_adjusted,
     )
     fd_arr = compute_fd(confound=motion_df, head_radius=head_radius)
     return np.sum(fd_arr <= fd_thresh) * TR
diff --git a/xcp_d/workflows/bold.py b/xcp_d/workflows/bold.py
index 4930e78..96aea8e 100644
--- a/xcp_d/workflows/bold.py
+++ b/xcp_d/workflows/bold.py
@@ -190,10 +190,12 @@ def init_postprocess_nifti_wf(
         run_data["confounds"],
     )
 
-    workflow.__desc__ = (
-        f"For each of the {num2words(n_runs)} BOLD runs found per subject "
-        "(across all tasks and sessions), the following post-processing was performed."
-    )
+    workflow.__desc__ = f"""\
+
+For each of the {num2words(n_runs)} BOLD runs found per subject (across all tasks and sessions),
+the following post-processing was performed.
+
+"""
 
     outputnode = pe.Node(
         niu.IdentityInterface(
diff --git a/xcp_d/workflows/cifti.py b/xcp_d/workflows/cifti.py
index d27fa62..a2ce0cb 100644
--- a/xcp_d/workflows/cifti.py
+++ b/xcp_d/workflows/cifti.py
@@ -179,10 +179,12 @@ def init_postprocess_cifti_wf(
 
     workflow = Workflow(name=name)
 
-    workflow.__desc__ = (
-        f"For each of the {num2words(n_runs)} BOLD runs found per subject "
-        "(across all tasks and sessions), the following post-processing was performed."
-    )
+    workflow.__desc__ = f"""\
+
+For each of the {num2words(n_runs)} BOLD runs found per subject (across all tasks and sessions),
+the following post-processing was performed.
+
+"""
 
     outputnode = pe.Node(
         niu.IdentityInterface(
diff --git a/xcp_d/workflows/postprocessing.py b/xcp_d/workflows/postprocessing.py
index e976f68..4a4bbfc 100644
--- a/xcp_d/workflows/postprocessing.py
+++ b/xcp_d/workflows/postprocessing.py
@@ -20,7 +20,11 @@ from xcp_d.interfaces.nilearn import DenoiseCifti, DenoiseNifti, Smooth
 from xcp_d.interfaces.plotting import CensoringPlot
 from xcp_d.interfaces.restingstate import DespikePatch
 from xcp_d.interfaces.workbench import CiftiConvert, FixCiftiIntent
-from xcp_d.utils.confounds import describe_censoring, describe_regression
+from xcp_d.utils.confounds import (
+    describe_censoring,
+    describe_motion_parameters,
+    describe_regression,
+)
 from xcp_d.utils.doc import fill_doc
 from xcp_d.utils.plotting import plot_design_matrix as _plot_design_matrix
 from xcp_d.utils.utils import fwhm2sigma
@@ -123,13 +127,18 @@ def init_prepare_confounds_wf(
             "regressors were discarded as non-steady-state volumes, or 'dummy scans'. "
         )
 
+    motion_description = describe_motion_parameters(
+        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,
+        TR=TR,
+    )
+
     if (fd_thresh > 0) or exact_scans:
         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,
             exact_scans=exact_scans,
         )
@@ -142,7 +151,9 @@ def init_prepare_confounds_wf(
         motion_filter_type=motion_filter_type,
     )
 
-    workflow.__desc__ = f" {dummy_scans_str}{censoring_description} {confounds_description}"
+    workflow.__desc__ = (
+        f" {dummy_scans_str}{motion_description} {censoring_description} {confounds_description}"
+    )
 
     inputnode = pe.Node(
         niu.IdentityInterface(
@@ -533,6 +544,7 @@ def init_denoise_bold_wf(TR, name="denoise_bold_wf"):
     """
     workflow = Workflow(name=name)
 
+    fd_thresh = config.workflow.fd_thresh
     low_pass = config.workflow.low_pass
     high_pass = config.workflow.high_pass
     bpf_order = config.workflow.bpf_order
@@ -542,10 +554,20 @@ def init_denoise_bold_wf(TR, name="denoise_bold_wf"):
     mem_gb = config.nipype.memory_gb
     omp_nthreads = config.nipype.omp_nthreads
 
-    workflow.__desc__ = (
-        "Nuisance regressors were regressed from the BOLD data using linear regression, "
-        "as implemented in *Nilearn*."
-    )
+    workflow.__desc__ = """\
+
+Nuisance regressors were regressed from the BOLD data using a denoising method based on *Nilearn*'s
+approach.
+"""
+    if fd_thresh > 0:
+        workflow.__desc__ += (
+            "Any volumes censored earlier in the workflow were first cubic spline interpolated in "
+            "the BOLD data. "
+            "Outlier volumes at the beginning or end of the time series were replaced with the "
+            "closest low-motion volume's values, "
+            "as cubic spline interpolation can produce extreme extrapolations. "
+        )
+
     if bandpass_filter:
         if low_pass > 0 and high_pass > 0:
             btype = "band-pass"
@@ -561,13 +583,23 @@ def init_denoise_bold_wf(TR, name="denoise_bold_wf"):
             filt_input = f"{low_pass}"
 
         workflow.__desc__ += (
-            " Any volumes censored earlier in the workflow were then interpolated in the residual "
-            "time series produced by the regression. "
-            f"The interpolated timeseries were then {btype} filtered using a(n) "
+            f"The timeseries were {btype} filtered using a(n) "
             f"{num2words(bpf_order, ordinal=True)}-order Butterworth filter, "
             f"in order to retain signals {preposition} {filt_input} Hz. "
-            "The filtered, interpolated time series were then re-censored to remove high-motion "
-            "outlier volumes."
+            "The same filter was applied to the confounds."
+        )
+
+    if fd_thresh > 0:
+        workflow.__desc__ += (
+            " The resulting time series were then denoised via linear regression, "
+            "in which the low-motion volumes from the BOLD time series and confounds were used to "
+            "calculate parameter estimates, and then the interpolated time series were denoised "
+            "using the low-motion parameter estimates. "
+            "The interpolated time series were then censored using the temporal mask."
+        )
+    else:
+        workflow.__desc__ += (
+            " The resulting time series were then denoised using linear regression. "
         )
 
     inputnode = pe.Node(
@@ -606,7 +638,6 @@ def init_denoise_bold_wf(TR, name="denoise_bold_wf"):
         n_procs=omp_nthreads,
     )
 
-    # fmt:off
     workflow.connect([
         (inputnode, regress_and_filter_bold, [
             ("preprocessed_bold", "preprocessed_bold"),
@@ -616,10 +647,9 @@ def init_denoise_bold_wf(TR, name="denoise_bold_wf"):
         (regress_and_filter_bold, outputnode, [
             ("denoised_interpolated_bold", "denoised_interpolated_bold"),
         ]),
-    ])
+    ])  # fmt:skip
     if not cifti:
         workflow.connect([(inputnode, regress_and_filter_bold, [("mask", "mask")])])
-    # fmt:on
 
     censor_interpolated_data = pe.Node(
         Censor(),
@@ -628,7 +658,6 @@ def init_denoise_bold_wf(TR, name="denoise_bold_wf"):
         omp_nthreads=omp_nthreads,
     )
 
-    # fmt:off
     workflow.connect([
         (inputnode, censor_interpolated_data, [("temporal_mask", "temporal_mask")]),
         (regress_and_filter_bold, censor_interpolated_data, [
@@ -637,13 +666,11 @@ def init_denoise_bold_wf(TR, name="denoise_bold_wf"):
         (censor_interpolated_data, outputnode, [
             ("censored_denoised_bold", "censored_denoised_bold"),
         ]),
-    ])
-    # fmt:on
+    ])  # fmt:skip
 
     if smoothing:
         resd_smoothing_wf = init_resd_smoothing_wf()
 
-        # fmt:off
         workflow.connect([
             (censor_interpolated_data, resd_smoothing_wf, [
                 ("censored_denoised_bold", "inputnode.bold_file"),
@@ -651,8 +678,7 @@ def init_denoise_bold_wf(TR, name="denoise_bold_wf"):
             (resd_smoothing_wf, outputnode, [
                 ("outputnode.smoothed_bold", "smoothed_denoised_bold"),
             ]),
-        ])
-        # fmt:on
+        ])  # fmt:skip
 
     return workflow
 
diff --git a/xcp_d/workflows/restingstate.py b/xcp_d/workflows/restingstate.py
index 8d2c8e6..6de9856 100644
--- a/xcp_d/workflows/restingstate.py
+++ b/xcp_d/workflows/restingstate.py
@@ -103,6 +103,7 @@ def init_alff_wf(
         )
 
     workflow.__desc__ = f""" \
+
 The amplitude of low-frequency fluctuation (ALFF) [@alff] was computed by transforming
 the mean-centered, standard deviation-normalized, denoised BOLD time series to the frequency
 domain{periodogram_desc}.


Test patch:
diff --git a/xcp_d/tests/test_utils_confounds.py b/xcp_d/tests/test_utils_confounds.py
index cd8be4b..f531bd0 100644
--- a/xcp_d/tests/test_utils_confounds.py
+++ b/xcp_d/tests/test_utils_confounds.py
@@ -9,6 +9,72 @@ from scipy import signal
 from xcp_d.utils import confounds
 
 
+def test_modify_motion_filter():
+    """Run a simple test of the motion filter modification function."""
+    band_stop_min = 6
+    TR = 0.8
+
+    with pytest.warns(match="The parameter 'band_stop_max' will be ignored."):
+        band_stop_min2, _, is_modified = confounds._modify_motion_filter(
+            motion_filter_type="lp",
+            band_stop_min=band_stop_min,
+            band_stop_max=18,
+            TR=TR,
+        )
+    assert band_stop_min2 == band_stop_min
+    assert is_modified is False
+
+    # 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 (37.5 BPM), "
+            "so it has been changed (42 --> 33.0 BPM)."
+        ),
+    ):
+        band_stop_min2, _, is_modified = confounds._modify_motion_filter(
+            TR=TR,  # 1.25 Hz
+            motion_filter_type="lp",
+            band_stop_min=42,  # 0.7 Hz > (1.25 / 2)
+            band_stop_max=None,
+        )
+    assert band_stop_min2 == 33.0
+    assert is_modified is True
+
+    # Now test band-stop filter
+    # 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 (37.5 BPM), "
+            "so they have been changed (42 --> 33.0, 45 --> 30.0 BPM)."
+        ),
+    ):
+        band_stop_min2, band_stop_max2, is_modified = confounds._modify_motion_filter(
+            TR=TR,
+            motion_filter_type="notch",
+            band_stop_min=42,
+            band_stop_max=45,  # 0.7 Hz > (1.25 / 2)
+        )
+
+    assert band_stop_min2 == 33.0
+    assert band_stop_max2 == 30.0
+    assert is_modified is True
+
+    # Notch without modification
+    band_stop_min2, band_stop_max2, is_modified = confounds._modify_motion_filter(
+        TR=TR,
+        motion_filter_type="notch",
+        band_stop_min=30,
+        band_stop_max=33,
+    )
+
+    assert band_stop_min2 == 30
+    assert band_stop_max2 == 33
+    assert is_modified is False
+
+
 def test_motion_filtering_lp():
     """Run low-pass filter on toy data, compare to simplified results."""
     raw_data = np.random.random(500)
@@ -34,32 +100,25 @@ def test_motion_filtering_lp():
 
     # 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 = confounds.motion_regression_filter(
-            raw_data,
-            TR=TR,
-            motion_filter_type="lp",
-            band_stop_min=band_stop_min,
-            band_stop_max=18,
-            motion_filter_order=2,
-        )
+    lowpass_data_test = confounds.motion_regression_filter(
+        raw_data,
+        TR=TR,
+        motion_filter_type="lp",
+        band_stop_min=band_stop_min,
+        band_stop_max=None,
+        motion_filter_order=2,
+    )
 
     # What's the difference from the verified data?
     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)."
-        ),
-    ):
+    # Using a filter type other than notch or lp should raise an exception.
+    with pytest.raises(ValueError, match="Motion filter type 'fail' not supported."):
         confounds.motion_regression_filter(
             raw_data,
-            TR=TR,  # 1.25 Hz
-            motion_filter_type="lp",
-            band_stop_min=42,  # 0.7 Hz > (1.25 / 2)
+            TR=TR,
+            motion_filter_type="fail",
+            band_stop_min=band_stop_min,
             band_stop_max=None,
             motion_filter_order=2,
         )
@@ -100,31 +159,92 @@ def test_motion_filtering_notch():
     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)."
-        ),
-    ):
-        confounds.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."):
-        confounds.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_describe_motion_parameters():
+    """Test confounds.describe_motion_parameters."""
+    # notch with no modification
+    desc = confounds.describe_motion_parameters(
+        motion_filter_type="notch",
+        motion_filter_order=2,
+        band_stop_min=12,
+        band_stop_max=20,
+        head_radius=50,
+        TR=0.8,
+    )
+    assert isinstance(desc, str)
+
+    # notch with modification
+    desc = confounds.describe_motion_parameters(
+        motion_filter_type="notch",
+        motion_filter_order=2,
+        band_stop_min=42,
+        band_stop_max=45,
+        head_radius=50,
+        TR=0.8,
+    )
+    assert isinstance(desc, str)
+
+    # lowpass with no modification
+    desc = confounds.describe_motion_parameters(
+        motion_filter_type="lp",
+        motion_filter_order=2,
+        band_stop_min=12,
+        band_stop_max=None,
+        head_radius=50,
+        TR=0.8,
+    )
+    assert isinstance(desc, str)
+
+    # notch with modification
+    desc = confounds.describe_motion_parameters(
+        motion_filter_type="lp",
+        motion_filter_order=2,
+        band_stop_min=42,
+        band_stop_max=None,
+        head_radius=50,
+        TR=0.8,
+    )
+    assert isinstance(desc, str)
+
+
+def test_describe_censoring():
+    """Test confounds.describe_censoring."""
+    # notch filter, no censoring, no exact-scan
+    desc = confounds.describe_censoring(
+        motion_filter_type="notch",
+        fd_thresh=0,
+        exact_scans=[],
+    )
+    assert isinstance(desc, str)
+
+    # notch filter, censoring, no exact-scan
+    desc = confounds.describe_censoring(
+        motion_filter_type="notch",
+        fd_thresh=0.1,
+        exact_scans=[],
+    )
+    assert isinstance(desc, str)
+
+    # notch filter, censoring, exact-scan
+    desc = confounds.describe_censoring(
+        motion_filter_type="notch",
+        fd_thresh=0.1,
+        exact_scans=[100, 150],
+    )
+    assert isinstance(desc, str)
+
+    # no filter, no censoring, no exact-scan
+    desc = confounds.describe_censoring(
+        motion_filter_type=None,
+        fd_thresh=0,
+        exact_scans=[],
+    )
+    assert isinstance(desc, str)
+
+    # no filter, no censoring, exact-scan
+    desc = confounds.describe_censoring(
+        motion_filter_type=None,
+        fd_thresh=0,
+        exact_scans=[100, 150],
+    )
+    assert isinstance(desc, str)