[MRG] Filter Report

neurodsp/filt/checks.py

@@ -1,6 +1,8 @@
 """Checker functions for filtering."""
 
 from warnings import warn
+import os
+import json
 
 import numpy as np
 
@@ -89,8 +91,8 @@ def check_filter_definition(pass_type, f_range):
     return f_lo, f_hi
 
 
-def check_filter_properties(filter_coefs, a_vals, fs, pass_type, f_range,
-                            transitions=(-20, -3), verbose=True):
+def check_filter_properties(filter_coefs, a_vals, fs, pass_type, f_range, transitions=(-20, -3),
+                            filt_type=None, verbose=True, save_properties=None):
     """Check a filters properties, including pass band and transition band.
 
     Parameters
@@ -117,8 +119,12 @@ def check_filter_properties(filter_coefs, a_vals, fs, pass_type, f_range,
         a tuple and is assumed to be (None, f_hi) for 'lowpass', and (f_lo, None) for 'highpass'.
     transitions : tuple of (float, float), optional, default: (-20, -3)
         Cutoffs, in dB, that define the transition band.
+    filt_type : str, optional, {'FIR', 'IIR'}
+        The type of filter being applied.
     verbose : bool, optional, default: True
-        Whether to print out transition and pass bands.
+        Whether to print out filter properties.
+    save_properties : str
+        Path, including file name, to save filter properites to as a json.
 
     Returns
     -------
@@ -138,8 +144,8 @@ def check_filter_properties(filter_coefs, a_vals, fs, pass_type, f_range,
     """
 
     # Import utility functions inside function to avoid circular imports
-    from neurodsp.filt.utils import (compute_frequency_response,
-                                     compute_pass_band, compute_transition_band)
+    from neurodsp.filt.utils import (compute_frequency_response, compute_pass_band,
+                                     compute_transition_band, gen_filt_report, save_filt_report)
 
     # Initialize variable to keep track if all checks pass
     passes = True
@@ -164,18 +170,25 @@ 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:
         passes = False
         warn('Transition bandwidth is  {:.1f}  Hz. This is greater than the desired'\
              'pass/stop bandwidth of  {:.1f} Hz'.format(transition_bw, pass_bw))
 
-    # Print out transition bandwidth and pass bandwidth to the user
+    # Report filter properties
+    if verbose or save_properties:
+        filt_report = gen_filt_report(pass_type, filt_type, fs, f_db, db, pass_bw,
+                                      transition_bw, f_range, f_range_trans)
+
     if verbose:
-        print('Transition bandwidth is {:.1f} Hz.'.format(transition_bw))
-        print('Pass/stop bandwidth is {:.1f} Hz.'.format(pass_bw))
+        print('\n'.join('{} : {}'.format(key, value) for key, value in filt_report.items()))
+
+    if save_properties is not None:
+        save_filt_report(save_properties, filt_report)
 
     return passes
 

neurodsp/filt/filter.py

@@ -8,9 +8,9 @@
 ###################################################################################################
 ###################################################################################################
 
-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):
+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, save_properties=None, return_filter=False, verbose=False):
     """Apply a bandpass, bandstop, highpass, or lowpass filter to a neural signal.
 
     Parameters
@@ -50,8 +50,12 @@ def filter_signal(sig, fs, pass_type, f_range, filter_type='fir',
         If True, print out the transition and pass bandwidths.
     plot_properties : bool, optional, default: False
         If True, plot the properties of the filter, including frequency response and/or kernel.
+    save_properties : str, optional, default: None
+        Path, including file name, to save filter properites to as a json.
     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
     -------
@@ -73,13 +77,13 @@ 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)
+                                 remove_edges, print_transitions, plot_properties,
+                                 save_properties, return_filter, 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)
+                                 print_transitions, plot_properties, save_properties,
+                                 return_filter, 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, save_properties=None,
+                      return_filter=False, verbose=False, file_path=None, file_name=None):
     """Apply an FIR filter to a signal.
 
     Parameters
@@ -46,8 +47,12 @@ def filter_signal_fir(sig, fs, pass_type, f_range, n_cycles=3, n_seconds=None, r
         If True, print out the transition and pass bandwidths.
     plot_properties : bool, optional, default: False
         If True, plot the properties of the filter, including frequency response and/or kernel.
+    save_properties : str
+        Path, including file name, to save filter properites to as a json.
     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 +83,9 @@ 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]),
+                            save_properties=save_properties)
 
     # 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,9 @@
 ###################################################################################################
 ###################################################################################################
 
-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, save_properties=None, return_filter=False,
+                      verbose=False):
     """Apply an IIR filter to a signal.
 
     Parameters
@@ -39,8 +41,12 @@ def filter_signal_iir(sig, fs, pass_type, f_range, butterworth_order,
         If True, print out the transition and pass bandwidths.
     plot_properties : bool, optional, default: False
         If True, plot the properties of the filter, including frequency response and/or kernel.
+    save_properties : str
+        Path, including file name, to save filter properites to as a json.
     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 +71,9 @@ 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]),
+                            save_properties=save_properties)
 
     # Remove any NaN on the edges of 'sig'
     sig, sig_nans = remove_nans(sig)

neurodsp/filt/utils.py

@@ -1,5 +1,8 @@
 """Utility functions for filtering."""
 
+import os
+import json
+
 import numpy as np
 from scipy.signal import freqz, sosfreqz
 
@@ -136,7 +139,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 +152,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 +185,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
 
@@ -241,3 +256,96 @@ def remove_filter_edges(sig, filt_len):
     sig[-n_rmv:] = np.nan
 
     return sig
+
+
+def gen_filt_report(pass_type, filt_type, fs, f_db, db, pass_bw,
+                    transition_bw, f_range, f_range_trans):
+    """Create a filter report.
+
+    Parameters
+    ----------
+    pass_type : {'bandpass', 'bandstop', 'lowpass', 'highpass'}
+        Which type of filter was applied.
+    filt_type : str, {'FIR', 'IIR'}
+        The type of filter being applied.
+    fs : float
+        Sampling rate, in Hz.
+    f_db : 1d array
+        Frequency vector corresponding to attenuation decibels, in Hz.
+    db : 1d array
+        Degree of attenuation for each frequency specified in `f_db`, in dB.
+    pass_bw : float
+        The pass bandwidth of the filter.
+    transition_band : float
+        The transition bandwidth of the filter.
+    f_range : tuple of (float, float) or float
+        Cutoff frequency(ies) used for filter, specified as f_lo & f_hi.
+    f_range_trans : tuple of (float, float)
+        The lower and upper frequencies of the transition band.
+
+    Returns
+    -------
+    filt_report : dict
+        A dicionary of filter parameter keys and corresponding values.
+    """
+    filt_report = {}
+
+    # Filter type (high-pass, low-pass, band-pass, band-stop, FIR, IIR)
+    filt_report['Pass Type'] = '{pass_type}'.format(pass_type=pass_type)
+
+    # Cutoff frequenc(ies) (including definition)
+    filt_report['Cutoff (Half-Amplitude)'] = '{cutoff} Hz'.format(cutoff=f_range)
+
+    # Filter order (or length)
+    filt_report['Filter Order'] = '{order}'.format(order=len(f_db)-1)
+
+    # Roll-off or transition bandwidth
+    filt_report['Transition Bandwidth'] = '{:.1f} Hz'.format(transition_bw)
+    filt_report['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]:])
+    filt_report['Passband Ripple'] = '{pb_ripple} db'.format(pb_ripple=pb_ripple)
+    filt_report['Stopband Attenuation'] = '{sb_atten} db'.format(sb_atten=sb_atten)
+
+    # Filter delay (zero-phase, linear-phase, non-linear phase)
+    filt_report['Filter Type'] = filt_type
+
+    if filt_type == 'FIR':
+
+        filt_report['Phase'] = '{filt_class}'.format(filt_class='linear-phase')
+        filt_report['Group Delay'] = '{delay}s'.format(delay=(len(f_db)-1) / (2 * fs))
+        filt_report['Direction'] = 'one-pass reverse'
+
+    elif filt_type == 'IIR':
+
+        # Group delay isn't reported for IIR since it varies from sample to sample
+        filt_report['Phase'] = '{filt_class}'.format(filt_class='non-linear-phase')
+        filt_report['Direction'] = 'two-pass forward and reverse'
+
+    return filt_report
+
+
+def save_filt_report(save_properties, filt_report):
+    """Save filter properties as a json file.
+
+    Parameters
+    ----------
+    save_properties : str
+        Path, including file name, to save filter properites to as a json.
+    filt_report : dict
+        Contains filter report info.
+    """
+
+    # Ensure parents exists
+    if not os.path.isdir(os.path.dirname(save_properties)):
+        raise ValueError("Unable to save properties. Parent directory does not exist.")
+
+    # Enforce file extension
+    if not save_properties.endswith('.json'):
+        save_properties = save_properties + '.json'
+
+    # Save
+    with open(save_properties, 'w') as file_path:
+        json.dump(filt_report, file_path)

neurodsp/tests/filt/test_checks.py

@@ -1,5 +1,6 @@
 """Tests for filter check functions."""
 
+import tempfile
 from pytest import raises
 
 from neurodsp.tests.settings import FS
@@ -46,7 +47,7 @@ def test_check_filter_definition():
 def test_check_filter_properties():
 
     filter_coefs = design_fir_filter(FS, 'bandpass', (8, 12))
-    
+
     passes = check_filter_properties(filter_coefs, 1, FS, 'bandpass', (8, 12))
     assert passes is True
 
@@ -58,6 +59,10 @@ def test_check_filter_properties():
     passes = check_filter_properties(filter_coefs, 1, FS, 'bandpass', (8, 12))
     assert passes is False
 
+    temp_path = tempfile.NamedTemporaryFile()
+    check_filter_properties(filter_coefs, 1, FS, 'bandpass', (8, 12),
+                            verbose=True, save_properties=temp_path.name)
+    temp_path.close()
 
 def test_check_filter_length():