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)

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

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

1. Normalise HRIRs to the frontal direction.

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)}\).

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.

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.

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.

Caching

bayesian_listener.utils.cache_load_target(cache_dir, sofa_file)

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

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

Load a cached template for a SOFA file and interpolation method, or None.

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

Save target to the HRTF pickle, creating the cache entry if needed.

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

Add template to the HRTF pickle under the given interpolation key.

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.