From ac84a222c7e80330a577564915c3e06e38928546 Mon Sep 17 00:00:00 2001
From: smoia <s.moia.research@gmail.com>
Date: Sat, 12 Aug 2023 23:41:02 +0200
Subject: [PATCH] Fix bandpass filter order to get better bandpass, add option
 to CLI

---
 phys2cvr/cli/run.py  |  8 ++++++++
 phys2cvr/phys2cvr.py | 18 ++++++++++++------
 phys2cvr/signal.py   | 10 ++++------
 3 files changed, 24 insertions(+), 12 deletions(-)

diff --git a/phys2cvr/cli/run.py b/phys2cvr/cli/run.py
index 087621c..c700eb0 100644
--- 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
index dab319f..a1e7230 100755
--- 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
index c8c509e..9520e3c 100644
--- 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