Repository: PennLINC/xcp_d

Issue:
Possible bug in DCANBoldProc that may or may not be relevant to XCP-D
@cindyhfls noticed a possible bug in our code base for DCANBoldProc. Our team is going to check it out soon. I wanted to share in case the issue is present in XCP-D as well! 

See: https://github.com/DCAN-Labs/dcan_bold_processing/issues/28

Discussion:
tsalo: @cindyhfls would you mind reviewing my changes? Also, given that you detected the problem and proposed the fix, would you like me to add you to our CITATION.cff file? If so, can you provide your information to match this?

https://github.com/PennLINC/xcp_d/blob/16b16f72817b79b5c98382e5c4e4abd5da1be183/CITATION.cff#L13-L17

Patch:
diff --git a/xcp_d/cli/parser.py b/xcp_d/cli/parser.py
index b229afe..1c97449 100644
--- a/xcp_d/cli/parser.py
+++ b/xcp_d/cli/parser.py
@@ -373,7 +373,15 @@ This parameter is used in conjunction with ``motion-filter-order`` and ``band-st
         dest="motion_filter_order",
         default=4,
         type=int,
-        help="Number of filter coeffecients for the motion parameter filter.",
+        help=(
+            "Number of filter coefficients for the motion parameter filter. "
+            "If the motion filter type is 'lp', the order is divided by 2 as filtfilt applies "
+            "the filter twice. "
+            "If the motion filter type is 'notch', the order is divided by 4 as iirnotch is a "
+            "second-order filter and filtfilt applies the filter twice. "
+            "Make sure to set this parameter to a multiple of 2 if you choose the 'lp' filter and "
+            "a multiple of 4 if you choose the 'notch' filter."
+        ),
     )
 
     g_censor = parser.add_argument_group("Censoring and scrubbing options")
diff --git a/xcp_d/interfaces/censoring.py b/xcp_d/interfaces/censoring.py
index 5a5e003..92d54a0 100644
--- a/xcp_d/interfaces/censoring.py
+++ b/xcp_d/interfaces/censoring.py
@@ -523,7 +523,7 @@ class ProcessMotion(SimpleInterface):
                 filters = col_metadata.get("SoftwareFilters", {})
                 filters["Butterworth low-pass filter"] = {
                     "cutoff": band_stop_min_adjusted / 60,
-                    "order": self.inputs.motion_filter_order,
+                    "order": int(np.floor(self.inputs.motion_filter_order / 2)),
                     "cutoff units": "Hz",
                     "function": "scipy.signal.filtfilt",
                 }
@@ -536,7 +536,7 @@ class ProcessMotion(SimpleInterface):
                         band_stop_max_adjusted / 60,
                         band_stop_min_adjusted / 60,
                     ],
-                    "order": self.inputs.motion_filter_order,
+                    "order": int(np.floor(self.inputs.motion_filter_order / 4)),
                     "cutoff units": "Hz",
                     "function": "scipy.signal.filtfilt",
                 }
@@ -879,7 +879,7 @@ class GenerateConfounds(SimpleInterface):
                         filters = col_metadata.get("SoftwareFilters", {})
                         filters["Butterworth low-pass filter"] = {
                             "cutoff": band_stop_min_adjusted / 60,
-                            "order": self.inputs.motion_filter_order,
+                            "order": int(np.floor(self.inputs.motion_filter_order / 2)),
                             "cutoff units": "Hz",
                             "function": "scipy.signal.filtfilt",
                         }
@@ -892,7 +892,7 @@ class GenerateConfounds(SimpleInterface):
                                 band_stop_max_adjusted / 60,
                                 band_stop_min_adjusted / 60,
                             ],
-                            "order": self.inputs.motion_filter_order,
+                            "order": int(np.floor(self.inputs.motion_filter_order / 4)),
                             "cutoff units": "Hz",
                             "function": "scipy.signal.filtfilt",
                         }
diff --git a/xcp_d/interfaces/connectivity.py b/xcp_d/interfaces/connectivity.py
index 2717b5f..5bce9ec 100644
--- a/xcp_d/interfaces/connectivity.py
+++ b/xcp_d/interfaces/connectivity.py
@@ -239,6 +239,7 @@ def correlate_timeseries(timeseries, temporal_mask):
         if censoring_df.shape[0] == timeseries_df.shape[0]:
             # The time series is not censored
             timeseries_df = timeseries_df.loc[censoring_df["framewise_displacement"] == 0]
+            timeseries_df.reset_index(drop=True, inplace=True)
 
         # Now create correlation matrices limited to exact scan numbers
         censored_censoring_df = censoring_df.loc[censoring_df["framewise_displacement"] == 0]
diff --git a/xcp_d/utils/boilerplate.py b/xcp_d/utils/boilerplate.py
index 4661fec..b25122d 100644
--- a/xcp_d/utils/boilerplate.py
+++ b/xcp_d/utils/boilerplate.py
@@ -29,6 +29,7 @@ def describe_motion_parameters(
     desc : :obj:`str`
         A text description of the motion parameters.
     """
+    import numpy as np
     from num2words import num2words
 
     desc = ""
@@ -40,6 +41,7 @@ def describe_motion_parameters(
             TR=TR,
         )
         if motion_filter_type == "notch":
+            n_filter_applications = int(np.floor(motion_filter_order / 4))
             if is_modified:
                 desc = (
                     "The six translation and rotation head motion traces were "
@@ -47,7 +49,7 @@ def describe_motion_parameters(
                     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, "
+                    f"{num2words(n_filter_applications, ordinal=True)}-order notch filter, "
                     "based on @fair2020correction. "
                 )
             else:
@@ -55,17 +57,18 @@ def describe_motion_parameters(
                     "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, "
+                    f"{num2words(n_filter_applications, ordinal=True)}-order notch filter, "
                     "based on @fair2020correction. "
                 )
         else:  # lp
+            n_filter_applications = int(np.floor(motion_filter_order / 2))
             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, "
+                    f"{num2words(n_filter_applications, ordinal=True)}-order Butterworth filter, "
                     "based on @gratton2020removal. "
                 )
             else:
@@ -73,7 +76,7 @@ def describe_motion_parameters(
                     "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, "
+                    f"{num2words(n_filter_applications, ordinal=True)}-order Butterworth filter, "
                     "based on @gratton2020removal. "
                 )
 
diff --git a/xcp_d/utils/confounds.py b/xcp_d/utils/confounds.py
index 42900bd..f02e5e5 100644
--- a/xcp_d/utils/confounds.py
+++ b/xcp_d/utils/confounds.py
@@ -139,6 +139,8 @@ def filter_motion(
     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.
+    Since filtfilt applies the filter twice, motion_filter_order is divided by 2 before applying
+    the 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`.
@@ -147,6 +149,8 @@ def filter_motion(
     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.
+    Because iirnotch is a second-order filter and filtfilt applies the filter twice,
+    motion_filter_order is divided by 4 before applying the filter.
     The filter is applied with constant-padding, using the default padding size determined by
     :func:`scipy.signal.filtfilt`.
 
@@ -162,8 +166,9 @@ def filter_motion(
     sampling_frequency = 1 / TR
 
     if motion_filter_type == "lp":  # low-pass filter
+        n_filter_applications = int(np.floor(motion_filter_order / 2))
         b, a = butter(
-            motion_filter_order / 2,
+            n_filter_applications,
             lowpass_hz,
             btype="lowpass",
             output="ba",
@@ -180,7 +185,9 @@ def filter_motion(
 
         # 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))
+        # Both iirnotch and filtfilt are second-order,
+        # so we need to divide the motion_filter_order by 4.
+        n_filter_applications = int(np.floor(motion_filter_order / 4))
 
         filtered_data = data.copy()
         for _ in range(n_filter_applications):
diff --git a/xcp_d/utils/doc.py b/xcp_d/utils/doc.py
index 7e0f849..b3a3f0f 100644
--- a/xcp_d/utils/doc.py
+++ b/xcp_d/utils/doc.py
@@ -247,6 +247,9 @@ motion_filter_order : :obj:`int`
     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``.
+    For "lp" filters, the order is divided by 2 as filtfilt applies the filter twice.
+    For "notch" filters, the order is divided by 4 as iirnotch is a second-order filter and
+    filtfilt applies the filter twice.
 """
 
 docdict[


Test patch:
diff --git a/xcp_d/tests/test_cli.py b/xcp_d/tests/test_cli.py
index 16003d0..939f715 100644
--- a/xcp_d/tests/test_cli.py
+++ b/xcp_d/tests/test_cli.py
@@ -103,6 +103,7 @@ def test_ds001419_cifti(data_dir, output_dir, working_dir):
         "--warp_surfaces_native2std=n",
         "--head_radius=40",
         "--motion-filter-type=notch",
+        "--motion-filter-order=4",
         "--band-stop-min=12",
         "--band-stop-max=18",
         "--dummy-scans=auto",
diff --git a/xcp_d/tests/test_utils_confounds.py b/xcp_d/tests/test_utils_confounds.py
index 5399ce4..9a73414 100644
--- a/xcp_d/tests/test_utils_confounds.py
+++ b/xcp_d/tests/test_utils_confounds.py
@@ -154,7 +154,7 @@ def test_motion_filtering_notch():
         motion_filter_type="notch",
         band_stop_min=band_stop_min,
         band_stop_max=band_stop_max,
-        motion_filter_order=2,
+        motion_filter_order=4,
     )
     notch_data_test = np.squeeze(notch_data_test)
     assert np.allclose(notch_data_test, notch_data_true)