Gammatone filter bank: runnable experiments

This notebook accompanies the gammatone filters topic. It builds an auditory filter bank, listens to a single tone with it, and compares a cochleagram against an ordinary spectrogram on a synthetic vowel.

Everything imports from the clean module gammatone.py; the notebook only drives it.

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

from gammatone import erb, erb_space, gammatone_sos, gammatone_filterbank, cochleagram

fs = 16000

1. One channel = four biquads

A single 4th-order gammatone channel factors into four second-order sections. We design the 1 kHz channel and confirm its magnitude response peaks where it should.

fc = 1000.0
sos = gammatone_sos(fc, fs)
print('SOS shape:', sos.shape, '->', sos.shape[0], 'biquad sections')

w, H = signal.sosfreqz(sos, worN=8192, fs=fs)
print('peak at %.0f Hz, ERB = %.1f Hz' % (w[np.argmax(np.abs(H))], erb(fc)))

plt.figure(figsize=(9, 3.5))
plt.semilogx(w[1:], 20*np.log10(np.maximum(np.abs(H[1:]), 1e-6)))
plt.axvline(fc, color='k', ls=':', alpha=0.4)
plt.xlim(100, fs/2); plt.ylim(-60, 5)
plt.xlabel('Frequency [Hz]'); plt.ylabel('Magnitude [dB]')
plt.title('1 kHz gammatone channel'); plt.grid(True, alpha=0.3, which='both')
plt.tight_layout(); plt.show()

2. A tone excites the nearest place

Send a pure tone through a 32-channel bank and see which channel rings loudest. This is the tonotopic map in miniature: frequency maps to place.

cfs = erb_space(100, 7000, 32)
t = np.arange(fs) / fs
f0 = 1200.0
x = np.sin(2*np.pi*f0*t)

out = gammatone_filterbank(x, fs, cfs)
energy = np.mean(out[:, fs//2:]**2, axis=1)  # steady-state energy per channel
winner = cfs[np.argmax(energy)]
print('tone at %.0f Hz lights up the %.0f Hz channel' % (f0, winner))

plt.figure(figsize=(9, 3.5))
plt.semilogx(cfs, energy / energy.max(), 'o-')
plt.axvline(f0, color='C3', ls='--', label='tone frequency')
plt.xlabel('Channel centre frequency [Hz]'); plt.ylabel('Normalised energy')
plt.title('Which channel responds to a %.0f Hz tone' % f0)
plt.legend(); plt.grid(True, alpha=0.3, which='both')
plt.tight_layout(); plt.show()

3. Cochleagram vs spectrogram of a synthetic vowel

We synthesise a steady vowel: a 140 Hz glottal pulse train shaped by three formant resonances (a stylised /a/). The cochleagram resolves the low harmonics and formants the way the ear does, devoting its resolution to the perceptually important low end.

def synth_vowel(f0=140.0, formants=(700, 1220, 2600), dur=0.7, fs=16000):
    n = int(dur * fs)
    t = np.arange(n) / fs
    # Impulse train at the pitch period.
    period = int(round(fs / f0))
    glottal = np.zeros(n)
    glottal[::period] = 1.0
    # Cascade of resonant biquads = vocal-tract formants.
    y = glottal
    for fF in formants:
        bw = 80.0
        r = np.exp(-np.pi * bw / fs)
        theta = 2 * np.pi * fF / fs
        b = [1 - r]
        a = [1, -2 * r * np.cos(theta), r * r]
        y = signal.lfilter(b, a, y)
    return y / np.max(np.abs(y))

x = synth_vowel(fs=fs)
cfs = erb_space(120, 7000, 64)
coch, times = cochleagram(x, fs, cfs, frame_ms=25, hop_ms=10)

f_stft, t_stft, Z = signal.stft(x, fs=fs, nperseg=512, noverlap=384)
S_db = 20*np.log10(np.maximum(np.abs(Z)/np.abs(Z).max(), 1e-4))

fig, (a1, a2) = plt.subplots(1, 2, figsize=(11, 4))
a1.pcolormesh(t_stft, f_stft, S_db, shading='auto', vmin=-60, vmax=0, cmap='magma')
a1.set_ylim(0, 4000); a1.set_title('STFT spectrogram')
a1.set_xlabel('Time [s]'); a1.set_ylabel('Frequency [Hz]')
a2.pcolormesh(times, np.arange(len(cfs)), coch, shading='auto', vmin=-60, vmax=0, cmap='magma')
yt = np.linspace(0, len(cfs)-1, 6).astype(int)
a2.set_yticks(yt, ['%.0f' % cfs[i] for i in yt])
a2.set_title('Gammatone cochleagram')
a2.set_xlabel('Time [s]'); a2.set_ylabel('Channel centre [Hz]')
fig.tight_layout(); plt.show()

The horizontal bands in the cochleagram are the resolved harmonics of the 140 Hz pitch; the brighter clusters are the formants. Because the channels are ERB-spaced, the low harmonics that a linear spectrogram crams together are pulled apart here, which is exactly the representation a speech recogniser or pitch tracker wants to start from. See the embedded page for running this bank in real time.