code: def bilinear_zpk(zeros, poles, gain, fs=1.0, unit='Hz'):
      
      zeros = numpy.array(zeros, dtype=float, copy=False)
      zeros = zeros[numpy.isfinite(zeros)]
      poles = numpy.array(poles, dtype=float, copy=False)
      gain = gain
  
      # convert from Hz to rad/s if needed
      unit = Unit(unit)
      if unit == Unit('Hz'):
          zeros *= -2 * pi
          poles *= -2 * pi
      elif unit != Unit('rad/s'):
          raise ValueError("zpk can only be given with unit='Hz' "
                           "or 'rad/s'")
  
      # convert to Z-domain via bilinear transform
      fs = 2 * Quantity(fs, 'Hz').value
      dpoles = (1 + poles/fs) / (1 - poles/fs)
      dzeros = (1 + zeros/fs) / (1 - zeros/fs)
      dzeros = numpy.concatenate((
          dzeros, -numpy.ones(len(dpoles) - len(dzeros)),
      ))
      dgain = gain * numpy.prod(fs - zeros)/numpy.prod(fs - poles)
      return dzeros, dpoles, dgain
docs: Convert an analogue ZPK filter to digital using a bilinear transform
  
      Parameters
      ----------
      zeros : array-like
          list of zeros
  
      poles : array-like
          list of poles
  
      gain : `float`
          filter gain
  
      fs : `float`, `~astropy.units.Quantity`
          sampling rate at which to evaluate bilinear transform, default: 1.
  
      unit : `str`, `~astropy.units.Unit`
          unit of inputs, one or 'Hz' or 'rad/s', default: ``'Hz'``
  
      Returns
      -------
      zpk : `tuple`
          digital version of input zpk