APGF (all pole gammatone filterbank) implementation in Python

gammatone_apgf.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""
/*
** Copyright (C) 2020 Ricard Marxer
**
** This program is free software; you can redistribute it and/or modify
** it under the terms of the GNU General Public License as published by
** the Free Software Foundation; either version 3 of the License, or
** (at your option) any later version.
**
** This program is distributed in the hope that it will be useful,
** but WITHOUT ANY WARRANTY; without even the implied warranty of
** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
** GNU General Public License for more details.
**
** You should have received a copy of the GNU General Public License
** along with this program; if not, write to the Free Software
** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA.
*/

"""

import numpy as np
from scipy.special import factorial
import matplotlib.pyplot as plt
import scipy.signal

"""Create a Gammatone AGPF in Python.

 * Based on:
 *
 * http://www.dicklyon.com/tech/Hearing/APGF_Lyon_1996.pdf
 *
 * Implementing a GammaTone Filter Bank
 * Holdsworth et al. 1988
 *
"""


def linear_to_erb(x):
    return 21.4 * np.log10(4.37e-3 * x + 1.0)


def erb_to_linear(x):
    return (10.0 ** (x / 21.4) - 1.0) / 4.37e-3


def linear_to_erb_bandwidth(x):
    return (x * 4.37e-3 + 1.0) * 24.7


def gammatonepy(sample_rate,
                low_frequency=50.,
                high_frequency=3500.,
                channel_count=32,
                peak_phase=True,
                ):
    # Implemented using SOS filters so order is forced to be 4
    order = 4

    # Compute the center frequencies
    low_erb = linear_to_erb(low_frequency)
    high_erb = linear_to_erb(high_frequency)

    center_erbs = np.linspace(low_erb, high_erb, channel_count)
    _centerFrequencies = erb_to_linear(center_erbs)

    two_pi = 2 * np.pi

    # Prepare filters parameters
    an_inverse = pow(factorial(order - 1), 2) / \
                 (np.pi * factorial(2 * order - 2) * pow(2, -(2 * order - 2)))

    _bandwidths = linear_to_erb_bandwidth(_centerFrequencies)
    _bandwidths *= an_inverse
    delay = 3 / (two_pi * _bandwidths)

    # Compute the poles positions and the phases
    phi = _bandwidths * (two_pi / sample_rate)
    theta = _centerFrequencies * (two_pi / sample_rate)
    atilde = np.exp(-phi + theta * 1j)

    btmp = (1 - np.exp(-phi))
    b2 = (btmp ** order).astype(np.complex128)
    if peak_phase:
        b2 *= np.exp(_centerFrequencies.astype(np.complex128) * delay.astype(np.complex128) * (two_pi * 1j))

    poles = np.repeat(atilde[:, np.newaxis], order / 2, axis=1)
    a = np.apply_along_axis(np.poly, arr=poles, axis=1)

    b = np.c_[b2, np.zeros_like(b2), np.zeros_like(b2)]

    sos = np.c_[b, a]
    sos = sos[:, np.newaxis, :]
    sos = sos.repeat(2, axis=1)

    # Set the numerator of the second filter to 1s
    sos[:, 1, 0] = 1

    # Filter an impulse response
    x = scipy.signal.unit_impulse(700, dtype=np.complex128)
    y = np.zeros((sos.shape[0], x.shape[0]), dtype=np.complex128)
    for i in range(sos.shape[0]):
        y[i] = scipy.signal.sosfilt(sos[i], x)

    plt.figure()
    plt.plot(y[[0, 15, 31], :].T)

    # Plot frequency response
    plt.figure()
    for i in range(sos.shape[0]):
        w, h = scipy.signal.sosfreqz(sos[i], worN=1500)
        plt.subplot(2, 1, 1)
        db = 20 * np.log10(np.maximum(np.abs(h), 1e-5))
        plt.plot(w / np.pi, db)

        plt.subplot(2, 1, 2)
        plt.plot(w / np.pi, np.angle(h))

    plt.subplot(2, 1, 1)
    plt.ylim(-80, 5)
    plt.grid(True)
    plt.yticks([0, -20, -40, -60])
    plt.ylabel('Gain [dB]')
    plt.title('Frequency Response')

    plt.subplot(2, 1, 2)
    plt.grid(True)
    plt.yticks([-np.pi, -0.5 * np.pi, 0, 0.5 * np.pi, np.pi], [r'$-\pi$', r'$-\pi/2$', '0', r'$\pi/2$', r'$\pi$'])
    plt.ylabel('Phase [rad]')
    plt.xlabel('Normalized frequency (1.0 = Nyquist)')

    plt.show()


def main():
    gammatonepy(16000)


if __name__ == "__main__":
    main()