Utils

Feature extraction, spherical grid utilities, caching helpers, and sampling. The recommended entry point for feature extraction is compute_features.

Feature extraction

bayesian_listener.utils.compute_features(hrir, coords, fs, spectral_range=[700.0, 18000.0], halfwave_rectifier=True, reference='frontal')

Compute ITD, ILD, and monaural spectral cues from binaural HRIRs.

Implements the feature extraction of Eq. 1 of Barumerli et al.[1]:

1. Normalise the whole HRIR set by a single level reference (reference).

2. Estimate ITD via itdestimator and apply the signed-log perceptual warp \(\mathrm{sgn}(t)\,(\log(a + b\,|t|) - \log a)/b\) with \(a = 32.5\,\mu\mathrm{s}\) and \(b = 0.095\).

3. Compute ILD as the dB ratio of broadband RMS energies.

4. Filter both ears with an ERB-spaced gammatone bank (gammatone) and convert per-band amplitude to dB.

Gammatone filtering uses a Numba-compiled fused kernel (_gammatone_rms_numba) that processes all bands concurrently across CPU cores without storing any intermediate filtered signal. Each band is zero-padded by its minimum impulse response length (minimum_ir_length) rather than a fixed 50 ms, reducing peak memory from ~3.4 GB to ~0.4 MB for a 793-direction HRTF at 48 kHz. The first call incurs a one-time JIT compilation cost of a few seconds.

Parameters:

• hrir (numpy.ndarray) – Head-related impulse responses, shape (n_dirs, 2, n_samples).

• coords (pyfar.Coordinates) – Source positions, one per HRIR row.

• fs (int) – Sampling rate in Hz.

• spectral_range (list of float, default=[700.0, 18000.0]) – [low_Hz, high_Hz] bracket for the gammatone filterbank.

• halfwave_rectifier (bool, default=True) – If True, apply half-wave rectification before computing the per-band mean: \(\sqrt{\mathrm{mean}(\max(x,0))}\). If False, compute the full-wave RMS instead: \(\sqrt{\mathrm{mean}(x^2)}\).

• reference ({'frontal', 'global', 'none'}, default='frontal') –

  Level reference \(k\) by which the whole HRIR set is divided before any cue is extracted.

  • 'frontal' — peak sample of the impulse response nearest to the frontal direction (the historical behaviour, kept as the default).

  • 'global' — peak sample over the whole set, both ears and all directions.

  • 'none' — no normalisation, \(k = 1\).

  ITD and ILD are unaffected by this choice: ITD is a timing readout and ILD is a level ratio, so \(k\) cancels. It shifts spectral_cues by a constant, identical for every direction and band.

  The choice therefore only matters when target and template come from different HRTF sets: their references \(k_A \neq k_B\) then leave a residual offset in every band that no later stage can tell apart from a genuine spectral difference. 'frontal' derives \(k\) from a single impulse response, whose peak is precisely what spectral smoothing moves most; 'global' references the whole set and so moves only when the set’s overall level moves, which is what makes cross-file comparison valid.

Raises:

ValueError – If reference is not one of 'frontal', 'global', 'none'.

Returns:

• itd (numpy.ndarray) – Warped ITD, shape (n_dirs, 1) (dimensionless after warping).

• ild (numpy.ndarray) – ILD in dB, shape (n_dirs, 1).

• spectral_cues (numpy.ndarray) – Monaural spectral amplitudes in dB, shape (n_dirs, n_freqs, 2).

• freqs (numpy.ndarray) – Filterbank centre frequencies in Hz, shape (n_freqs,).

bayesian_listener.utils.gammatone(fc, fs, n=4, betamul=None, scale='0dBforall', phase='causalphase')

Complex-valued, all-pole gammatone filter coefficients as in Lyon, 1997. The function has been taken from the AMT/gammatone.m.

Parameters:

• fc (array-like or scalar) – Center frequency(ies) in Hz. Must be within (0, fs/2].

• fs (float) – Sampling rate (Hz), positive.

• n (int, default=4) – Filter order (positive integer).

• betamul (float or None) –

  Multiplier for ERB bandwidth. If None, uses the MATLAB formula:

  betamul = (factorial(n-1)^2) / (pi*factorial(2n-2) * 2^(-(2n-2)))

• scale ({"0dBforall", "6dBperoctave"}) – Amplitude scaling mode.

• phase ({"causalphase", "peakphase", "exppeakphase"}) – Phase option. “exppeakphase” additionally aligns maxima using an impulse response simulation (requires SciPy).

Returns:

• b (np.ndarray, shape (n_channels,)) – Numerator scalar for each channel (complex).

• a (np.ndarray, shape (n_channels, n+1)) – Denominator polynomial coefficients (complex), leading 1 per row.

• delay (np.ndarray, shape (n_channels,)) – Peak time (seconds) of the envelope: 3 / (2*pi*beta).

• z (list) – Zeros per channel (always empty lists here, all-pole).

• p (list) – Poles per channel (length n per channel; all identical).

• k (list) – Gain per channel (always 1 here, to mirror the zpk form in MATLAB).

bayesian_listener.utils.minimum_ir_length(fc, fs, n=4, tolerance=5e-05)

Minimum impulse response length for Lyon 1997 gammatone band(s).

Adapts the pole-decay algorithm from pyfar / MATLAB impzlength to the Lyon 1997 all-pole gammatone filter. The dominant pole of each band has magnitude |z| = exp(-2π·β/fs); the impulse response of an n-th order filter with that repeated pole decays as |z|^k / k^(n-1), so the number of samples required to fall below tolerance is

n · log10(tolerance) / log10(|z|)

which is the pyfar formula for the “no-oscillation” IIR case with pole multiplicity n.

Parameters:

• fc (float or array-like) – Centre frequency (Hz) of each gammatone band. Must satisfy 0 < fc <= fs/2.

• fs (float) – Sampling rate (Hz).

• n (int, default=4) – Filter order (must match the n passed to gammatone).

• tolerance (float, default=5e-5) – Amplitude threshold below which the impulse response is considered negligible. The default 5e-5 matches pyfar’s default.

Returns:

lengths – Minimum number of samples for the impulse response of each band to decay to tolerance. Always at least n + 1.

Return type:

ndarray of int, shape (n_bands,)

bayesian_listener.utils.itdestimator(signals, fs=None)

Estimate the interaural time difference (ITD) by Hilbert-envelope MaxIACCe.

For each direction, low-pass filters the binaural HRIRs at 3 kHz with a 10th-order Butterworth filter, computes the analytic envelope of each ear via the Hilbert transform, and locates the peak of their cross- correlation (MaxIACCe mode of AMT 1.x itdestimator.m).

Parameters:

• signals (numpy.ndarray) – Binaural HRIRs of shape (n_dirs, 2, n_samples).

• fs (int) – Sampling rate in Hz. Required.

Returns:

ITD in seconds, shape (n_dirs, 1).

Return type:

numpy.ndarray

Raises:

ValueError – If fs is None.

Spherical grid

bayesian_listener.utils.load_n_design(degree)

Load a spherical t-design (n-design) grid of the given degree.

Points are loaded from a bundled .mat file covering degrees 1–124. These are Chebyshev-type quadrature rules on the unit sphere, equivalent to modern t-designs.

Parameters:

degree (int) – Degree of exactness, between 1 and 124.

Returns:

vecs – Cartesian coordinates of the grid points on the unit sphere.

Return type:

ndarray, shape (M, 3)

References

The grid data (n_designs_1_124.mat) was originally published by Manuel Gräf at https://homepage.univie.ac.at/manuel.graef/quadrature.php and redistributed by spaudiopy (MIT License, Copyright 2019 Chris Hold). The upstream data source does not carry an explicit license.

bayesian_listener.utils.vbap_interpolate(src, grid, norm=1)

Compute VBAP interpolation weights on the unit sphere.

For each source direction, finds the enclosing triangle on the convex hull of grid and returns the panning gains normalised according to norm.

Parameters:

• src (ndarray, shape (n_src, 3)) – Cartesian coordinates of target directions.

• grid (ndarray, shape (n_grid, 3)) – Cartesian coordinates of the source grid (unit sphere).

• norm ({1, 2}, default=1) – Gain normalisation: 1 — gains sum to 1 (anechoic, equivalent to barycentric interpolation); 2 — sum of squared gains equals 1 (energy-preserving / reverberant).

Returns:

weights – Sparse-like weight matrix. Each row has at most 3 non-zero entries normalised according to norm.

Return type:

ndarray, shape (n_src, n_grid)

References

Algorithm adapted from spaudiopy (MIT License, Copyright 2019 Chris Hold, https://github.com/chris-hold/spaudiopy), based on: Pulkki, V. (1997). Virtual Sound Source Positioning Using Vector Base Amplitude Panning. JAES, 45(6), 456–466.

Sampling

bayesian_listener.utils.scatter_von_mises(dirs, kappa, seed=None)

Perturb unit-direction vectors with von Mises–Fisher noise.

Implements Eq. 7 of Barumerli et al.[1]: each input direction is replaced by a sample from \(\mathrm{vMF}(\boldsymbol{\mu}_i, \kappa)\). The output preserves the input shape.

The scalar, per-direction reference implementation is randvmf below; it is kept as the readable specification of the algorithm (Rubinstein 1981, Fisher et al. 1987). This function is its vectorised equivalent and must produce the same distribution.

U and psi are drawn for all rows up front from a single generator, and only the deterministic sample-and-rotate step is chunked (_VMF_CHUNK rows at a time). That keeps peak memory roughly flat in n while making the result independent of the chunk size.

Warning

The seed is consumed once per call. Within a call every direction receives an independent draw; two calls with the same seed and the same input return identical output (intended reproducibility).

Therefore: use one seeded call per simulation. Calling this function repeatedly with a fixed seed — for example once per repetition inside a loop — gives every batch the same scatter and silently removes trial-to-trial variability, which is a subtle way to under-estimate response noise. If several calls are unavoidable, derive independent substreams with np.random.SeedSequence(seed).spawn(k), or pass seed=None.

Parameters:

• dirs (numpy.ndarray) – Direction vectors in Cartesian coordinates, shape (n, 3) or (3,) (each row should be unit-norm).

• kappa (float) – Von Mises–Fisher concentration; higher values yield tighter samples. Must be positive.

• seed (int or None, default=None) – Seed forwarded to numpy.random.default_rng. See the warning above for the semantics.

Returns:

Perturbed direction vectors, same shape as dirs.

Return type:

numpy.ndarray

Raises:

ValueError – If dirs does not have 3 components in its last dimension, or if kappa is not positive.

bayesian_listener.utils.randvmf(kappa, mu, seed=None)

Draw a single sample from a 3-D von Mises–Fisher distribution.

Uses the z-axis tangent-rotation algorithm of Rubinstein (1981) and Fisher et al. (1987): sample on a vMF aligned with the north pole, then rotate to align with mu via a Rodrigues rotation.

This is the readable, single-direction reference implementation. For arrays of directions use scatter_von_mises, which is the vectorised equivalent.

Warning

A fresh generator is constructed on every call, so calling this in a loop with a fixed seed makes every iteration return the same sample — a constant rotation rather than independent noise. Prefer scatter_von_mises for more than one direction; if you must loop here, pass seed=None or distinct seeds (np.random.SeedSequence(seed).spawn(k)).

Parameters:

• kappa (float) – Concentration parameter (positive).

• mu (numpy.ndarray) – Mean direction (unit vector), shape (3,).

• seed (int or None, default=None) – Seed forwarded to numpy.random.default_rng. See the warning above.

Returns:

Single direction vector, shape (3,).

Return type:

numpy.ndarray

Caching

bayesian_listener.utils.feature_cache_key(convention, spectral_range, halfwave_rectifier, reference)

Build the cache key identifying one set of feature-extraction parameters.

Every argument that changes the output of compute_features must appear here, otherwise two differently-computed targets collide in the cache and one is silently served for the other.

Parameters:

• convention (str) – Auditory representation name, e.g. 'Barumerli2023'.

• spectral_range (list of float) – [low_Hz, high_Hz] limits of the gammatone filterbank.

• halfwave_rectifier (bool) – Whether half-wave rectification is applied before the per-band mean.

• reference ({'frontal', 'global', 'none'}) – Level reference used to normalise the HRIRs.

Returns:

Stable key, e.g. 'Barumerli2023|hw=True|ref=frontal|sr=700.0-18000.0'.

Return type:

str

bayesian_listener.utils.cache_load_target(cache_dir, sofa_file, feature_key)

Load the cached target for a SOFA file and feature key, or None if not found.

bayesian_listener.utils.cache_load_template(cache_dir, sofa_file, feature_key, interpolation)

Load a cached template for a SOFA file, feature key, and interpolation method.

bayesian_listener.utils.cache_save_target(cache_dir, sofa_file, feature_key, target)

Save target under feature_key, creating the cache entry if needed.

bayesian_listener.utils.cache_save_template(cache_dir, sofa_file, feature_key, interpolation, template)

Add template to the HRTF pickle under the given feature and interpolation keys.

bayesian_listener.utils.clear_cache(sofa_file=None)

Remove cached preprocessed data.

Parameters:

sofa_file (str or None) – If provided, only removes cache for this specific file. If None, removes all cached data.