mi@0: # Copyright 2014 Jason Heeris, jason.heeris@gmail.com
mi@0: # 
mi@0: # This file is part of the gammatone toolkit, and is licensed under the 3-clause
mi@0: # BSD license: https://github.com/detly/gammatone/blob/master/COPYING
mi@0: """
mi@0: This module contains functions for constructing sets of equivalent rectangular
mi@0: bandwidth gammatone filters.
mi@0: """
mi@0: from __future__ import division
mi@0: from collections import namedtuple
mi@0: 
mi@0: import numpy as np
mi@0: import scipy as sp
mi@0: from scipy import signal as sgn
mi@0: 
mi@0: DEFAULT_FILTER_NUM = 100
mi@0: DEFAULT_LOW_FREQ = 100
mi@0: DEFAULT_HIGH_FREQ = 44100/4
mi@0: 
mi@0: 
mi@0: def erb_point(low_freq, high_freq, fraction):
mi@0:     """
mi@0:     Calculates a single point on an ERB scale between ``low_freq`` and
mi@0:     ``high_freq``, determined by ``fraction``. When ``fraction`` is ``1``,
mi@0:     ``low_freq`` will be returned. When ``fraction`` is ``0``, ``high_freq``
mi@0:     will be returned.
mi@0:     
mi@0:     ``fraction`` can actually be outside the range ``[0, 1]``, which in general
mi@0:     isn't very meaningful, but might be useful when ``fraction`` is rounded a
mi@0:     little above or below ``[0, 1]`` (eg. for plot axis labels).
mi@0:     """
mi@0:     # Change the following three parameters if you wish to use a different ERB
mi@0:     # scale. Must change in MakeERBCoeffs too.
mi@0:     # TODO: Factor these parameters out
mi@0:     ear_q = 9.26449 # Glasberg and Moore Parameters
mi@0:     min_bw = 24.7
mi@0:     order = 1
mi@0: 
mi@0:     # All of the following expressions are derived in Apple TR #35, "An
mi@0:     # Efficient Implementation of the Patterson-Holdsworth Cochlear Filter
mi@0:     # Bank." See pages 33-34.
mi@0:     erb_point = (
mi@0:         -ear_q*min_bw
mi@0:         + np.exp(
mi@0:             fraction * (
mi@0:                 -np.log(high_freq + ear_q*min_bw)
mi@0:                 + np.log(low_freq + ear_q*min_bw)
mi@0:                 )
mi@0:         ) *
mi@0:         (high_freq + ear_q*min_bw)
mi@0:     )
mi@0:     
mi@0:     return erb_point
mi@0: 
mi@0: 
mi@0: def erb_space(
mi@0:     low_freq=DEFAULT_LOW_FREQ,
mi@0:     high_freq=DEFAULT_HIGH_FREQ,
mi@0:     num=DEFAULT_FILTER_NUM):
mi@0:     """
mi@0:     This function computes an array of ``num`` frequencies uniformly spaced
mi@0:     between ``high_freq`` and ``low_freq`` on an ERB scale.
mi@0:     
mi@0:     For a definition of ERB, see Moore, B. C. J., and Glasberg, B. R. (1983).
mi@0:     "Suggested formulae for calculating auditory-filter bandwidths and
mi@0:     excitation patterns," J. Acoust. Soc. Am. 74, 750-753.
mi@0:     """
mi@0:     return erb_point(
mi@0:         low_freq,
mi@0:         high_freq,
mi@0:         np.arange(1, num+1)/num
mi@0:         )
mi@0: 
mi@0: 
mi@0: def centre_freqs(high_freq, num_freqs, cutoff):
mi@0:     """
mi@0:     Calculates an array of centre frequencies (for :func:`make_erb_filters`)
mi@0:     from a sampling frequency, lower cutoff frequency and the desired number of
mi@0:     filters.
mi@0:     
mi@0:     :param fs: sampling rate
mi@0:     :param num_freqs: number of centre frequencies to calculate
mi@0:     :type num_freqs: int
mi@0:     :param cutoff: lower cutoff frequency
mi@0:     :return: same as :func:`erb_space`
mi@0:     """
mi@0:     # print 'center freq', erb_space(cutoff, fs/2, num_freqs)
mi@0:     return erb_space(cutoff, high_freq, num_freqs)
mi@0: 
mi@0: 
mi@0: def make_erb_filters(fs, centre_freqs, width=1.0):
mi@0:     """
mi@0:     This function computes the filter coefficients for a bank of 
mi@0:     Gammatone filters. These filters were defined by Patterson and Holdworth for
mi@0:     simulating the cochlea. 
mi@0:     
mi@0:     The result is returned as a :class:`ERBCoeffArray`. Each row of the
mi@0:     filter arrays contains the coefficients for four second order filters. The
mi@0:     transfer function for these four filters share the same denominator (poles)
mi@0:     but have different numerators (zeros). All of these coefficients are
mi@0:     assembled into one vector that the ERBFilterBank can take apart to implement
mi@0:     the filter.
mi@0:     
mi@0:     The filter bank contains "numChannels" channels that extend from
mi@0:     half the sampling rate (fs) to "lowFreq". Alternatively, if the numChannels
mi@0:     input argument is a vector, then the values of this vector are taken to be
mi@0:     the center frequency of each desired filter. (The lowFreq argument is
mi@0:     ignored in this case.)
mi@0:     
mi@0:     Note this implementation fixes a problem in the original code by
mi@0:     computing four separate second order filters. This avoids a big problem with
mi@0:     round off errors in cases of very small cfs (100Hz) and large sample rates
mi@0:     (44kHz). The problem is caused by roundoff error when a number of poles are
mi@0:     combined, all very close to the unit circle. Small errors in the eigth order
mi@0:     coefficient, are multiplied when the eigth root is taken to give the pole
mi@0:     location. These small errors lead to poles outside the unit circle and
mi@0:     instability. Thanks to Julius Smith for leading me to the proper
mi@0:     explanation.
mi@0:     
mi@0:     Execute the following code to evaluate the frequency response of a 10
mi@0:     channel filterbank::
mi@0:     
mi@0:         fcoefs = MakeERBFilters(16000,10,100);
mi@0:         y = ERBFilterBank([1 zeros(1,511)], fcoefs);
mi@0:         resp = 20*log10(abs(fft(y')));
mi@0:         freqScale = (0:511)/512*16000;
mi@0:         semilogx(freqScale(1:255),resp(1:255,:));
mi@0:         axis([100 16000 -60 0])
mi@0:         xlabel('Frequency (Hz)'); ylabel('Filter Response (dB)');
mi@0:     
mi@0:     | Rewritten by Malcolm Slaney@Interval.  June 11, 1998.
mi@0:     | (c) 1998 Interval Research Corporation
mi@0:     |
mi@0:     | (c) 2012 Jason Heeris (Python implementation)
mi@0:     """
mi@0:     T = 1/fs
mi@0:     # Change the followFreqing three parameters if you wish to use a different
mi@0:     # ERB scale. Must change in ERBSpace too.
mi@0:     # TODO: factor these out
mi@0:     ear_q = 9.26449 # Glasberg and Moore Parameters
mi@0:     min_bw = 24.7
mi@0:     order = 1
mi@0: 
mi@0:     erb = width*((centre_freqs/ear_q)**order + min_bw**order)**(1/order)
mi@0:     B = 1.019*2*np.pi*erb
mi@0: 
mi@0:     arg = 2*centre_freqs*np.pi*T
mi@0:     vec = np.exp(2j*arg)
mi@0: 
mi@0:     A0 = T
mi@0:     A2 = 0
mi@0:     B0 = 1
mi@0:     B1 = -2*np.cos(arg)/np.exp(B*T)
mi@0:     B2 = np.exp(-2*B*T)
mi@0:     
mi@0:     rt_pos = np.sqrt(3 + 2**1.5)
mi@0:     rt_neg = np.sqrt(3 - 2**1.5)
mi@0:     
mi@0:     common = -T * np.exp(-(B * T))
mi@0:     
mi@0:     # TODO: This could be simplified to a matrix calculation involving the
mi@0:     # constant first term and the alternating rt_pos/rt_neg and +/-1 second
mi@0:     # terms
mi@0:     k11 = np.cos(arg) + rt_pos * np.sin(arg)
mi@0:     k12 = np.cos(arg) - rt_pos * np.sin(arg)
mi@0:     k13 = np.cos(arg) + rt_neg * np.sin(arg)
mi@0:     k14 = np.cos(arg) - rt_neg * np.sin(arg)
mi@0: 
mi@0:     A11 = common * k11
mi@0:     A12 = common * k12
mi@0:     A13 = common * k13
mi@0:     A14 = common * k14
mi@0: 
mi@0:     gain_arg = np.exp(1j * arg - B * T)
mi@0: 
mi@0:     gain = np.abs(
mi@0:             (vec - gain_arg * k11)
mi@0:           * (vec - gain_arg * k12)
mi@0:           * (vec - gain_arg * k13)
mi@0:           * (vec - gain_arg * k14)
mi@0:           * (  T * np.exp(B*T)
mi@0:              / (-1 / np.exp(B*T) + 1 + vec * (1 - np.exp(B*T)))
mi@0:             )**4
mi@0:         )
mi@0: 
mi@0:     allfilts = np.ones_like(centre_freqs)
mi@0:     
mi@0:     fcoefs = np.column_stack([
mi@0:         A0*allfilts, A11, A12, A13, A14, A2*allfilts,
mi@0:         B0*allfilts, B1, B2,
mi@0:         gain
mi@0:     ])
mi@0:     
mi@0:     return fcoefs
mi@0: 
mi@0: 
mi@0: def erb_filterbank(wave, coefs):
mi@0:     """
mi@0:     :param wave: input data (one dimensional sequence)
mi@0:     :param coefs: gammatone filter coefficients
mi@0:     
mi@0:     Process an input waveform with a gammatone filter bank. This function takes
mi@0:     a single sound vector, and returns an array of filter outputs, one channel
mi@0:     per row.
mi@0:     
mi@0:     The fcoefs parameter, which completely specifies the Gammatone filterbank,
mi@0:     should be designed with the :func:`make_erb_filters` function.
mi@0:     
mi@0:     | Malcolm Slaney @ Interval, June 11, 1998.
mi@0:     | (c) 1998 Interval Research Corporation
mi@0:     | Thanks to Alain de Cheveigne' for his suggestions and improvements.
mi@0:     |
mi@0:     | (c) 2013 Jason Heeris (Python implementation)
mi@0:     """
mi@0:     output = np.zeros((coefs[:,9].shape[0], wave.shape[0]))
mi@0:     
mi@0:     gain = coefs[:, 9]
mi@0:     # A0, A11, A2
mi@0:     As1 = coefs[:, (0, 1, 5)]
mi@0:     # A0, A12, A2
mi@0:     As2 = coefs[:, (0, 2, 5)]
mi@0:     # A0, A13, A2
mi@0:     As3 = coefs[:, (0, 3, 5)]
mi@0:     # A0, A14, A2
mi@0:     As4 = coefs[:, (0, 4, 5)]
mi@0:     # B0, B1, B2
mi@0:     Bs = coefs[:, 6:9]
mi@0:     
mi@0:     # Loop over channels
mi@0:     for idx in range(0, coefs.shape[0]):
mi@0:         # These seem to be reversed (in the sense of A/B order), but that's what
mi@0:         # the original code did...
mi@0:         # Replacing these with polynomial multiplications reduces both accuracy
mi@0:         # and speed.
mi@0:         y1 = sgn.lfilter(As1[idx], Bs[idx], wave)
mi@0:         y2 = sgn.lfilter(As2[idx], Bs[idx], y1)
mi@0:         y3 = sgn.lfilter(As3[idx], Bs[idx], y2)
mi@0:         y4 = sgn.lfilter(As4[idx], Bs[idx], y3)
mi@0:         output[idx, :] = y4/gain[idx]
mi@0:         
mi@0:     return output