[BUG] - Technical Audit: Fix numerical instabilities for IIR filters

neurodsp/filt/checks.py

@@ -133,26 +133,103 @@ def check_filter_properties(b_vals, a_vals, fs, pass_type, f_range,
     >>> passes = check_filter_properties(filter_coefs, 1, fs, pass_type, f_range)
     Transition bandwidth is 0.5 Hz.
     Pass/stop bandwidth is 24.0 Hz.
+    """
+
+    # Import utility functions inside function to avoid circular imports
+    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
+
+    # Compute the frequency response
+    f_db, db = compute_frequency_response(b_vals, a_vals, fs)
+
+    # Check that frequency response goes below transition level (has significant attenuation)
+    if np.min(db) >= transitions[0]:
+        passes = False
+        warn('The filter attenuation never goes below {} dB.'\
+             'Increase filter length.'.format(transitions[0]))
+        # If there is no attenuation, cannot calculate bands, so return here
+        return passes
 
+    # Check that both sides of a bandpass have significant attenuation
+    if pass_type == 'bandpass':
+        if db[0] >= transitions[0] or db[-1] >= transitions[0]:
+            passes = False
+            warn('The low or high frequency stopband never gets attenuated by'\
+                 'more than {} dB. Increase filter length.'.format(abs(transitions[0])))
+
+    # 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])
+
+    # 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
+    if verbose:
+        print('Transition bandwidth is {:.1f} Hz.'.format(transition_bw))
+        print('Pass/stop bandwidth is {:.1f} Hz.'.format(pass_bw))
+
+    return passes
+
+def sos_check_filter_properties(sos, fs, pass_type, f_range,
+                                transitions=(-20, -3), verbose=True):
+    """Check a filters properties, including pass band and transition band.
+
+    Parameters
+    ----------
+    sos : 2d array of shape (n, 6)
+        Second order series coefficients for an IIR filter.
+    fs : float
+        Sampling rate, in Hz.
+    pass_type : {'bandpass', 'bandstop', 'lowpass', 'highpass'}
+        Which kind of filter to apply:
+
+        * 'bandpass': apply a bandpass filter
+        * 'bandstop': apply a bandstop (notch) filter
+        * 'lowpass': apply a lowpass filter
+        * 'highpass' : apply a highpass filter
+    f_range : tuple of (float, float) or float
+        Cutoff frequency(ies) used for filter, specified as f_lo & f_hi.
+        For 'bandpass' & 'bandstop', must be a tuple.
+        For 'lowpass' or 'highpass', can be a float that specifies pass frequency, or can be
+        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.
+    verbose : bool, optional, default: True
+        Whether to print out transition and pass bands.
+
+    Returns
+    -------
+    passes : bool
+        Whether all the checks pass. False if one or more checks fail.
+
+    Examples
+    --------
     Check the properties of an IIR filter:
 
     >>> from neurodsp.filt.iir import design_iir_filter
     >>> fs, pass_type, f_range, order = 500, 'bandstop', (55, 65), 7
-    >>> b_vals, a_vals = design_iir_filter(fs, pass_type, f_range, order)
-    >>> passes = check_filter_properties(b_vals, a_vals, fs, pass_type, f_range)
+    >>> sos = design_iir_filter(fs, pass_type, f_range, order)
+    >>> passes = sos_check_filter_properties(sos, fs, pass_type, f_range)
     Transition bandwidth is 1.5 Hz.
     Pass/stop bandwidth is 10.0 Hz.
     """
 
     # Import utility functions inside function to avoid circular imports
-    from neurodsp.filt.utils import (compute_frequency_response,
+    from neurodsp.filt.utils import (sos_compute_frequency_response,
                                      compute_pass_band, compute_transition_band)
 
     # Initialize variable to keep track if all checks pass
     passes = True
 
     # Compute the frequency response
-    f_db, db = compute_frequency_response(b_vals, a_vals, fs)
+    f_db, db = sos_compute_frequency_response(sos, fs)
 
     # Check that frequency response goes below transition level (has significant attenuation)
     if np.min(db) >= transitions[0]:

neurodsp/filt/fir.py

@@ -174,15 +174,14 @@ def design_fir_filter(fs, pass_type, f_range, n_cycles=3, n_seconds=None):
     f_lo, f_hi = check_filter_definition(pass_type, f_range)
     filt_len = compute_filter_length(fs, pass_type, f_lo, f_hi, n_cycles, n_seconds)
 
-    f_nyq = compute_nyquist(fs)
     if pass_type == 'bandpass':
-        filter_coefs = firwin(filt_len, (f_lo, f_hi), pass_zero=False, nyq=f_nyq)
+        filter_coefs = firwin(filt_len, (f_lo, f_hi), pass_zero=False, fs=fs)
     elif pass_type == 'bandstop':
-        filter_coefs = firwin(filt_len, (f_lo, f_hi), nyq=f_nyq)
+        filter_coefs = firwin(filt_len, (f_lo, f_hi), fs=fs)
     elif pass_type == 'highpass':
-        filter_coefs = firwin(filt_len, f_lo, pass_zero=False, nyq=f_nyq)
+        filter_coefs = firwin(filt_len, f_lo, pass_zero=False, fs=fs)
     elif pass_type == 'lowpass':
-        filter_coefs = firwin(filt_len, f_hi, nyq=f_nyq)
+        filter_coefs = firwin(filt_len, f_hi, fs=fs)
 
     return filter_coefs
 

neurodsp/filt/iir.py

@@ -2,11 +2,11 @@
 
 from warnings import warn
 
-from scipy.signal import butter, filtfilt
+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.checks import check_filter_definition, check_filter_properties
+from neurodsp.filt.utils import compute_nyquist, compute_frequency_response, sos_compute_frequency_response
+from neurodsp.filt.checks import check_filter_definition, check_filter_properties, sos_check_filter_properties
 from neurodsp.plts.filt import plot_frequency_response
 
 ###################################################################################################
@@ -42,14 +42,14 @@ def filter_signal_iir(sig, fs, pass_type, f_range, butterworth_order,
     plot_properties : bool, optional, default: False
         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 IIR filter.
+        If True, return the second order series coefficients of the IIR filter.
 
     Returns
     -------
     sig_filt : 1d array
         Filtered time series.
-    filter_coefs : tuple of (1d array, 1d array)
-        Filter coefficients of the IIR filter, as (b_vals, a_vals).
+    sos : 2d array of shape (n, 6)
+        Second order series coefficients of the IIR filter.
         Only returned if `return_filter` is True.
 
     Examples
@@ -64,42 +64,40 @@ def filter_signal_iir(sig, fs, pass_type, f_range, butterworth_order,
     """
 
     # Design filter
-    b_vals, a_vals = design_iir_filter(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(b_vals, a_vals, fs, pass_type, f_range, verbose=print_transitions)
+    sos_check_filter_properties(sos, fs, pass_type, f_range, verbose=print_transitions)
 
     # Remove any NaN on the edges of 'sig'
     sig, sig_nans = remove_nans(sig)
 
     # Apply filter
-    sig_filt = apply_iir_filter(sig, b_vals, a_vals)
+    sig_filt = apply_iir_filter(sig, sos)
 
     # Add NaN back on the edges of 'sig', if there were any at the beginning
     sig_filt = restore_nans(sig_filt, sig_nans)
 
     # Plot frequency response, if desired
     if plot_properties:
-        f_db, db = compute_frequency_response(b_vals, a_vals, fs)
+        f_db, db = sos_compute_frequency_response(sos, fs)
         plot_frequency_response(f_db, db)
 
     if return_filter:
-        return sig_filt, (b_vals, a_vals)
+        return sig_filt, sos
     else:
         return sig_filt
 
 
-def apply_iir_filter(sig, b_vals, a_vals):
+def apply_iir_filter(sig, sos):
     """Apply an IIR filter to a signal.
 
     Parameters
     ----------
     sig : array
         Time series to be filtered.
-    b_vals : 1d array
-        B value filter coefficients for an IIR filter.
-    a_vals : 1d array
-        A value filter coefficients for an IIR filter.
+    sos : 2d array of shape (n, 6)
+        Second order series coefficients for an IIR filter.
 
     Returns
     -------
@@ -113,12 +111,12 @@ def apply_iir_filter(sig, b_vals, a_vals):
     >>> from neurodsp.sim import sim_combined
     >>> sig = sim_combined(n_seconds=10, fs=500,
     ...                    components={'sim_powerlaw': {}, 'sim_oscillation' : {'freq': 10}})
-    >>> b_vals, a_vals = design_iir_filter(fs=500, pass_type='bandstop',
+    >>> sos = design_iir_filter(fs=500, pass_type='bandstop',
     ...                                    f_range=(55, 65), butterworth_order=7)
-    >>> filt_signal = apply_iir_filter(sig, b_vals, a_vals)
+    >>> filt_signal = apply_iir_filter(sig, sos)
     """
 
-    return filtfilt(b_vals, a_vals, sig)
+    return sosfiltfilt(sos, sig)
 
 
 def design_iir_filter(fs, pass_type, f_range, butterworth_order):
@@ -146,23 +144,17 @@ def design_iir_filter(fs, pass_type, f_range, butterworth_order):
 
     Returns
     -------
-    b_vals : 1d array
-        B value filter coefficients for an IIR filter.
-    a_vals : 1d array
-        A value filter coefficients for an IIR filter.
+    sos : 2d array of shape (n, 6)
+        Second order series coefficients for an IIR filter.
 
     Examples
     --------
     Compute coefficients for a bandstop IIR filter:
 
-    >>> b_vals, a_vals = design_iir_filter(fs=500, pass_type='bandstop',
-    ...                                    f_range=(55, 65), butterworth_order=7)
+    >>> sos = design_iir_filter(fs=500, pass_type='bandstop',
+    ...                         f_range=(55, 65), butterworth_order=7)
     """
 
-    # Warn about only recommending IIR for bandstop
-    if pass_type != 'bandstop':
-        warn('IIR filters are not recommended other than for notch filters.')
-
     # Check filter definition
     f_lo, f_hi = check_filter_definition(pass_type, f_range)
 
@@ -175,6 +167,6 @@ def design_iir_filter(fs, pass_type, f_range, butterworth_order):
         win = f_hi / f_nyq
 
     # Design filter
-    b_vals, a_vals = butter(butterworth_order, win, pass_type)
+    sos = butter(butterworth_order, win, pass_type, output='sos')
 
-    return b_vals, a_vals
+    return sos

neurodsp/filt/utils.py

@@ -1,7 +1,7 @@
 """Utility functions for filtering."""
 
 import numpy as np
-from scipy.signal import freqz
+from scipy.signal import freqz, sosfreqz
 
 from neurodsp.utils.decorators import multidim
 from neurodsp.filt.checks import check_filter_definition
@@ -66,16 +66,42 @@ def compute_frequency_response(b_vals, a_vals, fs):
     >>> from neurodsp.filt.fir import design_fir_filter
     >>> filter_coefs = design_fir_filter(fs=500, pass_type='bandpass', f_range=(8, 12))
     >>> f_db, db = compute_frequency_response(filter_coefs, 1, fs=500)
+    """
+
+    w_vals, h_vals = freqz(b_vals, a_vals, worN=int(fs * 2))
+    f_db = w_vals * fs / (2. * np.pi)
+    db = 20 * np.log10(abs(h_vals))
+
+    return f_db, db
+
+def sos_compute_frequency_response(sos, fs):
+    """Compute the frequency response of a filter.
+
+    Parameters
+    ----------
+    sos: 2d array of shape (n, 6)
+        Second-order filter coefficients.
+    fs : float
+        Sampling rate, in Hz.
+
+    Returns
+    -------
+    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.
 
-    Compute the frequency response for an IIR filter:
+    Examples
+    --------
+    Compute the frequency response for an IIR filter given second order series coefficients:
 
     >>> from neurodsp.filt.iir import design_iir_filter
-    >>> b_vals, a_vals = design_iir_filter(fs=500, pass_type='bandstop',
+    >>> sos = design_iir_filter(fs=500, pass_type='bandstop',
     ...                                    f_range=(55, 65), butterworth_order=7)
-    >>> f_db, db = compute_frequency_response(b_vals, a_vals, fs=500)
+    >>> f_db, db = sos_compute_frequency_response(sos, fs=500)
     """
 
-    w_vals, h_vals = freqz(b_vals, a_vals, worN=int(fs * 2))
+    w_vals, h_vals = sosfreqz(sos, worN=int(fs * 2))
     f_db = w_vals * fs / (2. * np.pi)
     db = 20 * np.log10(abs(h_vals))
 
@@ -158,9 +184,9 @@ def compute_transition_band(f_db, db, low=-20, high=-3):
     Compute the transition band of an IIR filter, using the computed frequency response:
 
     >>> from neurodsp.filt.iir import design_iir_filter
-    >>> b_vals, a_vals = design_iir_filter(fs=500, pass_type='bandstop',
+    >>> sos = design_iir_filter(fs=500, pass_type='bandstop',
     ...                                    f_range=(10, 20), butterworth_order=7)
-    >>> f_db, db = compute_frequency_response(b_vals, a_vals, fs=500)
+    >>> f_db, db = sos_compute_frequency_response(sos, fs=500)
     >>> compute_transition_band(f_db, db, low=-20, high=-3)
     2.0
     """

neurodsp/tests/test_filt_iir.py

@@ -22,7 +22,7 @@ def test_filter_signal_iir_2d(tsig2d):
 
 def test_apply_iir_filter(tsig):
 
-    out = apply_iir_filter(tsig, np.array([1, 1, 1, 1, 1]), np.array([1, 1, 1, 1, 1]))
+    out = apply_iir_filter(tsig, np.array([1, 1, 1, 1, 1, 1]))
     assert out.shape == tsig.shape
 
 def test_design_iir_filter():