code: def gammatone(freq, bandwidth, phase=0, eta=4):
    
    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))
docs: ``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)