A Physically Motivated Approach for Binaural Simulation of Moving Sound
Sources and Receivers
Christoph URBANIETZ(1) , Gerald ENZNER(2)(1) Ruhr University Bochum, Institute of Communication Acoustics, Germany, christoph.urbanietz@rub.de(2) Ruhr University Bochum, Institute of Communication Acoustics, Germany, gerald.enzner@rub.de

Abstract
In this paper we derive a physically motivated approach to handle the simulation of moving sound sources andreceivers including phenomena such as the Doppler effect in audio signal processing by using the concept ofretarded time in a very rigorous way. After recapitulating the physical basis of moving point-like sound sources
and receivers we derive an alternative formulation of the received signal based on superposition of Green’s
functions and transfer it to a system theory point of view. A first proof of the validity of this formulation is shownfor continuous-time systems. Afterwards we consider a discrete-time implementation and reveal the strengths
and the limits of this approach. Finally a binaural receiver is added to the model. The proposed formulation isanchored on a physical basis and therefore it handles physical effects inherently. Beside the correct reproductionof the Doppler frequency shift it also includes the physically attested amplitude shift and further phenomena.
Nevertheless the specific implementation for discrete-time signal processing can be challenging. We present a
possible implementation and verify some simulation examples against physically predictable features. It turns
out that the results of the presented approach and the physical prediction match with high precision.
Keywords: Dynamic Rendering, Binaural, Doppler Effect1

INTRODUCTION

Virtual acoustics is a key feature in virtual reality and multimedia applications [1] and has been investigated in
the past years in various ways. For the directional listening with headphones the use of head-related impulse
responses (HRIRs) is the state-of-the-art technique [2] although there is still much ongoing research in thisfield. The theory is well known for static scenes, where sound source and receiver have static positions andorientations and, therefore, the distance between sound source and receiver can be described as a delay-line andamplitude attenuation. Furthermore, many approaches have been presented how to handle motion of source andreceiver. The dynamic or time-varying acoustic environment, however, exhibits a number of sound features that
require particular care, e.g., Doppler shift, amplitude modulation, or spatial mislocalization of sound sources.
To address dynamic virtual acoustic environments there are phenomenological approaches which try to mimic
physical features without calculating the whole true physical process and, on the other side, physical simulationsof the whole acoustic field which suffer from a very high computational effort. The first kind calculates expectedfeatures such as Doppler shift or time delays and mimics it, e.g., by fractional delay lines that cover only a partof the physical truth [10, 11, 3, 5, 4]. The second kind uses comprehensive simulations of the acoustic field,e.g., by boundary elements methods [7] but is limited due to high computational effort.
The presented approach is positioned inbetween as it derives a point-to-point description for moving soundsources and receivers from a deep physical starting point. In contrast to work such as [8] where an analyticallycalculated transfer function based on the physical basis is used for the inverse process, our approach starts at amore basic level of physical description. The desired physical phenomena, such as the Doppler shift, amplitude
shift, or further small nonlinearities, are thus generated inherently. Basically the derived description is grounded
on a rigorous utilization of the concept of retarded time. In [12] it is shown that already pure azimuthal binaural
rendering with constant distances benefits from this treatment of retarded time. Nevertheless, in this paper, theconcept is derived from a more fundamental physical basis and also radial velocities are considered.

4833

~xr (t)R = c · (t − τ )θ
~xs (τ ) = vτ e~1

Source~xs (t) = vt~
e1v · (t − τ ) = vR/c~e1

xr,1 (t) − vt

R cos θ = xr,1 (t) − vt + MR

Figure 1. Sound source moving along the ~e1 axis, according to [9].

2

SOUND FIELD OF MOVING SOUND SOURCE AT MOVING RECEIVING POINT

In a first step we review the physical description of a moving sound source. In particular we will constrain
ourselves to point-like sound sources. We start with static source and receiver positions and follow [9] defining
a velocity potential ψ by ∂∂ ψxi = −ρ0 νi where ρ0 is the static pressure, and νi is the sound particle velocity inthe i-th spatial direction. Thus, according to linearized acoustic Euler equation,
∂ψ
p=
(1)
∂t
is the sound pressure. We define the sound source by a rate of mass flux q(t). In a scene with a statically
positioned source q(t) at the origin of the coordinate system we get the non-homogeneous wave equation1 ∂2∆ − 2 2 ψ = −q(t)δ (x1 )δ (x2 )δ (x3 )(2)c ∂twhere c is the speed of sound and δ (xi ) is the spatial delta distribution. The solution is given byqq(t − r/c)ψ=
, r = x12 + x22 + x32 .
4πr
Now we move the source with the velocity v along the ~e1 -axis as visualized in Fig. 1. Therefore we get1 ∂2∆ − 2 2 ψ = −q(t)δ (x1 − vt)δ (x2 )δ (x3 ).
c ∂t

(3)

(4)

Morse and Ingard [9] obtained, using Lorentz transform and multiple substitutions, the solutionψ=q(t − R/c)
4πR(1 − M cos θ )(5)where R is the distance between emission point ~xs (τ) at retarded time τ = t − R/c and receiving point ~xr (t) of
the sound received at time t. The angle θ is visualized in Fig. 1 and R has been found to beqM(x1 − vt) + (x1 − vt)2 + (x22 + x32 )(1 − M 2 )v, M= .
(6)R=
1 − M2c
2.1 Alternative way of the solution
For a more general trajectory of the source (and receiver) we now choose another description. The solution to
the non-homogeneous wave equation for acoustic waves


1 ∂2
−
∆
ψ = δ (~x)δ (t)(7)
c2 ∂t 24834

with a given point-like source in space and time is Green’s function
ψ = G(~x,t) =1 δ (t − |~x|/c)
.
4π
|~x|

(8)

For a source distribution with only a single point-like source q(t) positioned at ~xs (t)Q(~x,t) = q(t)δ (~x −~xs (t))(9)

we can get the velocity potential at receiving point ~xr as superposition of elementary solutions byq(τ)δ (~x −~xs (τ))
1δ (t − τ − |~xr −~x|/c) d3~x dτ
ψ(~xr ,t) =G(~xr −~x,t − τ)Q(~x, τ) d ~x dτ =4π|~xr −~x|
Z
q(τ)
1
=
δ (t − τ − |~xr −~xs (τ)|/c) dτ .
4π |~xr −~xs (τ)|
ZZ

ZZ

3

(10)
(11)

This new description potentially can be used with arbitrary trajectories. For first verification
we apply ~xs (τ) =qvτ~e1 and compare the result to (5). The argument of the delta distribution is t − τ −2 + x2 /c(xr,1 − vτ)2 + xr,2r,3which has its root at τ = t − Rc with R as defined in (6). Due to the derivative of the argument of the deltaq
2 + x2 /c) = 1 − M cos(θ ) we get during the integration an additional
distribution ∂∂τ (t − τ − (xr,1 − vτ)2 + xr,2r,3factor of 1/(1 − M cos θ ) and therefore the result

ψ(~xr ,t) =1 q(t − R/c)
4π R(1 − M cos θ )(12)is indeed equal to (5).

3

SYSTEM THEORY POINT OF VIEW

In the next step we transfer (11) into a system theory context. We absorb all constants in r0 and therefore
ψ(t) =Zr0
δ (t − τ − |~xr (t) −~xs (τ)|/c)q(τ) dτ
|~xr (t) −~xs (τ)|(13)

where we consider now q(t) to be an input signal to our system and ψ(t) to be the output signal at position~xr (t). Thus, we can identify a time varying, i.e., t-dependent, transfer function ash(~xr ,~xs ,t, τ) =r0
δ (τ − |~xr (t) −~xs (t − τ)|/c) .
|~xr (t) −~xs (t − τ)|(14)

When we start with a sound pressure signal s(t) which was measured with a microphone in a stationary caseand, therefore, was achieved by differentiating a velocity potential by time, we have vice versa to integrates(t) over time to get a mass flux equivalent q(t). The resulting sound pressure p(t) is achieved by a final
differentiation by
q(t) =Zt0s(τ)dτ , p(t) =

dψ
.
dt

(15)

3.1 Proof on specific example
We shall confirm that the proposed approach leads to reasonable results, especially that this approach reproduces
Doppler effect and amplitude modulation inherently. We therefore do an explicit calculation of a moving source

4835

vs

Receivervr

xs,0 = 0

Sourceyr,0Figure 2. Visualization of the tested scene.
and receiver along the first axis of the global coordinate system. The sound source is given by s(t) = sin(ω0t).
The receiver moves with velocity vr starting at xr,0 and the source moves with velocity vs starting at xs,0 = 0 as
visualized in Fig. 2. We further restrict to xs (t) = vst < xr (t) = xr,0 + vr t. The output signal received at xr isdp(t) =
dt

Z


 Zτ
r0
|xr (t) − xs (τ)|
s(τ 0 ) dτ 0 dτ
δ t −τ −|xr (t) − xs (τ)|
c

(16)

0


d c
r0
1

=
(1 − cos(ω0 τt ))x 
dt c − vs |xr (t) − xs (τt )| ω0
τt = t c−vr − r,0
c

c


c
c−vs



(17)

c(c − vr )r0
c2 (vr − vs )r0
1
sin(ω0 τt ) −(1 − cos(ω0 τt )).
(18)2(c − vs ) |xr (t) − xs (τt )|(c − vs )2 |xr (t) − xs (τt )|2 ω0
xr,0  cr
The root of the delta distribution argument τt = t c−v
c − cc−vs , called “retarded time”, is the emission timeof the sound received at time t and therefore |xr (t) − xs (τt )| is the distance between emission and receiving
point of this sound part. The first of the two terms in (18) equals the results in [9] if we choose vr = 0. The
second term is decreasing with 1/R2 and rather small and therefore not further considered. Nevertheless thepresented approach also recovers these non-linear effects. This second term differs a bit from [9] which can beq (cos θ −M)V
explained by a mistake found in [9]. In [9] the second term of 11.2.15 has to be 4π
due to theR2 (1−M cos θ )3
outer derivation, and this coincides with (18).
Having a deeper look at the first term of (18), we see the sound pressure is by a factor (c + vs )2 /(c − vs )2 larger
for a sound source moving towards the receiver compared to a sound source moving away from the receiver.
This amplitude factor is known in literature [9] and reproduced by the presented calculation. An even morecommonly known phenomenon is the Doppler-shift. If we differentiate the phase of the sine in (18) we get=which is the Doppler shift as expected.

dc − vrω0 τt = ω0dt
c − vs(19)

3.2 With binaural receiver
As the human ears are no point-like receivers we have to consider their characteristics. Usually this is done by
an HRIR ĥl,r (θ , φ , τ) depending on the relative angular position, given by azimuth φ and elevation θ , of thesource in the head coordinate system. For a scene with the position ~xr and orientation (e.g, given by a triple ofEuler angles ~αr ) of the receiving head and the position ~xs of the point-like source we can calculate the relative
angles φ and θ . Instead of denoting the HRIR with relative angles we can also denote the relative orientationof the source by absolute positions and orientations with ĥl,r (θ , φ , τ) = ȟl,r (~xr , ~αr ,~xs , τ). As the HRIR describes
the difference between a fictional point-like microphone in the center of the head and the signals at the ears we
can describe the whole transfer function from source to binaural receiver ashl,r (τ) = ĥl,r (θ , φ , τ) ∗ h(~xr ,~xs , τ) = ȟl,r (~xr , ~αr ,~xs , τ) ∗ h(~xr ,~xs , τ) .

(20)

For treating time-varying scenes we associate the position and orientation of the receiver (i.e., ~xr , ~αr ) withthe receiving time t and the position of the source (i.e., ~xs ) with the emission time τ, in analogy to (13).

4836This approach has been presented and validated in [12], although we have to admit that the difference between
velocity potential and sound pressure was not taken into account in [12]. The binaural received velocity potential
is then given byψl,r (t) =ZZh(~xr (t),~xs (τ), τ1 )ȟl,r (~xr (t), ~αr (t),~xs (τ), (t − τ) − τ1 )q(τ) dτ1 dτ

1
δ (τ1 − |~xr (t) −~xs (τ)|/c)ȟl,r (~xr (t), ~αr (t),~xs (τ), (t − τ) − τ1 )q(τ) dτ1 dτ
|~xr (t) −~xs (τ)|/r0
Z1ȟl,r (~xr (t), ~αr (t),~xs (τ), (t − τ) − |~xr (t) −~xs (τ)|/c)q(τ) dτ.
=|~xr (t) −~xs (τ)|/r0

=

ZZ

(21)

By choosing a unity transfer function for the HRIR, i.e., ĥl,r (θ , φ , τ) = δ (τ), it reverts to (13).

4DISCRETE SOUND SCENE RENDERINGThe next step is to transform (21) into a discrete form to use it for digital sound rendering. We denote discrete
sampled signals with square brackets like ψ[k] = ψ(kT ) with discrete-time index k and sampling time T . Thusr0ȟl,r (~xr [k], ~αr [k],~ys [κ], (k − κ)T − |~xr [k] −~ys [κ]|/c)q[κ].
xr [k] −~ys [κ]|κ |~

ψ[k] = ∑(22)

There are now some obstacles to overcome:• The limits for κ range from −∞ to ∞, which means infinite computational effort.
• (k − κ)T − |~xr [k] −~ys [κ]|/c is in general not an integer multiple of T , so that we cannot simply make useof a discrete sampled version of ȟl,r .
To overcome the second issue we have to temporally interpolate the discrete HRIR for a given orientationȟl,r [~xr [k], ~αr [k],~ys [κ], · ], shortly denoted by hl/r [· ] or hl/r (· ), respectively. As an optimal interpolation we use
N−1

hl/r (τ) = ∑ hl/r [κ]si(κ − τ/T ) with si(x) =κ=0sin(πx)πx(23)

where N is the order of the discrete HRIR-filter. As a result we get
N−1



h
i 
r0
|~xr [k] −~ys [κ]|~
ψ[k] = ∑ ∑ȟl,r ~xr [k], αr [k],~ys [κ], m si m − (k − κ) −q[κ]
xr [k] −~ys [κ]|Tc
κ m=0 |~

(24)

At this stage we would have to sum over all possible κ, meaning over the whole input signal. Especially wewould have to do this N times for every output sample. This could be very time consuming and would result in
an infeasible problem for continuous audio rendering. One approach to overcome this is to chose an interpolator
with finite size which we will later denote by I(· ) instead of si(· ).
4.1 More efficient concatenated processing
For more tangible processing we rely on (21) and have a closer look at (14). The impulse response is only
non-zero if the argument of the delta distribution is zero and, thus, we can find a unique mapping from t to τand vise versa as long as the velocities are lower than the speed of sound. In other words, this is a mapping
between emission time τ and receiving time t of a distinct sound event. In (17) this mapping from t to τ forexample is done by the definition of τt . We denote the general mapping from t to τ as τ = o(t) where o stands

4837

10 2

8

Received relative amplitude

Received frequency (kHz)10

Moving Receiver, simualtion
Moving Receiver, expectation
Moving Source, simualtion
Moving Source, expectation

6

4

2

0.2

0.4

0.6

0.8

10 1

10 0

0
0

Moving receiver, simualtion
Moving receiver, expectation
Moving source, simualtion
Moving source, expectation

1

0

0.2

0.4

0.6

M

M

(a) Doppler frequency shift(b) Amplitude shift

0.8

1

Figure 3. Measured amplitude and frequency of received signals, simulation vs. expectation, w = 100.
for originating. Therefore we can write (21) as
ψl,r (t) ===ZZh(~xr (t),~ys (τ), (t − τ) − τ1 )ȟl,r (~xr (t), ~αr (t),~ys (τ), τ1 )q(τ) dτ1 dτh(~xr (t),~ys (τ), (t − τ1 ) − τ)q(τ) dτ ȟl,r (~xr (t), ~αr (t),~ys (o(t − τ1 )), τ1 ) dτ1

Z ZZψm (t − τ1 )ȟl,r (~xr (t), ~αr (t),~ys (o(t − τ1 )), τ1 ) dτ1 .

(25)This describes a subsequent binaural rendering of a mono distance-rendered signal with a retarded source position. Thus, in a first step we render the mono sound and simultaneously save the emission time for each output
sample and do a simple binaural rendering of the mono sound using a source position that belongs to this saved
retarded time in a second step. The binaural rendering part in the discrete case is then given byih
N−1ψl,r [k] = ∑ ȟl,r ~xr [k], ~αr [k],~ys [o[k − `]], ` ψm [k − `] .

(26)

`=0

h
i
For the radial (distance) rendered part we can set ȟl,r ~xr [k], ~αr [k],~ys [κ], m = δ [m] in (24) and get
r0
I
|~
x
[k]
−~
ys [κ]|rκψm [k] = ∑|~xr [k] −~ys [κ]|(k − κ) −
q[κ]
Tc

(27)

where we use a finite interpolator I(· ) with one-sided width w. Thus, the interpolator has to be evaluated onlyfor κ in the range−w ≤ κ − k −|~xr [k] −~ys [κ]|
≤ w.
Tc

(28)

As an additional, external step, we need a discrete integration to get q[k] from s[k] and a discrete differentiation
to get p[k] from ψ[k]. In a first test this is done by first order integration and differentiation asq[k + 1] = q[k] + s[k]/T and p[k + 1] = (ψ[k + 1] − ψ[k])/T.
For continuously operating systems an additional anti-windup in the integration can be useful.

4838

(29)

0

-20

-20

-40

-40

SNR [dB]

SNR [dB]0-60

w = 10, full processing
w = 10, only inner step
w = 100, full processing
w = 100, only inner step
w = 1000, full processing
w = 1000, only inner step

-80

-100

-60

w = 10, full processing
w = 10, only inner step
w = 100, full processing
w = 100, only inner step
w = 1000, full processing
w = 1000, only inner step

-80

-100

-120

-120
0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

M

M

(a) 1st order integration and differentiation(b) higher order integration and differentiation

Figure 4. SNR between simulated and predicted signal for various source velocities and window-width w,
performance of full processing and inner time-retardation step compared.

5

EXPERIMENTAL RESULTS

The evaluation of the whole system consisting of radial and binaural part is quite complicated due to the absenceof a ground-truth reference signal. We implemented the whole processing including the binaural receiver using acontinuous-azimuth representation of the HRIR as described in [13] for motions in a 2D plane. Some examples
were rendered, e.g., a passing ambulance that sounded plausible.
A way to evaluate the proposed method with objective measures is the comparison of physically expectedfeatures with the simulation results. Although it tackles only the part in (11) without the integration and differentiation steps, for the angular part we have done such comparisons in [12]. Thus, below we concentrate on theevaluation of the radial part. A verification of the continuous-time calculation was already done in the sections
before. Therefore the numerical implementation is left for evaluation in this section. For the radial part weknow the results for frequency and amplitude shift in theory. We consider a scene as shown in Fig. 2 with ar
1 kHz sine as source signal. For this scene we know the frequency shift is given by ωDoppler = ω0 c−v
c−vs and thec(c−vr )amplitude shift is given by Adyn = A0 (c−v
2 . The results for various M = v/c for independent source or receivers)motion is shown in Fig. 3. We see that the numerical simulation meets the physical expectation very well.
For the radial part we can furthermore predict the received signal for the given scene by (18). We use amoving source in this case with various M = vs /c and compare the numerically simulated signal with thepredicted one. Therefore we calculate an waveform distance defined by SNR =∑(usimulated −upredicted )2∑ u2predicted

and show

the results in Fig. 4a. The simulation was done for the window-widths w = 10; 100; 1000. We see that thefull processing suffers from higher velocities. Therefore we separated the time retarding part (11) from theintegration and differentiation and plotted the results of inner step (27) additionally in Fig. 4a. This part delivers
constant performance over all velocities and therefore the deteriorating effect for larger velocities is causedby the implementation of integration and differentiation. When we use a more sophisticated integration anddifferentiation Fig. 4b shows an improvement in performance. Here we used q[k] = q[k − 1] + (s[k] + s[k − 1])/2
and an 8th-order central finite difference for differentiation. It should be noted that the SNR is a very rigorous
measure as also small phase-shifts would dramatically affect it, even if they are non-audible.

4839

6

CONCLUSIONSIn this paper we introduced a new approach to handle moving sound sources and receivers within virtual acoustic environments. The presented approach is positioned between phenomenological approaches which try to
mimic physical features without calculating the whole true physical process and, on the other side, comprehensive physical simulations of the whole acoustic field which suffer from very high computational effort. In thispaper we derived a description for moving sound sources and receivers from a deep physical starting point.
Therefore, our approach handles many physical phenomena, such as Doppler shift, amplitude shift, or further
small nonlinearities inherently. We showed the reproduction of these effects in theory and in numerical simulations. This paper is meant to introduce this kind of description for advanced binaural rendering in a generalway. Further improvements in the numerical implementation are possible, e.g., by more sophisticated interpolation algorithms [6] or better integration and differentiation techniques.

ACKNOWLEDGEMENTS
This work is supported by the European Regional Development Fund EFRE-0800372 NRW grant.

REFERENCES[1] D. R. Begault. 3D Sound for Virtual Reality and Multimedia. Academic Press, Boston, Oct. 1994.
[2] J. Blauert, editor. The Technology of Binaural Listening. Springer, 2013.
[3] J. M. Chowning. The Simulation of Moving Sound Sources. JAES, 19(1):2–6, Jan. 1971.
[4] Y. Iwaya and Y. Suzuki. Rendering moving sound with the Doppler effect in sound space. Applied
Acoustics, 68(8):916–922, Aug. 2007.
[5] R. Kronland-Martinet and T. Voinier. Real-Time Perceptual Simulation of Moving Sources: Application to
the Leslie Cabinet and 3D Sound Immersion. J. Audio, Speech, Music Process., 2008(1):849696, June 2008.
[6] T. I. Laakso, V. Valimaki, M. Karjalainen, and U. K. Laine. Splitting the unit delay [FIR/all pass filters
design]. IEEE Signal Processing Magazine, 13(1):30–60, Jan. 1996.
[7] P. L. Lee and J. H. Wang. The simulation of binaural hearing caused by a moving sound source. Computers & Structures, 87(17):1102–1110, Sept. 2009.
[8] F. Meng, G. Behler, and M. Vorländer. A Synthesis Model for a Moving Sound Source Based on Beamforming. Acta acustica united with acustica, 104(2):351–362, 2018.
[9] P. M. Morse and K. U. Ingard. Theoretical Acoustics. Princeton University Press, 1968.
[10] J. O. Smith and D. P. Berners. Doppler Simulation and the Leslie. In COST-G6 Conf. Digital Audio Effects(DAFx-02), pages 13–20, Sept. 2002.
[11] H. Strauss. Implementing Doppler Shifts for Virtual Auditory Environments. In Audio Eng. Soc. Conv. 104,
May 1998.
[12] C. Urbanietz and G. Enzner. Binaural Rendering of Dynamic Head and Sound Source Orientation Using
High-Resolution HRTF and Retarded Time. In IEEE Int. Conf. Acoustics, Speech, Signal Process. (ICASSP),
pages 566–570, Apr. 2018.
[13] C. Urbanietz and G. Enzner. Spatial-Fourier Retrieval of Head-related Impulse Responses from Fast
Continuous-Azimuth Recordings in the Time-Domain. In IEEE Int. Conf. Acoustics, Speech and Signal
Process. (ICASSP), pages 950–954, May 2019.

4840