Apply notch filter appropriate number of times based on `--motion-fil…

…ter-order` (PennLINC#1300)

xcp_d/cli/parser.py

@@ -373,7 +373,15 @@ def _build_parser():
         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")

xcp_d/interfaces/censoring.py

@@ -523,7 +523,7 @@ def _run_interface(self, runtime):
                 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 @@ def _run_interface(self, runtime):
                         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 @@ def _run_interface(self, runtime):
                         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 @@ def _run_interface(self, runtime):
                                 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",
                         }

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]

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",

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)

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,40 +41,42 @@ 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 "
                     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, "
+                    f"{num2words(n_filter_applications, 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, "
+                    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:
                 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, "
+                    f"{num2words(n_filter_applications, ordinal=True)}-order Butterworth filter, "
                     "based on @gratton2020removal. "
                 )
 

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):

xcp_d/utils/doc.py

@@ -247,6 +247,9 @@
     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[