MAINT: signal: adjust digital filter design warnings

scipy/signal/filter_design.py

@@ -28,7 +28,7 @@
            'BadCoefficients', 'freqs_zpk', 'freqz_zpk',
            'tf2sos', 'sos2tf', 'zpk2sos', 'sos2zpk', 'group_delay',
            'sosfreqz', 'iirnotch', 'iirpeak', 'bilinear_zpk',
-           'lp2lp_zpk', 'lp2hp_zpk', 'lp2bp_zpk', 'lp2bs_zpk']
+           'lp2lp_zpk', 'lp2hp_zpk', 'lp2bp_zpk', 'lp2bs_zpk', 'is_stable']
 
 
 class BadCoefficients(UserWarning):
@@ -1022,8 +1022,7 @@ def tf2zpk(b, a):
 
     Notes
     -----
-    If some values of `b` are too close to 0, they are removed. In that case,
-    a BadCoefficients warning is emitted.
+    If leading values of `b` are zeros, they are removed.
 
     The `b` and `a` arrays are interpreted as coefficients for positive,
     descending powers of the transfer function variable.  So the inputs
@@ -1062,8 +1061,7 @@ def tf2zpk(b, a):
 
     """
     b, a = normalize(b, a)
-    b = (b + 0.0) / a[0]
-    a = (a + 0.0) / a[0]
+    b = np.trim_zeros(b, 'f')
     k = b[0]
     b /= b[0]
     z = roots(b)
@@ -1564,9 +1562,6 @@ def _align_nums(nums):
 def normalize(b, a):
     """Normalize numerator/denominator of a continuous-time transfer function.
 
-    If values of `b` are too close to 0, they are removed. In that case, a
-    BadCoefficients warning is emitted.
-
     Parameters
     ----------
     b: array_like
@@ -1608,30 +1603,66 @@ def normalize(b, a):
     # Normalize transfer function
     num, den = num / den[0], den / den[0]
 
-    # Count numerator columns that are all zero
-    leading_zeros = 0
-    for col in num.T:
-        if np.allclose(col, 0, atol=1e-14):
-            leading_zeros += 1
-        else:
-            break
-
-    # Trim leading zeros of numerator
-    if leading_zeros > 0:
-        warnings.warn("Badly conditioned filter coefficients (numerator): the "
-                      "results may be meaningless", BadCoefficients)
-        # Make sure at least one column remains
-        if leading_zeros == num.shape[1]:
-            leading_zeros -= 1
-        num = num[:, leading_zeros:]
-
     # Squeeze first dimension if singular
     if num.shape[0] == 1:
         num = num[0, :]
 
     return num, den
 
 
+def is_stable(den):
+    """Determine digital filter stability
+
+    Parameters
+    ----------
+    b: array_like
+        Numerator of the transfer function. Can be a 2d array to normalize
+        multiple transfer functions.
+    a: array_like
+        Denominator of the transfer function. At most 1d.
+
+    Returns
+    -------
+    stable: bool
+        Indicates whether all poles of the filter are stable.
+
+    Notes
+    -----
+    This function checks if all poles of the filter lie within the z-plane unit circle.
+    As direct root finding is numerically unreliable, this function uses step-down,
+    reverse Levinson recursion to determine the denominator polynomial's reflection
+    coefficients, as shown in [1]_. For a stable filter, all reflection coefficients
+    have magnitude less than unity.
+
+    References
+    ----------
+    .. [1] S. M. Kay, "Modern Spectral Estimation: Theory and Application",
+           O.P. EngleWood Cliffs, N.J., Ch. 6, p.172, Eq. 6.51
+
+    """
+
+    den = np.atleast_1d(den)
+
+    if den.ndim != 1:
+        raise ValueError("Denominator polynomial must be rank-1 array.")
+
+    if np.all(den == 0):
+        raise ValueError("Denominator must have at least on nonzero element.")
+
+    den = np.trim_zeros(den, 'f')
+    a = den / den[0]  # Normalize to first element
+
+    stable = True
+    for ii in reversed(range(1, len(a))):
+        k = a[ii]
+        if abs(k) >= 1:
+            stable = False
+            break
+        a = (a[:ii] - k * a[ii:0:-1]) / (1 - k ** 2)
+
+    return stable
+
+
 def lp2lp(b, a, wo=1.0):
     r"""
     Transform a lowpass filter prototype to a different frequency.

scipy/signal/tests/test_dltisys.py

@@ -11,7 +11,7 @@
 from scipy._lib._numpy_compat import suppress_warnings
 from scipy.signal import (dlsim, dstep, dimpulse, tf2zpk, lti, dlti,
                           StateSpace, TransferFunction, ZerosPolesGain,
-                          dfreqresp, dbode, BadCoefficients)
+                          dfreqresp, dbode)
 
 
 class TestDLTI(object):
@@ -470,10 +470,8 @@ def test_from_state_space(self):
 
         system_SS = dlti(A, B, C, D)
         w = 10.0**np.arange(-3,0,.5)
-        with suppress_warnings() as sup:
-            sup.filter(BadCoefficients)
-            w1, H1 = dfreqresp(system_TF, w=w)
-            w2, H2 = dfreqresp(system_SS, w=w)
+        w1, H1 = dfreqresp(system_TF, w=w)
+        w2, H2 = dfreqresp(system_SS, w=w)
 
         assert_almost_equal(H1, H2)
 

scipy/signal/tests/test_filter_design.py

@@ -1,7 +1,5 @@
 from __future__ import division, print_function, absolute_import
 
-import warnings
-
 from distutils.version import LooseVersion
 import numpy as np
 from numpy.testing import (assert_array_almost_equal,
@@ -10,10 +8,9 @@
                            assert_allclose, assert_warns)
 import pytest
 from pytest import raises as assert_raises
-from scipy._lib._numpy_compat import suppress_warnings
 
 from numpy import array, spacing, sin, pi, sort, sqrt
-from scipy.signal import (BadCoefficients, bessel, besselap, bilinear, buttap,
+from scipy.signal import (bessel, besselap, bilinear, buttap,
                           butter, buttord, cheb1ap, cheb1ord, cheb2ap,
                           cheb2ord, cheby1, cheby2, ellip, ellipap, ellipord,
                           firwin, freqs_zpk, freqs, freqz, freqz_zpk,
@@ -182,13 +179,6 @@ def test_simple(self, dt):
         assert_array_almost_equal(k, 1.)
         assert k.dtype == dt
 
-    def test_bad_filter(self):
-        # Regression test for #651: better handling of badly conditioned
-        # filter coefficients.
-        with suppress_warnings():
-            warnings.simplefilter("error", BadCoefficients)
-            assert_raises(BadCoefficients, tf2zpk, [1e-15], [1.0, 1.0])
-
 
 class TestZpk2Tf(object):
 
@@ -255,8 +245,7 @@ def test_fewer_zeros(self):
 
         sos = butter(12, [5., 30.], 'bandpass', fs=1200., analog=False,
                     output='sos')
-        with pytest.warns(BadCoefficients, match='Badly conditioned'):
-            z, p, k = sos2zpk(sos)
+        z, p, k = sos2zpk(sos)
         assert len(z) == 24
         assert len(p) == 24
 

scipy/signal/tests/test_ltisys.py

@@ -7,12 +7,10 @@
                            assert_)
 from pytest import raises as assert_raises
 
-from scipy._lib._numpy_compat import suppress_warnings
 from scipy.signal import (ss2tf, tf2ss, lsim2, impulse2, step2, lti,
                           dlti, bode, freqresp, lsim, impulse, step,
                           abcd_normalize, place_poles,
                           TransferFunction, StateSpace, ZerosPolesGain)
-from scipy.signal.filter_design import BadCoefficients
 import scipy.linalg as linalg
 from scipy.sparse.sputils import matrix
 
@@ -394,16 +392,10 @@ def test_multioutput(self):
 
 
 class TestLsim(object):
-    def lti_nowarn(self, *args):
-        with suppress_warnings() as sup:
-            sup.filter(BadCoefficients)
-            system = lti(*args)
-        return system
-
     def test_first_order(self):
         # y' = -y
         # exact solution is y(t) = exp(-t)
-        system = self.lti_nowarn(-1.,1.,1.,0.)
+        system = lti(-1.,1.,1.,0.)
         t = np.linspace(0,5)
         u = np.zeros_like(t)
         tout, y, x = lsim(system, u, t, X0=[1.0])
@@ -413,7 +405,7 @@ def test_first_order(self):
 
     def test_integrator(self):
         # integrator: y' = u
-        system = self.lti_nowarn(0., 1., 1., 0.)
+        system = lti(0., 1., 1., 0.)
         t = np.linspace(0,5)
         u = t
         tout, y, x = lsim(system, u, t)
@@ -426,7 +418,7 @@ def test_double_integrator(self):
         A = matrix("0. 1.; 0. 0.")
         B = matrix("0.; 1.")
         C = matrix("2. 0.")
-        system = self.lti_nowarn(A, B, C, 0.)
+        system = lti(A, B, C, 0.)
         t = np.linspace(0,5)
         u = np.ones_like(t)
         tout, y, x = lsim(system, u, t)
@@ -444,7 +436,7 @@ def test_jordan_block(self):
         A = matrix("-1. 1.; 0. -1.")
         B = matrix("0.; 1.")
         C = matrix("1. 0.")
-        system = self.lti_nowarn(A, B, C, 0.)
+        system = lti(A, B, C, 0.)
         t = np.linspace(0,5)
         u = np.zeros_like(t)
         tout, y, x = lsim(system, u, t, X0=[0.0, 1.0])
@@ -457,7 +449,7 @@ def test_miso(self):
         B = np.array([[1.0, 0.0], [0.0, 1.0]])
         C = np.array([1.0, 0.0])
         D = np.zeros((1,2))
-        system = self.lti_nowarn(A, B, C, D)
+        system = lti(A, B, C, D)
 
         t = np.linspace(0, 5.0, 101)
         u = np.zeros_like(t)
@@ -470,7 +462,7 @@ def test_miso(self):
         assert_almost_equal(x[:,1], expected_x1)
 
     def test_nonzero_initial_time(self):
-        system = self.lti_nowarn(-1.,1.,1.,0.)
+        system = lti(-1.,1.,1.,0.)
         t = np.linspace(1,2)
         u = np.zeros_like(t)
         tout, y, x = lsim(system, u, t, X0=[1.0])
@@ -520,21 +512,14 @@ def test_04(self):
         assert_almost_equal(x[:,0], expected_x)
 
     def test_05(self):
-        # The call to lsim2 triggers a "BadCoefficients" warning from
-        # scipy.signal.filter_design, but the test passes.  I think the warning
-        # is related to the incomplete handling of multi-input systems in
-        # scipy.signal.
-
         # A system with two state variables, two inputs, and one output.
         A = np.array([[-1.0, 0.0], [0.0, -2.0]])
         B = np.array([[1.0, 0.0], [0.0, 1.0]])
         C = np.array([1.0, 0.0])
         D = np.zeros((1, 2))
 
         t = np.linspace(0, 10.0, 101)
-        with suppress_warnings() as sup:
-            sup.filter(BadCoefficients)
-            tout, y, x = lsim2((A,B,C,D), T=t, X0=[1.0, 1.0])
+        tout, y, x = lsim2((A,B,C,D), T=t, X0=[1.0, 1.0])
         expected_y = np.exp(-tout)
         expected_x0 = np.exp(-tout)
         expected_x1 = np.exp(-2.0 * tout)
@@ -1177,10 +1162,8 @@ def test_from_state_space(self):
         B = np.array([[0.0], [0.0], [1.0]])
         C = np.array([[1.0, 0.0, 0.0]])
         D = np.array([[0.0]])
-        with suppress_warnings() as sup:
-            sup.filter(BadCoefficients)
-            system = lti(A, B, C, D)
-            w, mag, phase = bode(system, n=100)
+        system = lti(A, B, C, D)
+        w, mag, phase = bode(system, n=100)
 
         expected_magnitude = 20 * np.log10(np.sqrt(1.0 / (1.0 + w**6)))
         assert_almost_equal(mag, expected_magnitude)
@@ -1242,10 +1225,8 @@ def test_from_state_space(self):
         B = np.array([[0.0],[0.0],[1.0]])
         C = np.array([[1.0, 0.0, 0.0]])
         D = np.array([[0.0]])
-        with suppress_warnings() as sup:
-            sup.filter(BadCoefficients)
-            system = lti(A, B, C, D)
-            w, H = freqresp(system, n=100)
+        system = lti(A, B, C, D)
+        w, H = freqresp(system, n=100)
         s = w * 1j
         expected = (1.0 / (1.0 + 2*s + 2*s**2 + s**3))
         assert_almost_equal(H.real, expected.real)

scipy/signal/tests/test_signaltools.py

@@ -25,7 +25,7 @@
     fftconvolve, oaconvolve, choose_conv_method,
     hilbert, hilbert2, lfilter, lfilter_zi, filtfilt, butter, zpk2tf, zpk2sos,
     invres, invresz, vectorstrength, lfiltic, tf2sos, sosfilt, sosfiltfilt,
-    sosfilt_zi, tf2zpk, BadCoefficients, detrend, unique_roots, residue,
+    sosfilt_zi, tf2zpk, detrend, unique_roots, residue,
     residuez)
 from scipy.signal.windows import hann
 from scipy.signal.signaltools import (_filtfilt_gust, _compute_factors,
@@ -2341,14 +2341,10 @@ def test_shape(self):
         assert_equal(d1.shape, (30, 15))
 
     def test_phaseshift_FIR(self):
-        with suppress_warnings() as sup:
-            sup.filter(BadCoefficients, "Badly conditioned filter")
-            self._test_phaseshift(method='fir', zero_phase=False)
+        self._test_phaseshift(method='fir', zero_phase=False)
 
     def test_zero_phase_FIR(self):
-        with suppress_warnings() as sup:
-            sup.filter(BadCoefficients, "Badly conditioned filter")
-            self._test_phaseshift(method='fir', zero_phase=True)
+        self._test_phaseshift(method='fir', zero_phase=True)
 
     def test_phaseshift_IIR(self):
         self._test_phaseshift(method='iir', zero_phase=False)