diff --git a/neurodsp/filt/checks.py b/neurodsp/filt/checks.py index 9dfcefa8..3b74b2d3 100644 --- a/neurodsp/filt/checks.py +++ b/neurodsp/filt/checks.py @@ -89,8 +89,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), + return_properties=False, verbose=True): """Check a filters properties, including pass band and transition band. Parameters @@ -117,6 +117,8 @@ 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. + return_properties : bool, optional, default: False + Returns the frequency response, pass band and transition band. verbose : bool, optional, default: True Whether to print out transition and pass bands. @@ -124,6 +126,9 @@ def check_filter_properties(filter_coefs, a_vals, fs, pass_type, f_range, ------- 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 -------- @@ -138,8 +143,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) # Initialize variable to keep track if all checks pass passes = True @@ -164,7 +169,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: @@ -177,6 +183,14 @@ def check_filter_properties(filter_coefs, a_vals, fs, pass_type, f_range, 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} + + return passes, properties + return passes diff --git a/neurodsp/filt/filter.py b/neurodsp/filt/filter.py index 1b56f11c..685c5f8d 100644 --- a/neurodsp/filt/filter.py +++ b/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, return_filter=False, save_report=None): """Apply a bandpass, bandstop, highpass, or lowpass filter to a neural signal. Parameters @@ -52,6 +52,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. + save_report : str, optional, default: None + Path, including file name, to save a filter report to as a pdf. Returns ------- @@ -73,13 +75,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, + 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, - return_filter) + 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 f3384015..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,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, + save_report=None): """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. + save_report : str, optional, default: None + Path, including file name, to save a filter report to as a pdf. 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) + _, 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) @@ -93,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 70161caf..a34b07a8 100644 --- a/neurodsp/filt/iir.py +++ b/neurodsp/filt/iir.py @@ -3,15 +3,15 @@ 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 ################################################################################################### ################################################################################################### -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, save_report=None): """Apply an IIR filter to a signal. Parameters @@ -41,6 +41,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. + save_report : str, optional, default: None + Path, including file name, to save a filter report to as a pdf. Returns ------- @@ -65,7 +67,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) + _, 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) @@ -76,11 +79,26 @@ def filter_signal_iir(sig, fs, pass_type, f_range, butterworth_order, # 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 a4355d0a..46810dc5 100644 --- a/neurodsp/filt/utils.py +++ b/neurodsp/filt/utils.py @@ -1,10 +1,14 @@ """Utility functions for filtering.""" +import os + +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 ################################################################################################### ################################################################################################### @@ -136,7 +140,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 +153,16 @@ 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, in Hz. + The transition bandwidth of the filter. + f_range : tuple of (float, float) + The lower and upper frequencies of the transition band. + Only returned is return_freqs is True. Examples -------- @@ -177,8 +186,17 @@ 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]])]) + + # 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] + 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 +259,174 @@ def remove_filter_edges(sig, filt_len): sig[-n_rmv:] = np.nan return sig + + +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 + ---------- + 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. + + Returns + ------- + filt_str : str + Filter properties as a string that is ready to embed into a pdf report. + """ + + filt_str = [] + + # Filter type (high-pass, low-pass, band-pass, band-stop, FIR, IIR) + filt_str.append('Pass Type: {pass_type}'.format(pass_type=pass_type)) + + # Cutoff frequenc(ies) (including definition) + filt_str.append('Cutoff (Half-Amplitude): {cutoff} Hz'.format(cutoff=f_range)) + + # 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_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_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_str.append('Filter Type: {filt_type}'.format(filt_type=filt_type)) + + if filt_type == 'FIR': + + filt_str.append('Phase: linear-phase') + filt_str.append('Group Delay: 0s') + filt_str.append('Direction: one-pass') + + elif filt_type == 'IIR': + + # Group delay isn't reported for IIR since it varies from sample to sample + 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', + '', + + # Settings + *filt_str, + + # Footer + '', + '=' + ] + + 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 + ---------- + 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 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)): + raise ValueError("Unable to save properties. Parent directory does not exist.") + + # Enforce file extension + 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 + plt.savefig(pdf_path) + plt.close() diff --git a/neurodsp/tests/filt/test_checks.py b/neurodsp/tests/filt/test_checks.py index 5c9d260b..77720978 100644 --- a/neurodsp/tests/filt/test_checks.py +++ b/neurodsp/tests/filt/test_checks.py @@ -58,6 +58,9 @@ def test_check_filter_properties(): assert passes is False + check_filter_properties(filter_coefs, 1, FS, 'bandpass', (8, 12), + verbose=True, return_properties=True) + def test_check_filter_length(): check_filter_length(1000, 500) diff --git a/neurodsp/tests/filt/test_fir.py b/neurodsp/tests/filt/test_fir.py index 9da9a83f..c9e94273 100644 --- a/neurodsp/tests/filt/test_fir.py +++ b/neurodsp/tests/filt/test_fir.py @@ -1,5 +1,6 @@ """Tests for neurodsp.filt.fir.""" +import tempfile from pytest import raises import numpy as np @@ -13,7 +14,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, which should should give little attenuation diff --git a/neurodsp/tests/filt/test_iir.py b/neurodsp/tests/filt/test_iir.py index ac17c438..0c39a07e 100644 --- a/neurodsp/tests/filt/test_iir.py +++ b/neurodsp/tests/filt/test_iir.py @@ -1,5 +1,7 @@ """Tests for neurodsp.filt.iir.""" +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): diff --git a/neurodsp/tests/filt/test_utils.py b/neurodsp/tests/filt/test_utils.py index 44539d10..882779e1 100644 --- a/neurodsp/tests/filt/test_utils.py +++ b/neurodsp/tests/filt/test_utils.py @@ -1,10 +1,15 @@ """Tests for neurodsp.filt.utils.""" -from pytest import raises +import tempfile +from pytest import raises, mark, param -from neurodsp.tests.settings import FS +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 +from neurodsp.filt.iir import design_iir_filter +from neurodsp.filt.checks import check_filter_definition, check_filter_properties from neurodsp.filt.utils import * @@ -55,3 +60,56 @@ 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_str(pass_type, filt_type): + + 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_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_str + assert filt_type in report_str + + +@mark.parametrize("dir_exists", [True, param(False, marks=mark.xfail)]) +@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 + + temp_path = tempfile.NamedTemporaryFile() + + if not dir_exists: + save_filt_report('/bad/path/', pass_type, filt_type, FS, f_db, db, pass_bw, + transition_bw, f_range, f_range_trans, order, filter_coefs=filter_coefs) + else: + print(temp_path.name) + save_filt_report(temp_path.name, pass_type, filt_type, FS, f_db, db, pass_bw, + transition_bw, f_range, f_range_trans, order, filter_coefs=filter_coefs) + + temp_path.close()