Most current simulation techniques rely on drastic simplifications, such as image source methods and ray tracing [1], and do not easily allow diffraction effects, or are non-physical, such as feedback delay methods [2]. FDTD methods, proposed in room acoustics in the 1990s [3], however, are now gaining in popularity [4, 5], but there remains a great deal of work in adapting such methods (originally designed for electromagnetic simulation) to the peculiarities of room acoustics.

**numerical dispersion****air viscosity****moving sources/listening positions****numerical modelling of frequency dependent boundary conditions****modelling of irregular boundaries and obstacles**

The first two features are characteristics of the numerical scheme itself, and must be studied in a first step. Numerical dispersion in simple FDTD schemes leads to very unpleasant artefacts, including bandwidth limitation, and unnatural clustering of modal frequencies, remedied with recourse to more involved schemes and, perhaps, a good model of viscous air damping, which has not been included in most previous time-domain room simulation work. Moving source/listener positions need careful attention in an FDTD framework, due to the possibility of audible artefacts (hearing the grid). The last two features are related to the room boundary (walls and obstructions)—more realistic impedance boundary conditions lead to spectral shaping and frequency-dependent loss in resulting room responses. Specialization to irregular boundaries/obstacles constitutes the final step towards a virtual room emulation.

**References**

[1] U. Zölzer, Ed., *DAFX: Digital Audio Effects*, John Wiley and Sons, Chichester, UK, 2002.

[2] D. Rocchesso. *Maximally Diffusive Yet Efficient Feedback Delay Networks for Artificial Reverberation*, IEEE Signal Processing Letters, 4(9):252—255, 1997.

[3] D. Botteldooren. *Finite-difference time-domain simulation of low-frequency room acoustic problems*, Journal of the Acoustical Society of America, 98(6):3302—3308, 1995.

[4] K. Kowalczyk and M. van Walstijn. *Modeling Frequency-Dependent Boundaries as Digital Impedance Filters in FDTD and K-DWM Room Acoustics Simulations*, J. Audio Eng. Soc., 56(7):569—583, 2008.

[5] M. Beeson and D. Murphy. *RoomWeaver: A digital waveguide mesh based room acoustics research tool*, Proceedings of the 7th International Digital Audio Effects Conference, Naples, Italy, 2004.

### Numerical Dispersion

In simple room acoustics models, sound waves of different frequencies should travel through the air at the same speed. In finite difference (FDTD) methods, “numerical dispersion” causes waves of different frequencies to travel at different speeds. Numerical dispersion manifests as a smearing of transients and a mistuning of room modes. These are artefacts that may be audible. Below is a video of how numerical dispersion acts on a two-dimensional circular wavefront.

Without dispersion:

With dispersion:

In the above video, numerical dispersion causes ripples to be left in the wake of the propagating wavefront.

Below are some sound examples where artefacts due to numerical dispersion may be heard.

The clean sound (no artefacts):

The same sound after travelling one metre at 340 metres per second (sound speed in air) in a typical finite difference simulation running at a 44.1 kHz sample rate:

… after travelling 5 metres:

… after travelling 10 metres:

… after travelling 20 metres:

… and after travelling 50 metres:

Ideally, the travelling sounds should not differ from the original sound (Huygens principle), besides a change in amplitude.

The severity of the numerical dispersion in a numerical scheme may be affected by the way in which space is discretised. There are generally three ways in which to regularly distribute points in three-dimensional space. These are the cubic, the face-centered cubic (FCC), and the body-centered cubic (BCC) grids. These grids (lattices) are also found in the crystal structures of salt (cubic), gold (FCC), and tungsten (BCC), to name a few examples.

Of these three choices, the FCC grid is a better starting point for minimising for numerical dispersion [1]. Below is a simplified display of some dispersion error surfaces at a fixed wavenumber (spatial frequency magnitude), in different spherical directions. As can be seen, the error is smaller for the FCC grid. See [1] for more information.

There is still much work to be done in tackling numerical dispersion while keeping computational costs to a minimum. Concerns of stability at the boundaries [2] must be balanced with available methods for reducing numerical dispersion [3]. See this page for new work using implicit methods.

# References

[1] B. Hamilton and S. Bilbao. *On Finite Difference Schemes for the 3-D Wave Equation Using Non-Cartesian Grids. *Proceedings of the Sound and Music Computing Conference, Stockholm, July/Aug. 2013.

[2] S. Bilbao. *Modeling of Complex Geometries and Boundary Conditions in Finite Difference/Finite Volume Time Domain Room Acoustics Simulation.* IEEE Transactions on Audio Speech and Language Processing, July, 2013.

[3] S. Bilbao. *Optimized FDTD Schemes for 3-D Acoustic Wave Propagation .* IEEE Transactions on Audio Speech and Language Processing, April, 2012.

### Modelling Irregular Boundaries

It is important to be able to accurately represent irregular geometries when discretising space for room acoustics modelling. The simple “shoe-box” geometry can be easily reproduced with a grid of cubic cells, but in most cases (concert halls, recording studios), the room geometry is not an ideal shoe-box. In real listening spaces, sound waves can encounter curved walls, slanted walls, and re-entrant corners (non-convex geometries). These details must be reproduced to convincingly recreate the acoustics of the space.

Finite difference methods sometimes fail to capture the room geometry because they are usually based on regular grids of points (or cells), which can sometimes lead to a “staircase” approximation at the boundaries. The following example will illustrate this problem.

Below and left is a square domain discretised with a grid of square cells. Notice that the square cells fit the outer square perfectly. On the other hand (below right), when the square domain undergoes a rotation with respect to the square cells, a staircase approximation is obtained.

What is needed is an “unstructured grid”, where cells can vary throughout the space. Below is an unstructured grid that conforms to the rotated square domain.

This is also possible in 3D. Below is an analogous example of a rotated cubic domain meshed with cubic cells (left), and on the right, modified at the boundaries to conform to the domain.

And the two examples rotating in 3D:

We can also employ the face-centered cubic (FCC) grid, in which case the volumetric cell is a rhombic dodecahedron. The same rotated cube is now meshed with a staircase approximation and with fitted cells at the boundaries.

Finite element and finite volume are alternative methods which are able to employ unstructured grids. Finite volume methods are closely related to finite difference methods and the two methods can be easily combined. This hybrid modelling results in computationally efficient regular updates on the interior of the domain with finite differences and tailored updates at the boundaries with finite volumes. See [1] for more details.

See [2] for preliminary room acoustic studies with FCC grids and finite volume boundaries.

# References

[1] S. Bilbao. *Modeling of Complex Geometries and Boundary Conditions in Finite Difference/Finite Volume Time Domain Room Acoustics Simulation.* IEEE Transactions on Audio Speech and Language Processing, July, 2013.

[2] B. Hamilton and C. Webb. *Room Acoustics Modeling Using GPU-Accelerated Finite Difference and Finite Volume Methods on a Face-Centred Cubic Grid.*Proceedings of the 16th International Digital Audio Effects Conference, Maynooth, Ireland, September, 2013.

### Implicit Methods

Implicit time integration may offer improved accuracy in finite difference schemes, which can lead to improved numerical dispersion (an ongoing problem) over explicit time integration. For an example, compare the following videos.

First, a two-dimensional slice of an exact solution (in 3-D):

Next, an approximation to the above solution, obtained from the standard explicit scheme:

Note the trailing ripples; this is caused by numerical dispersion.

A more accurate approximation can be obtained from an implicit scheme:

This case is more similar to the exact solution.

While implicit schemes can offer improvements over explicit schemes, they do come at a higher computational cost; a linear system of equations must be solved at each time-step.

Simple iterative methods, such as the Jacobi method, can be used in order to solve the linear system of equations. This an approach that is commonly employed in computational fluid dynamics. Jacobi iterations reduce to sparse matrix-vector multiplications that are easily parallelised on a GPU.

Using the Jacobi method, favourable implicit schemes can be designed that, for a given dispersion error tolerance (below 1%), are more efficient than explicit counterparts, even after taking into account the additional computational costs. Such schemes are also amenable to GPU implementations. See [1] for more information.

Some comparative sound examples are provided below.

The following is the `clean’ sound that has travelled roughly 7 metres in air without numerical dispersion effects.

The following three examples were obtained from GPU simulations modelling a 525 cubic metre room (11.7 m x 8.0 m x 5.61 m) with a specific acoustic admittance of 0.9 at the walls (very absorptive).

This is a sound obtained from the simplest (7-point) scheme running at 32 kHz, after a distance of 7.7 metres travelled

and from an explicit 19-point scheme running at 32 kHz

and from an implicit 19-point scheme running at 32 kHz

The last two examples are better than the first, and the third is incrementally better than the second, but work remains to be done.

It is important to keep in mind that these sounds are based on a lossless wave equation model. In reality, viscothermal effects may play a large role in higher frequencies. For more information, see this page on modelling viscothermal effects in air.

[1] **B. Hamilton, S. Bilbao and C. Webb**. *Revisiting Implicit Finite Difference Schemes for Three-Dimensional Room Acoustics Simulations on GPU.*Proceedings of the International Conference on Digital Audio Effects, Erlangen, Germany, 2014. download paper

### Viscothermal Effects

The lossless wave equation is the usual starting model for room acoustics modelling. A more detailed model should take into account sound absorption in air due to viscosity and thermal effects [1] and possibly other relaxation effects [2].

The simplest model of sound attenuation in air is due to Stokes [3], and this classic model predicts an attenuation that varies as a the inverse-square of the frequency. To illustrate these attenuation effects, consider the following two videos.

The first example is obtained from a lossless wave equation scheme:

The next example is obtained from a viscothermal wave equation scheme:

Note the smoothing effect; this is due to the presence of the viscothermal component.

In even a reasonably-sized room, these additional details can have significant effects on the modelled sound. The following are two sound examples from GPU simulations for a modelled 284 cubic metre room (9.5m x 6.5m x 4.6m) with a specific acoustic admittance of 0.02 at the walls (slightly absorptive).

The first example was obtained from an explicit scheme on the FCC grid, without viscothermal effects.

The next example was obtained from the same scheme, but with viscothermal effects that correspond, roughly, to standard indoor conditions.

[1] Morse, P.M. and Ingard, K.U. “Theoretical Acoustics”, Princeton University Press (1986)

[2] A.D. Pierce. Acoustics: An Introduction to Its Physical Principles and Applications, Acoustical Society of America (1989)

[3] Stokes, G.G. “On the theories of the internal friction in fluids in motion, and of the equilibrium and motion of elastic solids”, Transaction of the Cambridge Philosophical Society, vol.8, 22, pp. 287-342 (1845)