The variation of interaural level difference (ILD) with direction and frequency is particularly complex and convoluted. The purpose of this work was to determine a set of parametric equations that can be used to calculate ILDs continuously at any value of frequency and azimuth in the horizontal plane. They were derived by fitting equations to ILDs derived from the azimuthal-dependence data tabulated by Shaw and Vaillancourt [(1985). J. Acoust. Soc Am. 78, 1120–1123] and assuming left-right symmetry. The equations are shown to fit those data to an overall RMS error less than 0.5 dB.

To complement some current projects on spatial hearing, it was necessary to have a set of closed-form equations from which interaural time and level differences (ITDs and ILDs, respectively) could be calculated continuously for any value of azimuth or frequency in the horizontal plane. For ITDs, it is typical to use the expression (rθ + r sin θ)/c, where r is the radius of the head, θ is the azimuth of the sound source, in radians, and c is the speed of sound (Woodworth, 1938; Blauert, 1997). This expression is derived from a simple geometrical model in which it is assumed that the head is spherical in shape (Duda and Martens, 1998), the two ears are diametrically opposite each other, and the source of sound is sufficiently far away for the wavefronts to be planes. It has proved remarkably popular, even though the equation is not an exact fit to any particular individual's ITDs nor does it take into account the frequency dependence of ITDs. However, the acoustics of ILDs are far more complicated, and to the best of our knowledge there is no published set of expressions for ILDs that allow accurate calculations. The best simple approximation that we know of is that published by van Opstal (2016), which is $ILD=0.18fsin(θ)$, where f is the frequency in Hz. This equation captures two overall phenomena pertaining to the shape of the function of ILD across azimuth and frequency: the dependence on azimuth essentially follows the first half (0-π) of a sinusoid, and the overall magnitude of the ILDs increases with frequency.

An ILD is the across-ear difference in spectral magnitude of the head-related transfer function (HRTF). As HRTFs are generally measured experimentally only at discrete angles [see Li and Peissig (2020) for a recent review of measurement methods], there has been substantial work on computational methods to generate continuous values at any angle. These methods include numerical interpolation [e.g., Brinkmann et al. (2015), Gamper (2013), and Grijalva et al. (2017)], acoustic simulations [e.g., Katz (2001), Kreuzer et al. (2009), and Mokhtari et al. (2011)], neural networks (Qi and Tao, 2017), or weighted sums of various forms with, in general, considerable analytic complexity, including diffraction [e.g., Blauert (1997)], Fourier-Bessel series [e.g., Zhang et al. (2009)], principal components [e.g., Kistler and Wightman (1992), Middlebrooks and Green (1992), and Mokhtari et al. (2019)], Slepian functions [e.g., Bates et al. (2015)], spherical harmonics [e.g., Evans et al. (1998)], spherical wavelets (Hu et al., 2019), and from the analysis of a spherical head [Duda and Martens (1998), who also include pseudocode]. Nevertheless, none of these cited papers has published any tables of coefficients or weights that are suitable for direct numerical calculation by their readers.

To resolve this, we set out to develop continuous parametric equations that could be fitted to published ILD data. We derived these from tabulated data reported by Shaw and Vaillancourt (1985) [henceforth S-V], which quantify “self-consistent” families of curves that best fitted the azimuthal-dependence (see below) data from twelve previous studies (Shaw, 1974). These curves were described as “syntheses,” were visually fitted and are internally consistent with smooth progressions, and represent the earlier data as a whole, though being neither true averages across listeners nor the functions of any particular individual or manikin. Shaw (1974) noted that individuals would generally be within 1 dB below 500 Hz but 5 dB or more at about 5 kHz. Though they are not any individual's ILDs, for our purposes this was much outweighed by the convenience of using smoothed, tabulated, and published data. Overall, we set ourselves the pragmatic target of fitting these within 1 dB or less within the range of 200 to 10 000 Hz. We avoided higher frequencies as they can be particularly idiosyncratic, and we could not find a simple set of expressions that were suitable.

We derived ILDs from S-V's data on the “azimuthal dependence” of at-ear level at 33 frequencies from 200 to 10 000 Hz (spaced at 50, 100, 200, and 500 Hz) and 24 azimuths from 0 to 3450 (spaced at 150; see Fig. 1(A)]; we excluded data at the irregularly spaced frequencies of 320, 630, 1250, 2900, 6300 Hz and, as noted, above 10 500 Hz. The azimuthal dependence is the difference in level at the left eardrum between a sound presented at an azimuth θ and a sound presented at 0° (i.e., directly in front of the listener). We then calculated the ILD θ,f as the difference between the azimuthal dependence at θ and at −θ (or, in terms of S-V's tabulated data, at 360°-θ); note that we assumed left-right symmetry in ILDs.

Fig. 1.

(A) ILDs derived from the data reported by Shaw and Vaillancourt (1985); (B) computed ILDs from the equations reported here; (C) the error between the data and the computations, plotted using the same y axis limits in dB as panels A and B; (D) the same error, plotted using much finer y axis limits.

Fig. 1.

(A) ILDs derived from the data reported by Shaw and Vaillancourt (1985); (B) computed ILDs from the equations reported here; (C) the error between the data and the computations, plotted using the same y axis limits in dB as panels A and B; (D) the same error, plotted using much finer y axis limits.

Close modal

To derive the fitting equations, we started with a “principal” sinusoid dependent on azimuth θ, capturing the overall dependence on azimuth that essentially follows the first half (0-π) of a sinusoid. Next, for some frequencies, the ILD function is front-back asymmetric, in that the ILD at an azimuth of θ is considerably different to the ILD at 180-θ. We represented this by varying asymmetrically the magnitude of the principal sinusoid by a factor dependent on the argument 2θ. Next, for many frequencies, there is a large perturbation, or “dip,” at azimuths around 90°. These were represented by adding a normal distribution dependent on azimuth θ, whose parameters of magnitude, mean and standard deviation were themselves functions of frequency. A second perturbation around 5000 Hz required another normal distribution. In sum, our general equation for ILD as a function of frequency and azimuth was

$ILDθ,f=PfsinθAfsin2θ+1+D90,fNθ,μ90,f,σ90,f+D5k,fNθ,μ5k,σ5k,$
(1)

where Pf is the magnitude of the principal sinusoid at frequency f; Af is the magnitude of the front-back asymmetry; D90,f, μ90,f, and σ90,f are the parameters of the normal distribution accounting for the “dip”; D5k,f, μ5k, and σ5k are the parameters for the normal distribution accounting for the 5-kHz perturbation; and N is the standard normal distribution. We found the best-fitting values of each parameter separately for each value of frequency, using the generalised reduced gradient (GRG) algorithm within Excel's “Solver” tool. Figure 2 shows these values (circles). Their dependencies on frequency are mostly smooth but in some cases are of considerable extent and complexity.

Fig. 2.

(symbols) the values of the six key parameters found from within-frequency fits to the data; (lines) the computed values from the equations reported here.

Fig. 2.

(symbols) the values of the six key parameters found from within-frequency fits to the data; (lines) the computed values from the equations reported here.

Close modal

We next needed to determine further equations that would give continuous values of these parameters as a function of frequency. Trial-and-error work demonstrated that we could fit the functions of these parameters across frequency by summing a linear function of f * = log10(f) with four perturbing normal distributions [see Eqs. (2)–(8) below]. The best-fitting values of the parameters of these sums were found by again running the GRG nonlinear algorithm. We set limits of –10 ≤  a  ≤ 10, σ ≤ 0.03, and the various values of μ were at least 0.1 log(kHz) apart. The lines in Fig. 2 show these fits.

The parameters that we found are listed in Eqs. (2)–(8) below. The various perturbations to the main functions (p1600, etc.) are indexed by the centre-frequency of their peak, in Hz. These seven equations, together with the model set out in Eq. (1), form our continuous parametric equations for ILDs in the horizontal plane. The units are decibels for Pf, D90,f, D5k,f, and the magnitudes of their normal-distribution perturbations, and the units are log(kHz) for the means and standard-deviations of those normal-distributions, while Af is a dimensionless number. Computational programs are available as supplementary data.1

$Pf=13.062f*+11.696+p1600+p2010+p4880+p6140p1600=+2.865 Nf*,+0.203, 0.351p2010=−1.772 Nf*,+0.303, 0.114p4880=−4.022 Nf*,+0.688, 0.222p6140=+0.619 Nf*,+0.788, 0.049$
(2)
$Af=−0.578f*+0.504+a220+a1080+a1970+a4510a220=−0.959 Nf*,−0.649, 0.374a1080=−0.144 Nf*,+0.035, 0.168a1970=−0.044 Nf*,+0.295, 0.068a4510=+0.039 Nf*,+0.654, 0.064$
(3)
$D90,f=+0.647f*+1.118+d1130+d1850+d3180+d21400d1130=−7.516 Nf*,+0.052, 0.352d1850=+1.242 Nf*,+0.266, 0.107d3180=−0.128 Nf*,+0.503, 0.030d21400=−9.038 Nf*,+1.330, 0.220$
(4)
$μ90,f=−0.224f*+1.605+m100+m2980+m6720+m21890m100=+2.868 Nf*,−1.010, 0.135m6720=−0.041 Nf*,+0.827, 0.030m2980=−0.028 Nf*,+0.475, 0.054m21890=−6.000 Nf*,+1.340, 0.030$
(5)
$σ90,f=−0.063f*+0.326+s500+s1760+s4450+s5610s500=+0.040 Nf*,−0.302, 0.149s4450=+0.160 Nf*,+0.649, 0.584s5610=−0.115 Nf*,+0.749, 0.126s1760=−0.055 Nf*,+0.245, 0.092$
(6)
$D5k,f=−0.266 Nf*,+0.707, 0.076$
(7)
$μ5k=0.750; σ5k=0.200$
(8)

Figure 1(B) shows the resulting ILDs calculated from these equations. Note that we used a much finer “mesh” of frequency by azimuth in order to check there were no issues or instabilities at intermediate values. Figure 1(C) shows the difference between our calculations and the ILDs derived from the S-V data, using the same y-scale as panels (A) and (B). Figure 1(D) shows the same differences using a finer y-scale of ±1 dB. It can be seen that the difference is mostly less than half a decibel. The overall RMS error, across the entire grid of 13 azimuths by 33 frequencies, was 0.40 dB, with 97% of the values being less than 1 dB. The RMS error was generally slightly smaller at low frequencies than high ones (RMS = 0.3 dB for 0.2–2 kHz but 0.5 dB for 8–10 kHz).

We report here a set of equations that enable ILDs to be calculated for any angle of azimuth up to 10 kHz. The equations are derived from tabulated data reported by Shaw and Vaillancourt (1985) and assume left-right symmetry.

To recap, ILDs are represented first as a sinusoid dependent on azimuth θ, whose magnitude is front-back asymmetric and varies with frequency f. This front-back asymmetry depends on 2θ and also varies with frequency. Two major corrections (or “perturbations”) are then added to the ILDs to represent effects at azimuths around 90° and 5 kHz, which are both modelled as normal distributions dependent on azimuth and frequency. The magnitudes of all of these effects are modelled as linear functions of log(f), corrected by four normal distributions. We appreciate that the equations are somewhat daunting at first sight, and are free of any acoustical theory. Further, though the assumption of left-right symmetry is often made in experimental work—for instance, the popular Gardner and Martin (1995) dataset is symmetric—it is known that HRTFs are not perfectly symmetric [e.g., Zhong et al. (2013)]. An analysis of the CIPIC database of individual HRTFs (Algazi et al., 2001) shows that, in terms of ILD, the asymmetry is on the order of 2 dB up to about 4 kHz and then gradually increases to about 6 dB (also, the standard deviation across listeners of the mean ILDs across azimuths follows essentially the same pattern). On the other hand, these equations conveniently give a value for ILD at any desired azimuth and frequency due to the continuous, rather than discrete, numerical modelling of the data.

We tested various reduced equations to determine their effectiveness, but we have not found any simpler way to represent the complexity in the frequential and azimuthal dependence of the ILD to the accuracy we desired. Excluding the 5-kHz perturbation barely affected the overall RMS error (0.5 dB)—though it doubled the RMS error at 5-kHz to 1 dB—but additionally excluding the 90° dip increased the overall RMS error to 1.6 dB. Further excluding the front-back asymmetry, so the ILD was modelled by the primary sinusoid only, increased the error to 1.9 dB. Excluding the various normal corrections in Eq. (2)–(8), i.e., modelling each parameter as a linear function only, increased the RMS error to 3 dB, and with the mean ILD being about 2 dB greater than the data. In this regard, the equation $ILD=0.18fsin(θ)$ of Van Opstal (2016) does slightly better than this much-simplified reduction of our equations, as it gives a RMS error of 2.7 dB on our data, and with a mean ILD about 2 dB less than the data. Thus, that equation can be recommended if a front-back symmetric ILD is sufficient. However, if the front-back asymmetry matters (or if the major perturbations that we identified are required), then the equations reported here should be suitable for the purpose of calculating ILDs for any value of azimuth or frequency in the horizontal plane for use in computational research.

This work was supported by the Medical Research Council (Grant No. MR/S002898/1). We thank the Editor and the two expert reviewers for their comments, which improved the manuscript greatly.

1

See supplementary material at https://www.scitation.org/doi/suppl/10.1121/10.0004261 or Graetzer (2021) for matlab and Python functions implementing the equations.

1.
Algazi
,
V. R.
,
Duda
,
R. O.
,
Thompson
,
D. M.
, and
Avendano
,
C.
(
2001
). “
The CIPIC HRTF database
,” in
Proceedings of the 2001 IEEE ASSP Workshop on Applications of Signal Processing to Audio and Acoustics
, New Paltz, NY.
2.
Bates
,
A. P.
,
Khalid
,
Z.
, and
Kennedy
,
R. A.
(
2015
). “
On the use of Slepian functions for the reconstruction of the head-related transfer function on the sphere
,” in
9th International Conference on Signal Processing and Communication Systems (ICSPCS)
,
QLD
,
Cairns
, pp.
1
7
.
3.
Blauert
,
J.
(
1997
).
Spatial Hearing: The Psychophysics of Human Sound Localization
(
MIT Press
,
Cambridge, MA
).
4.
Brinkmann
,
F.
,
Roden
,
R.
,
Lindau
,
A.
, and
Weinzierl
,
S.
(
2015
). “
Audibility and Interpolation of Head-Above-Torso Orientation in Binaural Technology
,”
IEEE J. Sel. Top. Sign. Proc.
9
(
5
),
931
942
.
5.
Duda
,
R. O.
, and
Martens
,
W. L.
(
1998
). “
Range dependence of the response of a spherical head model
,”
J. Acoust. Soc. Am.
104
,
3048
3058
.
6.
Evans
,
M. J.
,
Angus
,
J. A. S.
, and
Tew
,
A. I.
(
1998
). “
Analyzing head-related transfer function measurements using surface spherical harmonics
,”
J. Acoust. Soc. Am.
104
,
2400
2411
.
7.
Gamper
,
H.
(
2013
). “
Head-related transfer function interpolation in azimuth, elevation, and distance
,”
J. Acoust. Soc. Am.
134
,
EL547
EL553
.
8.
Gardner
,
W. G.
, and
Martin
,
K. D.
(
1995
). “
HRTF measurements of a KEMAR
,”
J. Acoust. Soc. Am.
97
,
3907
3908
.
9.
Graetzer
,
S.
(
2021
). “
Akeroyd_parametricilds
,” https://github.com/sgraetzer/Akeroyd_parametricilds (Last viewed 4/19/2021).
10.
Grijalva
,
F.
,
Martini
,
L. C.
,
Florencio
,
D.
, and
Goldenstein
,
S.
(
2017
). “
Interpolation of head-related transfer functions using manifold learning
,”
IEEE Sign. Proc. Lett.
24
,
221
225
.
11.
Hu
,
S.
,
Trevino
,
J.
,
,
C.
,
Sakamoto
,
S.
, and
Suzuki
,
Y.
(
2019
). “
Modeling head-related transfer functions with spherical wavelets
,”
Appl. Acoust.
146
,
81
88
.
12.
Katz
,
B. F. G.
(
2001
). “
Boundary element method calculation of individual head-related transfer function. II. Impedance effects and comparisons to real measurements
,”
J. Acoust. Soc. Am.
110
,
2449
2455
.
13.
Kistler
,
D. J.
, and
Wightman
,
F. L.
(
1992
). “
A model of head-related transfer functions based on principal components analysis and minimum-phase reconstruction
,”
J. Acoust. Soc. Am.
91
,
1637
1647
.
14.
Kreuzer
,
W.
,
Majdak
,
P.
, and
Chen
,
Z.
(
2009
). “
Fast multipole boundary element method to calculate head-related transfer functions for a wide frequency range
,”
J. Acoust. Soc. Am.
126
,
1280
1290
.
15.
Li
,
S.
, and
Pessig
,
J.
(
2020
). “
Measurement of head-related transfer functions: A review
,”
Appl. Sci.
10
,
5014
.
16.
Middlebrooks
,
J. C.
, and
Green
,
D. M.
(
1992
). “
Observations on a principal components analysis of head-related transfer functions
,”
J. Acoust. Soc. Am.
92
,
597
599
.
17.
Mokhtari
,
P.
,
Kato
,
H.
,
Takemoto
,
H.
, et al (
2019
). “
Further observations on a principal components analysis of head-related transfer functions
,”
Sci. Rep.
9
,
7477
.
18.
Mokhtari
,
P.
,
Takemoto
,
H.
,
Nishimura
,
R.
, and
Kato
,
H.
(
2011
). “
Computer simulation of KEMAR's head-related transfer functions: Verification with measurements and acoustic effects of modifying head shape and pinna concavity
,” in
Principles and Applications of Spatial Hearing
, edited by
Y.
Suzuki
,
D.
Brungart
,
Y.
Iwaya
,
K.
Iida
,
D.
Cabrera
, and
H.
Kato
(
World Scientific
,
Singapore
).
19.
Qi
,
X.
, and
Tao
,
J.
(
2017
). “
A domain knowledge-assisted nonlinear model for head-related transfer functions based on bottleneck deep neural network
,” in
Proceedings of Interspeech 2017
,
3058
3062
.
20.
Shaw
,
E. A. G.
(
1974
). “
Transformation of sound pressure level from the free field to the eardrum in the horizontal plane
,”
J. Acoust. Soc. Am.
56
,
1848
1861
.
21.
Shaw
,
E. A. G.
, and
Vaillancourt
,
M. M.
(
1985
). “
Transformation of sound‐pressure level from the free field to the eardrum presented in numerical form
,”
J. Acoust. Soc Am.
78
,
1120
1123
.
22.
Woodworth
,
R. S.
(
1938
).
Experimental Psychology
(
Holt
,
New York
).
23.
Van Opstal
,
J.
(
2016
).
The Auditory System and Human Sound-Localisation Behavior
(
,
London
).
24.
Zhang
,
W.
,
Kennedy
,
R. A.
, and
Abhayapala
,
T. D.
(
2009
). “
Efficient continuous HRTF model using data independent basis functions: Experimentally guided approach
,”
IEEE Trans. Audio Speech Lang. Proc.
17
(
4
),
819
829
.
25.
Zhong
,
X.
,
Zhang
,
F.
, and
Xie
,
N.
(
2013
). “
On the spatial symmetry of head-related transfer functions
,”
Appl. Acoust.
74
(
6
),
856
864
.