Source code for gammatone

"""This module implements gammatone filters a filtering routine and a filterbank class.

The filterbank is based on [Hohmann2002]_.

.. plot::

    import pyfilterbank.gammatone

.. plot::

    import pyfilterbank.gammatone


    - Tests,
    - nice introduction with example,


.. [Hohmann2002]
   Hohmann, V., Frequency analysis and synthesis using a Gammatone filterbank,
   Acta Acustica, Vol 88 (2002), 433--442


import numpy as np
from numpy.fft import rfft, rfftfreq
from numpy import (arange, array, pi, cos, exp, log10, ones_like, sqrt, zeros)
# from scipy.misc import factorial
from scipy.signal import lfilter

# ERB means "Equivalent retangular band(-width)"
# Constants:
_ERB_L = 24.7
_ERB_Q = 9.265

[docs]def erb_count(centerfrequency): """Returns the equivalent rectangular band count up to centerfrequency. Parameters ---------- centerfrequency : scalar /Hz The center frequency in Hertz of the desired auditory filter. Returns ------- count : scalar Number of equivalent bandwidths below `centerfrequency`. """ return 21.4 * log10(4.37 * 0.001 * centerfrequency + 1)
[docs]def erb_aud(centerfrequency): """Retrurns equivalent rectangular band width of an auditory filter. Implements Equation 13 in [Hohmann2002]_. Parameters ---------- centerfrequency : scalar /Hz The center frequency in Hertz of the desired auditory filter. Returns ------- erb : scalar Equivalent rectangular bandwidth of an auditory filter at `centerfrequency`. """ return _ERB_L + centerfrequency / _ERB_Q
[docs]def hertz_to_erbscale(frequency): """Returns ERB-frequency from frequency in Hz. Implements Equation 16 in [Hohmann2002]_. Parameters ---------- frequency : scalar The Frequency in Hertz. Returns ------- erb : scalar The corresponding value on the ERB-Scale. """ return _ERB_Q * np.log(1 + frequency / (_ERB_L * _ERB_Q))
[docs]def erbscale_to_hertz(erb): """Returns frequency in Hertz from ERB value. Implements Equation 17 in [Hohmann2002]_. Parameters ---------- erb : scalar The corresponding value on the ERB-Scale. Returns ------- frequency : scalar The Frequency in Hertz. """ return (exp(erb/_ERB_Q) - 1) * _ERB_L * _ERB_Q
[docs]def frequencies_gammatone_bank(start_band, end_band, norm_freq, density): """Returns centerfrequencies and auditory Bandwidths for a range of gamatone filters. Parameters ---------- start_band : int Erb counts below norm_freq. end_band : int Erb counts over norm_freq. norm_freq : scalar The reference frequency where all filters are around density : scalar ERB density 1would be `erb_aud`. Returns ------- centerfrequency_array : ndarray """ norm_erb = hertz_to_erbscale(norm_freq) centerfrequencies = erbscale_to_hertz( arange(start_band, end_band, density) + norm_erb) return centerfrequencies
[docs]def design_filter( sample_rate=44100, order=4, centerfrequency=1000.0, band_width=None, band_width_factor=1.0, attenuation_half_bandwidth_db=-3): """Returns filter coefficient of a gammatone filter [Hohmann2002]_. Parameters ---------- sample_rate : int/scalar order : int centerfrequency : scalar band_width : scalar band_width_factor : scalar attenuation_half_bandwidth_db : scalar Returns ------- b, a : ndarray, ndarray """ if band_width: phi = pi * band_width / sample_rate # alpha = 10**(0.1 * attenuation_half_bandwidth_db / order) # p = (-2 + 2 * alpha * cos(phi)) / (1 - alpha) # lambda_ = -p/2 - sqrt(p*p/4 - 1) elif band_width_factor: erb_audiological = band_width_factor * erb_aud(centerfrequency) phi = pi * erb_audiological / sample_rate # a_gamma = ((factorial(pi * (2*order - 2)) * # 2**(-(2*order - 2))) / (factorial(order - 1)**2)) # b = erb_audiological / a_gamma # lambda_ = exp(-2 * pi * b / sample_rate) else: raise ValueError( 'You need to specify either `band_width` or `band_width_factor!`') alpha = 10**(0.1 * attenuation_half_bandwidth_db / order) p = (-2 + 2 * alpha * cos(phi)) / (1 - alpha) lambda_ = -p/2 - sqrt(p*p/4 - 1) beta = 2*pi * centerfrequency / sample_rate coef = lambda_ * exp(1j*beta) factor = 2 * (1 - abs(coef))**order b, a = array([factor]), array([1., -coef]) return b, a
[docs]def fosfilter(b, a, order, signal, states=None): """Return signal filtered with `b` and `a` (first order section) by filtering the signal `order` times. This Function was created for filtering signals by first order section cascaded complex gammatone filters. Parameters ---------- b, a : ndarray, ndarray Filter coefficients of a first order section filter. Can be complex valued. order : int Order of the filter to be applied. This will be the count of refiltering the signal order times with the given coefficients. signal : ndarray Input signal to be filtered. states : ndarray, default None Array with the filter states of length `order`. Initial you can set it to None. Returns ------- signal : ndarray Output signal, that is filtered and complex valued (analytical signal). states : ndarray Array with the filter states of length `order`. You need to loop it back into this function when block processing. """ if not states: states = zeros(order, dtype=np.complex128) for i in range(order): state = [states[i]] signal, state = lfilter(b, a, signal, zi=state) states[i] = state[0] b = ones_like(b) return signal, states
[docs]def freqz_fos(b, a, order, nfft, plotfun=None): """Returns Frequency Response and impulse response of first order section coeffs. Parameters ---------- b : arraylike a : arraylike order : int nfft : int plotfun : function(frequencies, response) Returns ------- freqresponse : arraylike frequencies : arraylike response : arraylike """ impulse = _create_impulse(nfft) response, states = fosfilter(b, a, order, impulse) freqresponse = rfft(np.real(response)) frequencies = rfftfreq(nfft) if plotfun: plotfun(frequencies, freqresponse) return freqresponse, frequencies, response
[docs]def design_filtbank_coeffs( samplerate, order, centerfrequencies, bandwidths=None, bandwidth_factor=None, attenuation_half_bandwidth_db=-3): """Designs complex coefficients for a gammatone filter bank. Parameters ---------- samplerate : int order : int centerfrequencies : arraylike bandwidths : arraylike (optional) bandwidth_factor : scalar (1.0 is auditory erb) attenuation_half_bandwidth_db : scalar Returns ------- generator of b, a coefficients """ for i, cf in enumerate(centerfrequencies): if bandwidths: bw = bandwidths[i] bwf = None else: bw = None bwf = bandwidth_factor yield design_filter( samplerate, order, cf, band_width=bw, band_width_factor=bwf, attenuation_half_bandwidth_db=attenuation_half_bandwidth_db)
[docs]class GammatoneFilterbank: """Returns a GammatoneFilterbank instance for filtering signals. Parameters ---------- samplerate : scalar Default: 44100. order : scalar Default: 4. startband : scalar Default: -12, endband : scalar Default: 12, normfreq : scalar Default: 1000.0, density : scalar Default: 1.0, bandwidth_factor : scalar Default: 1.0, desired_delay_sec : scalar Default: 0.02 Attributes ---------- """ def __init__( self, samplerate=44100, order=4, startband=-12, endband=12, normfreq=1000.0, density=1.0, bandwidth_factor=1.0, desired_delay_sec=0.02): self.samplerate = samplerate self.order = order self.centerfrequencies = frequencies_gammatone_bank( startband, endband, normfreq, density) self._coeffs = tuple(design_filtbank_coeffs( samplerate, order, self.centerfrequencies, bandwidth_factor=bandwidth_factor)) self.init_delay(desired_delay_sec) self.init_gains()
[docs] def init_delay(self, desired_delay_sec): """Initializes delay estimation and variables for delay bands to achieve a optimized impulse response. """ self.desired_delay_sec = desired_delay_sec self.desired_delay_samples = int(self.samplerate*desired_delay_sec) self.max_indices, self.slopes = self.estimate_max_indices_and_slopes( delay_samples=self.desired_delay_samples) self.delay_samples = self.desired_delay_samples - self.max_indices self.delay_memory = np.zeros((len(self.centerfrequencies), np.max(self.delay_samples)))
[docs] def init_gains(self): """Initializes gains for weighting bands leading to a equalized freqresponse when summing (synthesizing) bands. """ self.gains = np.ones(len(self.centerfrequencies)) # not correct until now: # x, s = list(zip(*self.analyze(_create_impulse(self.samplerate/10)))) # rss = [np.sqrt(np.sum(np.real(b)**2)) for b in x] # self.gains = 1/np.array(rss)
[docs] def analyze(self, signal, states=None): """Returns a generator yielding filtered signal bands and states. """ for i, (b, a) in enumerate(self._coeffs): st = None if not states else states[i] yield fosfilter(b, a, self.order, signal, states=st)
[docs] def reanalyze(self, bands, states=None): """Refilters already filtered signal bands if bandwidening effect of a used algorithm needs to be reduced. """ for i, ((b, a), band) in enumerate(zip(self._coeffs, bands)): st = None if not states else states[i] yield fosfilter(b, a, self.order, band, states=st)
[docs] def synthesize(self, bands): """Returns summed dleayed and weighted bands. """ return np.array(list(self.delay( [b*g for b, g in zip(bands, self.gains)]))).sum(axis=0)
[docs] def delay(self, bands): """Returns delayed bands for an optimized impulseresponse at synthesis. """ self.phase_factors = np.abs(self.slopes)*1j / self.slopes for i, band in enumerate(bands): phase_factor = self.phase_factors[i] delay_samples = self.delay_samples[i] if delay_samples == 0: yield np.real(band) * phase_factor else: yield np.concatenate( (self.delay_memory[i, :delay_samples], np.real(band[:-delay_samples])), axis=0) self.delay_memory[i, :delay_samples] = np.real( band[-delay_samples:])
[docs] def estimate_max_indices_and_slopes(self, delay_samples=None): """Returns Estimate of maximum index and slopes at max. Used for estimation of phase-factors for delaying bands. """ if not delay_samples: delay_samples = self.samplerate/10 sig = _create_impulse(delay_samples) bands = list(zip(*self.analyze(sig)))[0] ibandmax = [np.argmax(np.abs(b[:delay_samples])) for b in bands] slopes = [b[i+1]-b[i-1] for (b, i) in zip(bands, ibandmax)] return np.array(ibandmax), np.array(slopes)
[docs] def freqz(self, nfft=4096, plotfun=None): """Returns frequency and impulse responses. Accepts `plotfun` function for cusomized plotting. """ def gen_freqz(): for b, a in self._coeffs: yield freqz_fos(b, a, self.order, nfft, plotfun) return list(gen_freqz())
def _create_impulse(num_samples): sig = zeros(num_samples) + 0j sig[0] = 1.0 return sig def example_filterbank(): from pylab import plt import numpy as np x = _create_impulse(2000) gfb = GammatoneFilterbank(density=1) analyse = gfb.analyze(x) imax, slopes = gfb.estimate_max_indices_and_slopes() fig, axs = plt.subplots(len(gfb.centerfrequencies), 1) for (band, state), imx, ax in zip(analyse, imax, axs): ax.plot(np.real(band)) ax.plot(np.imag(band)) ax.plot(np.abs(band)) ax.plot(imx, 0, 'o') ax.set_yticklabels([]) [ax.set_xticklabels([]) for ax in axs[:-1]] axs[0].set_title('Impulse responses of gammatone bands') fig, ax = plt.subplots() def plotfun(x, y): ax.semilogx(x, 20*np.log10(np.abs(y)**2)) plt.hold(True) gfb.freqz(nfft=2*4096, plotfun=plotfun) plt.grid(True) plt.title('Absolute spectra of gammatone bands.') plt.xlabel('Normalized Frequency (log)') plt.ylabel('Attenuation /dB(FS)') plt.axis('Tight') plt.ylim([-90, 1]) plt.xlim([0.001, 0.1]) # return gfb def example_filter(): from pylab import plt, np sample_rate = 44100 order = 4 b, a = design_filter( sample_rate=sample_rate, order=order, centerfrequency=1000.0, attenuation_half_bandwidth_db=-3, band_width_factor=1.0) x = _create_impulse(1000) y, states = fosfilter(b, a, order, x) y = y[:800] plt.plot(np.real(y), label='Re(z)') plt.plot(np.imag(y), label='Im(z)') plt.plot(np.abs(y), label='|z|') plt.legend() # return y, b, a if __name__ == '__main__': from pylab import plt gfb = example_filterbank() y = example_filter()