Repository: scipy/scipy
Pull request: BUG: signal.bilinear handles complex input, and strips leading zeros
#### Reference issue
Closes gh-13823
Closes gh-6606.

#### What does this implement/fix?
This handles complex numerators and denominators in `signal.bilinear`.  While searching for this issue I spotted gh-6606, which is simple to fix.  I'll split this into a separate PR if that's desired.

#### Additional information
It would be useful to have a reference for the algorithm, which I found non-obvious at first look.  If we can't find a reference, maybe a hint in the comments, say  "after substituting in the bilinear transform, multiply by $(z+1)^k/(z+1)^k$, where k is the maximum of the degree of numerator and denominator; then apply the binomial expansion to $(z-1)^i$ and $(z+1)^j$ terms".

I implemented a Sympy version of the algorithm to generate reference values (see below).  Since the operations involved are "basic" arithmetic (addition, subtraction, multiplication, division), computing with rational numbers for the non-complex case gives exact results, though I've left the allowed error in the test at the default 1 ULP.  (Previously these tests were to 5 decimal places.)  (Please challenge this paragraph!  I'm not a numerical analyst.)

For the complex test I added, the ULP difference gets to about 2.2 (a single metric is computed for the complex difference, so one gets non-integer ULP).  I'm not sure why it's not also 0 - I think all operations should still be "basic" arithmetic.

#### Motivating example

This applies the bilinear transform to a filter with passband at +50Hz.  The example is somewhat artificial: in this case one should use `signal.bilinear_zpk`, but in general the S-plane numerator and denominator could be computed by other means, or could represent a physical process (e.g., three-wire three-phase circuit).

```python
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal

def complex_filt():
    # centre frequency for filter [Hz]
    fcentre = 50
    # filter passband width [Hz]
    fwidth = 5

    # generate via ZPK
    z, p, k = signal.ellip(2, 0.5, 40, 2*np.pi*fwidth/2, output='zpk', analog=True)

    z = z.astype(complex)
    p = p.astype(complex)

    z += 2j * np.pi * fcentre
    p += 2j * np.pi * fcentre

    b, a = signal.zpk2tf(z, p, k)

    return b, a


def plot_filt():
    # sample frequency [Hz]
    fs = 1000

    b, a = complex_filt()
    bz, az = signal.bilinear(b, a, fs)

    fpos = np.geomspace(5, 5e2, 801)

    # show suppression at -50Hz, passband at +50Hz
    f = np.concat((-fpos[::-1], fpos))

    resps = signal.freqs(b, a, 2*np.pi*f)[1]
    respz = signal.freqz(bz, az, f, fs=fs)[1]

    fig, axs = plt.subplots(2, sharex=True, constrained_layout=True)

    axs[0].plot(f, abs(resps), label='s', lw=5)
    axs[0].plot(f, abs(respz), label='z')
    axs[0].set_xscale('symlog')
    axs[0].set_yscale('log')
    axs[0].set_ylabel('mag [1]')
    axs[0].grid()

    axs[1].plot(f, np.rad2deg(np.angle(resps)), lw=5)
    axs[1].plot(f, np.rad2deg(np.angle(respz)))
    axs[1].set_ylabel('phase [deg]')
    axs[1].grid()

    axs[1].set_xlabel('freq [hz]')

    fig.legend(loc='outside right')
    fig.suptitle('bilinear mapping of +50Hz-centered cheby. filter')


if __name__ == '__main__':
    plot_filt()
    plt.show()
```

![image](https://github.com/user-attachments/assets/4e16bb34-7581-406e-a66b-6e57b9546b4b)


#### Sympy code

Run with Sympy v1.12.1.

```python
import sympy

def bilinear_poly(coeffs, k, fs):
    """
    coeffs: sequence of coefficients in decreasing powers
    k: max degree of num, den; must have len(coeffs) <= k+1
    fs: scalar
    """
    if not (len(coeffs) <= k + 1):
        raise ValueError(f'require {len(coeffs)=} <= {k+1=}')
    
    z = sympy.Symbol('z')
    dz = sympy.Number(0)
    for j, coeff in enumerate(coeffs[::-1]):
        pre = coeff * 2**j * fs**j
        for p in range(j+1):
            for q in range(k-j+1):
                # comment edit: deleted a debug line here
                dz += pre * sympy.binomial(j, p) * sympy.binomial(k-j, q) * (-1)**p * z**(k-p-q)

    return dz.as_poly(z).all_coeffs()


def bilinear(nscoeffs, dscoeffs, fs=None):
    """Bilinear transform of rational function"""
    if fs is None:
        fs = sympy.Symbol('f_s')

    k = max(len(nscoeffs)-1, len(dscoeffs)-1)

    nzcoeffs = bilinear_poly(nscoeffs, k, fs)
    dzcoeffs = bilinear_poly(dscoeffs, k, fs)

    norm_nzcoeffs = [c/dzcoeffs[0] for c in nzcoeffs]
    norm_dzcoeffs = [c/dzcoeffs[0] for c in dzcoeffs]

    # leave normalization for caller
    return norm_nzcoeffs, norm_dzcoeffs


def torat(x):
    """Converter float or complex to sympy rational (or 'complex rational')"""
    return sympy.Rational(x.real) + sympy.I * sympy.Rational(x.imag)


def tonumeric(x):
    """Converter sympy result to float or complex"""
    if x.is_real:
        return float(x)
    else:
        return complex(x)


def bilinear_numeric(b, a, fs, to_float=True):
    """Compute bilinear transformation with rational numbers"""
    b = [torat(b_) for b_ in b]
    a = [torat(a_) for a_ in a]
    fs = sympy.Rational(fs)

    bzs, azs = bilinear(b, a, fs)

    if to_float:
        return ([tonumeric(bz) for bz in bzs],
                [tonumeric(az) for az in azs])
    else:
        return bzs, azs


def demo_bilinear_poly(m=2, k=None):
    if k is None:
        k = m
    coeffs = sympy.symbols(f'a:{m}')[::-1]
    fs = sympy.Symbol('fs')
    print(bilinear_poly(coeffs, k, fs))


if __name__ == '__main__':
    demo_bilinear_poly()
    demo_bilinear_poly(3, 4)

    b= [0.14879732743343033]
    a= [1, 0.54552236880522209, 0.14879732743343033]
    fs = 0.5

    print(f'{bilinear_numeric(b, a, fs, False)=}')
    print(f'{bilinear_numeric(b, a, fs, True)=}')
```

Output:
```
[a0 + 2*a1*fs, 2*a0, a0 - 2*a1*fs]
[a0 + 2*a1*fs + 4*a2*fs**2, 4*a0 + 4*a1*fs, 6*a0 - 8*a2*fs**2, 4*a0 - 4*a1*fs, a0 - 2*a1*fs + 4*a2*fs**2]
bilinear_numeric(b, a, fs, False)=([1340247176765845/15261075105253774, 1340247176765845/7630537552626887, 1340247176765845/15261075105253774], [1, -7666952077975147/7630537552626887, 2716908878879950/7630537552626887])
bilinear_numeric(b, a, fs, True)=([0.08782128175913713, 0.17564256351827426, 0.08782128175913713], [1.0, -1.0047722097030667, 0.3560573367396151])
```

Problem statement:
BUG: signal.bilinear: can't handle leading zeros
It looks to me that bilinear(b, a) can't handle cases where there is leading zeros in arrays a and b. Here is some examples:

```
In [57]: signal.bilinear(1, [2], 1)
Out[57]: (array([ 0.5]), array([ 1.]))

In [58]: signal.bilinear(1, [0, 2], 1)
Out[58]: (array([ 0.5,  0.5]), array([ 1.,  1.]))

In [59]: signal.bilinear([0,1], [0, 2], 1)
Out[59]: (array([ 0.5,  0.5]), array([ 1.,  1.]))

In [60]: signal.bilinear([0,1], [2], 1)
Out[60]: (array([ 0.5,  0.5]), array([ 1.,  1.]))
```

What is expected that all return the same answer as first one. One may not review it as an issue. However, in cases that b and a are computed on-the-fly, one can't predict when we have leading zeros because of limited computational precision. IMO, it's better removing the leading zeros inside bilinear()


BUG: signal.bilinear: doesn't work for complex valued arrays
scipy.signal.bilinear(b,a,fs) gives incorrect ouput when input arrays (b,a) are complex and doesn't give any warning about ignoring the imaginary part.
I am also providing the function to fix this and implement bilinear with complex values.
```
def poly_power(a,n):
    result = np.array([1])
    for i in range(n):
        result = np.convolve(result,a)
    return result
def my_bilinear(b,a,fs):
    B = np.zeros((b.shape),dtype = complex)
    A = np.zeros((a.shape),dtype = complex)
    s1 = np.array([1,-1])
    s2 = np.array([1,1])
    num = np.zeros(b.shape)
    den = np.zeros(a.shape)
    n = len(b) - 1
    for i in range(len(b)):
        B[i] = b[i]*((2*fs)**(n-i))
        num = num + B[i]*np.convolve(poly_power(s1,n-i),poly_power(s2,i))
    n = len(a) -1
    for i in range(len(a)):
        A[i] = a[i]*((2*fs)**(n-i))
        den = den + A[i]*np.convolve(poly_power(s1,n-i),poly_power(s2,i))
    if(len(b) > len(a)):
        den = np.convolve(poly_power(s2,len(b) - len(a)),den)
    if(len(a) > len(b)):
        num = np.convolve(poly_power(s2,len(a) - len(b)),num)
    num = num/den[0]
    den = den/den[0]
    return num,den
```

Hints from issue discussion:
The result is not strictly wrong it is just not a minimal realization. So I don't think it actually promises a result with minimal realization anyways.

But until the developers chime in, you can also pass them through [`numpy.trim_zeros`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.trim_zeros.html)

```
sig.bilinear(*map(np.trim_zeros,([0,1], [0, 2])), 1)
(array([ 0.5]), array([ 1.]))
```

There exists the function `scipy.signal.normalize`, which trims leading zeros from a `(b, a)` transfer function; it is automatically used in the construction of a `TransferFunction` object. If you'd like `bilinear` to call this automatically on its input, a PR is welcome!

Commit message: Reduce the width of the longest equation in the  docstr of `signal.blinear` a little more.

Patch to review:
diff --git a/scipy/signal/_filter_design.py b/scipy/signal/_filter_design.py
index 459d7087c786..b6a49b29c792 100644
--- a/scipy/signal/_filter_design.py
+++ b/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 denominatory 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_P(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
+    compensate 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):
diff --git a/scipy/signal/tests/test_filter_design.py b/scipy/signal/tests/test_filter_design.py
index 863aa2bf2e50..0e8d31f4dbba 100644
--- a/scipy/signal/tests/test_filter_design.py
+++ b/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]


Code review comments:
diff_hunk: @@ -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 denominatory by
path: scipy/signal/_filter_design.py
text: @author:
```suggestion
    This expression can be simplified by multiplying numerator and denominator by
```

diff_hunk: @@ -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
path: scipy/signal/_filter_design.py
text: @author:
```suggestion
        Sample rate, as ordinary frequency (e.g., hertz). No pre-warping is
```

see, e.g., https://www.nist.gov/pml/owm/writing-si-metric-system-units .

diff_hunk: @@ -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 denominatory 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_P(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
+    compensate those deviations.
path: scipy/signal/_filter_design.py
text: @author:
```suggestion
    compensate for those deviations.
```

"reduce those deviations" ?

@user1:
"reduce" is better