From 17593b7764991ba257c0b2e73c77f47b96976a5f Mon Sep 17 00:00:00 2001 From: ryanhammonds Date: Tue, 12 Jan 2021 19:38:52 -0800 Subject: [PATCH 1/9] filter report added --- neurodsp/filt/checks.py | 51 +++++++++++++++++++++++++++++++++++++---- neurodsp/filt/filter.py | 9 +++++--- neurodsp/filt/fir.py | 8 +++++-- neurodsp/filt/iir.py | 10 +++++--- neurodsp/filt/utils.py | 16 +++++++++++-- 5 files changed, 80 insertions(+), 14 deletions(-) diff --git a/neurodsp/filt/checks.py b/neurodsp/filt/checks.py index d6cef1b5..f1e76e56 100644 --- a/neurodsp/filt/checks.py +++ b/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 diff --git a/neurodsp/filt/filter.py b/neurodsp/filt/filter.py index 1b56f11c..4f716db3 100644 --- a/neurodsp/filt/filter.py +++ b/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.') diff --git a/neurodsp/filt/fir.py b/neurodsp/filt/fir.py index f3384015..6f0e1f29 100644 --- a/neurodsp/filt/fir.py +++ b/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) diff --git a/neurodsp/filt/iir.py b/neurodsp/filt/iir.py index 70161caf..aeddb6f7 100644 --- a/neurodsp/filt/iir.py +++ b/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) diff --git a/neurodsp/filt/utils.py b/neurodsp/filt/utils.py index 93629bfe..dcc3aef8 100644 --- a/neurodsp/filt/utils.py +++ b/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 From 6a5695ef881a9635538a840958c696b2c94c8ff5 Mon Sep 17 00:00:00 2001 From: ryanhammonds Date: Wed, 13 Jan 2021 17:59:20 -0800 Subject: [PATCH 2/9] save_properties of a filt --- neurodsp/filt/checks.py | 68 ++++++------------- neurodsp/filt/filter.py | 17 ++--- neurodsp/filt/fir.py | 9 ++- neurodsp/filt/iir.py | 8 ++- neurodsp/filt/utils.py | 105 +++++++++++++++++++++++++++++ neurodsp/tests/filt/test_checks.py | 7 +- neurodsp/tests/filt/test_utils.py | 38 ++++++++++- 7 files changed, 187 insertions(+), 65 deletions(-) diff --git a/neurodsp/filt/checks.py b/neurodsp/filt/checks.py index f1e76e56..29748896 100644 --- a/neurodsp/filt/checks.py +++ b/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), filt_type=None, 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 @@ -173,52 +179,16 @@ def check_filter_properties(filter_coefs, a_vals, fs, pass_type, f_range, 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('\n'.join('{} : {}'.format(key, value) for key, value in filt_report.items())) - # 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') + if save_properties is not None: + save_filt_report(save_properties, filt_report) return passes diff --git a/neurodsp/filt/filter.py b/neurodsp/filt/filter.py index 4f716db3..cc61229b 100644 --- a/neurodsp/filt/filter.py +++ b/neurodsp/filt/filter.py @@ -8,10 +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, - verbose=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 @@ -51,6 +50,8 @@ 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 @@ -76,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, verbose=verbose) + 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, verbose=verbose) + print_transitions, plot_properties, save_properties, + return_filter, verbose) else: raise ValueError('Filter type not understood.') diff --git a/neurodsp/filt/fir.py b/neurodsp/filt/fir.py index 6f0e1f29..6e2042ae 100644 --- a/neurodsp/filt/fir.py +++ b/neurodsp/filt/fir.py @@ -14,8 +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, - verbose=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 @@ -47,6 +47,8 @@ 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 @@ -82,7 +84,8 @@ def filter_signal_fir(sig, fs, pass_type, f_range, n_cycles=3, n_seconds=None, r # Check filter properties: compute transition bandwidth & run checks check_filter_properties(filter_coefs, 1, fs, pass_type, f_range, filt_type="FIR", - verbose=np.any([print_transitions, verbose])) + verbose=np.any([print_transitions, verbose]), + save_properties=save_properties) # Remove any NaN on the edges of 'sig' sig, sig_nans = remove_nans(sig) diff --git a/neurodsp/filt/iir.py b/neurodsp/filt/iir.py index aeddb6f7..08c8fb56 100644 --- a/neurodsp/filt/iir.py +++ b/neurodsp/filt/iir.py @@ -12,7 +12,8 @@ ################################################################################################### def filter_signal_iir(sig, fs, pass_type, f_range, butterworth_order, print_transitions=False, - plot_properties=False, return_filter=False, verbose=False): + plot_properties=False, save_properties=None, return_filter=False, + verbose=False): """Apply an IIR filter to a signal. Parameters @@ -40,6 +41,8 @@ def filter_signal_iir(sig, fs, pass_type, f_range, butterworth_order, print_tran 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 @@ -69,7 +72,8 @@ def filter_signal_iir(sig, fs, pass_type, f_range, butterworth_order, print_tran # Check filter properties: compute transition bandwidth & run checks check_filter_properties(sos, None, fs, pass_type, f_range, filt_type="IIR", - verbose=np.any([print_transitions, verbose])) + verbose=np.any([print_transitions, verbose]), + save_properties=save_properties) # Remove any NaN on the edges of 'sig' sig, sig_nans = remove_nans(sig) diff --git a/neurodsp/filt/utils.py b/neurodsp/filt/utils.py index dcc3aef8..6ba90e1a 100644 --- a/neurodsp/filt/utils.py +++ b/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 @@ -253,3 +256,105 @@ 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' and pass_type in ['bandstop', 'lowpass']: + + filt_report['Filter Class'] = '{filt_class}'.format(filt_class='linear-phase') + filt_report['Group Delay'] = '{delay}s'.format(delay=(len(f_db)-1) / 2 * fs) + + elif filt_type == 'FIR' and pass_type in ['bandpass', 'highpass']: + + filt_report['Filter Class'] = '{filt_class}'.format(filt_class='zero-phase') + filt_report['Group Delay'] = '0s' + + elif filt_type == 'IIR': + + # Group delay isn't reported for IIR since it varies from sample to sample + filt_report['Filter Class'] = '{filt_class}'.format(filt_class='non-linear-phase') + + # Direction of computation (one-pass forward/reverse, or two-pass forward and reverse) + if filt_type == 'FIR': + filt_report['Direction'] = 'one-pass reverse' + else: + 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) diff --git a/neurodsp/tests/filt/test_checks.py b/neurodsp/tests/filt/test_checks.py index 32b5776e..d5b07324 100644 --- a/neurodsp/tests/filt/test_checks.py +++ b/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(): diff --git a/neurodsp/tests/filt/test_utils.py b/neurodsp/tests/filt/test_utils.py index c839a75d..9b8f338d 100644 --- a/neurodsp/tests/filt/test_utils.py +++ b/neurodsp/tests/filt/test_utils.py @@ -1,8 +1,11 @@ """Tests for filter utilities.""" -from pytest import raises -from neurodsp.tests.settings import FS +import tempfile +from pytest import raises, mark, param + +import numpy as np +from neurodsp.tests.settings import FS from neurodsp.filt.utils import * from neurodsp.filt.fir import design_fir_filter, compute_filter_length @@ -53,3 +56,34 @@ def test_remove_filter_edges(): assert np.all(np.isnan(dropped_sig[:n_rmv])) assert np.all(np.isnan(dropped_sig[-n_rmv:])) assert np.all(~np.isnan(dropped_sig[n_rmv:-n_rmv])) + + +@mark.parametrize("pass_type", ['bandpass', 'bandstop', 'lowpass', 'highpass']) +@mark.parametrize("filt_type", ['IIR', 'FIR']) +def test_gen_filt_report(pass_type, filt_type): + + fs = 1000 + f_db = np.arange(0, 50) + db = np.random.rand(50) + pass_bw = 10 + transition_bw = 4 + f_range = (10, 40) + f_range_trans = (40, 44) + + report = gen_filt_report(pass_type, filt_type, fs, f_db, db, pass_bw, + transition_bw, f_range, f_range_trans) + + assert pass_type in report.values() + assert filt_type in report.values() + + +@mark.parametrize("dir_exists", [True, param(False, marks=mark.xfail)]) +def test_save_filt_report(dir_exists): + + filt_report = {'Pass Type': 'bandpass', 'Cutoff (half-amplitude)': 50} + temp_path = tempfile.NamedTemporaryFile() + if not dir_exists: + save_filt_report('/bad/path/', filt_report) + else: + save_filt_report(temp_path.name, filt_report) + temp_path.close() From 04e5de7c44cc2a4c83290121e7f3d2f725bf07f3 Mon Sep 17 00:00:00 2001 From: ryanhammonds Date: Thu, 14 Jan 2021 11:42:02 -0800 Subject: [PATCH 3/9] report update - all fir is linear-phase --- neurodsp/filt/utils.py | 25 ++++++++----------------- 1 file changed, 8 insertions(+), 17 deletions(-) diff --git a/neurodsp/filt/utils.py b/neurodsp/filt/utils.py index 6ba90e1a..30e7fa49 100644 --- a/neurodsp/filt/utils.py +++ b/neurodsp/filt/utils.py @@ -294,14 +294,14 @@ def gen_filt_report(pass_type, filt_type, fs, f_db, db, pass_bw, 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) + 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) + 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) + 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]]) @@ -312,25 +312,16 @@ def gen_filt_report(pass_type, filt_type, fs, f_db, db, pass_bw, # Filter delay (zero-phase, linear-phase, non-linear phase) filt_report['Filter Type'] = filt_type - if filt_type == 'FIR' and pass_type in ['bandstop', 'lowpass']: + if filt_type == 'FIR': - filt_report['Filter Class'] = '{filt_class}'.format(filt_class='linear-phase') + filt_report['Phase'] = '{filt_class}'.format(filt_class='linear-phase') filt_report['Group Delay'] = '{delay}s'.format(delay=(len(f_db)-1) / 2 * fs) - - elif filt_type == 'FIR' and pass_type in ['bandpass', 'highpass']: - - filt_report['Filter Class'] = '{filt_class}'.format(filt_class='zero-phase') - filt_report['Group Delay'] = '0s' + 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['Filter Class'] = '{filt_class}'.format(filt_class='non-linear-phase') - - # Direction of computation (one-pass forward/reverse, or two-pass forward and reverse) - if filt_type == 'FIR': - filt_report['Direction'] = 'one-pass reverse' - else: + filt_report['Phase'] = '{filt_class}'.format(filt_class='non-linear-phase') filt_report['Direction'] = 'two-pass forward and reverse' return filt_report From fdd109c4c32863669cc323c344474398d23aaabb Mon Sep 17 00:00:00 2001 From: ryanhammonds Date: Thu, 14 Jan 2021 17:44:00 -0800 Subject: [PATCH 4/9] group delay fix --- neurodsp/filt/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/neurodsp/filt/utils.py b/neurodsp/filt/utils.py index 30e7fa49..73c2205b 100644 --- a/neurodsp/filt/utils.py +++ b/neurodsp/filt/utils.py @@ -315,7 +315,7 @@ def gen_filt_report(pass_type, filt_type, fs, f_db, db, pass_bw, 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['Group Delay'] = '{delay}s'.format(delay=(len(f_db)-1) / (2 * fs)) filt_report['Direction'] = 'one-pass reverse' elif filt_type == 'IIR': From b98d6693bea253f1614170211f852168664e8384 Mon Sep 17 00:00:00 2001 From: ryanhammonds Date: Thu, 21 Jan 2021 14:27:11 -0800 Subject: [PATCH 5/9] review updates --- neurodsp/filt/checks.py | 35 +++---- neurodsp/filt/filter.py | 14 ++- neurodsp/filt/fir.py | 33 ++++--- neurodsp/filt/iir.py | 34 ++++--- neurodsp/filt/utils.py | 152 ++++++++++++++++++++++------- neurodsp/tests/filt/test_checks.py | 6 +- neurodsp/tests/filt/test_utils.py | 46 ++++++--- 7 files changed, 222 insertions(+), 98 deletions(-) diff --git a/neurodsp/filt/checks.py b/neurodsp/filt/checks.py index 29748896..8a189b14 100644 --- a/neurodsp/filt/checks.py +++ b/neurodsp/filt/checks.py @@ -1,8 +1,6 @@ """Checker functions for filtering.""" from warnings import warn -import os -import json import numpy as np @@ -92,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), - filt_type=None, verbose=True, save_properties=None): + return_properties=False, verbose=True): """Check a filters properties, including pass band and transition band. Parameters @@ -119,17 +117,18 @@ def check_filter_properties(filter_coefs, a_vals, fs, pass_type, f_range, transi 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. + return_properties : bool, optional, default: False + Returns the frequency response, pass band and transition band. verbose : bool, optional, default: True - Whether to print out filter properties. - save_properties : str - Path, including file name, to save filter properites to as a json. + Whether to print out transition and pass bands. Returns ------- passes : bool Whether all the checks pass. False if one or more checks fail. + properties : dict + The frequency response, pass band and transition band. + Only returned if return_properties is True. Examples -------- @@ -145,7 +144,7 @@ def check_filter_properties(filter_coefs, a_vals, fs, pass_type, f_range, transi # Import utility functions inside function to avoid circular imports from neurodsp.filt.utils import (compute_frequency_response, compute_pass_band, - compute_transition_band, gen_filt_report, save_filt_report) + compute_transition_band) # Initialize variable to keep track if all checks pass passes = True @@ -179,16 +178,18 @@ def check_filter_properties(filter_coefs, a_vals, fs, pass_type, f_range, transi warn('Transition bandwidth is {:.1f} Hz. This is greater than the desired'\ 'pass/stop bandwidth of {:.1f} Hz'.format(transition_bw, pass_bw)) - # 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) - + # Print out transition bandwidth and pass bandwidth to the user if verbose: - print('\n'.join('{} : {}'.format(key, value) for key, value in filt_report.items())) + print('Transition bandwidth is {:.1f} Hz.'.format(transition_bw)) + print('Pass/stop bandwidth is {:.1f} Hz.'.format(pass_bw)) + + # Return the filter properties + if return_properties: + + properties = {'f_db': f_db, 'db': db, 'pass_bw': pass_bw, 'transition_bw': transition_bw, + 'f_range_trans': f_range_trans} - if save_properties is not None: - save_filt_report(save_properties, filt_report) + return passes, properties return passes diff --git a/neurodsp/filt/filter.py b/neurodsp/filt/filter.py index cc61229b..685c5f8d 100644 --- a/neurodsp/filt/filter.py +++ b/neurodsp/filt/filter.py @@ -10,7 +10,7 @@ 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): + plot_properties=False, return_filter=False, save_report=None): """Apply a bandpass, bandstop, highpass, or lowpass filter to a neural signal. Parameters @@ -50,12 +50,10 @@ def filter_signal(sig, fs, pass_type, f_range, filter_type='fir', n_cycles=3, n_ 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. + save_report : str, optional, default: None + Path, including file name, to save a filter report to as a pdf. Returns ------- @@ -78,12 +76,12 @@ def filter_signal(sig, fs, pass_type, f_range, filter_type='fir', n_cycles=3, n_ 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, - save_properties, return_filter, verbose) + return_filter, save_report) 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, save_properties, - return_filter, verbose) + print_transitions, plot_properties, return_filter, + save_report) else: raise ValueError('Filter type not understood.') diff --git a/neurodsp/filt/fir.py b/neurodsp/filt/fir.py index 6e2042ae..e139e09e 100644 --- a/neurodsp/filt/fir.py +++ b/neurodsp/filt/fir.py @@ -6,7 +6,7 @@ from neurodsp.utils import remove_nans, restore_nans from neurodsp.utils.decorators import multidim from neurodsp.plts.filt import plot_filter_properties -from neurodsp.filt.utils import compute_frequency_response, remove_filter_edges +from neurodsp.filt.utils import compute_frequency_response, remove_filter_edges, save_filt_report from neurodsp.filt.checks import (check_filter_definition, check_filter_properties, check_filter_length) @@ -14,8 +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, save_properties=None, - return_filter=False, verbose=False, file_path=None, file_name=None): + print_transitions=False, plot_properties=False, return_filter=False, + save_report=None): """Apply an FIR filter to a signal. Parameters @@ -47,12 +47,10 @@ 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. + save_report : str, optional, default: None + Path, including file name, to save a filter report to as a pdf. Returns ------- @@ -83,9 +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, filt_type="FIR", - verbose=np.any([print_transitions, verbose]), - save_properties=save_properties) + _, properties = check_filter_properties(filter_coefs, 1, fs, pass_type, f_range, + return_properties=True, verbose=print_transitions) # Remove any NaN on the edges of 'sig' sig, sig_nans = remove_nans(sig) @@ -100,11 +97,25 @@ def filter_signal_fir(sig, fs, pass_type, f_range, n_cycles=3, n_seconds=None, r # Add NaN back on the edges of 'sig', if there were any at the beginning sig_filt = restore_nans(sig_filt, sig_nans) + # Unpack filter properties + if plot_properties or save_report is not None: + f_db = properties['f_db'] + db = properties['db'] + + if save_report is not None: + pass_bw = properties['pass_bw'] + transition_bw = properties['transition_bw'] + f_range_trans = properties['f_range_trans'] + # Plot filter properties, if specified if plot_properties: - f_db, db = compute_frequency_response(filter_coefs, 1, fs) plot_filter_properties(f_db, db, fs, filter_coefs) + # Save a pdf filter report containing plots and parameters + if save_report is not None: + save_filt_report(save_report, pass_type, 'IIR', fs, f_db, db, pass_bw, transition_bw, + f_range, f_range_trans, len(f_db)-1, filter_coefs=filter_coefs) + if return_filter: return sig_filt, filter_coefs else: diff --git a/neurodsp/filt/iir.py b/neurodsp/filt/iir.py index 08c8fb56..a34b07a8 100644 --- a/neurodsp/filt/iir.py +++ b/neurodsp/filt/iir.py @@ -1,10 +1,9 @@ """Filter signals with IIR filters.""" -import numpy as np from scipy.signal import butter, sosfiltfilt from neurodsp.utils import remove_nans, restore_nans -from neurodsp.filt.utils import compute_nyquist, compute_frequency_response +from neurodsp.filt.utils import compute_nyquist, compute_frequency_response, save_filt_report from neurodsp.filt.checks import check_filter_definition, check_filter_properties from neurodsp.plts.filt import plot_frequency_response @@ -12,8 +11,7 @@ ################################################################################################### 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): + plot_properties=False, return_filter=False, save_report=None): """Apply an IIR filter to a signal. Parameters @@ -41,12 +39,10 @@ def filter_signal_iir(sig, fs, pass_type, f_range, butterworth_order, print_tran 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. + save_report : str, optional, default: None + Path, including file name, to save a filter report to as a pdf. Returns ------- @@ -71,9 +67,8 @@ def filter_signal_iir(sig, fs, pass_type, f_range, butterworth_order, print_tran 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, filt_type="IIR", - verbose=np.any([print_transitions, verbose]), - save_properties=save_properties) + _, properties = check_filter_properties(sos, None, fs, pass_type, f_range, + return_properties=True, verbose=print_transitions) # Remove any NaN on the edges of 'sig' sig, sig_nans = remove_nans(sig) @@ -84,11 +79,26 @@ def filter_signal_iir(sig, fs, pass_type, f_range, butterworth_order, print_tran # Add NaN back on the edges of 'sig', if there were any at the beginning sig_filt = restore_nans(sig_filt, sig_nans) + # Unpack filter properties + if plot_properties or save_report is not None: + f_db = properties['f_db'] + db = properties['db'] + + if save_report is not None: + pass_bw = properties['pass_bw'] + transition_bw = properties['transition_bw'] + f_range_trans = properties['f_range_trans'] + # Plot frequency response, if desired if plot_properties: - f_db, db = compute_frequency_response(sos, None, fs) plot_frequency_response(f_db, db) + # Save a pdf filter report containing plots and parameters + if save_report is not None: + + save_filt_report(save_report, pass_type, 'IIR', fs, f_db, db, pass_bw, + transition_bw, f_range, f_range_trans, butterworth_order) + if return_filter: return sig_filt, sos else: diff --git a/neurodsp/filt/utils.py b/neurodsp/filt/utils.py index 73c2205b..088472ec 100644 --- a/neurodsp/filt/utils.py +++ b/neurodsp/filt/utils.py @@ -1,13 +1,14 @@ """Utility functions for filtering.""" import os -import json +import matplotlib.pyplot as plt import numpy as np from scipy.signal import freqz, sosfreqz from neurodsp.utils.decorators import multidim from neurodsp.filt.checks import check_filter_definition +from neurodsp.plts.filt import plot_frequency_response, plot_impulse_response ################################################################################################### ################################################################################################### @@ -159,8 +160,9 @@ def compute_transition_band(f_db, db, low=-20, high=-3, return_freqs=False): ------- transition_band : float The transition bandwidth of the filter. - f_range : tuple of (float, float), optional, default: False + f_range : tuple of (float, float) The lower and upper frequencies of the transition band. + Only returned is return_freqs is True. Examples -------- @@ -184,7 +186,8 @@ def compute_transition_band(f_db, db, low=-20, high=-3, return_freqs=False): # 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 + + # This determines at which frequencies the transition band occurs 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] @@ -258,8 +261,8 @@ def remove_filter_edges(sig, filt_len): return sig -def gen_filt_report(pass_type, filt_type, fs, f_db, db, pass_bw, - transition_bw, f_range, f_range_trans): +def gen_filt_str(pass_type, filt_type, fs, f_db, db, pass_bw, + transition_bw, f_range, f_range_trans, order): """Create a filter report. Parameters @@ -282,70 +285,149 @@ def gen_filt_report(pass_type, filt_type, fs, f_db, db, pass_bw, 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. + order : int + The filter length for FIR filter or butterworth order for IIR filters. Returns ------- - filt_report : dict - A dicionary of filter parameter keys and corresponding values. + filt_str : str + Filter properties as a string that is ready to embed into a pdf report. """ - filt_report = {} + + filt_str = [] # Filter type (high-pass, low-pass, band-pass, band-stop, FIR, IIR) - filt_report['Pass Type'] = '{pass_type}'.format(pass_type=pass_type) + filt_str.append('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) + filt_str.append('Cutoff (Half-Amplitude): {cutoff} Hz'.format(cutoff=f_range)) - # Filter order (or length) - filt_report['Filter Order'] = '{order}'.format(order=len(f_db)-1) + # Filter order (or length-1) for FIR or butterworth order for IIR + filt_str.append('Filter Order: {order}'.format(order=order)) # Roll-off or transition bandwidth - filt_report['Transition Bandwidth'] = '{:.1f} Hz'.format(transition_bw) - filt_report['Pass/Stop Bandwidth'] = '{:.1f} Hz'.format(pass_bw) + filt_str.append('Transition Bandwidth: {:.1f} Hz'.format(transition_bw)) + filt_str.append('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) + filt_str.append('Passband Ripple: {:1.4f} db'.format(pb_ripple)) + filt_str.append('Stopband Attenuation: {:1.4f} db'.format(sb_atten)) # Filter delay (zero-phase, linear-phase, non-linear phase) - filt_report['Filter Type'] = filt_type + filt_str.append('Filter Type: {filt_type}'.format(filt_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' + filt_str.append('Phase: linear-phase') + filt_str.append('Group Delay: {:1.4f} s'.format((len(f_db)-1) / (2 * fs))) + filt_str.append('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' + filt_str.append('Phase: non-linear-phase') + filt_str.append('Direction: two-pass forward and reverse') + + # Format the list into a string + filt_str = [ + + # Header + '=', + '', + 'FILTER REPORT', + '', - return filt_report + # Settings + *filt_str, + # Footer + '', + '=' + ] -def save_filt_report(save_properties, filt_report): + str_len = 50 + filt_str [0] = filt_str [0] * str_len + filt_str [-1] = filt_str [-1] * str_len + + filt_str = '\n'.join([string.center(str_len) for string in filt_str]) + + return filt_str + + +def save_filt_report(pdf_path, pass_type, filt_type, fs, f_db, db, pass_bw, transition_bw, + f_range, f_range_trans, order, filter_coefs=None): """Save filter properties as a json file. - Parameters + Parameters ---------- - save_properties : str - Path, including file name, to save filter properites to as a json. - filt_report : dict - Contains filter report info. + pdf_path: str + Path, including file name, to save a filter report to as a pdf. + 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. + order : int + The filter length for FIR filter or butterworth order for IIR filters. + filter_coefs : 1d array, optional, default: None + Filter coefficients of the FIR filter. """ - # Ensure parents exists - if not os.path.isdir(os.path.dirname(save_properties)): + # Ensure valid path + if not pdf_path.startswith('/') and not pdf_path.startswith('./'): + pdf_path = './' + pdf_path + + if not os.path.isdir(os.path.dirname(pdf_path)): + print(pdf_path) 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' + if not pdf_path.endswith('.pdf'): + pdf_path = pdf_path + '.pdf' + + # Create properties string + filt_str = gen_filt_str(pass_type, filt_type, fs, f_db, db, pass_bw, + transition_bw, f_range, f_range_trans, order) + + # Plot + if filter_coefs is not None: + + _, axes = plt.subplots(nrows=3, ncols=1, figsize=(8, 18), + gridspec_kw={'height_ratios': [1, 4, 4]}) + + # Plot impulse response for IIR filters + plot_impulse_response(fs, filter_coefs, ax=axes[2]) + + else: + + _, axes = plt.subplots(nrows=2, ncols=1, figsize=(8, 10), + gridspec_kw={'height_ratios': [1, 4]}) + + # Plot filter parameter string + font = {'family': 'monospace', 'weight': 'normal', 'size': 16} + axes[0].text(0.5, 0.7, filt_str, font, ha='center', va='center') + axes[0].set_frame_on(False) + axes[0].set_xticks([]) + axes[0].set_yticks([]) + + # Plot filter responses + plot_frequency_response(f_db, db, ax=axes[1]) # Save - with open(save_properties, 'w') as file_path: - json.dump(filt_report, file_path) + plt.savefig(pdf_path) + plt.close() diff --git a/neurodsp/tests/filt/test_checks.py b/neurodsp/tests/filt/test_checks.py index d5b07324..c47c639e 100644 --- a/neurodsp/tests/filt/test_checks.py +++ b/neurodsp/tests/filt/test_checks.py @@ -1,6 +1,5 @@ """Tests for filter check functions.""" -import tempfile from pytest import raises from neurodsp.tests.settings import FS @@ -59,10 +58,9 @@ 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() + verbose=True, return_properties=True) def test_check_filter_length(): diff --git a/neurodsp/tests/filt/test_utils.py b/neurodsp/tests/filt/test_utils.py index 9b8f338d..b89d9c2c 100644 --- a/neurodsp/tests/filt/test_utils.py +++ b/neurodsp/tests/filt/test_utils.py @@ -5,9 +5,11 @@ import numpy as np -from neurodsp.tests.settings import FS +from neurodsp.tests.settings import FS, FS_HIGH from neurodsp.filt.utils import * from neurodsp.filt.fir import design_fir_filter, compute_filter_length +from neurodsp.filt.iir import design_iir_filter +from neurodsp.filt.checks import check_filter_definition, check_filter_properties ################################################################################################### ################################################################################################### @@ -60,30 +62,52 @@ def test_remove_filter_edges(): @mark.parametrize("pass_type", ['bandpass', 'bandstop', 'lowpass', 'highpass']) @mark.parametrize("filt_type", ['IIR', 'FIR']) -def test_gen_filt_report(pass_type, filt_type): +def test_gen_filt_str(pass_type, filt_type): - fs = 1000 f_db = np.arange(0, 50) db = np.random.rand(50) pass_bw = 10 transition_bw = 4 f_range = (10, 40) f_range_trans = (40, 44) + order = 1 - report = gen_filt_report(pass_type, filt_type, fs, f_db, db, pass_bw, - transition_bw, f_range, f_range_trans) + report_str = gen_filt_str(pass_type, filt_type, FS, f_db, db, pass_bw, + transition_bw, f_range, f_range_trans, order) - assert pass_type in report.values() - assert filt_type in report.values() + assert pass_type in report_str + assert filt_type in report_str @mark.parametrize("dir_exists", [True, param(False, marks=mark.xfail)]) -def test_save_filt_report(dir_exists): +@mark.parametrize("filt_type", ['IIR', 'FIR']) +def test_save_filt_report(dir_exists, filt_type): + + pass_type = 'bandpass' + f_range = (10, 40) + + f_db = np.arange(1, 100) + db = np.random.rand(99) + + pass_bw = 10 + transition_bw = 4 + f_range_trans = (40, 44) + + order = 1 + + if pass_type == 'FIR': + filter_coefs = np.random.rand(10) + else: + filter_coefs = None - filt_report = {'Pass Type': 'bandpass', 'Cutoff (half-amplitude)': 50} temp_path = tempfile.NamedTemporaryFile() + if not dir_exists: - save_filt_report('/bad/path/', filt_report) + save_filt_report('/bad/path/', pass_type, filt_type, FS_HIGH, f_db, db, pass_bw, + transition_bw, f_range, f_range_trans, order, filter_coefs=filter_coefs) else: - save_filt_report(temp_path.name, filt_report) + print(temp_path.name) + save_filt_report(temp_path.name, pass_type, filt_type, FS_HIGH, f_db, db, pass_bw, + transition_bw, f_range, f_range_trans, order, filter_coefs=filter_coefs) + temp_path.close() From fcc4d93990bafea2e389325a5a1b8284a45bbbe4 Mon Sep 17 00:00:00 2001 From: ryanhammonds Date: Thu, 21 Jan 2021 14:51:16 -0800 Subject: [PATCH 6/9] increased test coverage --- neurodsp/tests/filt/test_fir.py | 8 +++++++- neurodsp/tests/filt/test_iir.py | 9 ++++++++- 2 files changed, 15 insertions(+), 2 deletions(-) diff --git a/neurodsp/tests/filt/test_fir.py b/neurodsp/tests/filt/test_fir.py index 55fd6488..e7949b3e 100644 --- a/neurodsp/tests/filt/test_fir.py +++ b/neurodsp/tests/filt/test_fir.py @@ -1,5 +1,6 @@ """Tests for FIR filters.""" +import tempfile from pytest import raises import numpy as np @@ -12,7 +13,12 @@ def test_filter_signal_fir(tsig, tsig_sine): - out = filter_signal_fir(tsig, FS, 'bandpass', (8, 12)) + temp_path = tempfile.NamedTemporaryFile() + + out = filter_signal_fir(tsig, FS, 'bandpass', (8, 12), save_report=temp_path.name) + + temp_path.close() + assert out.shape == tsig.shape # Apply lowpass to low-frequency sine. There should be little attenuation. diff --git a/neurodsp/tests/filt/test_iir.py b/neurodsp/tests/filt/test_iir.py index 4b209035..005ba39e 100644 --- a/neurodsp/tests/filt/test_iir.py +++ b/neurodsp/tests/filt/test_iir.py @@ -1,5 +1,7 @@ """Tests for IIR filters.""" +import tempfile + import numpy as np from neurodsp.tests.settings import FS @@ -11,7 +13,12 @@ def test_filter_signal_iir(tsig): - out = filter_signal_iir(tsig, FS, 'bandpass', (8, 12), 3) + temp_path = tempfile.NamedTemporaryFile() + + out = filter_signal_iir(tsig, FS, 'bandpass', (8, 12), 3, save_report=temp_path.name) + + temp_path.close() + assert out.shape == tsig.shape def test_filter_signal_iir_2d(tsig2d): From 1607d165373c9fd3eb657a50e9d732e83a1af251 Mon Sep 17 00:00:00 2001 From: ryanhammonds Date: Tue, 26 Jan 2021 16:17:11 -0800 Subject: [PATCH 7/9] remove unneeded print --- neurodsp/filt/utils.py | 1 - 1 file changed, 1 deletion(-) diff --git a/neurodsp/filt/utils.py b/neurodsp/filt/utils.py index 088472ec..6e5394dc 100644 --- a/neurodsp/filt/utils.py +++ b/neurodsp/filt/utils.py @@ -393,7 +393,6 @@ def save_filt_report(pdf_path, pass_type, filt_type, fs, f_db, db, pass_bw, tra pdf_path = './' + pdf_path if not os.path.isdir(os.path.dirname(pdf_path)): - print(pdf_path) raise ValueError("Unable to save properties. Parent directory does not exist.") # Enforce file extension From 59cd8768ae85adfc950fdbabd88d9587907b2af0 Mon Sep 17 00:00:00 2001 From: ryanhammonds Date: Thu, 3 Jun 2021 13:44:03 -0700 Subject: [PATCH 8/9] FIR direction updated to one-pass --- neurodsp/filt/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/neurodsp/filt/utils.py b/neurodsp/filt/utils.py index 6e5394dc..eaae4c53 100644 --- a/neurodsp/filt/utils.py +++ b/neurodsp/filt/utils.py @@ -322,7 +322,7 @@ def gen_filt_str(pass_type, filt_type, fs, f_db, db, pass_bw, filt_str.append('Phase: linear-phase') filt_str.append('Group Delay: {:1.4f} s'.format((len(f_db)-1) / (2 * fs))) - filt_str.append('Direction: one-pass reverse') + filt_str.append('Direction: one-pass') elif filt_type == 'IIR': From 993893c7c5b85ed7791c706ec04478cd7966e40e Mon Sep 17 00:00:00 2001 From: ryanhammonds Date: Thu, 3 Jun 2021 13:52:11 -0700 Subject: [PATCH 9/9] FIR group delay to zero --- neurodsp/filt/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/neurodsp/filt/utils.py b/neurodsp/filt/utils.py index eaae4c53..825ac28b 100644 --- a/neurodsp/filt/utils.py +++ b/neurodsp/filt/utils.py @@ -321,7 +321,7 @@ def gen_filt_str(pass_type, filt_type, fs, f_db, db, pass_bw, if filt_type == 'FIR': filt_str.append('Phase: linear-phase') - filt_str.append('Group Delay: {:1.4f} s'.format((len(f_db)-1) / (2 * fs))) + filt_str.append('Group Delay: 0s') filt_str.append('Direction: one-pass') elif filt_type == 'IIR':