Study of DIF boundary model with rectilinear FDTD scheme in voice booths

Vito ROMANELLI TRICANICO(1) , Marcel BORIN(2) , Carolina MONTEIRO(3) , Rânnely ARAÚJO(4)(1) Harmonia Acústica, Brazil, vito.romanelli@harmoniaacustica.com.br(2) Harmonia Acústica, Brazil, marcel.borin@harmoniaacustica.com.br(3) Harmonia Acústica, Brazil, carolina.monteiro@harmoniaacustica.com.br(4) Harmonia Acústica, Brazil, ranny.araujo@harmoniaacustica.com.br

Abstract
Computational simulations are helpful when studying acoustical properties of rooms. This represents a challenge when very small rooms, such as office video/phone booths, are taken into account. The volume of theserooms generally does not exceed 15 m3 and the first modes occur at higher frequencies if compared to larger
rooms. Computational room acoustics considers two main techniques of modelling sound propagation: geometrical and wave-based. As geometrical modelling is not the most suitable approach for low frequencies since thephase performs an important part in the physics at this frequency range, in the present study, a rectilinear FDTD
scheme routine was implemented to obtain the impulse responses for low frequency bands of a voice booth,using proper boundary conditions derived from digital impedance filters (DIF). In order to obtain broadband
impulse responses, high frequency impulse responses were acquired using a commercial software based on geometrical methods. Both responses were later filtered with low-pass and high-pass filters based on Schroeder
frequency, and then assembled. Finally, the results extracted from the impulse responses of simulations werecompared with data obtained from field measurements.
Keywords: Sound, Computational, Boundary

1

INTRODUCTION

The need to perform computational simulations to obtain acoustic room parameters is evident even when there
are pre and post in situ measurements in the room to be treated. They can predict room configurations without
the need of moving to the field. In Brazil, rooms such as voice booths have been increasingly used in offices inorder for the user to obtain a higher level of privacy, to make small meetings or calls. These rooms generally
do not exceed 15 m3 in volume, and Schroeder’s frequency occurs at higher frequencies when compared tolarger rooms, implying that booth has modal behavior in a broader band of frequencies. The modal frequencies
end up reaching the speech bands, and there is still no consensus on how to measure room parameters in suchsmall volume. In [1], a method is presented to measure the frequency response of this type of room, obtaining
alike curves between several voice booths with different configurations and dimensions, performed by Harmonia
Acústica 1 . Because of this behavior, this type of room is ideal for studying wave-based methods, such as
the finite-difference time-domain method (FDTD) [2], finite volumes and finite elements. Ray-based methods
are not the most suitable for this type of room because they suppress phenomena such as diffraction at low
frequencies, not adapting the room’s response, which can lead to large inaccuracies in the values generated by
the simulation. The FDTD has its advantages by being easily implementable, impulse responses obtained are
already in time-domain, and it has a DIY aproach, whereas other methods need a few other software available
for creating meshes, and ABCs based on parameters other than the coefficient of absorption, such as impedance,
often difficult to achieve for specific materials.
In 2008, Kowalczyk [3] presented boundary conditions for the FDTD method based on digital impedance filters
1 www.harmoniaacustica.com.br

6004

(DIF), which is obtained from the absorption coefficient of the boundary, taking into account important physicalconsiderations, but presenting only numerical results. In 2017, Reuben [4] showed that the model has a good
application for computational simulations, but also no tests have been done to attest the use of the program incomparison to in situ measurements. An advantage of this method however is that there is no need to do thesimulation for each frequency band since the modelled boundary is independent of the frequency, decreasing
the execution time intensively.
Although the model has shown convergent solutions for some simulated room configurations, it is desired to
study the applicability of this boundary model, comparing results obtained from simulations and measurements.
The intention is to find a point of similarity of parameters obtained with the simulation and with the measurements, presenting a mutual validation.

2

OBJECTIVES

The goals of this study are to compare measurements and simulation results of 1/3 octave band frequency
responses and T30 of the voice booth, in order to assess the applicability of the DIF boundary model in an
rectilinear FDTD grid for very small rooms.

3

METHODOLOGY

The process for obtaining the impulse responses presented a certain level of complexity amidst the simplicities
presented a priori. The method of measurements is also presented in this section.
3.1 The wave equation and FDTD method
After doing a discretization of the undamped homogeneous wave-equation [2, 5], and some algebra manipulation, one may get that the update equation for the sound pressure p(x, y, z,t) in a rectinilear 3D FDTD grid
is
1 n
n−1
pi,n+1j,k = 3 Si, j,k − pi, j,k ,(1)

Si,n j,k = pni+1, j,k + pni−1, j,k + pni, j+1,k + pni, j−1,k + pni, j,k+1 + pni, j,k−1 ,

(2)

where

with i, j and k denoting spatial indices for x, y and z directions and n denoting the time-step index. The Courant
number, a constraint given by λ = c XT , where T is the time-step in seconds, X is grid spacing in meters andc is sound speed in air, was chosen to be λ = √13 to get to Equation (1) in its form. This value satisfies thecondition to minimize the dispersion error with the rectilinear grid topology [3].

6005Figure 1. Rectinilear topology for the grid

Figure 2. Ghost node for x+ boundary

Different topologies would require Equation (1) to be more complicated. The rectilinear topology for the 3D
FDTD grid can be described by interconnected nodes, presented in Figure 1.
3.2 Boundary model
The boundaries must be considered as locally reacting surfaces (LRS), i.e., the normal component of the particle
velocity at the surface of the boundary depends on the sound pressure in front of the boundary element and not
on pressure in front of neighbouring elements [3]. This holds for boundaries that do not propagate vibrationsin the direction parallel to the boundary surface. For a boundary in the right, i.e., in the x+ direction, it isnecessary that the ghost node - the node that is out of simulation bounds - in position i + 1 is eliminated inequation (1), as illustrated in Figure 2. The nodes of sound pressure in positions i, j, k − 1 and i, j, k + 1 aresuppressed in the figure. For a pressure-centered approximation, Kowalczyk [3] showed that the equation for
the right boundary is
√
3 n−1
n
(p
− pn+1i, j,k ) + pi−1, j,k ,ξω i, j,k

pni+1, j,k =(3)where ξω is the specific acoustic impedance of the boundary. This approach may cause solution growth duringexecution depending on the form of the modelled room [6]. To reverse this problem, the velocity-centeredapproximation is used:
pni+1, j,k =

√
3 n−1
n
(p
− pn+1i, j,k ) + pi, j,k .
2ξω i, j,k(4)

For normal incidence in the boundary (which is the case for locally reacting surfaces), ξω may be written asξω (z) =1 + R0 (z),1 − R0 (z)(5)√
where R0 (z) is the boundary reflection coefficient at normal incidence [3, 6], and |R0 (z)| = 1 − α, with α
being the sound absorption coefficient of the boundary. Therefore, ξω is expressed as an arbitrary magnitude
IIR digital filter with transfer function given byξω (z) =6006

b0 + b1 z−1 + b2 z−2 + ... + bN z−N,
a0 + a1 z−1 + a2 z−2 + ... + aN z−N(6)where each z−m represents a delay of m samples in a multiplied discrete variable, by the time shifting propertyof the Z-transform.
Finally, substituting (4) in (1), and using (2) and (6), the final update equation for the right boundary is"
pi,n+1j,k =1 N− pn−m+1+ Si,n−m+1∑ bm (pn−m+1
i+1, j,k ) + ...
j,ki, j,k
3 m=1
N

+∑

pn−m
i, j,km=1!
√3
am − bm + ...
6
N

− ∑

pi,n−m+2
j,k

m=1

!#
√
3
am + bm÷
6!
√3
a1 + b1 . (7)6To get to (7), the transformation m → m + 1 is needed so that the sums start at 1 because the chosen programming language has no zero-based arrays. For the designed digital filter coefficients happens that am → am−1 andbm → bm−1 in equation (6), e.g., the first coefficients in the sums, a1 and b1 will actually be a0 and b0 in (6),respectively.
The purpose of writing Equation (7) in three distinct parts is because they have different constraints over thesums, so that the for loops iterations have variables with time index greater or equal than 1. The constraints
are given by

n>m,

(8a)

n > m+1 ,(8b)n > m−1

∧

m>2,

(8c)

for first, second and third sum in (7), respectivelly.
Equations for other boundaries are obtained analogously, but the subtracted term in the first sum would be theghost node to be eliminated for that boundary, e.g., for a boundary in negative x direction (x− ), the first partof the update equation becomes
1 N
n−m+1
+ Si,n−m+1− pi−1,∑ bm (pn−m+1
i, j,kj,kj,k ) ,
3 m=1

(9)

and if the boundary material is different, the digital filter coefficients values am and bm would also be different.
The equation for a surface boundary in the form of (7) has a straightforward approach, but equations for boundaries such as edges and corners can get very complicated, so the applied technique was to use an arithmetic
sum over the filter coefficients of each direction. As an example, for a (x+ , y+ , z+ ) corner the first sum for the
adopted update equation was
1 N
n−m+1
n−m+1
+ Si,n−m+1− (pn−m+1∑ (bx m + by m + bz m )[3pi,n−m+1j,kj,ki+1, j,k + pi, j+1,k + pi, j,k+1 )] ,
3 m=1

(10)

and for an (x+ , y+ ) edge was
1 N
n−m+1
+ Si,n−m+1− (pn−m+1∑ (bx m + by m )[2pn−m+1
i, j,kj,ki+1, j,k + pi, j+1,k )] .
3 m=1

(11)

6007

3.3 Grid excitation and constraints
The chosen pulse injection for the grid was a simple unit impulse hard source [7] in a position above the table,
although hard sources can cause some distortion in the results for lower frequencies. This kind of pulse is
known to have a flat response and excite all frequencies in the grid.
In general, the highest valid frequency of the simulation must be half the sampling frequency, given by fN = T1 .
In the case of the rectilinear FDTD scheme, it needs to be even lower [3, 4]: 0.196 fN . Thus, in the routine
was adopted that the Schroeder frequency of the room,√ fS ≈ 500 Hz, must satisfy that fS 6 0.196 fN , so thegrid spacing X must be chosen with care, since X = 3cT . Since c = 345.1 m/s and adopted grid spacing was
X = 0.055 m, then fN ≈ 10868 Hz.

Figure 3. Measurements in the voice booth

Figure 4. Graphical model of the voice booth (inmeters) and the 7 chosen points

3.4 Chosen language and tools
In order to test the model of digital impedance filters for the voice booth, a routine was created using theTorch 2 framework, which has Lua 3 as a base programming language, using LuaJIT compiler. Lua was chosenbecause of its syntax simplicity and speed, and Torch has multidimensional tensors, making implementation a
bit easier.
The boundary conditions were imposed with conditionals for the space. No visual tools or CAD programs wereused to aid in the room modelling.
The arbitrary magnitude IIR digital filter coefficients can be obtained from a function with the same purposein Matlab R 4 (function fdesign.arbmag) or GNU/Octave 5 (function fir2 from the signal package). In orderto get the correct coefficients, it is needed that the filter is modeled in function of the reflection coefficients[6] for each frequency, and the filter order adjusted for the required precision, fitting the resulting curve to the2 http://torch.ch/
3 https://www.lua.org/
4 https://www.mathworks.com/products/matlab.html
5 https://www.gnu.org/software/octave/6008

reflection coefficient curve for the boundary, then using the resulting coefficients.
For both functions, there is a need to adjust the coefficients of the IIR filter modelled for the boundary. If thechosen function for modelling the filter is fir2, for example, there will only be numerator coefficients for R0
since fir2 returns an arbitrary magnitude FIR filter. So, if R0 (z) = r1 + r2 z−1 + ... + rN zN−1 , from (5), results thatb1 = 1 + r1 , a1 = 1 − r1 and so on.
3.5 Field measurements
The measurements for the studied voice booth are described by Araújo, [8], as it is written below.
In order to evaluate the energy decay of the booth, the reverberation time (RT) was measured in situ for eachone of the setups presented in [8]. These measurements followed the procedure from the standard ISO 33822:2012 [9].
The 1/3 octave band frequency response of the room was also measured for all setups using the sound power
level of the pink noise from 50 Hz to 20 kHz characterized according to ISO 3745:2012 [10] and emitted
from an omnidirectional sound source positioned in one of the room corners. The measurement followed themethodology for voice booths presented by Monteiro and Borin [1].

4

RESULTS AND DISCUSSION

The result of the FDTD routine execution is the sound pressure for every node in the modelled voice booth insetup G presentend in [8], for every time sample. The room impulse responses (RIR) of 7 points distributed in
the modelled room were extracted from this field - check Figure 4. The RIRs for the same points were also
obtained from a commercial software that uses geometrical acoustics.

Figure 5. FDTD RIR for position [x y z] = [0.74 0.39 1.14] m, before and after filtering

6009

After this process, RIRs obtained with FDTD and geometrical acoustics were filtered at 500Hz - a low-pass filter
was used for FDTD and a high-pass for geometrical acoustics. Then, by doing a Fast Fourier Transform (FFT)in the filtered impulse responses, and transforming the values to 1/3 octave bands, the frequency responses
were assembled.
For assembling the frequency responses, a scaling was needed in order to adjust magnitudes for the lowerand higher frequencies. This happens because the impulse response obtained with geometrical acoustics isnormalized from −1 to 1, and the FDTD results are in Pascals. The scaling factor was chosen so that the 1/3
octave band frequency response of the room was acceptable.

Figure 6. Comparison of frequency response of
the voice booth

Figure 7. Comparison of reverberation time of thevoice booth

Moreover, for the simulation values, a level of 135 dB was subtracted from the original 1/3 frequency response,
which is the approximate sound power level of the source in the FDTD routine. As it can be noted, theresponse is fitted from 100 Hz to 160 Hz if compared to the measurements, and suffer some distortion 160 Hz
onwards. However, up to Schroeder frequency, the curve for the simulation follows the general behavior of thecurve obtained from measurements. Other causes for the distortions may be the equations chosen for edges andcorners, or differently chosen absorption coefficients for the materials used in the voice booth, since not all thematerials had this values registered.
The reverberation times were obtained for the low frequencies by filtering the RIRs obtained with the FDTD
method with band-pass filters for each band. Thereafter, the results were obtained by calculating a backward
cumulative integral [11] curve for the magnitude of the filtered RIRs, an then fitting a line to this curve forfirst −30dB decay, in order to get the T30 . This time, results extracted from simulation are very near themeasurements, except for lower frequencies, as it can be seen in Figure 7.

5

CONCLUSIONSAs it was seen, the boundary model using digital impedance filters for the rectilinear FDTD grid showed acceptable results for the presented equations, despite differences in some frequency bands and in the equations
for edges and corners. It showed a general similarity with the results obtained from measurements. Also, themethod for measurements in the voice booth may or may not be the most suitable for this type of room, despite
the similar results with the simulation, still needing more investigations [1]. The results for frequencies lower
than 100 Hz were not studied because at the time of the simulations, there was no information of absorption coefficients at these bands for the materials used in the voice booth. Further investigations about the scaling factor

6010

[12] need to be carried out. This factor was chosen in order to make the magnitude of the results obtained bysimulations as close as possible to the ones presented by the measurements, because the RIRs obtained with the
FDTD simulations are in pascals (Pa), and the ones obtained with geometrical acoustics are normalized from
−1 to 1. The reverberation time values obtained from the simulated RIRs showed to be very near the values
obtained by measurements up to the Schroeder frequency, which is the interest region of study for this method.
Execution time of simulation was very slow in gerenal (approximately 2 day of calculations), and in the future
it will be investigated the possibility of paralelizing the code with CUDA R 6 [4, 13] for shortening execution
time.

REFERENCES[1] C. Monteiro and M. Borin, “Acoustic quality evaluation of voice booths using 1/3rd octave band frequency
response,” Inter-Noise, no. October, 2018.
[2] N. J. Giordano and H. Nakanishi, Computational Physics. Pearson/Prentice Hall, 2006.
[3] K. Kowalczyk, Boundary and medium modelling using compact finite difference schemes in simulations of room
acoustics for audio and architectural design applications. PhD thesis, 2008.
[4] M. R. Thomas, Wayverb: A Graphical Tool for Hybrid Room Acoustics Simulation. PhD thesis, 2017.
[5] C. A. B. Kinsler L E, Frey A R and S. J. V, Fundamentals of Acoustics 4th edn (New York: Wiley). 2000.
[6] S. Oxnard, Efficient Hybrid Virtual Room Acoustic Modelling. PhD thesis, 2016.
[7] D. T. Murphy, A. Southern, and L. Savioja, “Source excitation strategies for obtaining impulse responses
in finite difference time domain room acoustics simulation,” Applied Acoustics, vol. 82, aug 2014.
[8] R. Araújo, “Room acoustic design of voice booths in open plan offices,” Inter-Noise, no. June, 2019.
[9] “ISO 3382:2008 - Acoustics – Measurement of room acoustic parameters – Part 2: Reverberation time inordinary rooms.,” standard, International Organization for Standardization, 2008.
[10] “ISO 3745:2012 - Acoustics – Determination of sound power levels and sound energy levels of noise
sources using sound pressure – Precision methods for anechoic rooms and hemianechoic rooms.,” standard,International Organization for Standardization, 2012.
[11] E. Brandão, Acústica de Salas: Projeto e Modelagem. Edgard Blücher, 2016.
[12] Marc Aretz, René Nöthen, Michael Vorländer, Dirk Schröder, “Combined Broadband Impulse Responses
Using FEM and Hybrid Ray-Based Methods,” EAA Symposium on Auralization, no. March 2015, 2009.
[13] C. J. Webb and S. Bilbao, “Computing room acoustics with cuda - 3d fdtd schemes with boundary losses
and viscosity,” in 2011 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP),May 2011.

6 https://developer.nvidia.com/cuda-zone

6011