BUG: signal.bilinear handles complex input, and strips leading zeros (scipy#22313)

* Add test for signal.bilinear for complex input

* BUG: signal.bilinear handles complex input

Fixes scipygh-13823.

* Add leading-zeros test for signal.bilinear

* BUG: trim leading zeros in inputs in signal.bilinear

* Rework of `signal.blinear` function.

* Simplified algorithm by utilizing NumPy's Polynomial class.
* Improved docstr and example.
* Added unit test to check for exceptions.

* Improve docstr of `signal.blinear` function by making the formulas a narrower by using subscript array notation. Also, some typos were removed.

* Reduce the width of the longest equation in the  docstr of `signal.blinear` a little more.

* Remove typos in function signal.bilinear

Co-authored-by: Rory Yorke <rory.yorke@gmail.com>

* Correct two typos in docst of `signal.bilinear`.

---------

Co-authored-by: Dietrich Brunn <Dietrich.Brunn@web.de>

Author: roryyorke

scipy/signal/_filter_design.py

@@ -2150,88 +2150,141 @@ def lp2bs(b, a, wo=1.0, bw=1.0):
 
 
 def bilinear(b, a, fs=1.0):
-    r"""
-    Return a digital IIR filter from an analog one using a bilinear transform.
-
-    Transform a set of poles and zeros from the analog s-plane to the digital
-    z-plane using Tustin's method, which substitutes ``2*fs*(z-1) / (z+1)`` for
-    ``s``, maintaining the shape of the frequency response.
+    r"""Calculate a digital IIR filter from an analog transfer function by utilizing
+    the bilinear transform.
 
     Parameters
     ----------
     b : array_like
-        Numerator of the analog filter transfer function.
+        Coefficients of the numerator polynomial of the analog transfer function in
+        form of a complex- or real-valued 1d array.
     a : array_like
-        Denominator of the analog filter transfer function.
+        Coefficients of the denominator polynomial of the analog transfer function in
+        form of a complex- or real-valued 1d array.
     fs : float
-        Sample rate, as ordinary frequency (e.g., hertz). No prewarping is
+        Sample rate, as ordinary frequency (e.g., hertz). No pre-warping is
         done in this function.
 
     Returns
     -------
-    b : ndarray
-        Numerator of the transformed digital filter transfer function.
-    a : ndarray
-        Denominator of the transformed digital filter transfer function.
+    beta : ndarray
+        Coefficients of the numerator polynomial of the digital transfer function in
+        form of a complex- or real-valued 1d array.
+    alpha : ndarray
+        Coefficients of the denominator polynomial of the digital transfer function in
+        form of a complex- or real-valued 1d array.
+
+    Notes
+    -----
+    The parameters :math:`b = [b_0, \ldots, b_Q]` and :math:`a = [a_0, \ldots, a_P]`
+    are 1d arrays of length :math:`Q+1` and :math:`P+1`. They define the analog
+    transfer function
+
+    .. math::
+
+        H_a(s) = \frac{b_0 s^Q + b_1 s^{Q-1} + \cdots + b_Q}{
+                       a_0 s^P + a_1 s^{P-1} + \cdots + a_P}\ .
+
+    The bilinear transform [1]_ is applied by substituting
+
+    .. math::
+
+        s = \kappa \frac{z-1}{z+1}\ , \qquad \kappa := 2 f_s\ ,
+
+    into :math:`H_a(s)`, with :math:`f_s`  being the sampling rate.
+    This results in the digital transfer function in the :math:`z`-domain
+
+    .. math::
+
+        H_d(z) = \frac{b_0 \left(\kappa \frac{z-1}{z+1}\right)^Q +
+                       b_1 \left(\kappa \frac{z-1}{z+1}\right)^{Q-1} +
+                       \cdots + b_Q}{
+                       a_0 \left(\kappa \frac{z-1}{z+1}\right)^P +
+                       a_1 \left(\kappa \frac{z-1}{z+1}\right)^{P-1} +
+                       \cdots + a_P}\ .
+
+    This expression can be simplified by multiplying numerator and denominator by
+    :math:`(z+1)^N`, with :math:`N=\max(P, Q)`. This allows :math:`H_d(z)` to be
+    reformulated as
+
+    .. math::
+
+       & & \frac{b_0 \big(\kappa (z-1)\big)^Q     (z+1)^{N-Q} +
+                b_1 \big(\kappa (z-1)\big)^{Q-1} (z+1)^{N-Q+1} +
+                \cdots + b_Q(z+1)^N}{
+                a_0 \big(\kappa (z-1)\big)^P     (z+1)^{N-P} +
+                a_1 \big(\kappa (z-1)\big)^{P-1} (z+1)^{N-P+1} +
+                \cdots + a_P(z+1)^N}\\
+        &=:& \frac{\beta_0 + \beta_1  z^{-1} + \cdots + \beta_N  z^{-N}}{
+                 \alpha_0 + \alpha_1 z^{-1} + \cdots + \alpha_N z^{-N}}\ .
+
+
+    This is the equation implemented to perform the bilinear transform. Note that for
+    large :math:`f_s`, :math:`\kappa^Q` or :math:`\kappa^P` can cause a numeric
+    overflow for sufficiently large :math:`P` or :math:`Q`.
+
+    References
+    ----------
+    .. [1] "Bilinear Transform", Wikipedia,
+           https://en.wikipedia.org/wiki/Bilinear_transform
 
     See Also
     --------
-    lp2lp, lp2hp, lp2bp, lp2bs
-    bilinear_zpk
+    lp2lp, lp2hp, lp2bp, lp2bs, bilinear_zpk
 
     Examples
     --------
+    The following example shows the frequency response of an analog bandpass filter and
+    the corresponding digital filter derived by utilitzing the bilinear transform:
+
     >>> from scipy import signal
     >>> import matplotlib.pyplot as plt
     >>> import numpy as np
+    ...
+    >>> fs = 100  # sampling frequency
+    >>> om_c = 2 * np.pi * np.array([7, 13])  # corner frequencies
+    >>> bb_s, aa_s = signal.butter(4, om_c, btype='bandpass', analog=True, output='ba')
+    >>> bb_z, aa_z = signal.bilinear(bb_s, aa_s, fs)
+    ...
+    >>> w_z, H_z = signal.freqz(bb_z, aa_z)  # frequency response of digitial filter
+    >>> w_s, H_s = signal.freqs(bb_s, aa_s, worN=w_z*fs)  # analog filter response
+    ...
+    >>> f_z, f_s = w_z * fs / (2*np.pi), w_s / (2*np.pi)
+    >>> Hz_dB, Hs_dB = (20*np.log10(np.abs(H_).clip(1e-10)) for H_ in (H_z, H_s))
+    >>> fg0, ax0 = plt.subplots()
+    >>> ax0.set_title("Frequency Response of 4-th order Bandpass Filter")
+    >>> ax0.set(xlabel='Frequency $f$ in Hertz', ylabel='Magnitude in dB',
+    ...         xlim=[f_z[1], fs/2], ylim=[-200, 2])
+    >>> ax0.semilogx(f_z, Hz_dB, alpha=.5, label=r'$|H_z(e^{j 2 \pi f})|$')
+    >>> ax0.semilogx(f_s, Hs_dB, alpha=.5, label=r'$|H_s(j 2 \pi f)|$')
+    >>> ax0.legend()
+    >>> ax0.grid(which='both', axis='x')
+    >>> ax0.grid(which='major', axis='y')
+    >>> plt.show()
 
-    >>> fs = 100
-    >>> bf = 2 * np.pi * np.array([7, 13])
-    >>> filts = signal.lti(*signal.butter(4, bf, btype='bandpass',
-    ...                                   analog=True))
-    >>> filtz = signal.lti(*signal.bilinear(filts.num, filts.den, fs))
-    >>> wz, hz = signal.freqz(filtz.num, filtz.den)
-    >>> ws, hs = signal.freqs(filts.num, filts.den, worN=fs*wz)
-
-    >>> plt.semilogx(wz*fs/(2*np.pi), 20*np.log10(np.abs(hz).clip(1e-15)),
-    ...              label=r'$|H_z(e^{j \omega})|$')
-    >>> plt.semilogx(wz*fs/(2*np.pi), 20*np.log10(np.abs(hs).clip(1e-15)),
-    ...              label=r'$|H(j \omega)|$')
-    >>> plt.legend()
-    >>> plt.xlabel('Frequency [Hz]')
-    >>> plt.ylabel('Amplitude [dB]')
-    >>> plt.grid(True)
+    The difference in the higher frequencies shown in the plot is caused by an effect
+    called "frequency warping". [1]_ describes a method called "pre-warping" to
+    reduce those deviations.
     """
+    b, a = np.atleast_1d(b), np.atleast_1d(a)  # convert scalars, if needed
+    if not a.ndim == 1:
+        raise ValueError(f"Parameter a is not a 1d array since {a.shape=}")
+    if not b.ndim == 1:
+        raise ValueError(f"Parameter b is not a 1d array since {b.shape=}")
+    b, a = np.trim_zeros(b, 'f'), np.trim_zeros(a, 'f')  # remove leading zeros
     fs = _validate_fs(fs, allow_none=False)
-    a, b = map(atleast_1d, (a, b))
-    D = len(a) - 1
-    N = len(b) - 1
-    artype = float
-    M = max([N, D])
-    Np = M
-    Dp = M
-    bprime = np.empty(Np + 1, artype)
-    aprime = np.empty(Dp + 1, artype)
-    for j in range(Np + 1):
-        val = 0.0
-        for i in range(N + 1):
-            for k in range(i + 1):
-                for l in range(M - i + 1):
-                    if k + l == j:
-                        val += (comb(i, k) * comb(M - i, l) * b[N - i] *
-                                pow(2 * fs, i) * (-1) ** k)
-        bprime[j] = real(val)
-    for j in range(Dp + 1):
-        val = 0.0
-        for i in range(D + 1):
-            for k in range(i + 1):
-                for l in range(M - i + 1):
-                    if k + l == j:
-                        val += (comb(i, k) * comb(M - i, l) * a[D - i] *
-                                pow(2 * fs, i) * (-1) ** k)
-        aprime[j] = real(val)
 
-    return normalize(bprime, aprime)
+    # Splitting the factor fs*2 between numerator and denominator reduces the chance of
+    # numeric overflow for large fs and large N:
+    fac = np.sqrt(fs*2)
+    zp1 = np.polynomial.Polynomial((+1, 1)) / fac  # Polynomial (z + 1) / fac
+    zm1 = np.polynomial.Polynomial((-1, 1)) * fac  # Polynomial (z - 1) * fac
+    # Note that NumPy's Polynomial coefficient order is backward compared to a and b.
+
+    N = max(len(a), len(b)) - 1
+    numerator   = sum(b_ * zp1**(N-q) * zm1**q for q, b_ in enumerate(b[::-1]))
+    denominator = sum(a_ * zp1**(N-p) * zm1**p for p, a_ in enumerate(a[::-1]))
+    return normalize(numerator.coef[::-1], denominator.coef[::-1])
 
 
 def _validate_gpass_gstop(gpass, gstop):

scipy/signal/tests/test_filter_design.py

@@ -1,5 +1,7 @@
 import warnings
 
+from itertools import product
+
 from scipy._lib import _pep440
 import numpy as np
 from numpy.testing import (
@@ -1389,22 +1391,84 @@ def test_basic(self):
 
 
 class TestBilinear:
+    """Tests for function `signal.bilinear`. """
+
+    def test_exceptions(self):
+        """Raise all exceptions in `bilinear()`. """
+        with pytest.raises(ValueError, match="Parameter a is not .*"):
+            bilinear(1., np.array([[1, 2, 3]]))
+        with pytest.raises(ValueError, match="Parameter b is not .*"):
+            bilinear(np.ones((2,3)), 1. )
 
     def test_basic(self):
+        # reference output values computed with sympy
         b = [0.14879732743343033]
         a = [1, 0.54552236880522209, 0.14879732743343033]
+        b_zref = [0.08782128175913713, 0.17564256351827426, 0.08782128175913713]
+        a_zref = [1.0, -1.0047722097030667, 0.3560573367396151]
+
         b_z, a_z = bilinear(b, a, 0.5)
-        assert_array_almost_equal(b_z, [0.087821, 0.17564, 0.087821],
-                                  decimal=5)
-        assert_array_almost_equal(a_z, [1, -1.0048, 0.35606], decimal=4)
+
+        assert_array_almost_equal_nulp(b_z, b_zref)
+        assert_array_almost_equal_nulp(a_z, a_zref)
 
         b = [1, 0, 0.17407467530697837]
         a = [1, 0.18460575326152251, 0.17407467530697837]
+        b_zref = [0.8641286432189045, -1.2157757001964216, 0.8641286432189045]
+        a_zref = [1.0, -1.2157757001964216, 0.7282572864378091]
+
         b_z, a_z = bilinear(b, a, 0.5)
-        assert_array_almost_equal(b_z, [0.86413, -1.2158, 0.86413],
-                                  decimal=4)
-        assert_array_almost_equal(a_z, [1, -1.2158, 0.72826],
-                                  decimal=4)
+
+        assert_array_almost_equal_nulp(b_z, b_zref)
+        assert_array_almost_equal_nulp(a_z, a_zref)
+
+
+    def test_ignore_leading_zeros(self):
+        # regression for gh-6606
+        # results shouldn't change when leading zeros are added to
+        # input numerator or denominator
+        b = [0.14879732743343033]
+        a = [1, 0.54552236880522209, 0.14879732743343033]
+
+        b_zref = [0.08782128175913713, 0.17564256351827426, 0.08782128175913713]
+        a_zref = [1.0, -1.0047722097030667, 0.3560573367396151]
+
+        for lzn, lzd in product(range(4), range(4)):
+            b_z, a_z = bilinear(np.pad(b, (lzn, 0)),
+                                np.pad(a, (lzd, 0)),
+                                0.5)
+            assert_array_almost_equal_nulp(b_z, b_zref)
+            assert_array_almost_equal_nulp(a_z, a_zref)
+
+
+    def test_complex(self):
+        # reference output values computed with sympy
+        # this is an elliptical filter, 5Hz width, centered at +50Hz:
+        #     z, p, k = signal.ellip(2, 0.5, 20, 2*np.pi*5/2, output='zpk', analog=True)
+        #     z = z.astype(complex) + 2j * np.pi * 50
+        #     p = p.astype(complex) + 2j * np.pi * 50
+        #     b, a = signal.zpk2tf(z, p, k)
+        b = [(0.09999999999999991+0j),
+             -62.831853071795805j,
+             (-9505.857007071314+0j)]
+        a = [(1+0j),
+             (21.09511000343942-628.3185307179587j),
+             (-98310.74322875646-6627.2242613473845j)]
+        # sample frequency
+        fs = 1000
+        b_zref = [(0.09905575106715676-0.00013441423112828688j),
+                  (-0.18834281923181084-0.06032810039049478j),
+                  (0.08054306669414343+0.05766172295523972j)]
+        a_zref = [(1+0j),
+                  (-1.8839476369292854-0.606808151331815j),
+                  (0.7954687330018285+0.5717459398142481j)]
+
+        b_z, a_z = bilinear(b, a, fs)
+
+        # the 3 ulp difference determined from testing
+        assert_array_almost_equal_nulp(b_z, b_zref, 3)
+        assert_array_almost_equal_nulp(a_z, a_zref, 3)
+
 
     def test_fs_validation(self):
         b = [0.14879732743343033]