Resample

Functions for HRTF template interpolation and spherical grid management. The recommended entry point is resample.

Interpolation

bayesian_listener.resample.resample(cues, coordinates, template=None, method='SH', **kwargs)

Unified entry point for the four resampling methods.

Parameters:

• cues (numpy.ndarray or list of numpy.ndarray) – Single cue of shape (n_dirs, ...) or a list of cues whose first dimension matches.

• coordinates (pyfar.Coordinates) – Source coordinates of the input cues.

• template (pyfar.Coordinates or None, default=None) – Output directions. None selects a 64th-degree spherical t-design for every method.

• method ({'SH', 'SHMAX', 'barycentric', 'barumerli2023'}, default='SH') – Resampling method (case-insensitive).

• **kwargs – Forwarded to resample_two_step for 'SH', 'SHMAX', and 'barycentric': regularisation_coefficient, condition_threshold, norm.

Returns:

• result (numpy.ndarray or list of numpy.ndarray) – Resampled cues with first dimension equal to template.csize; same container type as cues.

• template_coords (pyfar.Coordinates) – Output directions.

Raises:

ValueError – If method is none of the four accepted values.

bayesian_listener.resample.resample_two_step(cues, coordinates, template, second_step, **kwargs)

Resample localisation cues using the two-step procedure of Ahrens et al.[1].

Stage 1: low-order SH extrapolation completes missing low-elevation directions by mirroring measured points across the horizontal plane (see complement_sampling). Stage 2: high-order interpolation of the complemented cues onto template.

Parameters:

• cues (numpy.ndarray or list of numpy.ndarray) – Cues as a single array or a list of arrays. For each array, shape[-2] must equal the number of source positions in coordinates. When a list is passed, the result is also a list in the same order.

• coordinates (pyfar.Coordinates) – Source coordinates corresponding to the cue rows.

• template (pyfar.Coordinates or None, default=None) – Output directions. None selects a 64th-degree spherical t-design (2,112 directions).

• second_step ({'SH', 'SHMAX', 'barycentric'}) –

  Stage-2 interpolator (case-insensitive).

  • 'SH': regularised SH interpolation at the maximum stable order for coordinates_complemented.

  • 'SHMAX': regularised SH interpolation at fixed order 44.

  • 'barycentric': VBAP weights on the convex hull.

• **kwargs –

  Forwarded options:

  • regularisation_coefficient (float, default 1e-2) — Tikhonov weight on the Bau damping matrix.

  • condition_threshold (float, default 12.25) — condition-number bound forwarded to find_max_order.

  • norm ({1, 2}, default 1) — gain normalisation for 'barycentric'; see vbap_interpolate.

Returns:

• cues (numpy.ndarray or list of numpy.ndarray) – Resampled cues. shape[-2] of each array equals template.csize. Same container type as the input.

• template_coords (pyfar.Coordinates) – Output directions.

Raises:

• TypeError – If coordinates or template is not pyfar.Coordinates.

• ValueError – If second_step is not one of the three accepted values.

bayesian_listener.resample.resample_barumerli2023(values, coords_in, template=None, flag_regularisation=True)

Resample with order-15 SH interpolation, as in Barumerli et al.[2].

Single-step SH interpolation at order \(N = 15\) with optional Tikhonov regularisation. Retained for backward compatibility with the original MATLAB implementation; assigns no probability mass to directions below the lowest measured elevation (see Known limitations).

Parameters:

• values (numpy.ndarray or list of numpy.ndarray) – Single cue of shape (n_dirs, ...) or a list of cues whose first dimension matches.

• coords_in (pyfar.Coordinates) – Source coordinates of the input cues.

• template (pyfar.Coordinates or None, default=None) – Output directions. None selects a 64th-degree spherical t-design.

• flag_regularisation (bool, default=True) – If True, apply a fixed Tikhonov regulariser (\(\\lambda = 4\)) ignoring the first three SH orders.

Returns:

• values_out (numpy.ndarray or list of numpy.ndarray) – Resampled cues; same container type as values.

• template_coords (pyfar.Coordinates) – Output directions.

bayesian_listener.resample.complement_sampling(coordinates)

Complement sampling grid.

The sampling grid is complemented by detecting the minimum elevation and adding points below. For example, if the minimum elevation is -30 degree, the grid is complemented by finding sampling points above 30 degree and mirroring them downward by flipping the sign of their elevation angles.

Note this method works for Gaussian-like sampling grids but might not work well in other cases.

Parameters:

coordinates (pyfar.Coordinates) – The incomplete sampling grid.

Returns:

• complemented (pyfar.Coordinates) – The complemented sampling grid.

• mask (numpy.ndarray) – Boolean array of shape (complemented.csize,); True at positions corresponding to mirrored (complemented) points and False for the original measurement directions.

bayesian_listener.resample.interpolate_HRTF(query_dirs, C, N)

Evaluate an SH expansion \(\mathbf{Y}(\mathbf{q})\mathbf{C}\) at query directions.

Parameters:

• query_dirs (numpy.ndarray) – Query directions of shape (n_query, 2) containing (azimuth, elevation) in radians.

• C (numpy.ndarray) – SH coefficient matrix of shape ((N + 1) ** 2, n_features).

• N (int) – SH order matching the columns of C.

Returns:

Interpolated cues of shape (n_query, n_features).

Return type:

numpy.ndarray

Spherical harmonics

bayesian_listener.resample.find_max_order(dirs, condition_threshold=12.25, N_max=35, regularised=True, regularisation_coefficient=0.01)

Return the largest SH order \(N \le N_{\max}\) with stable conditioning.

Iterates from \(N = 1\) upward and returns the highest order whose (optionally Tikhonov-regularised) Gram matrix \(\mathbf{Y}^\top\mathbf{Y}\) has condition number below condition_threshold.

Parameters:

• dirs (pyfar.Coordinates) – Sampling grid.

• condition_threshold (float, default=12.25) – Upper bound on \(\kappa(\mathbf{Y}^\top\mathbf{Y})\). The default follows Bau et al.[3], equivalent to \(\kappa(\mathbf{Y}) < 3.5\) (Ben-Hur et al.[4]).

• N_max (int, default=35) – Maximum SH order to test.

• regularised (bool, default=True) – If True, add the Bau damping matrix build_bau_damping scaled by regularisation_coefficient before computing the condition number.

• regularisation_coefficient (float, default=1e-2) – Tikhonov weight applied to the damping matrix.

Returns:

Largest admissible SH order. Returns 1 if no order satisfies the threshold.

Return type:

int

bayesian_listener.resample.solve_sh(Y, H)

Solve \(\mathbf{Y}\mathbf{C} \approx \mathbf{H}\) by Moore–Penrose pseudo-inverse.

Non-regularised least-squares. Use build_bau_damping and a direct solve when regularisation is required (as in resample_two_step).

Parameters:

• Y (numpy.ndarray) – SH basis matrix of shape (n_dirs, (N + 1) ** 2).

• H (numpy.ndarray) – Cue matrix of shape (n_dirs, n_features).

Returns:

Coefficient matrix of shape ((N + 1) ** 2, n_features).

Return type:

numpy.ndarray

bayesian_listener.resample.build_Y(dirs, N)

Build the matrix of real spherical-harmonic basis functions at each direction.

Parameters:

• dirs (numpy.ndarray) – Direction array of shape (n_dirs, 2) containing (azimuth, elevation) in radians.

• N (int) – SH order; the basis has \((N+1)^2\) columns.

Returns:

Real SH basis \(\mathbf{Y}\) of shape (n_dirs, (N + 1) ** 2).

Return type:

numpy.ndarray

bayesian_listener.resample.build_bau_damping(N)

Bau et al. damping matrix for Tikhonov-regularised SH inversion.

Builds the diagonal damping \(D_{ii} = 1 + n(n+1)\) indexed in SH order \(n\) (each \(n\) repeats \(2n+1\) times).

Parameters:

N (int) – Maximum SH order.

Returns:

Diagonal damping matrix of shape ((N + 1) ** 2, (N + 1) ** 2).

Return type:

numpy.ndarray

References

Bau et al.[3].

Visualisation

bayesian_listener.resample.plot_resampling_grid(coords_meas_cart, dirs_virt, missing_mask, z_min_meas)

Plot measured (black) and added (red) directions on a 3-D + 2-D figure.

Parameters:

• coords_meas_cart (numpy.ndarray) – Measured directions in Cartesian coordinates, shape (n_meas, 3).

• dirs_virt (numpy.ndarray) – Virtual grid directions in Cartesian coordinates, shape (n_virt, 3).

• missing_mask (numpy.ndarray) – Boolean mask of shape (n_virt,) selecting which virtual directions were added (i.e. fall in the unmeasured region).

• z_min_meas (float) – Minimum z-coordinate of the measured directions; the horizontal plane at this height is drawn as a translucent reference surface.

[1]

J. Ahrens, M. R. P. Thomas, and I. J. Tashev. HRTF magnitude modeling using a non-regularized least-squares fit of spherical harmonics coefficients on incomplete data. In APSIPA Annual Summit and Conference. Hollywood, USA, 2012.

[2]

R. Barumerli, P. Majdak, M. Geronazzo, D. Meijer, F. Avanzini, and R. Baumgartner. A Bayesian model for human directional localization of broadband static sound sources. Acta Acustica, 7:12, 2023. doi:10.1051/aacus/2023006.

[3] (1,2)

D. Bau, J. M. Arend, and C. Pörschmann. Estimation of the optimal spherical harmonics order for the interpolation of head-related transfer functions sampled on sparse irregular grids. Frontiers in Signal Processing, 2:884541, 2022. doi:10.3389/frsip.2022.884541.

[4]

Z. Ben-Hur, D. L. Alon, R. Mehra, and B. Rafaely. Loudness stability of binaural sound with spherical harmonic representation of sparse head-related transfer functions. EURASIP J. Audio Speech Music Process., 2019(1):5, 2019. doi:10.1186/s13636-019-0148-x.