filter report added

Author: ryanhammonds

neurodsp/filt/checks.py

@@ -90,7 +90,7 @@ def check_filter_definition(pass_type, f_range):
 
 
 def check_filter_properties(filter_coefs, a_vals, fs, pass_type, f_range,
-                            transitions=(-20, -3), verbose=True):
+                            transitions=(-20, -3), filt_type=None, verbose=True):
     """Check a filters properties, including pass band and transition band.
 
     Parameters
@@ -164,7 +164,8 @@ def check_filter_properties(filter_coefs, a_vals, fs, pass_type, f_range,
 
     # Compute pass & transition bandwidth
     pass_bw = compute_pass_band(fs, pass_type, f_range)
-    transition_bw = compute_transition_band(f_db, db, transitions[0], transitions[1])
+    transition_bw, f_range_trans = compute_transition_band(f_db, db, transitions[0], transitions[1],
+                                                           return_freqs=True)
 
     # Raise warning if transition bandwidth is too high
     if transition_bw > pass_bw:
@@ -174,8 +175,50 @@ def check_filter_properties(filter_coefs, a_vals, fs, pass_type, f_range,
 
     # Print out transition bandwidth and pass bandwidth to the user
     if verbose:
-        print('Transition bandwidth is {:.1f} Hz.'.format(transition_bw))
-        print('Pass/stop bandwidth is {:.1f} Hz.'.format(pass_bw))
+
+        # Filter type (high-pass, low-pass, band-pass, band-stop, FIR, IIR)
+        print('Pass Type: {pass_type}'.format(pass_type=pass_type))
+
+        # Cutoff frequency (including definition)
+        cutoff = round(np.min(f_range) + (0.5 * transition_bw), 3)
+        print('Cutoff (half-amplitude): {cutoff} Hz'.format(cutoff=cutoff))
+
+        # Filter order (or length)
+        print('Filter order: {order}'.format(order=len(f_db)-1))
+
+        # Roll-off or transition bandwidth
+        print('Transition bandwidth: {:.1f} Hz'.format(transition_bw))
+        print('Pass/stop bandwidth: {:.1f} Hz'.format(pass_bw))
+
+        # Passband ripple and stopband attenuation
+        pb_ripple = np.max(db[:np.where(f_db < f_range_trans[0])[0][-1]])
+        sb_atten = np.max(db[np.where(f_db > f_range_trans[1])[0][0]:])
+        print('Passband Ripple: {pb_ripple} db'.format(pb_ripple=pb_ripple))
+        print('Stopband Attenuation: {sb_atten} db'.format(sb_atten=sb_atten))
+
+        # Filter delay (zero-phase, linear-phase, non-linear phase)
+        if filt_type == 'FIR' and pass_type in ['bandstop', 'lowpass']:
+            filt_class = 'linear-phase'
+        elif filt_type == 'FIR' and pass_type in ['bandpass', 'highpass']:
+            filt_class = 'zero-phase'
+        elif filt_type == 'IIR':
+            filt_class = 'non-linear-phase'
+        else:
+            filt_class = None
+
+        if filt_type is not None:
+            print('Filter Class: {filt_class}'.format(filt_class=filt_class))
+
+        if filt_class == 'linear-phase':
+            print('Group Delay: {delay}s'.format(delay=(len(f_db)-1) / 2 * fs))
+        elif filt_class == 'zero-phase':
+            print('Group Delay: 0s')
+
+        # Direction of computation (one-pass forward/reverse, or two-pass forward and reverse)
+        if filt_type == 'FIR':
+            print('Direction: one-pass reverse.')
+        else:
+            print('Direction: two-pass forward and reverse')
 
     return passes
 

neurodsp/filt/filter.py

@@ -10,7 +10,8 @@
 
 def filter_signal(sig, fs, pass_type, f_range, filter_type='fir',
                   n_cycles=3, n_seconds=None, remove_edges=True, butterworth_order=None,
-                  print_transitions=False, plot_properties=False, return_filter=False):
+                  print_transitions=False, plot_properties=False, return_filter=False,
+                  verbose=False):
     """Apply a bandpass, bandstop, highpass, or lowpass filter to a neural signal.
 
     Parameters
@@ -52,6 +53,8 @@ def filter_signal(sig, fs, pass_type, f_range, filter_type='fir',
         If True, plot the properties of the filter, including frequency response and/or kernel.
     return_filter : bool, optional, default: False
         If True, return the filter coefficients.
+    verbose : bool, optional, default: False
+        If True, print out detailed filter information.
 
     Returns
     -------
@@ -74,12 +77,12 @@ def filter_signal(sig, fs, pass_type, f_range, filter_type='fir',
     if filter_type.lower() == 'fir':
         return filter_signal_fir(sig, fs, pass_type, f_range, n_cycles, n_seconds,
                                  remove_edges, print_transitions,
-                                 plot_properties, return_filter)
+                                 plot_properties, return_filter, verbose=verbose)
     elif filter_type.lower() == 'iir':
         _iir_checks(n_seconds, butterworth_order, remove_edges)
         return filter_signal_iir(sig, fs, pass_type, f_range, butterworth_order,
                                  print_transitions, plot_properties,
-                                 return_filter)
+                                 return_filter, verbose=verbose)
     else:
         raise ValueError('Filter type not understood.')
 

neurodsp/filt/fir.py

@@ -14,7 +14,8 @@
 ###################################################################################################
 
 def filter_signal_fir(sig, fs, pass_type, f_range, n_cycles=3, n_seconds=None, remove_edges=True,
-                      print_transitions=False, plot_properties=False, return_filter=False):
+                      print_transitions=False, plot_properties=False, return_filter=False,
+                      verbose=False):
     """Apply an FIR filter to a signal.
 
     Parameters
@@ -48,6 +49,8 @@ def filter_signal_fir(sig, fs, pass_type, f_range, n_cycles=3, n_seconds=None, r
         If True, plot the properties of the filter, including frequency response and/or kernel.
     return_filter : bool, optional, default: False
         If True, return the filter coefficients of the FIR filter.
+    verbose : bool, optional, default: False
+        If True, print out detailed filter information.
 
     Returns
     -------
@@ -78,7 +81,8 @@ def filter_signal_fir(sig, fs, pass_type, f_range, n_cycles=3, n_seconds=None, r
     check_filter_length(sig.shape[-1], len(filter_coefs))
 
     # Check filter properties: compute transition bandwidth & run checks
-    check_filter_properties(filter_coefs, 1, fs, pass_type, f_range, verbose=print_transitions)
+    check_filter_properties(filter_coefs, 1, fs, pass_type, f_range, filt_type="FIR",
+                            verbose=np.any([print_transitions, verbose]))
 
     # Remove any NaN on the edges of 'sig'
     sig, sig_nans = remove_nans(sig)

neurodsp/filt/iir.py

@@ -1,5 +1,6 @@
 """Filter signals with IIR filters."""
 
+import numpy as np
 from scipy.signal import butter, sosfiltfilt
 
 from neurodsp.utils import remove_nans, restore_nans
@@ -10,8 +11,8 @@
 ###################################################################################################
 ###################################################################################################
 
-def filter_signal_iir(sig, fs, pass_type, f_range, butterworth_order,
-                      print_transitions=False, plot_properties=False, return_filter=False):
+def filter_signal_iir(sig, fs, pass_type, f_range, butterworth_order, print_transitions=False,
+                       plot_properties=False, return_filter=False, verbose=False):
     """Apply an IIR filter to a signal.
 
     Parameters
@@ -41,6 +42,8 @@ def filter_signal_iir(sig, fs, pass_type, f_range, butterworth_order,
         If True, plot the properties of the filter, including frequency response and/or kernel.
     return_filter : bool, optional, default: False
         If True, return the second order series coefficients of the IIR filter.
+    verbose : bool, optional, default: False
+        If True, print out detailed filter information.
 
     Returns
     -------
@@ -65,7 +68,8 @@ def filter_signal_iir(sig, fs, pass_type, f_range, butterworth_order,
     sos = design_iir_filter(fs, pass_type, f_range, butterworth_order)
 
     # Check filter properties: compute transition bandwidth & run checks
-    check_filter_properties(sos, None, fs, pass_type, f_range, verbose=print_transitions)
+    check_filter_properties(sos, None, fs, pass_type, f_range, filt_type="IIR",
+                            verbose=np.any([print_transitions, verbose]))
 
     # Remove any NaN on the edges of 'sig'
     sig, sig_nans = remove_nans(sig)

neurodsp/filt/utils.py

@@ -136,7 +136,7 @@ def compute_pass_band(fs, pass_type, f_range):
     return pass_bw
 
 
-def compute_transition_band(f_db, db, low=-20, high=-3):
+def compute_transition_band(f_db, db, low=-20, high=-3, return_freqs=False):
     """Compute transition bandwidth of a filter.
 
     Parameters
@@ -149,11 +149,15 @@ def compute_transition_band(f_db, db, low=-20, high=-3):
         The lower limit that defines the transition band, in dB.
     high : float, optional, default: -3
         The upper limit that defines the transition band, in dB.
+    return_freqs : bool, optional, default: False
+        Whether to return a tuple of (lower, upper) frequency bounds for the transition band.
 
     Returns
     -------
     transition_band : float
         The transition bandwidth of the filter.
+    f_range : tuple of (float, float), optional, default: False
+        The lower and upper frequencies of the transition band.
 
     Examples
     --------
@@ -178,7 +182,15 @@ def compute_transition_band(f_db, db, low=-20, high=-3):
     # This gets the indices of transitions to the values in searched for range
     inds = np.where(np.diff(np.logical_and(db > low, db < high)))[0]
     # This steps through the indices, in pairs, selecting from the vector to select from
-    transition_band = np.max([(b - a) for a, b in zip(f_db[inds[0::2]], f_db[inds[1::2]])])
+    transition_pairs = [(a, b) for a, b in zip(f_db[inds[0::2]], f_db[inds[1::2]])]
+    pair_idx = np.argmax([(tran[1] - tran[0]) for tran in transition_pairs])
+    f_lo = transition_pairs[pair_idx][0]
+    f_hi = transition_pairs[pair_idx][1]
+    transition_band = f_hi - f_lo
+
+    if return_freqs:
+
+        return transition_band, (f_lo, f_hi)
 
     return transition_band