def gammatone(freq, bandwidth, phase=0, eta=4):
  """
    ``Bellini, D. J. S. "AudioLazy: Processamento digital de sinais
    expressivo e em tempo real", IME-USP, Mastership Thesis, 2013.``

  This implementation have the impulse response (for each sample ``n``,
  keeping the input parameter names):

  .. math::

    n^{{eta - 1}} e^{{- bandwidth \cdot n}} \cos(freq \cdot n + phase)
  """
  assert eta >= 1

  A = exp(-bandwidth)
  numerator = cos(phase) - A * cos(freq - phase) * z ** -1
  denominator = 1 - 2 * A * cos(freq) * z ** -1 + A ** 2 * z ** -2
  filt = (numerator / denominator).diff(n=eta-1, mul_after=-z)

  # Filter is done, but the denominator might have some numeric loss
  f0 = ZFilter(filt.numpoly) / denominator
  f0 /= abs(f0.freq_response(freq)) # Max gain == 1.0 (0 dB)
  fn = 1 / denominator
  fn /= abs(fn.freq_response(freq))
  return CascadeFilter([f0] + [fn] * (eta - 1))