code: def residue(b, a, tol=1e-3, rtype="avg"):
      """
      Compute partial-fraction expansion of b(s) / a(s).
  
      If `M` is the degree of numerator `b` and `N` the degree of denominator
      `a`::
  
                b(s)     b[0] s**(M) + b[1] s**(M-1) + ... + b[M]
        H(s) = ------ = ------------------------------------------
                a(s)     a[0] s**(N) + a[1] s**(N-1) + ... + a[N]
  
      then the partial-fraction expansion H(s) is defined as::
  
                 r[0]       r[1]             r[-1]
             = -------- + -------- + ... + --------- + k(s)
               (s-p[0])   (s-p[1])         (s-p[-1])
  
      If there are any repeated roots (closer together than `tol`), then H(s)
      has terms like::
  
            r[i]      r[i+1]              r[i+n-1]
          -------- + ----------- + ... + -----------
          (s-p[i])  (s-p[i])**2          (s-p[i])**n
  
      This function is used for polynomials in positive powers of s or z,
      such as analog filters or digital filters in controls engineering.  For
      negative powers of z (typical for digital filters in DSP), use `residuez`.
  
      Parameters
      ----------
      b : array_like
          Numerator polynomial coefficients.
      a : array_like
          Denominator polynomial coefficients.
  
      Returns
      -------
      r : ndarray
          Residues.
      p : ndarray
          Poles.
      k : ndarray
          Coefficients of the direct polynomial term.
  
      See Also
      --------
      invres, residuez, numpy.poly, unique_roots
  
      """
  
      b, a = map(asarray, (b, a))
      rscale = a[0]
      k, b = polydiv(b, a)
      p = roots(a)
      r = p * 0.0
      pout, mult = unique_roots(p, tol=tol, rtype=rtype)
      p = []
      for n in range(len(pout)):
          p.extend([pout[n]] * mult[n])
      p = asarray(p)
      # Compute the residue from the general formula
      indx = 0
      for n in range(len(pout)):
          bn = b.copy()
          pn = []
          for l in range(len(pout)):
              if l != n:
                  pn.extend([pout[l]] * mult[l])
          an = atleast_1d(poly(pn))
          # bn(s) / an(s) is (s-po[n])**Nn * b(s) / a(s) where Nn is
          # multiplicity of pole at po[n]
          sig = mult[n]
          for m in range(sig, 0, -1):
              if sig > m:
                  # compute next derivative of bn(s) / an(s)
                  term1 = polymul(polyder(bn, 1), an)
                  term2 = polymul(bn, polyder(an, 1))
                  bn = polysub(term1, term2)
                  an = polymul(an, an)
              r[indx + m - 1] = (
                  polyval(bn, pout[n])
                  / polyval(an, pout[n])
                  / factorial(sig - m)
              )
          indx += sig
      return r / rscale, p, k
  
imports: from __future__ import division, print_function, absolute_import
  import operator
  import threading
  import sys
  import timeit
  from scipy.spatial import cKDTree
  from .None import sigtools, dlti
  from ._upfirdn import upfirdn, _output_len
  from scipy._lib.six import callable
  from scipy._lib._version import NumpyVersion
  from scipy import fftpack, linalg
  from scipy.fftpack.helper import _init_nd_shape_and_axes_sorted
  from numpy import allclose, angle, arange, argsort, array, asarray, atleast_1d, atleast_2d, cast, dot, exp, expand_dims, iscomplexobj, mean, ndarray, newaxis, ones, pi, poly, polyadd, polyder, polydiv, polymul, polysub, polyval, product, r_, ravel, real_if_close, reshape, roots, sort, take, transpose, unique, where, zeros, zeros_like
  import numpy
  import math
  from scipy.special import factorial
  from .windows import get_window
  from ._arraytools import axis_slice, axis_reverse, odd_ext, even_ext, const_ext
  from .filter_design import cheby1, _validate_sos
  from .fir_filter_design import firwin
label_desc: Use image node_id 29 for calling a global function with example usage: residue(b, a, tol, rtype) and returns: unknown,  p,  k