[Local directional source modeling in wave-based acoustic simulation]
Stefan Bilbao(1) , Jens Ahrens(2) , Brian Hamilton(1)(1) Acoustics and Audio Group, University of Edinburgh, United Kingdom, sbilbao@ed.ac.uk, brian.hamilton@ed.ac.uk(2) Audio Technology Group, Division of Applied Acoustics, Chalmers University of Technology, Gothenburg, Sweden,

jens.ahrens@chalmers.se

Abstract
Time-domain wave-based simulation approaches such as the finite difference time domain (FDTD) method
allow for a complete solution to the problem of virtual acoustics over the entire frequency range, in contrast withthe methods of geometric acoustics which are valid in the limit of high frequencies. They also allow for flexible
modelling of sources and receivers, due to the inherently local nature of the computation, and complete access tothe computed acoustic field over an enclosure. In this paper, a method for the emulation of sources of arbitrary
directivity is presented, framed directly as an inhomogeneous wave equation. The additional terms in the wave
equation take the form of Dirac delta functions and their distributional derivatives, and collections of such terms
may be associated directly with an expansion of source directivity in terms of spherical harmonics. The local
nature of the model implies a locally-defined efficient computational approach for wave-based methods definedover a spatial grid. Numerical results are presented.
Keywords: virtual acoustics, room acoustics, computational acoustics, source modeling, finite differencetime domain

1

INTRODUCTION

Volumetric time-domain wave-based methods for room acoustics and auralisation purposes were proposed in
the mid 1990s [12, 22, 8, 9]. In such a simulation framework, the time evolution of acoustic field is modeled
in its entirety over the room volume, and the simulation operates as a recursive time-stepping method definedover a spatial grid. Many varieties are available, including the finite difference time domain (FDTD) method[30, 25], finite volume methods [4, 6], and pseudospectral methods [14, 17]. In particular, FDTD methods,
defined over regular grids, are equivalent to digital waveguide mesh methods, also proposed for room acoustics
simulation [18]. (The term “wave-based" is also used to refer to non-volumetric simulation methods such as theboundary element method (BEM) [29], which often operate in the frequency domain.) Time-domain wave-based
methods, which compute an approximate wideband solution to the wave equation over an enclosure, subject
to wall conditions and source excitation, can be viewed in contrast with so-called geometrical acoustics (GA)algorithms such as ray tracing [16] or the image source method [2], which represent an approximation to wavebased acoustics in the limit of high frequencies. For audio rate simulation, computational cost is necessarily
high, due to the grid density required (of at most 1 cm resolution), but simulation for reasonable-sized spaces
is now coming within range of commercially-available hardware such as GPUs [27].
Source modeling, in the setting of time-domain wave-based simulation, has been approached by various authors,usually focusing on omnidirectional, dipole or cardioid source directivities [11, 24, 10, 15, 23, 19]. In some approaches, there have been attempts at calibrating an FDTD scheme against measured source directivity [13, 26].
In almost all cases above, however, a continuous spatio-temporal model of the underlying system is unavailable;
sources are implemented in an FDTD scheme through the insertion of driving signals at grid locations. In theinterest of greater generality and flexibility, some models of source directivity do employ such a continuous
model, framed as the 3D wave equation accompanied by a driving term which has the form of a Dirac delta
function under repeated differentiation. A source model of this type corresponding to a multipole expansion has
been presented in [5, 3], and allows for a convenient starting point for discretisation, regardless of the particular
choice of method.

6377

A better choice, which matches well with typical representations of sources in terms of directivity, is a source
model framed directly in terms of spherical harmonics. Such a model is available in a full spatio-temporal form,
in which case source directivity is incorporated into the definition of spatial differential operators applied to a
Dirac delta function which appears as a driving term. Such a model is presented here, accompanied by simple
strategies for discretisation. Numerical results are presented.

2

A SPATIO-TEMPORAL MODEL OF SOURCE DIRECTIVITY

The standard model of wave propagation used in room acoustics is the 3D wave equation:
1 2∂ p − ∆p = 0 .
c2 t

(1)

Here, p(x,t) represents the acoustic pressure in an enclosure, as a function of coordinate r = [x, y, z] ∈ R3 , andt ∈ R. In this paper, concerned with localised models of sources, the problem is assumed defined in free space,
and wall conditions need not be supplied. c is the wave speed, in m/s, ∂t represents partial differentiation withrespect to time t, and ∆ is the 3D Laplacian, defined in terms of the 3D gradient operation∇ = [∂x , ∂y , ∂z ](2)where ∂x , ∂y and ∂z represent partial derivatives with respect to coordinate x, y and z, respectively.
A model including a source with aribitrary directivity has been presented in [3, 5]. Though presented in acompact tensor form in [5], it may be written in expanded form as
∞
∞
∞
1 2
η
∂
p
−
∆p
=
ψηx ,ηy ,ηz (t) ∂xηx ∂y y ∂zηz δ (3) (r)
∑
∑
∑
t
2
c
ηx =0 ηy =0 ηz =0(3)

Here, δ (3) (r) represents a 3D Dirac delta function selecting the source location at r = 0. Each term, indexed
by non-negative integers ηx , ηy , ηz represents the contribution of a single multipole to the corresponding source
directivity, and is accompanied by a distinct driving term ψηx ,ηy ,ηz (t). For example, the term with ηx = 0,ηy = 0, and ηz = 0 corresponds, in isolation, to a monopole [20]; that with ηx = 1, ηy = 0 and ηz = 0 to a
dipole oriented along the x coordinate, and that with ηx = 0, ηy = 0 and ηz = 2 to a longitudinal quadrupoleoriented along the z coordinate.
2.1 Spherical Harmonic Representation
The directivity pattern associated with a given term in (3) is, in general, frequency dependent in a nontrivial
way. Each term corresponds, in general, to a frequency-dependent combination of spherical harmonic directivity
patterns. To simplify the representation in (3), and render it more suitable for matching against known directivity
patterns, a spherical harmonic representation of the source term is of great utility.
To this end, consider the family of real-valued spherical harmonic functions Yl,m (γγ ), indexed by integers l ≥0 and −l ≤ m ≤ l. When written as a function of the components of theγ = r/|r|,p
p unit-length 3-vectorγ ) is a homogeneous polynomial of degree l. For example, Y0,0 = 1/4π, Y1,−1 = 3/4πγy , Y2,0 =
Y
(γ
l,m
p

5/16π 2γz2 − γx2 − γy2 . Define now the spatial differential operator Dl,m in terms of the gradient operator
as defined in (2) as [21, 7]
Dl,m = Yl,m (∇)

or

(ξξ )

Dl,m = ∑ σl,m
ξ ∈Vl

(ξξ )

ξ

∏ ∂ν ν

,Vl = {v ∈ Z3+ |kvk1 = l}

ν=x,y,zfor some coefficients σl,m . The differential operator Dl,m is homogeneous and of degree l. See Table 1.

6378

(4)

Table 1. Dl,m for l = 0, 1, 2, 3 in Cartesian form.
m\ l

0

1

-3

·

·

-2

·

·

-1
0
1

2
·
q

3
∂y
q 4π

1
4π

q

3
4π ∂x

2

·

·

3

·

·

5
2∂z2 − ∂x2 − ∂y2
16π q
15
16π

∂x2 − ∂y2



q

21
32π

q


2∂z3 − 3∂x2 ∂z − 3∂y2 ∂z

21
2
2
3
32π 4∂z ∂x − ∂y ∂x − ∂xq105∂x2 ∂z − 3∂y2 ∂z
16π
q

35
3
2
32π ∂x − 3∂y ∂x4∂z2 ∂y − ∂x2 ∂y − ∂y3



7
q16π



15
4π ∂x ∂zq

35
3∂ 2 ∂y − ∂y3
32π
q x
105
4π ∂x ∂y ∂z

∂x ∂y

15
4π ∂y ∂z

3
∂z
q 4π

·

15

q 4π

q

·
q

3
q



·

An alternative representation for (3) is
∞
l
1 2
∂
p
−
∆p
=
∑ ∑ cl fl,m (t) Dl,m δ (3) (r)c2 t
l=0 m=−lfl,m = al,m ∗ f(5)Here, the driving signals fl,m (t) are assumed related to an underlying scalar driving function f (t) through a filtering operation, where ∗ indicates convolution. The filter responses al,m (t) completely determine the directivityof the source model, as will be shown below.
2.2 Solutions
Consider the basic equation describing a monopole source:
1 2∂ p − ∆p = f (t)δ (3) (r)
c2 t

(6)

The well-known solution p(r,t) and its Fourier transform p̂(r, ω), for angular frequency ω (using the sign
convention of Williams [28]) arep (r,t) =f (t − r/c)
4πrp̂ (r, ω) =fˆeiωr/ciω ˆ (1)=f h (ωr/c)
4πr
4πc 0

(7)

(1)

where r = |r|, and where fˆ (ω) is the Fourier tranform of f (t). h0 is the zeroth order spherical Hankel functionof the first kind.
Through superposition, the solution to (5) is then∞lfl,m (t − r/c)iω ˆ ∞ l l(1)lp (r,t) = ∑ ∑ c fl,m (t) Dl,mp̂ (r, ω) =f ∑ ∑ c âl,m (ω) Dl,m h0 (ωr/c) (8)4πr4πcl=0 m=−ll=0 m=−lwhere fˆl,m (ω) = fˆ (ω) âl,m (ω). Using the identity (not proven here):
(1)

(1)

Dl,m h0 (ωr/c) = (−ω/c)l Yl,m (γγ ) hl (ωr/c)

(9)

(1)

where hl is the lth order spherical Hankel function of the first kind leads to the following frequency domain
expression for the solution p̂ (r, ω):p̂ =iω ˆ ∞ l(1)f ∑ ∑ (−ω)l âl,m (ω)Yl,m (γγ ) hl (ωr/c) .
4πc l=0
m=−l

(10)

6379

Dividing by the response fˆ of the source signal leaves an expression for the directivity which recovers thestandard form for the exterior field of a source—see, e.g., [1]. The model (5), however, is expressed entirely inthe spatio-temporal domain, and is a good candidate for discretisation through volumetric wave-based methodsof any type—note that, as yet, no discretisation has been employed.
2.3 Displaced Monopole
An interesting test case is that of the monopole source located at coordinates r = r0 , where r0 = r0 γ 0 , for adisplacement distance r0 and unit-length direction vector γ 0 :1 2∂ p − ∆p = f (t)δ (3) (r − r0 )
(11)
c2 t
In this case, an analytic time domain expression for the filters al,m (t) (necessarily non-causal) in model (5) isavailable:
l


t l−1
t2
1, |ξ | ≤ 1
whereu (ξ ) =t0 = r0 /cal,m (t) = 4π (−1)l Yl,m (γγ 0 ) 0
1 − 2 u (t/t0 )0, |ξ | > 1
2 · (2l)!!
t0
(12)

3

FDTD METHODS

Only the most basic FDTD scheme, sometimes referred to as the seven-point scheme, will be described here.
But the method described below can be applied to any volumetric time-domain method.
Consider a grid function pnq , indexed by integer q ∈ Z3 and n ∈ Z. pnq represents an approximation to p (r,t) at
t = nT and x = qX. Here, T is the time step, in s (and 1/T is the sample rate), and X is the grid spacing, inm. The simplest possible FDTD scheme for the general system (5) is
∞
l

2
n
n
2
n
pqn+1 = 2pnq − pn−1
+
λ
p
−
p
+
T
c2+l fl,mdl,m sq
∑
∑
∑
q
q+e
q
e∈Q

,

Q = {e ∈ Z3 |kek1 = 1},

λ = cT /X

l=0 m=−l(13)√
The constant λ is referred to as the Courant number for the scheme, and is bounded by λ ≤ 1/ 3 for numericalstability. sq is an approximation to the 3D Dirac delta function, and defined as
1/X 3 , q = 0
sq =
.
(14)
0,
otherwisen result from discrete time convolution between the driving function f n and the individual
The time series fl,mspherical harmonic channel filters anl,m , both assumed known a priori. The discrete operators dl,m are approximations to the partial differential operators Dl,m , and their construction is described below.

3.1 Discrete Spherical Harmonic Operators
Approximations to the operators Dl,m , as defined in (4) may be obtained by various means—see, e.g., [5].
A basic approach is through standard difference operations. Consider a grid function gq , defined for q ∈ Z3 .
Difference operators dν+ and dν− , approximating ∂ν , and an averaging operator µν− may be defined as1(gq+eν − gq )X
Here eν is a unit vector in direction ν.
dν+ gq =

ξ

dν− =

1
(gq − gq−eν )
X

µν− =

1
(gq + gq−eν ) .
2(15)

Consider now the approximation of ∂ν , a ξ th spatial derivative in coordinate ν; ξ may be decomposed uniquelyξas ξ = 2Mξ + αξ for integer Mξ ≥ 0 and αξ ∈ {0, 1}. A centered approximation dν follows as
α
M +α
M
ξ
ξ
dν , µν− ξ dν+ ξ ξ dν− ξ u ∂ν ,
(16)

6380

Figure 1. Stencils of the discrete operators dl,m , for l = 0, 1, 2..
Using the expanded representation of Dl,m from (4), one may then arrive at a centered approximation dl,m toDl,m :(ξξ )dl,m = ∑ σl,m
ξ ∈Vl

ξ

∏ dνν u Dl,m .

(17)

ν=x,y,zSee Figure 1 for a graphical representation of the operators dl,m .

4

SIMULATION RESULTS

In this section, some basic results are shown, using scheme (13), operating at 44.1 kHz, over a cube of side
length 2 m.
As a first step, it is useful to regenerate spherical harmonic directivity patterns by exciting only one of thechannels (i.e., setting all filters al,m to zero except for one in the scheme (13)). The scheme is excited by
setting f n to a unit impulse, and then drawing output from a set of receivers located on a sphere of radius 0.3
m around the source location. High-order separable interpolating approximants are used [5] for the receivers(60 azimuthal by 30 in inclination). The resulting outputs are Fourier-transformed, and directivity magnitude is
as shown in Figure 2 at 2 kHz, for all spherical harmonics up to l = 3.

Figure 2. Spherical harmonic directivity patterns generated using scheme (13) at 44.1 kHz.
As a second example, consider the case of the displaced monopole source, as described in Section 2.3. Here, theexcitation f n is chosen as a Gaussian function, of variance σ = 2 × 10−4 s, and the displaced source location ischosen as r0 = [−0.2, 0, 0]. The filters anl,m are sampled from the exact expressions given in (12). Snapshots of6381

the time evolution of the acoustic field are as shown in Figure 3, for different choices of the maximal sphericalharmonic order L.

L=3

L=6

L=9

Figure 3. Snapshots of the time evolution of the acoustic field in the case of a displaced monopole. The
driving term for the FDTD scheme is located within the white circle, and the location of the displaced source
is indicated by a × symbol. Results are shown for a maximal spherical harmonic order of l ≤ L, for L = 3 (toprow), L = 6 (middle row) and L = 9 (bottom row).

5

CONCLUDING REMARKS

A basic model of source directivity has been presented here, intended for use within a volumetric wave-based
simulation framework. As such, the model is expressed entirely in the spatio-temporal domain, as an extensionof the 3D wave equation. Such a model lends itself easily to a time stepping method defined over a regular
Cartesian grid, as presented here, but the model in (5) does not require this. Operation over unstructured grids isalso possible, requiring some approximation to the derivatives of the Dirac distribution. As with all wave-based
methods, such a model naturally reproduces near-field effects.
Regardless of the particular discretisation method, however, an important feature is that the approximation to
the source is local, and thus not computationally costly. Indeed, for the simulation scenario presented here, thecomputational cost associated with the incorporation of a source term in the scheme (13) is far less than thecost of the operation of the scheme over the problem interior. Such a localised source model is also flexible,
in the sense that it may be placed independently of any consideration of room boundary conditions, and also
that, once the appropriate filters anl,m and spatial difference operators have been determined, the source location

6382

may be displaced within the simulation volume without any additional computation. In addition, fitting of thefilters al,m against measured directivity patterns may be carried out independently of the particular choice oftime stepping scheme.

REFERENCES[1] J. Ahrens. Analytic Methods of Sound Field Synthesis. Springer, Heidelberg, Germany, 2012.
[2] J. Allen and D. Berkley. Image method for efficiently simulating small-room acoustics. J. Acoust. Soc. Am.,
66(4):943–950, 1979.
[3] M. Bencomo and W. Symes. Discretization of multipole sources in a finite difference setting for wave
propagation problems. J. Comp. Phys., 386:296 – 322, 2019.
[4] S. Bilbao. Modeling of complex geometries and boundary conditions in finite difference/finite volume time
domain room acoustics simulation. IEEE Trans. Audio Speech Language Proces., 21(7):1524–1533, 2013.
[5] S. Bilbao and B. Hamilton. Directional sources in wave-based acoustic simulation. IEEE Trans. Audio
Speech Language Proces., 27:415–428, 2019.
[6] S. Bilbao, B. Hamilton, J. Botts, and L. Savioja. Finite volume time domain room acoustics simulation under general impedance boundary conditions. IEEE/ACM Trans. Audio Speech Language Proces., 24(1):161–
173, 2016.
[7] S. Bilbao, A. Politis, and B. Hamilton. Local time-domain spherical harmonic spatial encoding for wavebased acoustic simulation. IEEE Signal Proces. Lett., 26(4):617–621, 2019.
[8] D. Botteldooren. Acoustical finite-difference time-domain simulation in a quasi-cartesian grid. J. Acoust.
Soc. Am., 95(5):2313–2319, 1994.
[9] D. Botteldooren. Finite-difference time-domain simulation of low-frequency room acoustic problems. J.
Acoust. Soc. Am., 98(6):3302–3308, 1995.
[10] J. Botts, A. Bockman, and N. Xiang. On the selection and implementation of sources for finite-difference
methods. In Proc. 20th Int. Congr. Acoust., Sydney, Australia, 2010.
[11] A. Celestinos and S. Nielsen. Low-frequency loudspeaker room simulation using finite differences in thetime domain part 1: Analysis. J. Audio Eng. Soc., 56(10):772–786, 2008.
[12] O. Chiba, T. Kashiwa, H. Shimoda, S. Kagami, and I. Fukai. Analysis of sound fields in three dimensional
space by the time-dependent finite-difference method based on the leap frog algorithm. J. Acoust. Soc.
Jpn.(J), 49:551–562, 1993.
[13] J. Escolano, J. Lopez, and B. Pueo. Directive sources in acoustic discrete-time domain simulations basedon directivity diagrams. JASA Express Lett., 121:256–262, 2007.
[14] M. Hornikx, W. De Roeck, and W. Desmet. A multi-domain Fourier pseudospectral time-domain method
for the linearized Euler equations. J. Comp. Phys., 231(14):4759–4774, 2012.
[15] H. Jeong and L. Y. Source implementation to eliminate low-frequency artifacts in finite difference time
domain room acoustic simulation. J. Acoust. Soc. Am., 258–268(1):1112–1118, 2012.
[16] A. Krokstad, S. Strom, and S. Sorsdal. Calculating the acoustical room response by the use of a ray
tracing technique. J. Sound Vib., 8(1):118–125, 1968.

6383

[17] R. Mehra, N. Raghuvanshi, L. Savioja, M. Lin, and D. Manocha. An efficient GPU-based time domain
solver for the acoustic wave equation. Appl. Acoust., 73(2):83–94, 2012.
[18] D. Murphy, A. Kelloniemi, J. Mullen, and S. Shelley. Acoustic modelling using the digital waveguidemesh. IEEE Sig. Proc. Mag., 24(2):55–66, 2007.
[19] D. Murphy, A. Southern, and L. Savioja. Source excitation strategies for obtaining impulse responses infinite difference time domain room acoustics simulation. Appl. Acoust., 82:6–14, 2014.
[20] A. Pierce. Acoustics: An introduction to its physical principles and applications. Acoustical Society ofAmerica, 1991.
[21] E. Rowe. Spherical delta functions and multipole expansions. J. Math. Phys., 19:1962–1968, 1978.
[22] L. Savioja, T. Rinne, and T. Takala. Simulation of room acoustics with a 3-D finite-difference mesh. In
Proc. Int. Comp. Music Conf., pages 463–466, Århus, Denmark, sep 1994.
[23] J. Sheaffer, M. van Walstijn, and B. Fazenda. Physical and numerical constraints in source modeling for
finite difference simulation of room acoustics. J. Acoust. Soc. Am., 135(1), 2014.
[24] A. Southern and D. Murphy. Low complexity directional sound sources for finite difference time domain
room acoustic models. In Proc. 126th Audio Eng. Soc. Conv., Munich, Germany, 2009.
[25] J. Strikwerda. Finite Difference Schemes and Partial Differential Equations. Wadsworth and Brooks/Cole
Advanced Books and Software, Pacific Grove, California, 1989.
[26] D. Takeuchi, K. yatabe, and Y. Oikawa. Source directivity approximation for finite-difference time-domain
simulation by estimating initial value. J. Acoust. Soc. Am., 145(4):2638–2649, April 2019.
[27] C. J. Webb and S. Bilbao. Computing room acoustics with CUDA - 3D FDTD schemes with boundary
losses and viscosity. In Proc. IEEE Int. Conf. Acoust., Speech Signal Proces., pages 317–320, 2011.
[28] E. Williams. Fourier Acoustics: Sound Radiation and Nearfield Acoustical Holography. Academic Press,
New York, 1999.
[29] Y. Yasuda and T. Sakuma. Boundary element method. In T. Sakuma, S. Sakamoto, and T. Otsuru, editors,
Computational Simulation in Architectural and Environmental Acoustics, page 80. Springer, 2014.
[30] K. Yee. Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic
media. IEEE Trans. Antennas Propagation, 14:302–307, 1966.

6384