From 6dd4a51ab229de41ca61fe93c8042caaa2774afd Mon Sep 17 00:00:00 2001
From: Eric Lybrand <elybrand@ucsd.edu>
Date: Thu, 25 Jun 2020 13:31:27 -0700
Subject: [PATCH 1/5] Remove nyquist frequency. Input sampling frequency to
 firwin to match scipy API.

---
 neurodsp/filt/fir.py | 9 ++++-----
 1 file changed, 4 insertions(+), 5 deletions(-)

diff --git a/neurodsp/filt/fir.py b/neurodsp/filt/fir.py
index cc773733..b6cfa0cb 100644
--- a/neurodsp/filt/fir.py
+++ b/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
 

From 86a51c9c2a1419809c6db45a18fd365f18be2be1 Mon Sep 17 00:00:00 2001
From: Eric Lybrand <elybrand@ucsd.edu>
Date: Thu, 25 Jun 2020 13:39:23 -0700
Subject: [PATCH 2/5] Remove warning for filters that are not bandpass.

---
 neurodsp/filt/iir.py | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/neurodsp/filt/iir.py b/neurodsp/filt/iir.py
index c55cb764..5bb6592e 100644
--- a/neurodsp/filt/iir.py
+++ b/neurodsp/filt/iir.py
@@ -159,10 +159,6 @@ def design_iir_filter(fs, pass_type, f_range, butterworth_order):
     ...                                    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)
 

From 69150a687f0e638648903b4e3e57b4320d8d288f Mon Sep 17 00:00:00 2001
From: Eric Lybrand <elybrand@ucsd.edu>
Date: Thu, 25 Jun 2020 15:17:13 -0700
Subject: [PATCH 3/5] Change all iir filter related functions to use second
 order series coefficients. Update test to take sos coefficients.

---
 neurodsp/filt/checks.py         | 85 +++++++++++++++++++++++++++++++--
 neurodsp/filt/iir.py            | 50 +++++++++----------
 neurodsp/filt/utils.py          | 40 +++++++++++++---
 neurodsp/tests/test_filt_iir.py |  2 +-
 4 files changed, 138 insertions(+), 39 deletions(-)

diff --git a/neurodsp/filt/checks.py b/neurodsp/filt/checks.py
index 3208ea4b..464df29b 100644
--- a/neurodsp/filt/checks.py
+++ b/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]:
diff --git a/neurodsp/filt/iir.py b/neurodsp/filt/iir.py
index 5bb6592e..2772c435 100644
--- a/neurodsp/filt/iir.py
+++ b/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,17 +144,15 @@ 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)
     """
 
     # Check filter definition
@@ -171,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
diff --git a/neurodsp/filt/utils.py b/neurodsp/filt/utils.py
index de670290..cb0a7514 100644
--- a/neurodsp/filt/utils.py
+++ b/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
     """
diff --git a/neurodsp/tests/test_filt_iir.py b/neurodsp/tests/test_filt_iir.py
index d6a75a48..bd37c149 100644
--- a/neurodsp/tests/test_filt_iir.py
+++ b/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():

From cbea58675fc49d3c5a22e090cf2264f9ebd1a988 Mon Sep 17 00:00:00 2001
From: TomDonoghue <tdonoghue@ucsd.edu>
Date: Fri, 26 Jun 2020 12:52:57 -0700
Subject: [PATCH 4/5] refactor filter updates

---
 neurodsp/filt/checks.py | 100 ++++------------------------------------
 neurodsp/filt/fir.py    |   2 +-
 neurodsp/filt/iir.py    |  22 ++++-----
 neurodsp/filt/utils.py  |  62 +++++++++----------------
 4 files changed, 42 insertions(+), 144 deletions(-)

diff --git a/neurodsp/filt/checks.py b/neurodsp/filt/checks.py
index 464df29b..c96b1e5a 100644
--- a/neurodsp/filt/checks.py
+++ b/neurodsp/filt/checks.py
@@ -89,16 +89,18 @@ def check_filter_definition(pass_type, f_range):
     return f_lo, f_hi
 
 
-def check_filter_properties(b_vals, a_vals, fs, pass_type, f_range,
+def check_filter_properties(filter_coefs, a_vals, fs, pass_type, f_range,
                             transitions=(-20, -3), verbose=True):
     """Check a filters properties, including pass band and transition band.
 
     Parameters
     ----------
-    b_vals : 1d array
-        B value filter coefficients for a filter.
-    a_vals : 1d array
-        A value filter coefficients for a filter.
+    filter_coefs : 1d or 2d array
+        If 1d, interpreted as the B-value filter coefficients.
+        If 2d, interpreted as the second-order (sos) filter coefficients.
+    a_vals : 1d array or None
+        The A-value filter coefficients for a filter.
+        If second-order filter coefficients are provided in `filter_coefs`, must be None.
     fs : float
         Sampling rate, in Hz.
     pass_type : {'bandpass', 'bandstop', 'lowpass', 'highpass'}
@@ -143,93 +145,7 @@ def check_filter_properties(b_vals, a_vals, fs, pass_type, f_range,
     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
-    >>> 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 (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 = sos_compute_frequency_response(sos, fs)
+    f_db, db = compute_frequency_response(filter_coefs, a_vals, fs)
 
     # Check that frequency response goes below transition level (has significant attenuation)
     if np.min(db) >= transitions[0]:
diff --git a/neurodsp/filt/fir.py b/neurodsp/filt/fir.py
index b6cfa0cb..5fff10d3 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_nyquist, compute_frequency_response, remove_filter_edges
+from neurodsp.filt.utils import compute_frequency_response, remove_filter_edges
 from neurodsp.filt.checks import (check_filter_definition, check_filter_properties,
                                   check_filter_length)
 
diff --git a/neurodsp/filt/iir.py b/neurodsp/filt/iir.py
index 2772c435..70161caf 100644
--- a/neurodsp/filt/iir.py
+++ b/neurodsp/filt/iir.py
@@ -1,12 +1,10 @@
 """Filter signals with IIR filters."""
 
-from warnings import warn
-
 from scipy.signal import butter, sosfiltfilt
 
 from neurodsp.utils import remove_nans, restore_nans
-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.filt.utils import compute_nyquist, compute_frequency_response
+from neurodsp.filt.checks import check_filter_definition, check_filter_properties
 from neurodsp.plts.filt import plot_frequency_response
 
 ###################################################################################################
@@ -48,8 +46,8 @@ def filter_signal_iir(sig, fs, pass_type, f_range, butterworth_order,
     -------
     sig_filt : 1d array
         Filtered time series.
-    sos : 2d array of shape (n, 6)
-        Second order series coefficients of the IIR filter.
+    sos : 2d array
+        Second order series coefficients of the IIR filter. Has shape of (n_sections, 6).
         Only returned if `return_filter` is True.
 
     Examples
@@ -67,7 +65,7 @@ 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
-    sos_check_filter_properties(sos, fs, pass_type, f_range, verbose=print_transitions)
+    check_filter_properties(sos, None, fs, pass_type, f_range, verbose=print_transitions)
 
     # Remove any NaN on the edges of 'sig'
     sig, sig_nans = remove_nans(sig)
@@ -80,7 +78,7 @@ def filter_signal_iir(sig, fs, pass_type, f_range, butterworth_order,
 
     # Plot frequency response, if desired
     if plot_properties:
-        f_db, db = sos_compute_frequency_response(sos, fs)
+        f_db, db = compute_frequency_response(sos, None, fs)
         plot_frequency_response(f_db, db)
 
     if return_filter:
@@ -96,8 +94,8 @@ def apply_iir_filter(sig, sos):
     ----------
     sig : array
         Time series to be filtered.
-    sos : 2d array of shape (n, 6)
-        Second order series coefficients for an IIR filter.
+    sos : 2d array
+        Second order series coefficients for an IIR filter. Has shape of (n_sections, 6).
 
     Returns
     -------
@@ -144,8 +142,8 @@ def design_iir_filter(fs, pass_type, f_range, butterworth_order):
 
     Returns
     -------
-    sos : 2d array of shape (n, 6)
-        Second order series coefficients for an IIR filter.
+    sos : 2d array
+        Second order series coefficients for an IIR filter. Has shape of (n_sections, 6).
 
     Examples
     --------
diff --git a/neurodsp/filt/utils.py b/neurodsp/filt/utils.py
index cb0a7514..93629bfe 100644
--- a/neurodsp/filt/utils.py
+++ b/neurodsp/filt/utils.py
@@ -40,15 +40,17 @@ def infer_passtype(f_range):
     return pass_type
 
 
-def compute_frequency_response(b_vals, a_vals, fs):
+def compute_frequency_response(filter_coefs, a_vals, fs):
     """Compute the frequency response of a filter.
 
     Parameters
     ----------
-    b_vals : 1d array
-        B value filter coefficients for a filter.
-    a_vals : 1d array
-        A value filter coefficients for a filter.
+    filter_coefs : 1d or 2d array
+        If 1d, interpreted as the B-value filter coefficients.
+        If 2d, interpreted as the second-order (sos) filter coefficients.
+    a_vals : 1d array or None
+        The A-value filter coefficients for a filter.
+        If second-order filter coefficients are provided in `filter_coefs`, must be None.
     fs : float
         Sampling rate, in Hz.
 
@@ -66,42 +68,24 @@ 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.
-
-    Examples
-    --------
-    Compute the frequency response for an IIR filter given second order series coefficients:
+    Compute the frequency response for an IIR filter, which uses SOS coefficients:
 
     >>> from neurodsp.filt.iir import design_iir_filter
-    >>> sos = design_iir_filter(fs=500, pass_type='bandstop',
-    ...                                    f_range=(55, 65), butterworth_order=7)
-    >>> f_db, db = sos_compute_frequency_response(sos, fs=500)
+    >>> sos_coefs = design_iir_filter(fs=500, pass_type='bandpass',
+    ...                               f_range=(8, 12), butterworth_order=3)
+    >>> f_db, db = compute_frequency_response(sos_coefs, None, fs=500)
     """
 
-    w_vals, h_vals = sosfreqz(sos, worN=int(fs * 2))
+    if filter_coefs.ndim == 1 and a_vals is not None:
+        # Compute response for B & A value filter coefficient inputs
+        w_vals, h_vals = freqz(filter_coefs, a_vals, worN=int(fs * 2))
+    elif filter_coefs.ndim == 2 and a_vals is None:
+        # Compute response for sos filter coefficient inputs
+        w_vals, h_vals = sosfreqz(filter_coefs, worN=int(fs * 2))
+    else:
+        raise ValueError("The organization of the filter coefficient inputs is not understood.")
+
     f_db = w_vals * fs / (2. * np.pi)
     db = 20 * np.log10(abs(h_vals))
 
@@ -135,7 +119,7 @@ def compute_pass_band(fs, pass_type, f_range):
 
     Examples
     --------
-    Compute the bandwith of a bandpass filter:
+    Compute the bandwidth of a bandpass filter:
 
     >>> compute_pass_band(fs=500, pass_type='bandpass', f_range=(5, 25))
     20.0
@@ -185,8 +169,8 @@ def compute_transition_band(f_db, db, low=-20, high=-3):
 
     >>> from neurodsp.filt.iir import design_iir_filter
     >>> sos = design_iir_filter(fs=500, pass_type='bandstop',
-    ...                                    f_range=(10, 20), butterworth_order=7)
-    >>> f_db, db = sos_compute_frequency_response(sos, fs=500)
+    ...                         f_range=(10, 20), butterworth_order=7)
+    >>> f_db, db = compute_frequency_response(sos, None, fs=500)
     >>> compute_transition_band(f_db, db, low=-20, high=-3)
     2.0
     """

From 922a4bfd7571d4237eb9f682ebbe3f5e55fcc870 Mon Sep 17 00:00:00 2001
From: ryanhammonds <Ryan.hammonds@utexas.edu>
Date: Fri, 26 Jun 2020 14:23:51 -0700
Subject: [PATCH 5/5] increased coverage for filters

---
 neurodsp/tests/test_filt_checks.py | 14 ++++++++++++--
 neurodsp/tests/test_filt_fir.py    |  5 +++++
 neurodsp/tests/test_filt_iir.py    |  3 +++
 neurodsp/tests/test_filt_utils.py  |  4 ++++
 4 files changed, 24 insertions(+), 2 deletions(-)

diff --git a/neurodsp/tests/test_filt_checks.py b/neurodsp/tests/test_filt_checks.py
index 144431fb..32b5776e 100644
--- a/neurodsp/tests/test_filt_checks.py
+++ b/neurodsp/tests/test_filt_checks.py
@@ -46,8 +46,18 @@ def test_check_filter_definition():
 def test_check_filter_properties():
 
     filter_coefs = design_fir_filter(FS, 'bandpass', (8, 12))
-    check_filter_properties(filter_coefs, 1, FS, 'bandpass', (8, 12))
-    assert True
+    
+    passes = check_filter_properties(filter_coefs, 1, FS, 'bandpass', (8, 12))
+    assert passes is True
+
+    filter_coefs = design_fir_filter(FS, 'bandstop', (8, 12))
+    passes = check_filter_properties(filter_coefs, 1, FS, 'bandpass', (8, 12))
+    assert passes is False
+
+    filter_coefs = design_fir_filter(FS, 'bandpass', (20, 21))
+    passes = check_filter_properties(filter_coefs, 1, FS, 'bandpass', (8, 12))
+    assert passes is False
+
 
 def test_check_filter_length():
 
diff --git a/neurodsp/tests/test_filt_fir.py b/neurodsp/tests/test_filt_fir.py
index e695457b..ddffdfd4 100644
--- a/neurodsp/tests/test_filt_fir.py
+++ b/neurodsp/tests/test_filt_fir.py
@@ -1,5 +1,6 @@
 """Tests for FIR filters."""
 
+from pytest import raises
 import numpy as np
 
 from neurodsp.tests.settings import FS
@@ -56,3 +57,7 @@ def test_compute_filter_length():
     expected_filt_len = int(np.ceil(fs * n_cycles / f_lo)) + 1
     filt_len = compute_filter_length(fs, 'bandpass', f_lo, f_hi, n_cycles=n_cycles, n_seconds=None)
     assert filt_len == expected_filt_len
+
+    with raises(ValueError):
+        filt_len = compute_filter_length(fs, 'bandpass', f_lo, f_hi)
+
diff --git a/neurodsp/tests/test_filt_iir.py b/neurodsp/tests/test_filt_iir.py
index bd37c149..4b209035 100644
--- a/neurodsp/tests/test_filt_iir.py
+++ b/neurodsp/tests/test_filt_iir.py
@@ -20,6 +20,9 @@ def test_filter_signal_iir_2d(tsig2d):
    assert out.shape == tsig2d.shape
    assert sum(~np.isnan(out[0, :])) > 0
 
+   out, sos = filter_signal_iir(tsig2d, FS, 'bandpass', (8, 12), 3, return_filter=True)
+   assert np.shape(sos)[1] == 6
+
 def test_apply_iir_filter(tsig):
 
     out = apply_iir_filter(tsig, np.array([1, 1, 1, 1, 1, 1]))
diff --git a/neurodsp/tests/test_filt_utils.py b/neurodsp/tests/test_filt_utils.py
index 01897d89..c839a75d 100644
--- a/neurodsp/tests/test_filt_utils.py
+++ b/neurodsp/tests/test_filt_utils.py
@@ -1,5 +1,6 @@
 """Tests for filter utilities."""
 
+from pytest import raises
 from neurodsp.tests.settings import FS
 
 from neurodsp.filt.utils import *
@@ -19,6 +20,9 @@ def test_compute_frequency_response():
     filter_coefs = design_fir_filter(FS, 'bandpass', (8, 12))
     f_db, db = compute_frequency_response(filter_coefs, 1, FS)
 
+    with raises(ValueError):
+        f_db, db = compute_frequency_response(filter_coefs, None, FS)
+
 def test_compute_pass_band():
 
     assert compute_pass_band(FS, 'bandpass', (4, 8)) == 4.