Merge pull request #90 from smoia/fix/bandpass

Fix bandpass Butterworth filter order to get a better filter, add option to change order in CLI
smoia · Aug 13, 2023 · ca5257e · ca5257e
2 parents 8bf8a68 + ac84a22
commit ca5257e
Show file tree

Hide file tree

Showing 3 changed files with 24 additions and 12 deletions.
diff --git a/phys2cvr/cli/run.py b/phys2cvr/cli/run.py
@@ -260,6 +260,14 @@ def _get_parser():
         ),
         default=None,
     )
+    opt_filt.add_argument(
+        "-bo",
+        "--butter-order",
+        dest="butter_order",
+        type=int,
+        help=("Order of Butterworth filter. Default is 9."),
+        default=None,
+    )
 
     opt_flow = parser.add_argument_group("Optional Arguments to modify the workflow")
     opt_flow.add_argument(

diff --git a/phys2cvr/phys2cvr.py b/phys2cvr/phys2cvr.py
@@ -62,8 +62,9 @@ def phys2cvr(
     n_trials=None,
     abs_xcorr=False,
     skip_xcorr=False,
-    highcut=None,
-    lowcut=None,
+    highcut=0.04,
+    lowcut=0.02,
+    butter_order=9,
     apply_filter=False,
     run_regression=False,
     lagged_regression=True,
@@ -140,11 +141,14 @@ def phys2cvr(
     highcut : str, int, or float, optional
         High frequency limit for filter.
         Required if applying a filter.
-        Default: empty
+        Default: 0.02
     lowcut : str, int, or float, optional
         Low frequency limit for filter.
         Required if applying a filter.
-        Default: empty
+        Default: 0.04
+    butter_order : int, optional
+        Butterworth filter order.
+        Default: 9
     apply_filter : bool, optional
         Apply a Butterworth filter to the functional input.
         Default: False
@@ -322,7 +326,9 @@ def phys2cvr(
             LGR.info(f"Loading {fname_func}")
             if apply_filter:
                 LGR.info("Applying butterworth filter to {fname_func}")
-                func_avg = signal.filter_signal(func_avg, tr, lowcut, highcut)
+                func_avg = signal.filter_signal(
+                    func_avg, tr, lowcut, highcut, butter_order
+                )
         else:
             raise NameError(
                 "Provided functional signal, but no TR specified! "
@@ -372,7 +378,7 @@ def phys2cvr(
 
         if apply_filter:
             LGR.info(f"Obtaining filtered average signal in {roiref}")
-            func_filt = signal.filter_signal(func, tr, lowcut, highcut)
+            func_filt = signal.filter_signal(func, tr, lowcut, highcut, butter_order)
             func_avg = func_filt[roi].mean(axis=0)
         else:
             LGR.info(f"Obtaining average signal in {roiref}")

diff --git a/phys2cvr/signal.py b/phys2cvr/signal.py
@@ -90,7 +90,7 @@ def create_hrf(freq=40):
     return hrf
 
 
-def filter_signal(data, tr, lowcut="", highcut=""):
+def filter_signal(data, tr, lowcut=0.02, highcut=0.04, order=9):
     """
     Create a bandpass filter given a lowcut and a highcut, then filter data.
 
@@ -104,20 +104,18 @@ def filter_signal(data, tr, lowcut="", highcut=""):
         Lower frequency in the bandpass
     highcut : float
         Higher frequency in the bandpass
+    order : int
+        The order of the butterworth filter
 
     Returns
     -------
     filt_data : np.ndarray
         Input `data`, but filtered.
     """
-    if not lowcut:
-        lowcut = 0.02
-    if not highcut:
-        highcut = 0.04
     nyq = (1 / tr) / 2
     low = lowcut / nyq
     high = highcut / nyq
-    a, b = butter(1, [low, high], btype="band")
+    a, b = butter(order, [low, high], btype="band")
     filt_data = filtfilt(a, b, data, axis=-1)
     return filt_data