Spatial upsampling of individual sparse head-related transfer function setsby directional equalization
Christoph PÖRSCHMANN(1)† , Johannes M. AREND(1)(2) , Fabian BRINKMANN(2)(1) Institute of Communications Engineering, TH Köln, D-50679 Cologne, Germany,(2) Audio Communication Group, TU Berlin, D-10587 Berlin, Germany† Corresponding author, E-mail: christoph.poerschmann@th-koeln.de

Abstract
Determining full-spherical individual sets of head-related transfer functions (HRTFs) based on sparse measurements is a prerequisite for various applications in virtual acoustics. However, when applying HRTF interpolationin the spatially continuous spherical harmonics (SH) domain, the number of measured HRTFs limits the maximal
accessible SH order. This results in a restricted spatial resolution and can cause perceptual artefacts like coloration or localization errors. In a previous publication we presented the SUpDEq method (Spatial Upsamplingby Directional Equalization), which reduces these artifacts by a directional equalization based on a sphericalhead model prior to the SH transform. This removes direction-dependent temporal and spectral components
and thus reduces the spatial complexity of the HRTF set enabling improved interpolation of HRTFs already at
low SH orders. A subsequent de-equalization recovers energy in higher spatial orders that was discarded in thesparse HRTF set. In this study we analyze 96 individual HRTF sets and investigate to what extent the performance of SUpDEq, which we already analyzed for dummy heads, can be transferred to individual HRTF sets.
The results show that the SUpDEq method clearly outperforms common SH interpolation of individual HRTFs
with respect to the spectral structure and to modeled localization performance.
Keywords: Binaural hearing, Localization, Head-related transfer functions, Virtual acoustic environments

1

INTRODUCTION

A spatial presentation of sound sources is a fundamental element of virtual acoustic environments (VAEs). For
this, monaural and binaural cues, which are mainly caused by the shape of the pinna and the head, need tobe considered. While spectral information serves as main cue to determine elevation, differences between thesignals reaching the left and the right ear allow lateral localization. These binaural cues manifest in interauraltime differences (ITDs) and interaural level differences (ILDs). In many headphone-based VAEs, head-relatedtransfer functions (HRTFs) are applied to describe the sound incidence from a source, which is typically inthe far-ﬁeld, to the left and right ear incorporating both, monaural and the binaural cues. Generally, the useof individual HRTFs is advantageous, for example regarding localization accuracy in the median plane [5].
However, a high number of HRTFs is required to adequately capture the relevant cues for all directions ofincidence which makes the measurements time-consuming and tedious.
To allow an optimized interpolation between the measured directions, complete sets of HRTFs can be measured
on a spherical grid and described in the spherical harmonics (SH) domain [14, 12]. In this case a decomposition
into spherical base functions of different spatial orders N is applied, where higher orders correspond to a higher
spatial resolution. A subsequent inverse spatial Fourier transform at arbitrary angles can be used to recover a
spatially upsampled HRTF set. However, describing sparse HRTF sets in the SH domain results in a limited
spatial order and incorporates an incomplete description of the spatial properties resulting in spatial aliasingor truncation errors. To avoid spatial alisaing, an order N ≥ kr with k = ω/c, and r being the head radius
is required [11, 4]. For the full audio bandwidth ( f ≤ 20 kHz) this leads to N = 32 requiring at least 1089
measured directions when assuming r = 8.75 cm and c = 343 m/s.
Different studies analyzed the artifacts of sparsely measured HRTF sets or examined methods to reduce them4848

Equalization

De-Equalization

Equalization
DatasetDe-Equalization
DatasetL

SH
Coefficients

Nhigh
SparseHRTF Set

SparseSampling
Grid

L

R

ISFT

Dense
Sampling
Grid

R

ISFT

Nhigh

EqualizationSH
Coefficients

Nhigh

Nhigh

De-Equalization

L

L

R

R

SFT

Dense
HRTF Set

ISFT
Nlow

Nlow

Figure 1. Block diagram of the SUpDEq method. Left panel: A sparse HRTF set is equalized on the corresponding sparse sampling grid before transformed to the SH domain with N = Nlow . Right panel: The equalized
set is de-equalized on a dense sampling grid. If required, the resulting dense HRTF set can again be transformed
to the SH domain with N = Nhigh .
(e.g. [4, 3, 15, 7]). In this scope we recently introduced the SUpDEq (Spatial Upsampling by Directional Equalization) method [10], which removes frequency-dependent ITDs and ILDs as well as head-related elevationdependent spectral features from the HRTFs. SUpDEq applies a spectral division (equalization) of the HRTF
with a corresponding equalization function prior to the SH transform. A directional rigid sphere transfer function can be used here as equalization function, resulting in a signiﬁcantly reduced spatial order N. After spatial
upsampling, a de-equalization by means of a spectral multiplication with the same equalization function recovers
a spatially upsampled HRTF set. In this paper we analyze the SUpDEq method for a large number of measured
and simulated datasets.

2

METHOD

The SUpDEq method has been described in detail in [10]. In the following we thus brieﬂy outline the basic
concept. The corresponding block diagram is given in Fig. 1. First, the sparse HRTF set HHRTF measured at S
sampling points Ωs = {(φ1 , θ1 ), . . . , (φS , θS )} is spatially equalized with an appropriate equalization dataset HEQ
HHRTF,EQ (ω, Ωs ) =HHRTF (ω, Ωs ).
HEQ (ω, Ωs )(1)

While generally different equalization datasets can be applied, in this study a rigid sphere transfer function
is used [14, p. 227]. The radius of the sphere corresponds to the physical dimensions of a human head, as
ear position φ = ±90◦ and θ = 0◦ is considered. The rigid sphere transfer function can thus be regarded asa simpliﬁed HRTF set featuring basic temporal and spectral components, but leaving out information on the
shape of the outer ears or the ﬁne structure of the head. Thus, by the equalization a time-alignment of the
HRTFs is performed and direction-dependent inﬂuences of the spherical shape of the head are compensated.
As a consequence, the equalization with the rigid sphere transfer function considerably reduces the directional
complexity of HHRTF,EQ and thus the required order for the SH transform. As the equalization dataset HEQ canbe calculated based on an analytical description, it can be determined at a freely chosen maximal order, typically
Nhigh ≥ 35. The SH coefﬁcients for the equalized sparse HRTF set are obtained by applying the SH transform4849

on the equalized HRTFs up to an appropriate low maximal order Nlow which corresponds to the maximal orderHRTF,EQ is calculated on a dense sampling gridthat can be resolved by Ωs . Then an upsampled HRTF set H
Ωd = {(φ1 , θ1 ), . . . , (φD , θD )}, with D  S by using the inverse SH transform. Finally, HRTFs are reconstructedby a subsequent de-equalization by means of spectral multiplication with a de-equalization dataset HDEQHRTF,DEQ (ω, Ωd ) = H
HRTF,EQ (ω, Ωd ) · HDEQ (ω, Ωd ) .
H(2)

For de-equalization, again the rigid sphere transfer function is used in the present study. This last step recovers energy at higher spatial orders that was transformed to lower orders within the equalization. Again,
HRTF,DEQ holds if Nlow and Nhigh are chosen appropriately. Energy which, after the equalization, still
HHRTF = H
is apparent at high modal orders N > Nlow results in spatial aliasing and truncation errors as it is irreversibly
HRTF,DEQ . The following section analyzes
mirrored to lower orders N ≤ Nlow [4]. Thus we obtain HHRTF ≈ H
the inﬂuence of these deviations for individual datasets and investigates which advantage the SUpDEq method
provides compared to common (order-limited) SH interpolation without any pre- or postprocessing.

3

EVALUATIONIn previous publications [10, 9] we investigated the performance of SUpDEq for different artiﬁcial heads. However, one of the target applications of the SUpDEq method is the reduction of the measurement effort of individual HRTF sets. Thus, in this study we analyze the performance of the SUpDEq method for the HUTUBS
database which is online available on http://dx.doi.org/10.14279/depositonce-8487. The database contains of 96 acoustically measured and 96 numerically simulated datasets of full-spherical HRTFs (94 subjects
plus 2 repeated measurements of a human subject and an artiﬁcial head). For more detailed information on
the database please refer to [6]. We apply the HRTF sets to compare the performance of the SUpDEq method(de-equalized HRTFs) to HRTFs obtained with with strictly order limited SH interpolation, i.e., without any
pre- or post-processing before or after the SH transform. For this we generated SH coefﬁcients from 15 sparse
sampling grids equaling (limited) orders of N = 1 − 15. Thus, both order-limited (OL) and de-equalized (DEQ)
sets are based on the same respective sparse grid. To generate various sparse HRTF sets which we used as
input data for the evaluation, we simply spatially subsampled each individual reference set in the SH domain
by means of the inverse SH transform at the required directions. We calculated the optimal radius for the rigid
sphere model for each of the sets according to Algazi et al. [1] based on the individual anthropometry resulting
in an average value over the complete set of r = 9.1 cm (SD = 0.23 cm).
3.1 Spectral differences
First we analyze the spectral deviations to the reference set as a function of N on various test sampling gridswith T sampling points Ωt = {(φ1 , θ1 ), . . . , (φT , θT )}. For this the frequency-dependent spectral differences per
sampling point were calculated in dB as
Δg(ω, Ωt ) = 20lg

| HHRTF,REF (ω, Ωt ) |,| HHRTF,TEST (ω, Ωt ) |(3)

where HHRTF,REF is the left ear HRTF extracted from the reference set and HHRTF,TEST the one extracted fromthe order-limited or the de-equalized datasets at the sampling point Ωt . Then, the absolute value of Δg(ω, Ωt )was averaged across the temporal frequency ω to obtain one value ΔGsp (Ωt ) (in dB) per sampling point
ΔGsp (Ωt ) =1 nω
∑ | Δg(ω, Ωt ) |,nω ω=1(4)

across all sampling points Ωt to obtain the frequency-dependent measure ΔG f (ω) (in dB)
nΔG f (ω) =

4850

1 Ωt
| Δg(ω, Ωt ) |,nΩt Ω∑
t =1

(5)

8

8

OL
DEQ

6

G [dB]

G [dB]

6
4
2

4

5

10

0
1

15

10

Order [N]

(a)

(b)

8

OL-N4
OL-N7
OL-N10
OL-N13
DEQ-N4
DEQ-N7
DEQ-N10
DEQ-N136205Order [N]Gf [dB]Gf [dB]

6

4
2

0
1

8

OL
DEQ

4

15

OL-N4
OL-N7
OL-N10
OL-N13
DEQ-N4
DEQ-N7
DEQ-N10
DEQ-N13

2

1k

0

10k

1k

10k

Frequency [Hz]Frequency [Hz]

(c)

(d)

Figure 2. Spectral differences in dB (left ear) between the reference HRTF sets and the order-limited (OL) or
de-equalized HRTF sets (DEQ), both based on the respective sparse set, averaged over all 96 datasets. Additionally the standard deviations are plotted (shaded). The left row (a,c) illustrates the results for the simulated
datasets, in (b,d) the ones for the measured datasets are given. In (a,b) the spectral differences ΔG averaged
over the full audio bandwidth across N for order-limited datasets (red) and the de-equalized datasets (blue) are
given, in (c,d) the frequency-dependent spectral differences ΔG f (ω) for N = 4, 7, 10, 13 (color saturation).
and across ω and Ωt , resulting in a single value ΔG (in dB) describing the spectral differencen

ΔG =1 1 Ωt nω
∑ | Δg(ω, Ωt ) | .
nΩt nω Ω∑
t =1 ω=1(6)

Finally, the average values and standard deviations over all 96 datasets were calculated for the simulated andthe measured datasets.
Fig. 2 (a,b) show the spectral differences ΔG across N for order-limited interpolation and the SUpDEq method(de-equalized datasets) over the full audio bandwidth using the reference Lebedev2702 grid as test sampling gridΩt . The SUpDEq method clearly outperforms the order-limited interpolation both for the simulated and themesaured HRTF sets. The spectral differences are about 2 − 3 dB lower than for order-limited interpolation.
Fig. 2 (c,d) show the frequency-dependent spectral differences ΔG f (ω) at N = 4, 7, 10, 13. Generally, the spectral
differences are quite small at low frequencies. For order-limited interpolation they suddenly rise within one
octave from about 2 dB up to about 5 dB or more above a speciﬁc alias frequency. For the SUpDEq method,
however, the spectral differences show a much more gentle rise. The differences exceed 2 dB for frequencies

4851

g

90

Elevation [°]

4
6

-90
180

270

0

90

180

2
4

0
6

-45

8

-90
180

10

8

270

Azimuth [°]

90

180

10

6

-45

8

90

180

10

[
0
2

45

Elevation [°]

4

0

g

90

2

45

Elevation [°]

[
0

0

90

(b)
g

270

0

Azimuth [°]

(a)

-90
180

[
0

45

0

-45

g

90

2

45

Elevation [°]

[
0

4

0
6

-45

-90
180

8

270

0

Azimuth [°]Azimuth [°]

(c)

(d)

90

180

10

Figure 3. Spectral differences ΔGsp (Ωt ) per sampling point for order-limited interpolation (a,b) and for the
SUpDEq method (c,d) at N = 4 and f ≤ 10 kHz averaged over all 96 datasets. The left row (a,c) shows theresults for the simulated datasets, the right row (b,d) the results for the measured ones.
above 3 kHz for N = 4, while differences stay below 2 dB for orders of N ≥ 10 up to 10 kHz (DEQ).
Fig. 3 concludes the spectral analysis and shows the spectral differences ΔGsp (Ωt ) per sampling point at N = 4,f ≤ 10 kHz, and a full spherical test sampling grid Ωt with a resolution of 1◦ in azimuth and elevation. As
depicted in Fig. 3 (a,b), the order-limited interpolation results in distinct spectral differences spread over theentire angular range. On the contrary, Fig. 3 (c,d) shows that for the SUpDEq method the spectral differences
are mainly located at contralateral directions. At frontal directions, where order-limited interpolation typically
performs badly, the SUpDEq method shows good results. The same can be observed for various ipsilateral
directions. The spectral differences are generally higher for order-limited interpolation, with a maximum of
about ΔGsp (Ωt ) = 10.4 dB at φ = 262◦ and θ = −10◦ averaged over all subjects for the simulated datasets.
For these datasets applying the SUpDEq method results in a maximal spectral difference ΔGsp (Ωt ) of 6.6 dB
at φ = 257◦ and θ = 2◦ . Finally, Fig. 3 (b,d) show the same trend for the measured datasets, but reveal large
deviations for the downward directions. This is caused by the acoustic shadowing of the measurement equipment
and is described in detail in [6].
3.2 Localization performance
To compare the localization performance of order-limited HRTFs and de-equalized HRTFs in the median sagittal
plane, we used the model from Baumgartner et al. [2] which compares the spectral structure of a reference
HRTF set to a set of test HRTFs. Based on a probabilistic estimate of the perceived sound source location,
the model determines the polar RMS error which describes the expected angular error between the actual and

4852

15

OL
DEQ

10

PE [°]

PE [°]

15

5

5

10

0
1

15

(a)

(b)

20

15

15

10

15

10
5

5

10

0
1

15

5

10

Order [N]Order [N]

(c)

(d)

15

10

10

LE[°]

15

5

0
1

10

Order [N]

20

0
1

5

Order [N]5

LE [°]

10

5

QE [%]

QE [%]

0
1

OL
DEQ

15

5

5

10

15

0
1

5

10

Order [N]Order [N]

(e)

(f)

15

Figure 4. Absolute polar error difference ΔPE (a,b), quadrant error difference ΔQE (c,d), and lateral error
difference ΔLE (e,f) over SH order N for order-limited interpolation (red) and the SUpDEq method (blue)
averaged over the 96 individual datasets. Additionally, the standard deviations are shown (shaded). In theleft row (a,c,e) the results for the simulated HRTF sets are shown, in the right row (b,d,f) the results for themeasured HRTF sets.
perceived source positions. Additionally, it determines the quadrant error rate which speciﬁes the front-backand up-down confusions. Regarding the localization performance in the horizontal plane, we used the modelfrom May et al. [8] which weighs the frequency-dependent binaural cues (ILDs, ITDs) to estimate the azimuthal4853

position of a sound source. A lateral error can be calculated by comparing the intended and the estimated source
position. For the analysis of both models we used the Auditory Modeling Toolbox (AMT) [13]. The procedurefor determining the errors has been described in detail in [10] and can be outlined as follows. To estimate
median sagittal plane localization performance, we used a test sampling grid Ωt with φ = {0◦ , 180◦ } and −30◦ ≤θ ≤ 90◦ in steps of 1◦ , and assumed a median listener sensitivity of S = 0.76 (according to Baumgartner et al.
[2]). For the horizontal plane localization performance, we used a test sampling grid with φ = ±90◦ in steps of5◦ . We determined the absolute polar error difference (PE in degree)
ΔPE =| PEREF − PETEST |,(7)

the absolute quadrant error difference (QE in percent)
ΔQE =| QEREF − QETEST |,(8)

as well as the absolute lateral error difference (LE in degree)
ΔLE =

1 T
∑ | LEREF (Ωt ) − LETEST (Ωt ) |,
T t=1

(9)

for each order N with the subscripts REF describing the referene dataset and TEST the dataset under test.
Again we calculated the averages and standard deviations over all datasets separated for the simulated and themeasured sets.
As plotted in Fig. 4 (a–d), in the median sagittal plane the order-limited interpolation leads both for the simulated and the measured datasets to higher errors than the SUpDEq method. High-frequency deviations of theorder-limited HRTFs affect spectral cues which are relevant for sagittal plane localization. For the de-equalized
datasets, ΔPE decreases with increasing order N, ΔPE ≤ 2◦ holds for N ≥ 4. Thus the spectral cues seem to
be mostly unimpaired here. The extent of the quadrant error ΔQE varies greatly between the measured andthe simulated sets and lies for the order-limited sets between 4 % (simulated) and 10 % (measured) at N ≥ 7.
However, for the de-equalized datasets, ΔQE is below 2 % at N ≥ 7. Generally, in the median sagittal plane
the average errors are much higher for the measured datasets than for the simulated ones. This is probably a
result of the measurement inaccuracies for downward directions, which as well have been observed in Sec. 3.
In Fig. 4 (e–f) the localization performance in the horizontal plane is shown. Here the order-limited interpolation
performs quite well, even though lateral errors are distinctly ampliﬁed at orders N ≤ 3. This might be causedby strong pre-ringing artifacts causing wrong ITDs, as already discussed in [10]. The SUpDEq method leads tohardly any increase in lateral error over the entire tested range of N.

4

CONCLUSION

In this paper we analyzed the performance of the SUpDEq method for spatial upsampling of individual sparse
HRTF sets. Regarding the spectral structure, the deviations from the reference HRTF set are signiﬁcantly smaller
for the SUpDEq method than for order-limited interpolation. The average difference is about 2 dB, both for thesimulated and the measured datasets. Furthermore, the analysis of the spectral diffences showed for the SUpDEq
methods a much more gentle rise over frequency than for the order-limited interpolation. Finally, the spectral
differences induced by the SUpDEq method are mainly at contralateral directions, while the differences due toorder-limited interpolation spread over the entire angular range, with distinct clusters at frontal and contralateral
directions. Regarding the modeled localization performance the SUpDEq method performed better in both planes
because spectral and binaural cues are less impaired in comparison to the order-limited interpolation.
Generally, the evaluation showed that the results found for dummy heads in [10] can be generalized to individually measured or simulated datasets. Thus, the SUpDEq approach can help closing the gap between a practicaland fast measurement procedure and sufﬁcient accuracy of the upsampled HRTF set. However, for such a simpliﬁed procedure other inﬂuencing factors like e.g. the elimination of room reﬂections [9] or the compensationof small displacements of the human head during the measurement need to be considered.

4854

The research presented in this paper has been funded by the German Federal Ministry of Education and Research. Support Code: BMBF 03FH014IX5-NarDasS. A Matlab-based implementation of the SUpDEq method
is available on https://github.com/AudioGroupCologne/SUpDEq.

REFERENCES[1] V. Algazi, C. Avendano, and R. O. Duda. Estimation of a Spherical-Head Model from Anthropometry. J.
Audio Eng. Soc., 49(6):472 – 479, 2001.
[2] R. Baumgartner, P. Majdak, and B. Laback. Modeling sound-source localization in sagittal planes forhuman listeners. J. Acous. Soc. Am., 136(2):791–802, 2014.
[3] Z. Ben-Hur, F. Brinkmann, J. Sheaffer, S. Weinzierl, and B. Rafaely. Spectral equalization in binauralsignals represented by order-truncated spherical harmonics. The Journal of the Acoustical Society of America,
141(6):4087–4096, 2017.
[4] B. Bernschütz, A. Vázquez Giner, C. Pörschmann, and J. M. Arend. Binaural reproduction of plane waves
with reduced modal order. Acta Acustica united with Acustica, 100(5):972–983, 2014.
[5] J. Blauert. Spatial Hearing - The Psychophysics of Human Sound Localization. MIT Press, Cambridge, MA,
revised edition, 1996.
[6] F. Brinkmann, M. Dinakaran, R. Pelzer, P. Grosche, D. Voss, and S. Weinzierl. A cross-evaluated databaseof measured and simulated HRTFs including 3D head meshes, anthropometric features, and headphone
impulse responses. Journal of the Audio Engineering Society, in press, 2019.
[7] F. Brinkmann and S. Weinzierl. Comparison of head-related transfer functions pre-processing techniques
for spherical harmonics decomposition. In Proceedings of the AES International Conference on Audio for
Virtual and Augmented Reality, pages 1–10, 2018.
[8] T. May, S. Van De Par, and A. Kohlrausch. A probabilistic model for robust localization based on abinaural auditory front-end. IEEE Trans. Audio, Speech, Lang. Process., 19(1):1–13, 2011.
[9] C. Pörschmann and J. M. Arend. Obtaining Dense HRTF Sets from Sparse Measurements in Reverberant
Environments. In Proceedings of the AES Conference on Immersive and Interactive Audio, 2019.
[10] C. Pörschmann, J. M. Arend, and F. Brinkmann. Directional Equalization of Sparse Head-Related Transfer
Function Sets for Spatial Upsampling. IEEE/ACM Trans. Audio, Speech, Lang. Process., 27(6):1060 – 1071,
2019.
[11] B. Rafaely. Analysis and Design of Spherical Microphone Arrays. IEEE Transaction on Speech and Audio
Processing, 13(1):135–143, 2005.
[12] B. Rafaely. Fundamentals of Spherical Array Processing. Springer-Verlag, Berlin Heidelberg, 2015.
[13] P. Søndergaard and P. Majdak. The Auditory Modeling Toolbox. In J. Blauert, editor, The Technology ofBinaural Listening, pages 33–56. Springer-Verlag, Berlin Heidelberg, 2013.
[14] E. G. Williams. Fourier Acoustics - Sound Radiation and Nearﬁeld Acoustical Holography. Academic Press,
London, UK, 1999.
[15] M. Zaunschirm, C. Schoerkhuber, and R. Hoeldrich. Binaural rendering of Ambisonic signals by HRIR
time alignment and a diffuseness constraint. J. Acous. Soc. Am., 143(6):3616 – 3627, 2018.

4855