A localization algorithm based on head-related transfer functions
Maike GERHARD(1) , Patrick SCHILLBERG(2) , Hermann WAGNER(2) , Hartmut FÜHR(3)(1) Lehrstuhl A für Mathematik, RWTH Aachen University, Germany, maike.gerhard@matha.rwth-aachen.de(2) Institute of Biology, RWTH Aachen University, Germany(1) Lehrstuhl A für Mathematik, RWTH Aachen University, Germany

Abstract
A localization algorithm based on head-related transfer functions (HRFTs) is introduced. The algorithm minimizes the Euclidean distance to direction-dependent binaural signal spaces, in order to find an estimator of theactual sound source position. Numerical experiments with different target sounds (clicks, white noise, speech,
rustling of leaves) and two signal-to-noise ratios (10 dB and 0 dB), where the masker is non-spatialized white
noise, show that this projection algorithm outperforms template matching and the cross-channel algorithm innearly all experimental conditions. In particular, the algorithm is robust to changes of the emitted signal’s phase
spectrum, unlike template matching.
For white noise maskers, it is possible to compute SNR dependent estimates for the error probability in the taskof discriminating two directions, based on the associated HRTFs. We present simulations that demonstrate theprecision of the estimates. We show how these probabilities can be employed as a means to mathematicallyanalyze HRTFs, in particular with the aim of predicting localization performance from the HRTF data set.
Keywords: Sound localization, HRTFs

1

INTRODUCTION

For the localization of sound source positions, humans and animals use various spatial cues contained in the twosound signals reaching the left and right ear. It is currently assumed that the interaural time difference (ITD),
the interaural level difference (ILD), and high frequency spectral notches are most important among these cues
[1, 2]. All these cues and perhaps even more spatial information can be extracted from head-related impulse
responses (HRIRs) [3], which encode the linear direction dependent acoustic filtering of sounds by the listener’souter ears, head, and body. This acoustic filtering models physical phenomena such as interference, reflection,and diffraction and depends on both the direction of arrival and the listener’s anthropometric parameters. The
frequency domain representation of a HRIR is the head-related transfer function (HRTF), obtained as the Fouriertransform of the HRIR.
In the past decades, HRTFs have been investigated from a variety of perspectives. A topic related to all directions of research is the mathematical analysis of HRTF data. The present paper is a further contribution to thisresearch by deepening the theoretical understanding of the acoustic information inherent in head-related transfer
functions. Starting from first principles, a sound localization algorithm based on head-related transfer functions
is developed. Numerical experiments with different target sounds show, that this algorithm performs highly efficiently, also in comparison to other algorithms. We also introduce a theoretical estimate for the error probability
in the task of discriminating two directions. This estimate depends on the signal-to-noise ratio (SNR) of theemitted signal as well as on the HRTFs associated to the two positions taking part in the discrimination task.
Section 2 introduces a new HRTF based sound localization algorithm. At the beginning of this section we
describe the model of acoustic signal transmission that forms the basis of the considerations throughout this paper. We also introduce the mathematical framework by giving the necessary definitions. Afterwards we present
a mathematical development of the mentioned localization algorithm. The algorithm performs orthogonal projections of the incoming binaural sound onto HRTF-dependent binaural subspaces and minimizes the Euclidean
distance between these projections and the original sound.

3379

Section 3 compares the projection algorithm to two other algorithms: The template matching algorithm and thecross-channel algorithm introduced by MacDonald [5]. First we give an overview of the operating modes of thetwo competing algorithms. Then we compare the performances of the three algorithms by conducting numericalexperiments implemented in MATLAB. Two different signal-to-noise ratios (10 dB and 0 dB) and several target
sounds were used in these experiments to simulate various localization scenarios.
In section 4 we consider 2-choice discrimination tasks, where the listener has to decide which of two given
positions is the actual sound source position. If white noise maskers are used in such paradigms, it is possible
to compute an upper bound for the error probability. This bound depends on the HRTFs associated to the two
positions as well as on the SNR of the incoming sound. We show, that this bound can be used to predict the
localization performance of individuals.
The paper’s conclusion in section 5 sums up the main results of our work.

2

A NOVEL HRTF-BASED LOCALIZATION ALGORITHM

In this section we introduce a new algorithm for sound localization, called projection algorithm. Prior to that,
we present the mathematical framework.
2.1 Model assumptions and mathematical basics
The basic problem of sound localization can be formulated as follows: An acoustic signal x is emitted fromdirection P. The binaural sound (xl , xr ) arriving at the left ear and the right ear is then described by theequationsxl = x ∗ Hl,P + εl , xr = x ∗ Hr,P + εr .
Here• Hl,P and Hr,P denote the left and right head-related impulse responses associated to position P,• the symbol ∗ denotes convolution, and• εl , εr is interfering noise that is not related to the emitted signal x.
The task of sound localization consists in estimating P from the pair (xl , xr ), usually without any prior knowledge of the emitted signal, or of the interfering noise.
When we study sound localization as a discrimination task, we ask whether two directions P and Q can be
sufficiently discriminated, for example in the setting of a behavioral experiment. Since our analysis is purely
based on the available physical acoustic information contained in the HRIRs, we do not take into account any
prior knowledge of the emitted signal or the masker. Thus, in view of the above equation, the pertinent question
becomes whetherx ∗ Hl,P ≈ y ∗ Hl,Qandx ∗ Hr,P ≈ y ∗ Hr,Q
hold with a suitably chosen signal y emitted from Q. If this is the case, it becomes very difficult to decide
whether the signal arrives from P or Q.
For the mathematical modeling we assume that all acoustic objects (signals, HRIRs, noise, etc.) have time domain representations of the same length N ∈ N and are measured with the same sampling rate S = ∆1t . Hence thespace of acoustic signals is given by l 2 (N) ∼= CN with N := {0, 1, . . . , N − 1}. The kth entry of each vector cor2
responds to the time k · ∆t . The energy of a signal x is given by its squared Euclidean norm kxk2 = ∑N−1
k=0 |x(k)|and the (cyclic) convolution of two such signals is given by
N−1(x ∗ y)(k) = ∑ x(n)y((k − n) mod N),
n=0

3380

k ∈ N.

The Fourier transform F converts signals from time domain to frequency domain:
N−1

F (x)(k) = x̂(k) = ∑ x(n)e−2πikn/N ,F : l 2 (N) → l 2 (N),k ∈ N.

n=0

Hence all frequency domain representations also have length N. Convolution in time domain corresponds tomultiplication in frequency domain and F conserves energy, at least up to a fixed factor:\(x∗ y) = x̂ŷ

and Nkxk2 = kx̂k2 .

Here and in the following multiplication and division is componentwise.
2.2 The projection algorithm
As mentioned before the task of sound localization consists in estimating the actual sound source position P0
from the binaural sound(xl , xr ) = (x ∗ Hl,P0 + εl , x ∗ Hr,P0 + εr ) ∈ l 2 (N) × l 2 (N).
Since we assume that the emitted signal x and the noise (εl , εr ) are unknown, the HRTFs have to be used forsolving this task. We utilize the fact, that in the case of a noise-free signal transmission all signals arrivingfrom the very same position have the same pattern:
Definition 2.1. Given a position P, the space composed of all binaural sounds arriving from that position P is denotedby HP , i.e.

HP = (x ∗ Hl,P , x ∗ Hr,P ) | x ∈ l 2 (N) = (x̂Ĥl,P , x̂Ĥr,P ) | x̂ ∈ l 2 (N) .
For every position P the set HP is a subspace of the space of all binaural sounds l 2 (N) × l 2 (N). In addition, fora second position Q different from P one expects HQ 6= HP , since the HRTFs associated to the two positions
will not be equal, or in more precisely: the quotients Ĥl,P /Ĥr,P and Ĥl,Q /Ĥr,Q will not be equal. That means in
the noise-free case, the HRTFs describe subspaces of the space of spatialized sounds. For correct localization of
an incoming sound it is then sufficient to determine which of these subspaces contains the sound. Unfortunately
the noisy incoming signal cannot be expect to lie in any of the subspaces HP . However, we can compute theEuclidean distances between the incoming sound and the different subspaces HP , and minimize of all available
positions P to obtain an estimator of the actual sound source position:
Algorithm 2.2 (Projection algorithm). Input: A binaural sound (xl , xr ) ∈ l 2 (N) × l 2 (N).
Output: The estimator Pest of the sound source position P0 given by
Pest = argminP dist ((xl , xr ), HP ) .
The theorem below shows, that the distance dist((xl , xr ), HP ) is computable in an explicit and efficient way.
The overall minimization task is solved subsequently by running a simple loop through all available directions
P.
Theorem 2.3. For a position P the quantity dist((xl , xr ), HP ) can be computed explicitly.
Proof. We have to solveminxP ∈l 2 (N)(xl , xr ) − (xP ∗ Ĥl,P , xP ∗ Ĥr,P )2.

Because of the energy conservation of the Fourier transform, we can solve this in the frequency domain, where
it becomesminx̂P ∈l 2 (N)(x̂l , x̂r ) − (x̂P Ĥl,P , x̂P Ĥr,P )

2

= min

∑ (x̂l (k), x̂r (k)) − x̂P (k)(Ĥl,P (k), Ĥr,P (k))

x̂P ∈l 2 (N) k∈N

2

.

3381

Now the minimization can be done for every term of the sum (i.e. for every frequency) individually. For everyk ∈ N we project the vector (x̂l (k), x̂r (k)) ∈ C2 orthogonal onto the closed subspace{x̂ · (Ĥl,P (k), Ĥr,P (k)) | x̂ ∈ C} ⊂ C2 .
The solution is given byx̂P (k) =(x̂l (k)Ĥl,P (k) + x̂r (k)Ĥr,P (k))|Ĥl,P (k)|2 + |Ĥr,P (k)|2for k ∈ N,

where the overline denotes the complex conjugate.
We note that in the noise-free case, directions P and Q can only be confused by the algorithm if there aresignals xP 6= 0 and xQ 6= 0 such that xP emitted from P coincides with xQ emitted from Q. Mathematically, thiscan be expressed via the following equivalent conditions:HP ∩ HQ ) {0}⇔∃ k∈N:

Ĥl,P (k) Ĥl,Q (k)=.
Ĥr,P (k) Ĥr,Q (k)Note that the right hand side of the equivalence revals that typical signals leading to nonempty intersections are
tones.

3

PERFORMANCE OF THE PROJECTION ALGORITHM

In this section we compare the performance of the projection algorithm to the performance of two other soundlocalization algorithms: Template matching and the cross-channel algorithm. Note that all of these algorithms
require the complete HRIR data set of a listener in order to localize the sound source.
3.1 Template matching
If the emitted signal x is a click, the binaural sound (xl , xr ) arriving at the ears is up to a fixed factor equal
to the HRIRs (Hl,P0 , Hr,P0 ) associated to the actual sound source position P0 . This observation leads to a simple
localization technique called template matching, which determines which pair of HRIRs (Hl,P , Hr,P ) is mostsimilar to the incoming sound. All correlation-based localization algorithms perform some sort of templatematching. If the HRIRs themselves are used as templates, this investigation can be conveniently performed in
the frequency domain. In general the energy of different HRTF pairs is also different. To ensure that this does
not influence the outcome, all HRTFs were normalized before comparison with the incoming sound. So thetemplate related to position P is given by(Tl,P , Tr,P ) = q(Ĥl,P , Ĥr,P )
Ĥl,P

2

+ Ĥr,P2.

The real part of a scalar product hv, wi is an often used measure for the similarity of vectors v and w, which isalso employed by the template matching algorithm.
Algorithm 3.1 (Template matching). Input: A binaural sound (xl , xr ) ∈ l 2 (N) × l 2 (N).
Output: The estimator Pest of the sound source position P0 given by
Pest = argmaxP ℜ hTl,P , x̂l i + hTr,P , x̂r i .
We are not aware of sources studying sound localization by directly using HRTFs as templates for sound localization. However, the localization algorithm presented in [4] can be interpreted as as a template matchingalgorithm applied to suitably preprocessed binaural signals.

3382

3.2 The cross-channel algorithm
The cross-channel algorithm was developed by Justin A. MacDonald [5]. This algorithm is based on a simple
yet elegant idea: If a signal x is emitted from position P0 the noise-free sounds arriving at the ears arexl = x ∗ Hl,P0

andxr = x ∗ Hr,P0 ,respectively. If we now convolve the left sound xl with the right ear HRIR associated to position P0 and theright sound xr with the left ear HRIR associated to P0 , we get equal results, because of the commutativity andassociativity of convolution:
xl ∗ Hr,P0 = x ∗ Hl,P0 ∗ Hr,P0 = x ∗ Hr,P0 ∗ Hl,P0 = xr ∗ Hl,P0 .
For any other position P 6= P0 , and general signals x, this cannot be expected:xl ∗ Hr,P = x ∗ Hl,P0 ∗ Hr,P 6= x ∗ Hr,P0 ∗ Hl,P = xr ∗ Hl,P .
In the presence of noise, the equation xl ∗ Hr,P0 = xr ∗ Hl,P0 is no longer true. Nevertheless it can be expectedthat the correlation between xl ∗ Hr,P0 and xr ∗ Hl,P0 is larger than the correlation between xl ∗ Hr,P and xr ∗ Hl,P .
Algorithm 3.2 (Cross-channel algorithm). Input: A binaural sound (xl , xr ) ∈ l 2 (N) × l 2 (N).
Output: The estimator Pest of the sound source position P0 given by
Pest = argmaxP r xl ∗ Hr,P , xr ∗ Hl,P ,where r denotes the Pearson correlation coefficient.
3.3 Comparison
We used a numerical localization experiment implemented in MATLAB (The MathWorks, Natick, MA) to compare the performance of the three algorithms. The spatial binaural sounds (xl , xr ) used as inputs were generated
by convolving a signal x with the HRIRs Hl,P0 and Hr,P0 associated to the simulated sound source position P0 .
Afterwards non-spatialized white noise was added to the convolution product, to achieve the desired SNR. We
used publically available HRIRs from the CIPIC database [6]. This database includes 1250 positions whose
azimuth and elevation angles are measured in a head-centered interaural poplar coordinate system. All positions
are located on a sphere with radius 1 m around the listeners head. The inputs (xl , xr ) differed from each otherwith respect to• the source position P0 : 35 sound source positions ranging in azimuth from -80 deg to 80 deg in steps of
40 deg (5 values for the azimuth) and in elevation from -45 deg to 225 deg in steps of 45 deg (7 values
for the elevation) take part in the experiment.
• the SNR: 10 dB or 0 dB.
• the emitted signal x: We used 8 different signals, each with a length of 100 ms.
1.) Signal with a phase spectrum equal to 0 and a flat amplitude spectrum over the frequency rangeΩlow = [0 Hz, 4000 Hz].
2.) The same as in 1.), but with a different frequency range Ωhigh = [4000 Hz, 22000 Hz].
3.) The same as in 1.), but with a different frequency range Ωbroad = [0 Hz, 22000 Hz].
4.) The same as in 1.), but with a random phase spectrum.
5.) The same as in 2.), but with a random phase spectrum.
6.) The same as in 3.), but with a random phase spectrum.

3383

7.) Speech (female voice).
8.) Rustling of leaves.
Each signal was combined with both SNRs and each signal/SNR-combination was emitted 50 times per position. The noise was random in each trial, but all algorithms were tested with exactly the same inputs.
As error measure we use the spherical error, which is the absolute value of the angular distance between
the actual sound source position P0 and the estimated position Pest in degrees. We performed the experimentdescribed above with the HRIR sets of six different subjects. The mean spherical errors presented in table 1
are averaged across subjects and source positions.
Table 1. Mean spherical errors for inputs with a SNR of 10 dB (left) and 0 dB (right)
Signal
1.)
2.)
3.)
4.)
5.)
6.)
7.)
8.)

Projection
0.02
0.02
0.02
0.05
0.01
0.01
1.08
0.01

10 dB
Template Matching
40.56
1.58
0.00
73.15
67.25
66.75
59.71
54.23

Cross-channel
7.96
2.01
0.32
15.84
1.53
0.31
49.26
20.45

Projection
3.19
0.51
0.46
8.27
1.45
0.51
18.60
5.34

0 dB
Template Matching
40.69
1.68
0.00
77.31
68.50
67.32
74.95
59.50

Cross-channel
26.35
11.43
4.62
35.19
13.39
5.42
58.43
45.15

The table shows, that the projection algorithm outperforms the other algorithms in nearly all conditions. The
only exception is signal 3.), which is best estimated by template matching. This is not surprising, since signal
3.) is a click and template matching is designed for such signals. A comparison of signals 1-3.) to signals 4-6.)
emphasizes, that the success of template matching highly depends on the phase spectrum of the emitted signal,
which is not the case for the projection algorithm. Note, that even at a low SNR of 0 dB the mean errors of
the projection algorithm are with one exception (speech) smaller than 10 degrees.

4

ERROR PROBABILITY IN DISCRIMINATION TASKSIn this section we consider 2-choice discrimination tasks, where the listener has to decide whether the signal
arrives from the actual sound source position P or from a second position, which we call Q. Assuming that themasker is white noise, we can use probability theory to compute an upper bound for the error probability insuch tasks. This bound depends on the positions’ HRTFs and the sound’s SNR.
In this context we need a stochastic rather than a deterministic model: In the frequency domain the binauralsound (x̂l , x̂r ) ∈ l 2 (N) × l 2 (N) is given byx̂l = Ĥl,P · x̂ + σ ε̂landx̂r = Ĥr,P · x̂ + σ ε̂r ,where σ ∈ R≥0 and ε̂l = (ε̂l,0 , . . . , ε̂l,N−1 ), ε̂r = (ε̂r,0 , . . . , ε̂r,N−1 ) are random vectors withℜ(ε̂l,i ), ℑ(ε̂l,i ), ℜ(ε̂r,i ), ℑ(ε̂r,i )iid∼

N (0, α 2 ) for all i ∈ N.

This means that the real and imaginary parts of the scalar components of ε̂l and ε̂r are independent and identically normal distributed with mean 0 and variance α 2 . We make two assumptions:k(Ĥl,P x̂, Ĥr,P x̂)k2 = 1

3384

and E(k(ε̂l , ε̂r )k2 ) = 1,

where E denotes the expected value. Because of the second assumption the SNR of the binaural sound isdetermined by the scalar σ . Let P(Q|P) denote the error probability, i.e.
P(Q|P) =b the probability that the listener chooses Q althoug the signal arrives from P.
Then the following theorem holds true.
Theorem 4.1. Let {u1 , . . . , uN } be an orthonormal basis of HQ⊥ and {v1 , . . . , vN } an orthonormal basis of HP⊥ . Then
we have
!
(


 )
1 N N
σ 4 1
σ 2
1
2
−
+
P(Q|P) ≤ UB(P, Q) := min 1,∑ ∑ |hui , v j i|
2N 2N 2 i=1
kzk
N kzk
j=1

with z := ΠH ⊥ Ĥl,P x̂, Ĥr,P x̂ and ΠH ⊥ denotes the orthogonal projection onto HQ⊥ .
Q

Q

For P = Q we define UB(P, Q) = 1. To examine whether the quantity UB(P, Q) can be used to predict theresults of localization experiments, we compare UB(P, Q) to the performance of the projection algorithm. Two
examples are shown in figures 1 and 2.
1

0.14

0.9

200

200

0.12

0.8
0.1

0.08

100
0.06

50

0.7

150

Elevation of Q

Elevation Pest

150

0.6

100

0.5
0.4

50
0.3

0.04

0.2

0

0.02

0
0.1

0

-80

-60

-40

-20

0

20

40

60

80

Azimuth Pest(a) Localization performance

0

-80

-60

-40

-20

0

20

40

60

80

Azimuth of Q(b) UB

Figure 1. Comparison between the results of a localization experiment (a) and UB(P0 , Q) (b). The sound source
position is P0 = (80, 0), the emitted signals are of type 1.), and the SNR is -5 dB.
Part (a) of figure 1 shows the results of a localization experiment implemented in MATLAB. In this experiment
the source position P0 = (80, 0) was fixed. The emitted signals were of type 1.) (see section 3.3). The masker
was non-spatialized white noise and the SNR was -5 dB. We performed 200 localization runs with random
noise, respectively. The color map displays the relative frequency with that a position, indicated by x- andy-axis, was estimated as sound source position by the projection algorithm. Part (b) of the same figure shows
the values of UB(P0 , Q) for a fixed P0 = (80, 0). Azimuth and elevation of position Q are indicated by the xand y-axis, respectively. The value of UB(P0 , Q) is indicated by the color. For the calculation of UB(P0 , Q) we
used the same SNR (-5 dB) and signal (type 1.)) as in the localization experiment related to part (a), to ensurethat the comparison is valid.
Figure 2 is of the same type as figure 1. Here the sound source position is given by P0 = (0, 0) and the signals
are of type 3.). Furthermore the SNR is even smaller, namely -11 dB.

3385

1
0.18
0.9

200

200

0.16

0.8
0.14
0.7

150
0.12
0.1

100

0.08

50

0.06

Elevation of Q

Elevation Pest

150

0.6

100

0.5
0.4

50
0.3

0.04

0

0.2

0
0.02

0.1

0

-80

-60

-40

-20

0

20

40

60

80

Azimuth Pest(a) Localization performance

0

-80

-60

-40

-20

0

20

40

60

80

Azimuth of Q(b) UB

Figure 2. Comparison between the results of a localization experiment (a) and UB(P0 , Q) (b). The sound source
position is P0 = (0, 0), the emitted signals are of type 3.), and the SNR is -11 dB.
Clearly, both figures show pairs of very similar plots. This can also be confirmed quantitatively: The correlation
coefficients between the maps in figure 1 and 2 are given by 0.92 and 0.91, respectively. Nevertheless, thealgorithm performs a little better than predicted by UB. For example, the uncertainty in elevation shown in part(b) of figure 2 is less pronounced in part (a) of the same figure. That was to be expected, since UB is an upperbound for the error probability.

5

CONCLUSIONS

The main purpose of the present paper is the introduction of a new HRTF based sound localization algorithm.
The presented results stress that this algorithm exhibits a strong performance in comparison to established algorithms, especially at low signal-to-noise ratios. Furthermore a stochastic method for the estimation of the error
probability in discrimination tasks is presented. Initial experiments suggest that this method provides a validmeans for the prediction of localization performance.

REFERENCES[1] Hebrank, J.; Wright, D. Spectral cues used in the localization of sound sources on the median plane. The
Journal of the Acoustical Society of America, 56 (6), 1974, 1829-1834.
[2] Wightman, F. L.; Kistler, D. J. The dominant role of low-frequency interaural time differences in soundlocalization. The Journal of the Acoustical Society of America, 91 (3), 1992, 1648-1661.
[3] Blauert, J. Spatial Hearing: The Psychophysics of Human Sound Localization, The MIT Press, Cambridge(USA), Revised edition, 1996 .
[4] Baumgartner, R.; Majdak, P.; Laback, B. Modeling sound-source localization in sagittal planes for humanlisteners. The Journal of the Acoustical Society of America, 136 (2), 2014, 791-802.
[5] MacDonald, J. A. A localization algorithm based on head-related transfer functions. The Journal of the
Acoustical Society of America, 123 (6), 2008, 4290-4296.
[6] Algazi, V. R.; Duda, R. O.; Thompson D. M.; Avendano, C. The CIPIC HRTF Database. Proceedings 2001
IEEE Workshop on Applications of Signal Processing to Audio and Electroacoustics, Mohonk Mountain
House, New Paltz, NY, October 21-24, 2001, 99-102.

3386