PROCEEDINGS of the
23rd International Congress on Acoustics
9 to 13 September 2019 in Aachen, Germany

Exploration of the pseudospectral method for auralizations
Jorge PETROSINO1; Lucas LANDINI2; Andrés BONINO RETA2; Georgina LIZASO2
1,2

Universidad Nacional de Lanús, Argentina

ABSTRACTSound propagation in auralizations can be simulated with geometric or wave-based models. Each one offers
particular advantages and disadvantages. Geometric models lead to the ray tracing method and the image
source method. Both methods achieve good results for high frequencies and are efficient in large rooms withcomplex structures, but are not able to represent in a simple manner specific wave phenomena, such as
diffraction. Typical wave-based models such as the finite element method, the boundary element method, orthe finite difference method are characterized for achieving very precise results for individual frequencies
applied to small and moderately sized rooms. A typical disadvantage is the need for high computing power,
mainly due to the necessity of running separate simulations for individual frequencies because of numericaldispersion. In this paper, auralizations with the k-space pseudospectral method are explored since it ispossible to reduce or eliminate the problem of numerical dispersion with this technique. For certain
conditions, complete range impulse response can be obtained from a single simulation execution. The room
response, including diffraction and diffusion is analyzed for different scenarios.
Keywords: Auralization, Pseudospectral Method, Numerical Methods

1. INTRODUCTION
Auralization is defined as the technique of creating audible signa ls from numerical methods basedon simulated or measured data that describes an acoustic space, sources and listener positions (1). A
first approach to auralization included physical scale models (2). The process is based on obtaining theimpulse response of an acoustic space. This requires establishing a point in the acoustical space where
the source will be. A mathematical model is also necessary to represent the acoustic space in which tosimulate the sound propagation as well as the information related to the location of the receiver and its
characteristics.The impulse response of the model is usually obtained from processes based o n either
geometrical or wave-based models, each with its own particular advantages and disadvantages.
Strictly speaking, the ecogram resulting from geometrical models does not include all the expectedcharacteristics of an impulse response, since specifically wave-based phenomena are notstraightforwardly taken into account.
Among wave-based models, the finite element method (FEM) and the boundary element method(BEM) are worth mentioning. These are characterized for achieving very precise results for single
frequencies on small or medium sized rooms (1). These methods solve differential equations by
discretizing continuous variables, but, in turn result in unwanted numerical dispersion. In practice, thishinders the study of cases in which the source radiates a complex audio wave (for example, spoken
word) or an impulse, therefore, they are limited to individual frequencies. To obtain the impulse
response, the complex task of independently solving individual frequencies or narrow frequency bands
is required, to later be combined into the total impulse response. This entails a considerable increase in
the total computation time.
This paper proposes the usage of a different method for solving wave propagation differential
equations in order to produce auralizations. The k-space pseudospectral method could make possible
simulations for both complex audio waves and an impulse in a single run, deeming it an interestingalternative to traditional wave-based methods.

1jorgepetrosino@gmail.com

5497

2. COMMONLY USED METHODS
2.1. Geometrical acoustics
For geometric models to be considered adequate wavelengths are required to be smaller than thesurfaces on which sound bounces. Therefore, the assumption that sound waves can be treated as energy
rays establishes considerable restrictions on the frequency range, since for the influence of small
irregularities to be taken into account such simulation would be limited to very high frequencies. The
turning point between wave and ray behavior is gradual and difficult to determinate for a single
frequency number but complementary processes to simulate phenomena such as diffraction anddiffusion can be additionally employed.
The assumption that a range close to Schroeder’s frequency can be considered as a lower limit forgeometric methods is widely accepted. Due to the great superposition of room modes and diffuseness,certain variables can be considered of a statistical nature (4, 5, 6). The image source method is notsuitable for curved surfaces, diffuse reflection, diffraction, or for complex rooms that represent hidden
image sources. Interference can be simulated using additional running time phase, nonetheless, this is
not precise on low frequencies where geometric methods are not valid.
For characteristics that depend on frequency (like absorption or diffusion) ray tracing must be donefor each frequency band of interest separately to then be convolved with any audio file. This produces
a signal representative of what a listener would perceive on the location of choice.
2.2. Wave-based models
On the other hand, wave-based models are not frequency limited, but difficulties arise whenmanaging low computational times and defining physical elements with accuracy. FEM and BEM
models are not limited on room geometry, diffraction, diffusion or interference but the elements
represented on the models are based on their physical properties, such as medium density, medium
propagation speed or compressibility. Simulation efficiency and accuracy is strongly related todiscretization resolution and equation solver choice. In regards to auralization, these methods are morecommonly used for simulating low frequencies in small and medium-sized rooms; due to thecomputational effort they require (7).
One of the main problems of numerical methods lays on the fact that differential equation solvers
exhibit numerical dispersion resulting from wavelength related discretization error (8). This, in turnproduces the unwanted and un-physical effect of propagation speed varying with frequency. This
forces simulations to be carried out for single frequencies at a time and later be combined into the final
room impulse response.
It is noteworthy that either, geometrical or the formerly mentioned wave bas ed models, usuallyrequire, for different reasons, multiple simulation runs to produce a complete impulse response . As
Aretz and Cols clearly state (9): "…simulation results were calculated for frequencies ranging from
200Hz to 2.5kHz in frequency steps of 1Hz…". Since computation time is directly proportional to thenumber of frequency steps and assuming a computation time of 5 minutes per frequency (1) this would
take 190 continuous computation hours.
From a computational point of view, a viable method that made this possible in a single run could be
an interesting alternative to already in use methods, even if such method did not offer any particularadvantage in a strictly computational sense.

3. K-SPACE PSEUDOSPECTRAL METHOD
Spectral methods are a class of numerical techniques that solve differential equations. They arebased on the principle that a solution can be obtained as the sum of certain basic functions ( for example,
sine functions) to then calculate suitable coefficients for the differential equation (10).
They are closely related to finite element and finite difference time domain methods (FDTD). The
principal difference is that spectral methods use basic functions that are solved for the whole domain,
while for other methods functions are solved across small subdomains. It could be argued that spectral
methods use a global approach instead of a local approach (11, 12).
Pseudospectral methods are a sub-class of spectral methods that have been adapted to represent
functions on a uniform rectangular grid, which simplifies how certain operators are evaluated andincreases computation speed by applying the Fast Fourier Transform. The level of precision between
spectral and pseudospectral methods is very similar (13, 14).

5498

For acoustic wave propagation, variations of pressure, density, temperature, and particle velocity
among others take place. Those variations can be described by coupled first order differential
equations related to conservation of medium mass, energy and momentum. In acoustics it is common
to combine these equations in relation to a variable (sound pressure or particle velocity) resulting in asecond order differential equation.
The k-space pseudospectral method can be applied to simulate sound wave propagation by solving
the coupled first order differential equations instead of the second order version (15). There are two
reasons behind this. Firstly, it makes possible the inclusion of mass and force sources quite
straightforwardly. In second place, because it enables the use of a special anisotropic Perfectly
Matched Layer (PML) in order to greatly reduce computation time by reducing grid size, incomparison to other similar methods (16).
The k-space pseudospectral method has been greatly utilized on the last few years on biomedical
applications related to ultrasonic wave propagation on tissue (17, 18). It has also been used to study
alterations on mechanical wave propagation produced by the variation on atmospheric temperature(19). Not many examples of this method applied to room acoustics are to be found (20).
By solving the coupled first order differential equations separately, different operators can be used
to reduce numerical dispersion enough for broadband simulation to be practical. Pseudospectral
methods replace the spatial gradient calculation by FDTD usually used in other methods with theFourier collocation method. This, in itself, reduces numerical dispersion for the spatial gradient but
not for the temporal gradient. Nonetheless, this can be corrected for a homogeneous medium up to the
Nyquist frequency employing a different operator (21).
About less than a decade ago Treeby, Cox (University College London) and Jaros (Brno Universityof Technology) developed a MATLAB toolbox named k-Wave that efficiently utilizes the k-space
pseudospectral method. The toolbox was designed for acoustic wave propagation on the time domain
for complex media. Specifically, simulates ultrasonic waves in biological tissue for medical imaging(22). An unlimited number of sound sources can be arbitrarily defined using any kind of MATLAB
function with relative ease. Simulations can be carried out in one, two and three dimensions,combining pressure and velocity sources interchangeably. Pressure and particle velocity sensors canbe arbitrarily placed in any and all grid points without noticeable computation impact, since themethod solves the equations for the entire grid anyway. Besides, medium density, propagation spee d
and attenuation coefficients can be defined for any and all grid points independently.
These aspects are sufficiently interesting to begin an exploratory analysis on pseudospectral
method oriented auralizations based on producing a complete impulse respo nse with a single
simulation.

4. SIMULATION PARAMETERS AND CONSTRAINTS
Similarly to FEM, pseudospectral methods solution stability is evaluated by the CFL (Courant Friedrichs - Lewy) coefficient shown in Eq. (1) that describes the rate at which a wave travels across
two contiguous grid points over two interleaved time points.
ο‫ݐ‬
(1)
ο‫ݔ‬
Different numerical methods require different CFL values. An example of CFL = 0.3 with 2 grid
points per wavelength on pseudospectral methods would require a Δt step five time less stringent than
FDTD of equal CFL value with 10 grid points per wavelength (23, 24). Simulations were carried out
with a value of CFL = 0.3 for managing a reasonable balance between precision and computation time
for homogeneous and slightly heterogeneous media. Internally K-Wave utilizes a Δt obtained by Eq.
(2).
‫ ܮܨܥ‬ൌ ܿ௠௔௫ ήο‫ ݐ‬ൌ

‫ ܮܨܥ‬ή ο‫ݔ‬
ܿ௠௔௫

(2)

Considering CFL = 0.3, Δx = 1 cm, ܿ௠௔௫ = 343 m/s, results in Δt = 8.75 μs. This Δt is not the same
sample period usually used in audio processing. The maximum valid frequency for temporal functions
corresponding to sources and simulated recordings is obtained by Eq. (3), where ܿ௠௜௡ is the minimum
medium propagation speed.

5499

݂௠௔௫ ൌܿ௠௜௡ʹǤ ο‫ݔ‬(3)

5. IMPULSE PROPAGATION IN 3D SPACES
A three-dimensional space of 128 x 128 x 128 grid points was defined for simulations with morethan 16 million nodes in total, representing a cubic acoustic space of 2.56 m. A spatial discretization
resolution of 2 cm for each spatial dimension was defined. The maximum valid frequency obtained by
Eq. (3) was 8575 Hz. The center of the space is considered the coordinate origin in all three
dimensions.
Propagation of an impulsive signal was simulated for a free field open acoustic space. An acousticpressure point source was positioned at the coordinates (0 m; -1 m; 0 m), and an omnidirectional point
pressure sensor located at (0.5 m; 0.5 m; 0.5 m). The distance between source and sensor was of 1.66
m. The source was defined in the k-Wave toolbox as a dirichlet type source. This means that the source
output is forced to have a time evolution that exactly matches the given signal, in this case, an impulse.

Figure 1 – a) Impulse propagation. Point sensor (black dot on the left) and point source (right); b) Time plotof the impulse response; c) Spectrum of the impulse response
Computation time on a personal computer (i7 series 3.4 GHz processor with 16 GHz of RAM) wasof 17 min. Figure 1b shows the frequency response of the pressure recorded at the sensor. Ideally therecording would have a constant spectral modulus for all frequencies. The former equations predicted
a maximum valid frequency of 8575 Hz for the simulation. Above this value simulation stability
lowers at a great rate as frequency increases. Filtering the emitter beforehand is recommended andresults in a very smooth frequency response with a marginal variation close to ݂௠௔௫ (Figure 1c). Since݂௠௔௫ is proportional to discretization level different choices of ο‫ ݔ‬and ο‫ ݐ‬directly relates to higher orlower values of ݂௠௔௫ . This must be taken into account as a simulation parameter and chosen
accordingly to what is being simulated since it has a direct impact on total computation time.
The computation time reported above only considers enough time for the simulated impulsive wave
to reach the sensor as proof of concept compared to FEM. For an actual auralization a higher
computation time is expected. More simulation time would be required for multiple reflections toreach the sensor from surfaces across the room. In any case, computation time is less than the expectedby individually simulating such a great range as 8.5 kHz with other methods.

6. SIMULATION
6.1. Complex wave propagation in 2D
As previously stated, the auralization process is based on obtaining the room response to thenconvolve with any audio signal. This produces an audio file with similar characteristics as to what alistener would perceive at the simulated receiver’s position.
Clearly, simulating just the propagation of an impulse is far more comp utationally efficient thansimulating the propagation of an entire musical or spoken program. For an impulse, it is only necessary

5500

to simulate enough time for reflections to decay to an acceptable level. This is independent of themusical program duration. Convoluting the desired signal with an impulse response generates similar
final results. For doing this, the entire system is required to be linear time-invariant (LTI) to the input
signal. System input is the signal emitted by the source, the output being the record and the system theentire space on which waves propagate.
With a k-space pseudospectral method it is possible to evaluate the propagation linearity of th esystem. This was tested by obtaining an auralization of the spoken word ‘hello’ previously recorded on
an audio file by simulating the propagation in 2D of the entire signal directly as the source’s input andthen comparing it to the auralization produced by impulse response convolution.
Propagation of complex audio waves is not possible with traditional methods but it is feasible withpseudospectral methods. Simulating the propagation of an entire complex program has an
unreasonably high computational cost, for this reason, it was decided to simulate only two dimensionsof propagation for the test. In the 2D simulations space resolution was set to 1 cm, rising the maximum
frequency for audio signals.
6.2. Direct propagation of an audio file
A homogenous medium (air) was defined with ୫୧୬ and ୫ୟ୶ being equally 343 m/s. The resulting୫ୟ୶ = 17150 Hz. The simulated space was of 2.56 m by 2.56 m. (256 by 256 elements). A source was
located at coordinates x 0 = 0 m; y0 = - 1 m and a pressure sensor at at x0 = 0 m; y0 = 1 m. The source wasdefined as type dirichlet (forcing the output to match as closely as possible the input signal). Only
propagation was tested, therefore no obstacles are present.

Figure 2 – a) Wave propagation for a given instant in the simulation; b) Time plot of original emitted audio
and wave recorded at the sensor
The source time evolution was defined by the spoken word file. The recording was comprised of abrief word (hello) with a time span of 317 ms and a sample rate of 48 kHz. Since k-Wave works witha Δt much lower than the audio, sample period resampling was required to match the simulation
parameters. The MATLAB function resample was used to sample it up to 114 kHz.
The simulation time required for the propagation of the entire program to reach the sensor 2 meters
away took 4 min of computation time for a personal dektop computer (i7 series 3.4 GHz processor with
16 GHz of RAM). Figure 2a shows shows wave propagation. Figure 2b displays the time plots of theemitted and received waves. It can be noted a brief delay at the beginning of the recorded wave which
agrees with the propagation time required to travel the 2m distance between source and sensor. The
resulting signal was played directly in MATLAB. Reproduction was audible clear with a slightincrease in the lower frequencies, as it is expected by the use of dirichlet type source in 2D spaces.
6.3. Audio convolution with the impulse response
For the second simulation the former setup was maintained but an impulse source was used instead.
A slightly longer than the necessary for propagation was chosen for simulation time (6.3 ms). In thiscase computation time was of 6 seconds. Figure 3 shows the time plot of the impulse response. The
amplitude raise at low frequencies is related to propagation in 2D spaces but not on 3D ones as seen in
the proof-of-concept simulation.

5501

Figure 3 – Time plot and spectrum of 2D impulse response
The impulse response was convolved with the original audio sampled up to 114 kHz (to maintainthe k-Wave internal Δt). The resulting signal was again clear and perfectly comprehensible audio withthe same slight increase on the lower frequencies. As expected, a short delay of about 6 ms identical tothe previous simulation is present.
6.4. Wave phenomena
Simulations were carried out to test the performance of k-Wave to simulate wave propagation in 3D
spaces. Two scenarios were used to simulate wave phenomena in non-suitable contexts for geometricalmodels. First, a rigid surface with a single hole, and second, the reflection pattern on a Schroeder
diffuser.
6.4.1. Scenario a - Diffraction
An ideally rigid surface with a square aperture was simulated. To achieve a compromise between
computation time and simulation accuracy a spatial resolution ο‫ = ݔ‬2 cm was defined for a grid of 96
x 96 x 64 points that represented a volume of 1.84 m by 1.84 m by 1.28 m. The resulting ݂௠௔௫ = 8575
Hz. The aperture of width W = 12 cm was 4 cm deep. Multiple cases for incident waves of different
frequency were simulated.
An incident plane wave of f = 4 kHz was generated. The diffraction pattern is shown in Figure 4.
The effects of increase λ in relation to width W are show in Figure 4c (on the rigth).

Figure 4 – Horizontal and vertical projection of the model.
6.4.2. Scenario b - Diffusion
A Schroeder diffuser of N = 7 for 800 Hz was simulated using the same parameters as in thediffraction scenario. The diffuser was defined as an ideally rigid body shaped. It was 24 cm wide by 12
cm deep and infinitely long to avoid diffraction on the edges other than the axis of diffusion. Multiple
cases for incident waves of different wavelength were simulated (λ/W=5.7; λ/W=2.8; λ/W=1.4;
λ/W=0.7). Incident energy is scattered in different directions as the wave reflects on the irregularsurface of the diffuser. Diffusion results are based entirely on the geometric distribution of wells, theirlength and shape. Figure 5 shows a horizontal projection of the waves and a polar plot of the energy
distribution obtained.

5502

Figure 5 – Wave reflection on a Schroeder diffuser

7. CONCLUSIONS AND FUTURE WORK
Results suggest that the k-space pseudospectral method may be suitable for use in auralizations.
The principal advantage relies on the use of impulsive sources with broad range of frequencies. The
tests carried out show that it is possible to obtain an impulse respo nse in just one computation run.
The propagation of a complete audio wave compared with the propagation of an impulse and later
convolution is consistently to what is expected of LTI system. The transfer function of the simulation
system is practically flat for three dimensions with a frequency range adequate for audio application.
For simulation on two dimensions an increase in level for low frequen cies is found but could becompensated by filters.
The method is capable of representing wave phenomena but a more detailed analysis must be
carried out to determinate how much agrees with physical behaviour under simulation constraints.
Preliminary tests have been carried out in larger room including a variety of reflecting surfaces
configurations that render interesting results. This must be verified and submitted to careful
examination before being reported.
Different materials could be simulated based on density, propagation speed and absorption tosimulate wave propagation in more complex conditions. A greater exploration of method limitations
must be carried out to verify if this does not hinder the advantage of simulating impulses on just one
simulation run.

ACKNOWLEDGEMENTS
Authors would like to gracefully thank and acknowledge the research team comprised of Damián
Fernández, Ianina Canalis, Ian Kuri y Nicolás Casais Dassie for their contributions.
This paper is part of a research project currently in progress at Universidad Nacional de Lanús

REFERENCES
1. Vorländer M. Auralization: fundamentals of acoustics, modelling, simulation, algorithms and acoustic
virtual reality, Berlin, Germany: Springer Science & Business Media, 2008
2. Meares DJ. The Use of Scale Models in the Acoustic Design of Studios. SMPTE J. 1978; 87 (10):684-687
3. Naylor GM. ODEON-Another hybrid room acoustical model. Appl. Acoust. 1993; 38 (2-4): 131-143
4. Schroeder M. Die statistischen Parameter der Frequenzkurven von groβen Räumen. Acta Acust. united
Ac. 1954; 4(5): 594-600
5. Schroeder MR, Kuttruff KH. On frequency response curves in rooms. Comparison of experimental,
theoretical, and Monte Carlo results for the average frequency spacing between maxima. J.Acoust. Soc.
Am. 1962; 34 (1): 76-80

5503

6. Fuchs HV, Applied Acoustics: Concepts, Absorbers, and Silencers for Acoustical Comfort and NoiseControl: Alternative Solutions- Innovative Tools-Practical Examples, Berlin, Germany: Springer
Science & Business Media, 2013
7. Vorländer M. Computer simulations in room acoustics: Concepts and uncertainties. J. Acoust. Soc. Am.,
2013; 133 (3): 1203- 1213
8. Smith JO, Physical audio signal processing: For virtual musical instruments and audio effects. San
Francisco, CA, USA: W3K Publishing, 2018
9. Aretz M, Nöthen R, Vorländer M, Schröder D. Combined broadband impulse responses using FEM andhybrid ray-based methods. Proc. EAA Symp. Auralization, Jun. 2009, p. 201-206.
10. Canuto C, Hussaini MY, Quarteroni A, Zang TA. Spectral methods: Fundamentals in single domains,
Berlin, Germany: Springer-Verlag, Berlin, 2006
11. Solomonoff A, Turkel E. Global properties of pseudospectral methods. J. Comput. Phys. 1989; 81 (2):239-276
12. Huntington GT, Rao AV. A comparison between global and local orthogonal collocation methods forsolving optimal control problems. Amer. Contr. Conf. IEEE. 2007. p.1950-1957.
13. Orszag SA. Comparison of pseudospectral and spectral approximation. Stud. Appl. Math. 1972; 51(3):253-259
14. Virieux J, Calandra H, Plessix RE. A review of the spectral, pseudo-spectral, finite-difference andfinite-element modelling techniques for geophysical imaging. Geophys. Prospect. 2011; 59 (5):(794-813)
15. Firouzi K, Cox B, Treeby B, Saffari N. A first-order k-space model for elastic wave propagation inheterogeneous media. J. Acoust. Soc. Am. 2011; 129 (4): 2611-2611
16. Liu QH. The PSTD algorithm: A time ‐ domain method combining the pseudospectral technique andperfectly matched layers. J. Acoust. Soc. Am. 1997; 101 (5): 3182
17. Jing Y, Meral FC, Clement GT. Time-reversal transcranial ultrasound beam focusing using a k-spacemethod. Phys. Med. Biol. 2012; 57 (4): 901
18. Treeby B, Tumen M, Cox B. Time domain simulation of harmonic ultrasound images and beam
patterns in 3D using the k-space pseudospectral method.,” Int. Conf. Medical Image Computing
Computer-Assisted Intervention, Berlin, 2011. p. 363- 370.
19. Hargreaves J, Kendrick P, Von Hünerbein S. Simulating acoustic scattering from atmospherictemperature fluctuations using a k-space method. J. Acoust. Soc. Am. 2014; 135 (1): 83-92
20. Pind F, Mejling MS, Engsig-Karup AP, Jeong CH, Strømann-Andersen J. Room Acoustic Simulations
using High-Order Spectral Element Methods. Euronoise. 2018. p. 2085-2092.
21. Treeby BE, Cox BT. k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic
wave fields. J. Biomed. Opt. 2010; 15 (2): 021314
22. Robertson JL, Cox BT, Treeby BE. Quantifying numerical errors in the simulation of transcranial
ultrasound using pseudospectral methods. IEEE Int. Ultrasonics Symposium (IUS) 2014: (2000-2003)
23. Liu QH. The pseudospectral time-domain (PSTD) algorithm for acoustic waves in absorptive media.
IEEE T. Ultrason. Ferr. 1998; 45 (4): 1044-1055
24. Zhao G, Liu QH. The unconditionally stable pseudospectral time-domain (PSTD) method. IEEE
Microw. Wirel. Co. 2003; 13 (11); 475-477

5504