Abstract

The human cochlea includes about 3000 inner hair cells which filter sounds at frequencies between 20 Hz and 20 kHz. This massively parallel frequency analysis is reflected in models of auditory processing, which are often based on banks of filters. However, existing implementations do not exploit this parallelism. Here we propose algorithms to simulate these models by vectorizing computation over frequency channels, which are implemented in "Brian Hears," a library for the spiking neural network simulator package "Brian." This approach allows us to use high-level programming languages such as Python, because with vectorized operations, the computational cost of interpretation represents a small fraction of the total cost. This makes it possible to define and simulate complex models in a simple way, while all previous implementations were model-specific. In addition, we show that these algorithms can be naturally parallelized using graphics processing units, yielding substantial speed improvements. We demonstrate these algorithms with several state-of-the-art cochlear models, and show that they compare favorably with existing, less flexible, implementations.

KEYWORDS:

Simple spiking model of the auditory periphery. The cochlea and inner hair cells are modeled using gammatone filtering followed by half-wave rectification and 1/3-power law compression. We model the auditory nerve fibers as leaky integrate-and-fire neurons defined by the stochastic differential equation where the current I is the output of the inner hair cell model and J(t) is physiological white noise (not the acoustic input noise). (A) Python implementation with the Brian Hears toolbox. Variable I in the neuron model is linked to the output of the filterbank. (B) Raster plot of the model output.

Pseudo code of the two algorithms for the direct form II transposed IIR filter, sequential channels (Algorithm 1) and vectorized channels and filter order (Algorithm 2). Bold faced variables are vectors of size the number of channels, subscripts j give elements of these vectors, and the * operation corresponds to element-wise multiplication. The input is X(s) and the output is Y(s), that is Xj(s) is the sample in channel j at time s and similarly for Y(s). Ai and Bi are the parameters of the filter (which can be different for each channel), and the Zi are a set of internal variables. Written without indices, A,B,Z refer to two-dimensional arrays, and the operations over these 2D arrays vectorize over both channels and the order of the filter. The notation Z⩲ refers to a shift with respect to the filter order index, so that The variable “order” is the order of the filter. For Algorithm 2, in Python, the code reads almost directly as above. In C⩲⩲ on the CPU, it reads as above but each line with a vector operation has a loop. On the GPU it reads as above, but is evaluated in parallel with one thread per channel and a loop over the filter order.

Auditory model with feedback. (A) Python program defining a time-varying filterbank with center frequency modulated by the output of a low-pass filter using the Brian Hears toolbox. (B) Corresponding box representation.

Performance of the dynamic compressive gammachirp model (DCGC) as a function of update interval. Solid line, left axis: computation time; Dashed lines, right axis: error-to-signal ratio for sounds at three different levels (30, 50, and 90 dB).

Computation time taken to simulate a gammatone filterbank as a function of the number of channels, with a 200-ms sound at 20 kHz [(B) is a magnified version of (A)]. Five different implementations are compared: Brian Hears in pure Python with vectorization over channels (BH Python), Brian Hears with C code generation, sequential channels (BH C1) or vectorized channels (BH C2), Brian Hears with GPU code generation, and Matlab with operations on time-indexed arrays.

Stereausis model of sound localization. (A) Schematic of the stereausis binaural network. Ipsilateral and contralateral auditory nerve fibers project to coincidence detectors (each pair of fibers projects to a specific neuron). (B) Python implementation with the Brian Hears toolbox. (C) Spike counts of all coincidence detectors (horizontal: characteristic frequency (CF) of the ipsilateral input, vertical: CF of the contralateral input) in response to a 500-ms tone at 200 Hz presented simultaneously at both ears. The inset shows spike counts along the short diagonal line, as a function of CF difference between the two inputs. (D) Same as (C) but with the tone shifted by 830 μs at the contralateral ear.