Combining image and equivalent sources for room acoustic simulations
Boris MONDET(1) , U. Peter SVENSSON(2) , Jonas BRUNSKOG(3) , Cheol-Ho JEONG(3) , Claus Lynge CHRISTENSEN(4) ,Jens Holger RINDEL(4)(1) Odeon A/S, Denmark, bojmo@elektro.dtu.dk(2) Dept. of Electronic Systems, Norges Teknisk-Naturvitenskaplige Universitet, Norway(3) Dept. of Electrical Engineering, Danmarks Tekniske Universitet, Denmark(4) Odeon A/S, Denmark

Abstract
Room acoustic simulation tools based on geometrical acoustics fail to predict accurately the behaviour of rooms
at low frequencies. This inaccuracy can lead to poor acoustic conditions especially in rooms of moderate volume,
such as classrooms and meeting rooms. The present paper tackles this problem by introducing a simulation
technique for room acoustics based on the Image Source Method and the Equivalent Source Method. Thesemethods are implemented in the time domain to achieve efficient calculations with reliable results for low to
medium frequencies. Given a source emitting a pulse signal in a room, image sources up to different orders aredefined to determine early reflections. With the direct sound field from the original source and the contributionsof the image sources, the remaining part of the sound field is solved by the Equivalent Source Method to accountfor late reflections and diffraction effects. The results of the test simulations carried out show that this combined
Image-Equivalent Source Method is a suitable option for room acoustics.
Keywords: Room Acoustic Simulation, Equivalent Source Method, Image Source Method

1

INTRODUCTION

The current state of room acoustic simulation is well-suited for mid to high frequencies or large volumes, suchas concert halls and auditoria, thanks to the principles of geometrical acoustics. Geometrical acoustics methods
are indeed reliable in these conditions and computationally efficient, but they fail at producing accurate results
at low frequencies. Wave-based methods are the most common solution in acoustics for simulating low frequencies, with the Finite Element Method and the Boundary Element Method being the most popular. However,they have the disadvantage of being computationally heavy and thus unpractical for room acoustics. The Equivalent Source Method (ESM), however, is a compromise solution. Similarly to the Boundary Element Method, it
simulates the sound field by placing sources around the boundaries of the domain, but the main difference liesin the approximations used to reduce the calculation load. It is therefore an interesting candidate for simulating room acoustics at low frequencies. Furthermore, it is possible to easily combine the ESM with the Image
Source Method (ISM) to enhance its accuracy at a low computational cost. Calculations performed in the time
domain also allow to limit computation time since one simulation is enough to cover the whole frequency range
considered.
The Equivalent Source Method is typically used to calculate sound radiation or scattering from objects. Kropp
and Svensson [1] laid the foundation of the ESM in the time domain for acoustic problems, and Lee, Brentnerand Morris [2] furthered its study. Recently, Lee [3] gathered the state-of-the-art knowledge on the topic andestablished guidelines to solve radiation and scattering problems. The present study is the first step towards
a different approach of solving interior problems with the ESM, where the sound field inside a room will besimulated rather than around an object.

5997

2

THEORYThe sound reflection on a single plane with rigid boundaries will be investigated in this paper. In all thesimulations considered, sound propagates freely from sources to receivers in a three-dimensional domain. This
propagation is described by the free field Green’s function in 3D:
G(r,t) =1 δ (t − r/c)
,
4π
r

(1)

with r the distance between a source and a receiver, c the speed of sound, and δ the Dirac function. Given asource with a prescribed volume velocity q(t), the sound pressure at a receiver is obtained after multiplication
with the Green’s function and integration over time:p(r,t) =1 q(t − r/c)
.
4π
r

(2)

When several sources emit at the same time, their individual sound pressures are summed up to obtain the totalsound pressure at the receiver.
The Equivalent Source Method allows to calculate the sound field inside or around an object by solving theboundary condition equation. To do so, the plane considered in this study is discretized with control points.
With the sound source on one side of the plane, a virtual surface with same dimensions is created on theother side and discretized with equivalent sources. The number of equivalent sources can be lower than thatof surface points to reduce computation time, leading to an overdetermined problem. The volume velocities ofthese sources are then determined so that the sum of the incident sound field from the original source and thereflected sound field from the equivalent sources satisfies the boundary conditions on the control points. In thispaper only rigid boundaries are considered. In terms of sound pressure, it is equivalent to
→
−
→
−
−
−
∇ pr · →n = − ∇ pi · →
n,

(3)

−
with pi the incident sound pressure, pr the reflected sound pressure, and →n the unit normal vector pointing
towards the original sound source. If the volume velocity time signal of the original source is known as ananalytical and derivable formula, the gradient of the incident pressure can be calculated directly. The reflectedsound pressure, however, is related to the volume velocities of the equivalent sources that need to be determined.
By noting qe the volume velocity of the equivalent source e and using equation (2), the reflected pressure
gradient is expressed as→−
1 ∂ qe1qe (t − re /c)∇ pr =−−(t − re /c) −
u→
(4)
re .
4π ∑
re2
cre ∂te
In this equation, r represents the distance and −u→ the unit directional vector from the equivalent source e to theeresurface point under consideration. Numerical methods, including the ESM, work in a discrete time domain. In
order to represent t −re /c in discrete time, the volume velocities qe are discretized assuming linear interpolation;
hence, their derivatives are approximated with first order. With τe = t −re /c, the neighbouring discrete time steps
are noted te− and te+ such that te− ≤ τe ≤ te+ . Applying the time discretization to equation (4) with a time step∆t, the boundary condition in equation (3) becomes


 +

→
−
1
τe − te−11
te − τe
1
→
→
−
→
−
→
−
+ −
−
−
qe (te )ure · n = − ∇ pi (τe ) · n −−+qe (te− )−u→(5)
re · n .
4π∆t ∑
re2cre
4π∆t ∑
re2creee
It appears that the volume velocities of the equivalent sources at a new time step can be determined from previous time steps, resulting in a marching algorithm. The system made of the equations for all surface points
can be written in matrix form and solved at once with the least squares method for each time step.

5998

-10

Volume velocity (dB re 1 m3/s)

Volume velocity (m3/s)

1

0.8

0.6

0.4

0.2

0

0

0.002

0.004

0.006

0.008

-15

-20

-25

-30

-35
101

0.01

Time (s)102103Frequency (Hz)Figure 1. Source signal in the time domain and its frequency content. Solid blue line: fs = 16000 Hz; Dashed
red line: fs = 1600 Hz.
The Image Source Method also relies on solving the boundary condition equation. Given a source in frontof a single plane, an image source is located in the symmetric position relative to the plane. The signal of
the image source is then adjusted to fulfill the boundary conditions. In the case of rigid boundaries, the image
source simply emits the same signal as the original source.
In order to combine image and equivalent sources in one method, the image source for the single plane isfirst implemented. Its contribution is added to that of the original source to form the incident sound pressure,which is then incorporated in the ESM algorithm. This will later on be called the Image-Equivalent Source
Method (IESM).
The Edge Diffraction Method [4] is based on the theory from Biot and Tolstoy [5] about reflection aroundan infinite rigid wedge. The sound field is decomposed between incident wave, specular reflection, and diffraction waves with an opposite polarity emanating from the edges of the scattering object.

3

NUMERICAL SETUP

The simulation setup consists of a single plane with rigid boundaries that extends from −2.5 m to 4.5 m in thex-direction and from −3 m to 4 m in the z-direction. A point source is located at the coordinates (−1 m, 1.5 m,1 m) and emits a Gaussian pulse which volume velocity isqi (t) = e−(t−t0 )2
2s2

.

(6)

t0 represents the center time of the pulse and s is the standard deviation controlling the width of the pulse in thetime domain. These parameters are set to t0 = 5 ms and s = 0.5 ms. This pulse was chosen in order to limit thefrequency bandwidth of the signal used in the simulations. The corresponding source signal and its frequency
content are shown in Figure 1 with two different sampling frequencies fs . By taking 3 dB below the maximum
amplitude, the cut-off frequency of the pulse is found to be 400 Hz. The sampling frequency is thereafter set to1600 Hz, limiting the simulations to frequencies below 800 Hz.
Firstly, the growth rate of computation time in the ESM relative to the size of the problem is investigated.
Different simulations are performed with varying spacing between surface points. The spacing between equivalent sources is set 4 times larger than between surface points, meaning that each source is centered between

5999

Figure 2. Close up view of the plane discretization. Blue dots: surface points; Red crosses: equivalent sources.
16 points, as seen in Figure 2. The distance between the surface plane and the virtual plane comprising theequivalent sources is 11.4 mm. The single plane reflection being a very simple case, the residual of equation(5) is nearly constant in all simulations. Thus, a convergence study is not performed here.
Secondly, the sound pressure at three receivers with coordinates (2 m, 2 m, 1 m), (−3.5 m, 1.5 m, 1 m) and(−5 m, 1 m, 1 m) is investigated. This placement is chosen so that the reflection points corresponding to thosereceivers are respectively located on the plane, on the edge of the plane, and outside the plane. In the ESM
the plane surface is discretized with control points spaced by 71.5 mm, corresponding to one sixth of the wavelength for the highest frequency. In the ISM only one image source is required for the single plane case, but
this method has the drawback of assuming an infinitely large plane where diffraction cannot be represented.
Concerning the IESM simulation, the same parameters are applied as with the ESM and ISM taken separately.
The Edge Diffraction Method calculates the response to a unit Dirac pulse. This impulse response is obtainedin the time domain with only the first diffraction order, and it is then convoluted with the Gaussian pulse inequation (6).

4

RESULTS AND DISCUSSION

Although the single plane reflection case with rigid boundaries does not allow to study convergence, the computation time of the Equivalent Source Method can be investigated. Of particular interest is the way computation
time increases with the size of the problem to be solved. Here, the size of the problem is varied by changing
the spacing between control points on the surface, hence modifying the density of surface points h. The ratio
between surface points and equivalent sources being fixed, the spacing of equivalent sources also varies accordingly. The result of this investigation is presented in Figure 3. The computation time of the ESM is seen to be
proportional to the square of the surface point density, therefore it grows with second order. This means thatrefining the discretization of the surface would be attractive if the residual in equation (5) converges with order
2 or higher. As previously mentioned, it is unfortunately not possible to verify this convergence rate with thepresent test case.

6000

Computation time (s)

102

101

100

10-1
101

102

Surface point density (m-2)

Figure 3. Evolution of the computation time of the ESM according to the surface point density h. Solid blue
line: computation time of the simulations; Dashed black line: O(h2 ).

The sound pressure at the receivers is now considered. For the first receiver, the reflection point is locatedon the plane. Consequently, the specular reflection should be the same as for an infinite plane, which is thecase simulated with the ISM. This is in agreement with the results in Figure 4. Indeed, the direct sound andspecular reflection is the same for all the methods used. However, discrepancies arise for the modelling of
the diffraction that reaches the receiver between 0.03 s and 0.04 s. The ESM and the Edge Diffraction Method
show a fair agreement in the sound pressure, even though the ESM seems to underestimate the amplitude of thediffraction. Besides, it appears clearly here that the waves emanating from the edges are in antiphase with the
specular reflection. With the Image-Equivalent Source Method, the result obtained is identical to the ISM. This
has a simple explanation: the sound field emitted by the image source satisfies the boundary condition in equation (3), therefore the equivalent sources remain inactive in the simulation. In order to account for diffraction
effects in this case, the single plane should be modelled with a finite thickness. The equivalent sources would
then be activated to fulfill the boundary conditions on all faces of the plate. The diffraction effect being absentof IESM simulations in this single plane case, this method will not be considered for the other two receivers.
The next receiver position was chosen so that the reflection point is located on the edge of the single plane. It
implies that the specular reflection and the diffraction wave coming from that edge reach the receiver simultaneously. As seen previously, the specular and diffraction reflections are in antiphase, thus desctructive interference
occurs. Given the size of the plane, this applies to the whole frequency range and the amplitude of the resulting
reflection is expected to be half of what would be observed for an infinite plane [6]. This is what is observedin Figure 5 where the specular reflection in the Edge Diffraction Method is half as strong as in the ISM. However, the reflection of the ESM is not as attenuated. This is a confirmation that the ESM underestimates theamplitude of the diffraction. The sound pressure at the third receiver yields the same observations. The locationof this receiver makes the reflection point to be situated outside the plane. Naturally, it is seen that the specularreflection is even more attenuated in this case, although not as much for the ESM as for the Edge Diffraction
Method.
In these examples it was seen that the ESM is not as accurate as the Edge Diffraction Method. However,

6001

0.03

0.025

0.025

Sound pressure (Pa)Sound pressure (Pa)

0.03

0.02
0.015
0.01
0.005
0
-0.005

0.02
0.015
0.01
0.005
0

0

0.01

0.02

0.03

0.04

-0.005

0.05

0

0.01

Time (s)

0.02

0.03

0.04

0.05

Time (s)

0.03

0.03

0.025

0.025

Sound pressure (Pa)Sound pressure (Pa)

Figure 4. Sound pressure at the receiver with reflection point on the plane. Solid black line: Image Source
Method; Dotted blue line: Edge Diffraction Method; Dashed red line: Equivalent Source Method (left) or
Image-Equivalent Source Method (right).

0.02
0.015
0.01
0.005
0
-0.005

0.02
0.015
0.01
0.005
0

0

0.01

0.02

0.03

Time (s)

0.04

0.05

-0.005

0

0.01

0.02

0.03

0.04

0.05

Time (s)Figure 5. Sound pressure at two different receivers with reflection point on the edge of the plane (left) or
outside the plane (right). Solid black line: Image Source Method; Dotted blue line: Edge diffraction Method;
Dashed red line: Equivalent Source Method.
the accuracy is expected to be improved in the IESM. With the help of an image source accounting for the
specular reflection, the equivalent sources will be devoted to model diffraction. Applying the IESM to a closed
room will also allow to avoid the problem encountered here where equivalent sources remain inactive. Another topic that was not covered in this paper and will require attention for room acoustic simulation is thedescription of boundary conditions. In realistic room cases, the boundaries have a finite impedance that isfrequency-dependent. This changes the boundary condition equation in the time domain with additional terms
depending on past events.

5

CONCLUSION

This paper investigated the performance of the Equivalent Source Method and the Image-Equivalent Source
Method for a single plane reflection with rigid boundaries. The ESM was seen to be accurate to model the
specular reflection but slightly underestimated the amplitude of diffraction effects. The computation time was6002

also found to grow with second order, although a convergence study could not be carried out to complement
this finding. The IESM, however, could not yield correct results due to the boundary conditions being fulfilledby the image source alone; hence, the diffraction could not be modelled. This problem can later be encounteredin room acoustic simulations with hanging reflectors for example. In order to activate equivalent sources andcorrectly model diffraction, such reflectors would need to be given a finite thickness. It is expected that the
IESM will be more accurate than the ESM with similar convergence rate and growth rate of computation time.

REFERENCES[1] Kropp, W.; Svensson U. P. Application of the time domain formulation of the method of equivalent sources
to radiation and scattering problems. Acustica, Vol 81, 1995, pp 528-543.
[2] Lee, S.; Brentner K. S.; Morris P. J. Acoustic scattering in the time domain using an equivalent source
method. American Institute of Aeronautics and Astronautics Journal, Vol 48 (12), 2010, pp 2772-2780.
[3] Lee, S. Review: The use of equivalent source method in computational acoustics. Journal of Computational
Acoustics, Vol 25 (1), 2017, pp 1630001 1-19.
[4] Torres R.; Svensson U. P.; Kleiner M. Computation of edge diffraction for more accurate room acoustics
auralization. Journal of the Acoustical Society of America, Vol 109 (2), 2001, pp 600–610.
[5] Biot M. A.; Tolstoy I. Formulation of wave propagation in infinite media by normal coordinates with anapplication to diffraction. Journal of the Acoustical Society of America, Vol 29 (3), 1957, pp 381–391.
[6] Rindel J. H. Attenuation of sound reflections due to diffraction. Proceedings of Nordic Acoustical Meeting,
Aalborg, Denmark, August 20-22, 1986.

6003