Repository: scipy/scipy
Pull request: ENH: add IIR comb filter to scipy.signal
#### Reference issue
Closes #11558 

#### What does this implement/fix?
Implements IIR notching and peaking comb filters in the scipy.signal module.

Problem statement:
ENH: IIR comb filter
**Describe the solution you'd like**
I think it would be very useful to support the design of IIR comb filters. 

**Describe alternatives you've considered**
The only alternative I could think of would be the use of multiple IIR notch filters, but doesn't seem like a great idea.

Hints from issue discussion:
I think this is a reasonable feature request. Would you be interested in putting together a PR for this?

Incidentally, I notice that the current peak and notch design functions only output TF form, these should get zpk and sos output options as well.
Happy to give it a go. Likely that I will need some help at some point.

Also happy to include zpk and sos support for peak and notch filters. Should this go into a separate PR?
Yes, a separate PR would be preferable. I'll be happy to review what you put together!
@agamemnonc @e-q @larsoner I don't see any PR on IIR comb filters. I just recently implemented gammatone filters, so let me know if I can work on this.
Sure, go for it!
@skresearcher Yes, please go for it -- I hardly have any free time these days to spend on this, unfortunately.

Commit message: ENH: add IIR comb filter to scipy.signal

Patch to review:
diff --git a/scipy/signal/__init__.py b/scipy/signal/__init__.py
index 97a49f73f6a3..66944e1c35a6 100644
--- a/scipy/signal/__init__.py
+++ b/scipy/signal/__init__.py
@@ -159,6 +159,7 @@
    bessel -- Bessel (no order selection available -- try butterod)
    iirnotch      -- Design second-order IIR notch digital filter.
    iirpeak       -- Design second-order IIR peak (resonant) digital filter.
+   iircomb       -- Design Nth-order IIR comb filter.
 
 Continuous-time linear systems
 ==============================
diff --git a/scipy/signal/filter_design.py b/scipy/signal/filter_design.py
index 4fcf5586507f..c5f5d716b94d 100644
--- a/scipy/signal/filter_design.py
+++ b/scipy/signal/filter_design.py
@@ -27,7 +27,7 @@
            'tf2sos', 'sos2tf', 'zpk2sos', 'sos2zpk', 'group_delay',
            'sosfreqz', 'iirnotch', 'iirpeak', 'bilinear_zpk',
            'lp2lp_zpk', 'lp2hp_zpk', 'lp2bp_zpk', 'lp2bs_zpk',
-           'gammatone']
+           'gammatone', 'iircomb']
 
 
 class BadCoefficients(UserWarning):
@@ -4837,7 +4837,7 @@ def _design_notch_peak_filter(w0, Q, ftype, fs=2.0):
     bw = bw*np.pi
     w0 = w0*np.pi
 
-    # Compute -3dB atenuation
+    # Compute -3dB attenuation
     gb = 1/np.sqrt(2)
 
     if ftype == "notch":
@@ -4864,6 +4864,191 @@ def _design_notch_peak_filter(w0, Q, ftype, fs=2.0):
     return b, a
 
 
+def iircomb(w0, Q, N, ftype='notch', fs=2.0):
+    """
+    Design Nth-order IIR notching or peaking digital comb filter.
+
+    A notching comb filter is a band-stop filter with a narrow bandwidth
+    (high quality factor). It rejects a narrow frequency band and
+    leaves the rest of the spectrum little changed.
+
+    A peaking comb filter is a band-pass filter with a narrow bandwidth
+    (high quality factor). It rejects components outside a narrow
+    frequency band.
+
+    Parameters
+    ----------
+    w0 : float
+        Frequency to remove or pass from a signal. If `fs` is specified, this
+        is in the same units as `fs`. By default, it is a normalized scalar
+        that must satisfy  ``0 < w0 < 1``, with ``w0 = 1`` corresponding to
+        half of the sampling frequency.
+    Q : float
+        Quality factor. Dimensionless parameter that characterizes
+        notch filter -3 dB bandwidth ``bw`` relative to its center
+        frequency, ``Q = w0/bw``.
+    N : int
+        The order of the filter. Must be between 0 and 24.
+    ftype : {'notch', 'peak'}
+        The type of comb filter generated by the function. Default is 'notch'.
+    fs : float, optional
+        The sampling frequency of the signal. Default is 2.0.
+
+    Returns
+    -------
+    b, a : ndarray, ndarray
+        Numerator (``b``) and denominator (``a``) polynomials
+        of the IIR filter.
+
+    Raises
+    ------
+    ValueError
+        If `w0` is less than or equal to 0 or greater than or equal to
+        ``fs/2``, if `ftype` is not 'notch' or 'peak', if `N` is less than
+        or equal to 0 or greater than 24.
+
+    See Also
+    --------
+    iirnotch
+    iirpeak
+
+    References
+    ----------
+    .. [1] Sophocles J. Orfanidis, "Introduction To Signal Processing",
+           Prentice-Hall, 1996
+
+    Examples
+    --------
+    Design and plot 10th order notching comb filter at 60 Hz for a
+    signal sampled at 200 Hz, using a quality factor Q = 30
+
+    >>> from scipy import signal
+    >>> import matplotlib.pyplot as plt
+
+    >>> fs = 200.0  # Sample frequency (Hz)
+    >>> f0 = 60.0  # Frequency to be removed from signal (Hz)
+    >>> Q = 30.0  # Quality factor
+    >>> # Design notching comb filter
+    >>> b, a = signal.iircomb(f0, Q, 10, ftype='notch', fs=fs)
+
+    >>> # Frequency response
+    >>> freq, h = signal.freqz(b, a, fs=fs)
+    >>> response = abs(h)
+    >>> # To avoid divide by zero when graphing
+    >>> response[response == 0] = 1e-20
+    >>> # Plot
+    >>> fig, ax = plt.subplots(2, 1, figsize=(8, 6))
+    >>> ax[0].plot(freq, 20*np.log10(abs(response)), color='blue')
+    >>> ax[0].set_title("Frequency Response")
+    >>> ax[0].set_ylabel("Amplitude (dB)", color='blue')
+    >>> ax[0].set_xlim([0, 100])
+    >>> ax[0].set_ylim([-30, 10])
+    >>> ax[0].grid()
+    >>> ax[1].plot(freq, np.unwrap(np.angle(h))*180/np.pi, color='green')
+    >>> ax[1].set_ylabel("Angle (degrees)", color='green')
+    >>> ax[1].set_xlabel("Frequency (Hz)")
+    >>> ax[1].set_xlim([0, 100])
+    >>> ax[1].set_yticks([-90, -60, -30, 0, 30, 60, 90])
+    >>> ax[1].set_ylim([-90, 90])
+    >>> ax[1].grid()
+    >>> plt.show()
+
+    Design and plot 10th order peaking comb filter at 300 Hz for a
+    signal sampled at 1000 Hz, using a quality factor Q = 30
+
+    >>> fs = 1000.0  # Sample frequency (Hz)
+    >>> f0 = 300.0  # Frequency to be retained (Hz)
+    >>> Q = 30.0  # Quality factor
+    >>> # Design peaking filter
+    >>> b, a = signal.iircomb(f0, Q, 10, ftype='peak', fs=fs)
+
+    >>> # Frequency response
+    >>> freq, h = signal.freqz(b, a, fs=fs)
+    >>> response = abs(h)
+    >>> # To avoid divide by zero when graphing
+    >>> response[response == 0] = 1e-20
+    >>> # Plot
+    >>> fig, ax = plt.subplots(2, 1, figsize=(8, 6))
+    >>> ax[0].plot(freq, 20*np.log10(np.maximum(abs(h), 1e-5)), color='blue')
+    >>> ax[0].set_title("Frequency Response")
+    >>> ax[0].set_ylabel("Amplitude (dB)", color='blue')
+    >>> ax[0].set_xlim([0, 500])
+    >>> ax[0].set_ylim([-80, 10])
+    >>> ax[0].grid()
+    >>> ax[1].plot(freq, np.unwrap(np.angle(h))*180/np.pi, color='green')
+    >>> ax[1].set_ylabel("Angle (degrees)", color='green')
+    >>> ax[1].set_xlabel("Frequency (Hz)")
+    >>> ax[1].set_xlim([0, 500])
+    >>> ax[1].set_yticks([-90, -60, -30, 0, 30, 60, 90])
+    >>> ax[1].set_ylim([-90, 90])
+    >>> ax[1].grid()
+    >>> plt.show()
+    """
+
+    # Convert w0, Q, and fs to float, N to int
+    w0 = float(w0)
+    Q = float(Q)
+    N = operator.index(N)
+    fs = float(fs)
+
+    # Check for invalid cutoff frequency or filter type
+    ftype = ftype.lower()
+    filter_types = ['notch', 'peak']
+    if not 0 < w0 < fs / 2:
+        raise ValueError("w0 must be between 0 and {}"
+                         " (nyquist), but given {}.".format(fs / 2, w0))
+    if ftype not in filter_types:
+        raise ValueError('ftype must be either notch or peak.')
+
+    # Check for invalid order
+    if not 0 < N <= 24:
+        raise ValueError("N must be > 0 and <= 24.")
+
+    # Compute frequency in radians and filter bandwith
+    # Eq. 11.3.1 (p. 574) from reference [1]
+    w0 = (2 * np.pi * w0) / fs
+    w_delta = w0 / Q
+
+    # Define base gain values depending on notch or peak filter
+    # Compute -3dB attenuation
+    # Eqs. 11.4.1 and 11.4.2 (p. 582) from reference [1]
+    if ftype == 'notch':
+        G0, G = [1, 0]
+    elif ftype == 'peak':
+        G0, G = [0, 1]
+    GB = 1 / np.sqrt(2)
+
+    # Compute beta
+    # Eq. 11.5.3 (p. 591) from reference [1]
+    beta = np.sqrt((GB**2 - G0**2) / (G**2 - GB**2)) * np.tan(N * w_delta / 4)
+
+    # Compute filter coefficients
+    # Eq 11.5.1 (p. 590) variables a, b, c from reference [1]
+    ax = (1 - beta) / (1 + beta)
+    bx = (G0 + G * beta) / (1 + beta)
+    cx = (G0 - G * beta) / (1 + beta)
+
+    # Compute numerator coefficients
+    # Eq 11.5.1 (p. 590) or Eq 11.5.4 (p. 591) from reference [1]
+    # b - cz^-N or b + cz^-N
+    b = np.zeros(N + 1)
+    b[0] = bx
+    b[-1] = cx
+    if ftype == 'notch':
+        b[-1] = -cx
+
+    # Compute denominator coefficients
+    # Eq 11.5.1 (p. 590) or Eq 11.5.4 (p. 591) from reference [1]
+    # 1 - az^-N or 1 + az^-N
+    a = np.zeros(N + 1)
+    a[0] = 1
+    a[-1] = ax
+    if ftype == 'notch':
+        a[-1] = -ax
+
+    return b, a
+
+
 def _hz_to_erb(hz):
     """
     Utility for converting from frequency (Hz) to the
@@ -4936,7 +5121,7 @@ def gammatone(freq, ftype, order=None, numtaps=None, fs=None):
     (array([ 0.00000000e+00,  2.22196719e-07,  1.64942101e-06,  4.99298227e-06,
         1.01993969e-05,  1.63125770e-05,  2.14648940e-05,  2.29947263e-05,
         1.76776931e-05,  2.04980537e-06, -2.72062858e-05, -7.28455299e-05,
-       -1.36651076e-04, -2.19066855e-04, -3.18905076e-04, -4.33156712e-04]), 
+       -1.36651076e-04, -2.19066855e-04, -3.18905076e-04, -4.33156712e-04]),
        [1.0])
 
     IIR Gammatone filter centered at 440 Hz
@@ -5032,8 +5217,8 @@ def gammatone(freq, ftype, order=None, numtaps=None, fs=None):
         g = np.abs(g)
 
         # Create empty filter coefficient lists
-        b = np.zeros((5))
-        a = np.zeros((9))
+        b = np.zeros(5)
+        a = np.zeros(9)
 
         # Calculate the numerator coefficients
         b[0] = (T ** 4) / g
diff --git a/scipy/signal/tests/test_filter_design.py b/scipy/signal/tests/test_filter_design.py
index cd65c8946f0e..4fef997ed637 100644
--- a/scipy/signal/tests/test_filter_design.py
+++ b/scipy/signal/tests/test_filter_design.py
@@ -10,14 +10,14 @@
 from pytest import raises as assert_raises
 
 from numpy import array, spacing, sin, pi, sort, sqrt
-from scipy.signal import (BadCoefficients, bessel, besselap, bilinear, buttap,
-                          butter, buttord, cheb1ap, cheb1ord, cheb2ap,
+from scipy.signal import (argrelextrema, BadCoefficients, bessel, besselap, bilinear,
+                          buttap, butter, buttord, cheb1ap, cheb1ord, cheb2ap,
                           cheb2ord, cheby1, cheby2, ellip, ellipap, ellipord,
                           firwin, freqs_zpk, freqs, freqz, freqz_zpk,
-                          gammatone, group_delay, iirdesign, iirfilter, iirnotch, 
-                          iirpeak, lp2bp, lp2bs, lp2hp, lp2lp, normalize, sos2tf, 
-                          sos2zpk, sosfreqz, tf2sos, tf2zpk, zpk2sos, zpk2tf,
-                          bilinear_zpk, lp2lp_zpk, lp2hp_zpk, lp2bp_zpk,
+                          gammatone, group_delay, iircomb, iirdesign, iirfilter, 
+                          iirnotch, iirpeak, lp2bp, lp2bs, lp2hp, lp2lp, normalize, 
+                          sos2tf, sos2zpk, sosfreqz, tf2sos, tf2zpk, zpk2sos, 
+                          zpk2tf, bilinear_zpk, lp2lp_zpk, lp2hp_zpk, lp2bp_zpk,
                           lp2bs_zpk)
 from scipy.signal.filter_design import (_cplxreal, _cplxpair, _norm_factor,
                                         _bessel_poly, _bessel_zeros)
@@ -3611,6 +3611,83 @@ def test_fs_param(self):
         assert_allclose(abs(hp[2]), 1, rtol=1e-10)
 
 
+class TestIIRComb(object):
+    # Test erroneus input cases
+    def test_invalid_input(self):
+        # w0 is <= 0 or >= fs / 2
+        fs = 1000
+        with pytest.raises(ValueError, match='w0 must be between '):
+            iircomb(-fs, 30, 4, fs=fs)
+            iircomb(0, 35, 5, fs=fs)
+            iircomb(fs / 2, 40, 6, fs=fs)
+            iircomb(fs, 35, 7, fs=fs)
+
+        # Filter type is not notch or peak
+        with pytest.raises(ValueError, match='ftype must be '):
+            iircomb(0.2, 30, 4, ftype='natch')
+            iircomb(0.5, 35, 5, ftype='comb')
+
+        # Order is <= 0 or > 24
+        with pytest.raises(ValueError, match='N must be '):
+            iircomb(440, 30, -50, ftype='notch', fs=16000)
+            iircomb(220, 35, 0, ftype='peak', fs=32000)
+            iircomb(110, 40, 25, ftype='peak', fs=44100)
+            iircomb(55, 35, 50, ftype='notch', fs=48000)
+
+    # Verify that the filter's frequency response contains a
+    # notch or peak at the cutoff frequency
+    def test_frequency_response(self):
+        ftypes = ['notch', 'peak']
+        for ftype in ftypes:
+            # Create a notching or peaking comb filter at 1000 Hz
+            b, a = iircomb(1000, 30, 10, ftype=ftype, fs=10000)
+
+            # Compute the frequency response
+            freqs, response = freqz(b, a, 1000, fs=10000)
+
+            # Find the notches/peaks using argrelextrema
+            if ftype == 'notch':
+                comb_points = argrelextrema(abs(response), np.less)
+            elif ftype == 'peak':
+                comb_points = argrelextrema(abs(response), np.greater)
+
+            # Verify that the first notch/peak sits at 1000 Hz
+            comb1 = comb_points[0]
+            freqs[comb1] == pytest.approx(1000)
+
+    # All built-in IIR filters are real, so should have perfectly
+    # symmetrical poles and zeros. Then ba representation (using
+    # numpy.poly) will be purely real instead of having negligible
+    # imaginary parts.
+    def test_iir_symmetry(self):
+        b, a = iircomb(440, 30, 4, fs=24000)
+        z, p, k = tf2zpk(b, a)
+        assert_array_equal(sorted(z), sorted(z.conj()))
+        assert_array_equal(sorted(p), sorted(p.conj()))
+        assert_equal(k, np.real(k))
+
+        assert_(issubclass(b.dtype.type, np.floating))
+        assert_(issubclass(a.dtype.type, np.floating))
+
+    # Verify filter coefficients with MATLAB's iircomb function
+    def test_ba_output(self):
+        b_notch, a_notch = iircomb(60, 35, 10, ftype='notch', fs=600)
+        b_notch2 = [0.957020174408697, 0.0, 0.0, 0.0, 0.0, 0.0,
+                    0.0, 0.0, 0.0, 0.0, -0.957020174408697]
+        a_notch2 = [1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+                    0.0, 0.0, 0.0, 0.0, -0.914040348817395]
+        assert_allclose(b_notch, b_notch2)
+        assert_allclose(a_notch, a_notch2)
+
+        b_peak, a_peak = iircomb(60, 35, 10, ftype='peak', fs=600)
+        b_peak2 = [0.0429798255913026, 0.0, 0.0, 0.0, 0.0, 0.0,
+                   0.0, 0.0, 0.0, 0.0, -0.0429798255913026]
+        a_peak2 = [1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+                   0.0, 0.0, 0.0, 0.0, 0.914040348817395]
+        assert_allclose(b_peak, b_peak2)
+        assert_allclose(a_peak, a_peak2)
+
+
 class TestIIRDesign(object):
 
     def test_exceptions(self):


Code review comments:
diff_hunk: @@ -3611,6 +3611,83 @@ def test_fs_param(self):
         assert_allclose(abs(hp[2]), 1, rtol=1e-10)
 
 
+class TestIIRComb(object):
+    # Test erroneus input cases
+    def test_invalid_input(self):
+        # w0 is <= 0 or >= fs / 2
+        fs = 1000
+        with pytest.raises(ValueError, match='w0 must be between '):
+            iircomb(-fs, 30, 4, fs=fs)
+            iircomb(0, 35, 5, fs=fs)
+            iircomb(fs / 2, 40, 6, fs=fs)
+            iircomb(fs, 35, 7, fs=fs)
+
+        # Filter type is not notch or peak
+        with pytest.raises(ValueError, match='ftype must be '):
+            iircomb(0.2, 30, 4, ftype='natch')
+            iircomb(0.5, 35, 5, ftype='comb')
+
+        # Order is <= 0 or > 24
+        with pytest.raises(ValueError, match='N must be '):
+            iircomb(440, 30, -50, ftype='notch', fs=16000)
+            iircomb(220, 35, 0, ftype='peak', fs=32000)
+            iircomb(110, 40, 25, ftype='peak', fs=44100)
+            iircomb(55, 35, 50, ftype='notch', fs=48000)
+
+    # Verify that the filter's frequency response contains a
+    # notch or peak at the cutoff frequency
+    def test_frequency_response(self):
+        ftypes = ['notch', 'peak']
path: scipy/signal/tests/test_filter_design.py
text: @user2:
nowadays we can `@user1.parametrize('ftypes', ('notch', 'peak')` the `test_frequency_response(self, ftypes)` function

@author:
Fixed.

diff_hunk: @@ -3611,6 +3611,83 @@ def test_fs_param(self):
         assert_allclose(abs(hp[2]), 1, rtol=1e-10)
 
 
+class TestIIRComb(object):
+    # Test erroneus input cases
+    def test_invalid_input(self):
+        # w0 is <= 0 or >= fs / 2
+        fs = 1000
+        with pytest.raises(ValueError, match='w0 must be between '):
+            iircomb(-fs, 30, 4, fs=fs)
+            iircomb(0, 35, 5, fs=fs)
path: scipy/signal/tests/test_filter_design.py
text: @user1:
These need to each be in their own context manager for them all to execute (the first failure will kick you out of the context without executing the subsequent lines). You could for example loop over
```
for args in [(-fs, 30, 4), ...]:
    with pytest.raises(...):
        iircomp(*args, fs=fs)
```

@author:
Fixed.

diff_hunk: @@ -3611,6 +3611,83 @@ def test_fs_param(self):
         assert_allclose(abs(hp[2]), 1, rtol=1e-10)
 
 
+class TestIIRComb(object):
+    # Test erroneus input cases
+    def test_invalid_input(self):
+        # w0 is <= 0 or >= fs / 2
+        fs = 1000
+        with pytest.raises(ValueError, match='w0 must be between '):
+            iircomb(-fs, 30, 4, fs=fs)
+            iircomb(0, 35, 5, fs=fs)
+            iircomb(fs / 2, 40, 6, fs=fs)
+            iircomb(fs, 35, 7, fs=fs)
+
+        # Filter type is not notch or peak
+        with pytest.raises(ValueError, match='ftype must be '):
+            iircomb(0.2, 30, 4, ftype='natch')
+            iircomb(0.5, 35, 5, ftype='comb')
+
+        # Order is <= 0 or > 24
+        with pytest.raises(ValueError, match='N must be '):
+            iircomb(440, 30, -50, ftype='notch', fs=16000)
+            iircomb(220, 35, 0, ftype='peak', fs=32000)
+            iircomb(110, 40, 25, ftype='peak', fs=44100)
+            iircomb(55, 35, 50, ftype='notch', fs=48000)
+
+    # Verify that the filter's frequency response contains a
+    # notch or peak at the cutoff frequency
+    def test_frequency_response(self):
+        ftypes = ['notch', 'peak']
+        for ftype in ftypes:
+            # Create a notching or peaking comb filter at 1000 Hz
+            b, a = iircomb(1000, 30, 10, ftype=ftype, fs=10000)
+
+            # Compute the frequency response
+            freqs, response = freqz(b, a, 1000, fs=10000)
+
+            # Find the notches/peaks using argrelextrema
+            if ftype == 'notch':
+                comb_points = argrelextrema(abs(response), np.less)
+            elif ftype == 'peak':
+                comb_points = argrelextrema(abs(response), np.greater)
+
+            # Verify that the first notch/peak sits at 1000 Hz
+            comb1 = comb_points[0]
+            freqs[comb1] == pytest.approx(1000)
+
+    # All built-in IIR filters are real, so should have perfectly
+    # symmetrical poles and zeros. Then ba representation (using
+    # numpy.poly) will be purely real instead of having negligible
+    # imaginary parts.
+    def test_iir_symmetry(self):
+        b, a = iircomb(440, 30, 4, fs=24000)
+        z, p, k = tf2zpk(b, a)
+        assert_array_equal(sorted(z), sorted(z.conj()))
+        assert_array_equal(sorted(p), sorted(p.conj()))
+        assert_equal(k, np.real(k))
+
+        assert_(issubclass(b.dtype.type, np.floating))
path: scipy/signal/tests/test_filter_design.py
text: @user1:
pytest directly handles `assert` with an informative traceback so you can just do `assert issubclass(...)`

@author:
Fixed.

diff_hunk: @@ -3611,6 +3611,83 @@ def test_fs_param(self):
         assert_allclose(abs(hp[2]), 1, rtol=1e-10)
 
 
+class TestIIRComb(object):
+    # Test erroneus input cases
+    def test_invalid_input(self):
+        # w0 is <= 0 or >= fs / 2
+        fs = 1000
+        with pytest.raises(ValueError, match='w0 must be between '):
+            iircomb(-fs, 30, 4, fs=fs)
+            iircomb(0, 35, 5, fs=fs)
+            iircomb(fs / 2, 40, 6, fs=fs)
+            iircomb(fs, 35, 7, fs=fs)
+
+        # Filter type is not notch or peak
+        with pytest.raises(ValueError, match='ftype must be '):
+            iircomb(0.2, 30, 4, ftype='natch')
+            iircomb(0.5, 35, 5, ftype='comb')
+
+        # Order is <= 0 or > 24
+        with pytest.raises(ValueError, match='N must be '):
+            iircomb(440, 30, -50, ftype='notch', fs=16000)
+            iircomb(220, 35, 0, ftype='peak', fs=32000)
+            iircomb(110, 40, 25, ftype='peak', fs=44100)
+            iircomb(55, 35, 50, ftype='notch', fs=48000)
+
+    # Verify that the filter's frequency response contains a
+    # notch or peak at the cutoff frequency
+    def test_frequency_response(self):
+        ftypes = ['notch', 'peak']
+        for ftype in ftypes:
+            # Create a notching or peaking comb filter at 1000 Hz
+            b, a = iircomb(1000, 30, 10, ftype=ftype, fs=10000)
+
+            # Compute the frequency response
+            freqs, response = freqz(b, a, 1000, fs=10000)
+
+            # Find the notches/peaks using argrelextrema
+            if ftype == 'notch':
+                comb_points = argrelextrema(abs(response), np.less)
+            elif ftype == 'peak':
+                comb_points = argrelextrema(abs(response), np.greater)
+
+            # Verify that the first notch/peak sits at 1000 Hz
+            comb1 = comb_points[0]
+            freqs[comb1] == pytest.approx(1000)
path: scipy/signal/tests/test_filter_design.py
text: @user1:
The SciPy codebase prefers using the numpy asserts like `assert_allclose` for stuff like this (no need / advantage to introduce a new convention for people to learn here AFAIK)

@author:
Fixed.