© A. Erraji et al., Published by EDP Sciences, 2021

1 Introduction

Sound level determination based on outdoor noise sources and sound propagation paradigms applied in urban or rural environments is a key task of environmental noise prediction [1]. The calculation of expected noise immission tables and the generation of geographical noise maps are a well established acoustic field that have seen stable development in research and formulate a legislative foundation for strategic planning, administration and noise policy [2]. Noise mapping techniques based on (weighted) energetic sound pressure levels, however, are not fully sufficient to study noise perception and noise impact on the society. The dynamic aspects of moving traffic in complex scenarios and its conscious perception by human cannot be covered by simplified metrics. Deeper understanding of urban soundscapes includes temporal and spatial features which finally enable the listener to experience the sound field as a whole and interact with sounding objects [3]. This is particularly relevant due to the apparent transformation of current mobility and transport systems towards electric drives, which could be studied virtually using auralization.

An anticipated long-term scientific achievement is a physics-based spatio-temporal sound propagation modelling for complex urban scenes. Especially in the context of urban environments, the feasibility of edge diffraction model integration into real-time implementation for outdoor scenarios must be proven.

In real-time auralization of virtual environments, geometrical methods for the simulation are commonly used, which theoretically incorporate the transfer path from sound source to sound receiver on a geometrical path along the direction of propagation of a plane sound wave.

The simulation of the diffraction of sound (wave-based phenomenon) within geometrical acoustic methods (wave-particle dualism) undoubtedly increases the required computational effort. The key challenge is to determine the geometrical propagation path emitted from a source point to a target point along one or more apex points (locations on edges that correspond to the shortest path of a diffracted wave front) fulfilling Keller’s law of edge diffraction in sufficient approximation to provide a basis for diffraction models, an undertaking that is the imperative subject of this article. Two paradigms can be found in literature, that are considered as a solution to this problem. The first approach is following Fermat’s principle stating that information travels at the shortest possible paths, which results in a mathematical formulation of a non-linear minimization problem. A second approach considers equality of incoming and outgoing angles with respect to the edge, comparable to a specular reflection with respect to the face’s normal [4, 5], which results in a system of linear equations – a problem that can be inverted.

The proposed image edge model follows the latter principle and solves effectively the equation system of equal angles. Furthermore, it integrates specular reflections and is able to deterministically find geometrical propagation paths with an arbitrary constellation of reflections and diffractions. The method is applicable in far field conditions for subsequent edge diffractions. It can be used to further explore path reduction strategies for real-time auralization of urban sound propagation. It is available in open source code.

2 State of the art

The scientific debate about the diffraction phenomenon has already started in the far past (Huygens, Fresnel and Kirchhoff, see [6]).

First mathematical formulations were initially promoted by the optical characteristics of electromagnetic waves, like those from MacDonald in the early 1900s [7, 8]. In 1896, Sommerfeld [9] already dealt with the question of describing the diffraction around thin wedges.

In 1957, the Geometric Theory of Diffraction (GTD) evolved, which was described by Keller [4, 10]. It is closely linked to the description of rays that refract close to edges of diffraction.

The publication of Biot and Tolstoy (BT) the same year can be counted as a milestone within the theory of diffraction because they are the first to describe the solution to a diffraction problem in the time domain by expansion of the wave equation solution into its angular modes in cylindrical coordinates [11]. This interpretation aligns well with the principle of Geometrical Acoustics that interprets acoustic waves as rays or particles. The achievement cleared the way for a series of comparisons and verifications of the approaches of the different theories with measurements [12–14]. In 1981, Medwin expedited the BT method further [15]. In his publication on the diffraction of sound in the acoustical shadow zone, the author could present a closed form of the Biot-Tolstoy formula, which can be applied to predict attenuation of sound barriers with arbitrary opening angles and which, in contrast to the original view, were finite in their dimension (BTM). Other views of the matter with different approaches (e.g. Huygens-Interpretation and multiple diffraction [16]) followed, as well as comparisons of simulation and measurement, which prove the weak points of the Kirchhoff-Approximation and which, once again, portray the relevance of the Biot-Tolstoy (BT) diffraction formula and the extension by Medwin.

Until the end of the 1980’s the aforementioned researchers were engaged in the mathematically correct specification of the problem, whereas, beginning in the 1940’s, the shadowing of sound sources through sound barriers was simultaneously examined by means of a pragmatic approach [17, 18]. Especially important is Maekawa’s publication from 1968, which determines the diffraction of sound at an infinite half-plane in approximation to the reduction of sound energy by experimental determination as a function of detour and wavelength [17]. The result is solely based on geometrical operations to identify the source’s way around the boundary to the receiver. The approach was then extended and applied to all kinds of transmission paths such as sets of houses or noise barriers, whose amount of energy is subsequently accumulated as issued by ISO-9613-2 [19]. To establish simplified procedures with reasonable mistakes of approximation, modifications that combine theory and empiricism also rely on the approach of the specification of shadowing through sound barriers, for example that by Pierce in 1974 [20].

While Maekawa’s by now well-established model, usually referred to as rubber band model, has found application in noise immission simulation software (e.g. via compliance with ISO-9613-2 [19] that make use of Maekawas approximation), more accurate methods remained subject to scientific investigation. In 1999, Svensson et al. could significantly ameliorate the mathematical description of the effects of diffraction in the range of time and frequency. His work on the analytic secondary source model of diffraction of sound through a closed form description – which is based on the BTM-Method and which includes finite edges through a line integral on its apex – is another milestone in the development of diffraction [21]. Up to today, this achievement can be seen as the exact and elegant solution for the description of directivity around a finite diffraction edge that abdicates discrete scanning and develops a continuous impulse response. Several publications followed the BTM-Model extended by Svensson (BTMS), which was reported to be successfully integrated into various room acoustic simulation methods such as mirror sound sources and ray-based approaches [22–26]. Calamia and Svensson developed a discrete procedure based on the continuous integral form that diminishes the computing time yet maintaining a similar quality [27]. Other works deal with the efficient combination of sound paths with the continuity problem at the transition area and different acoustical shadow and half shadow zones [28]. Also, sound diffraction created fields of application for auralization of interactive simulation with the help of various acceleration methods [29–31]. Furthermore, the operation of the BTMS-Method for diffusion of surfaces including acoustical reflectors was discussed by Svensson and Asheim [32, 33].

In 2009, Taylor et al. used techniques of computer graphics to combine the Uniform Theory of Diffraction (UTD) formulation with the ray-frustum procedure, which expands more frusta at defined edges and identifies covered sources in the shadow zone. As Taylor et al. [34] state, using this approach it is possible to achieve update rates sufficient to real-time applications using regular multicore processors.

A fundamentally different yet mention-worthy diffraction modeling approach by Stephenson and Peter Svensson [35] is inspired by Heisenberg’s uncertainty relation, that is proposed to be integrated into a Monte Carlo particle tracing algorithm. It is reported sufficient both for employment in room acoustics and noise immission prognosis [36–38]. The concept of deterministic image edges is mentioned in [39], but is overruled as computationally too complex to handle. Stephenson et al. argued that the exponential complexity class of the mirror image method is inferior compared to the scaling behavior of ray/particle tracing and beam tracing. In contrast to interior spaces, in urban environments the majority of boundaries are not redirecting rays back towards the space where the receiver resides, and therefore stochastic ray tracing suffers from high energy loss due to the diverging character of ray casting. Naturally, the most significant paths often require multiple combinations of specular reflections at surface boundaries as well as diffractions around edges, which results in an unmanageable number of rays to be emitted in order to hit a receiving area with the right amount of acoustic energy.

In 2012, Antani et al. [40] presented an approach to solve multiple reflections and diffraction for finite edges that use a discrete number of connecting rays between edges to be applied with the BTMS diffraction model. Scene geometry is considered static and a visibility tree is constructed that establishes links between scene primitives (triangles and edges). The extension of images of potential diffraction edges upon specular reflections is implemented and paths are validated by traversal of the tree for a given source/receiver pair taking occlusion into account. A particular focusis beeing placed on the visibility of regions from a location such as a source, a receiver or an image point as well as regions among edges, which are considered secondary (line) sources [21]. The additional treatment of diffraction further raises the exponentially increasing number of paths. In 2008, however, Calamia et al. demonstrated the huge impact of path reduction approaches by excluding diffraction components with moderate implication on the spectral error [31]. On the one hand they showed, that diffracted paths with low amplitudes can be ignored at large in a reverberate space. However, audible coloration changes can occur even at low decibel values, which requires closer examination from case to case, especially if diffracted paths are dominant.

We want to point out, that diffuse reflections play a role in the context of urban sound propagation and are reported to be succesfully integrated into interactive sound propagation rendering [41], that also incorporate specular reflections and diffractions. We believe, that our approach can be expanded to handle diffuse components, e.g. with the help of diffuse rain or the radiosity method, but this is out of the scope of this paper.

Therefore, we introduce the concept of the image edge model in the context of image sources and path identification. It should be noted that this approach does not solve double-diffraction problems of objects with edges that are not large compared with wavelength. In other words, the sequence of reflection and diffraction points (“interaction points”) contains far-field conditions between adjacent interaction points, where plane wave incidence is assumed for each edge.

3 Image edge model

When sound propagates through urban environments, there may be cases where the total sound field cannot be fully described considering reflections and diffractions separately. A homogeneous set of paths that only consists of reflected or diffracted energy, even for very high orders, is not sufficient to generate a continuous sound field. Algorithms determining homogeneous propagation paths may lead to fewer or even no audible paths, in contrast to heterogeneous propagation paths, that allow for arbitrary combinations of reflected and diffracted components. The implementation of a heterogeneous propagation path algorithm is significantly more challenging as it requires the integration of diffracting edges into an image model. Therefore, the Image Edge Model (IEM), depicted in Figure 1, is introduced. By mirroring edges along faces and applying the angle constraints for the diffraction coefficients (comparable to the reflected angle of a specular reflection), combined reflection and diffraction paths can be found. The classical Image Source Model (ISM) can be used to determine reflection paths as depicted in Figure 2a, which may lead to audible sound even if objects are occluding the line-of-sight. Homogeneous diffraction paths can be obtained by resolving the Equation System of Equal Angle (ESEA), which has been formulated by Tsingos et al. [42] in the context of a beam tracing approach. It states that the incoming and outgoing vectorial directions at interaction points I must maintain the same angle regarding the edge following:

$$\begin{array}{ll}\overrightarrow{I_{1}S}\cdot \overrightarrow{E_1}& = \overrightarrow{I_1{I}_2}\cdot (-\overrightarrow{E_1})\\ \overrightarrow{I_2{I}_1}\cdot \overrightarrow{E_2}& = \overrightarrow{I_2{I}_3}\cdot (-\overrightarrow{E_2})\\ & ⋮\\ \overrightarrow{I_n{I}_{n-1}}\cdot \overrightarrow{E_n}& = \overrightarrow{I_{n}R}\cdot (-\overrightarrow{E_n})\end{array},$$(1)where S and R represent source and receiver locations, I stands for interaction points, E is the edge unit direction vector and subscripts indicate diffraction orders. A visualization of the angle constraint and the corresponding vectors is shown in Figure 3. The vector notation $\overrightarrow{I_k{I}_{k-1}}$ describes the normalized direction from point I _k-1 to I _k of a segment. The entire segments of equation (1) are representing the shortest path from source to receiver via the considered edges.

Figure 1 Propagation path along three buildings determined by the image edge model.

Figure 2 Comparison of the construction of propagation paths without and with image edges. (a) By use of image source S 1, a first order reflection path can be found. One second order diffraction path is composed of two diffractions at edges E a and E b of the left building. (b) Neither source S nor receiver R are directly illuminating the right building. Irrespective of the order of reflection, neither image sources nor mirrored receivers directly lead to valid reflection paths. (c) By introducing image edges, a valid propagation paths consisisting of reflections and diffractions is found.

Figure 3 Constraint of equal incoming and outgoing diffraction angles and the relationship between direction vectors shown as top-down and side view of a two-building scene with an elevated receiver.

The solution of the ESEA can be obtained by transferring into a matrix and performing an inversion. The problem forms a symmetric tridiagonal matrix that can be solved by Gaussian elimination. However, the angle terms are non-linear due to the Euclidean distance normalization required to apply the scalar product (angle constraint). Hence, rewriting the matrix coefficients results in a system of non-linear equations. To circumvent the non-linearity and keep runtimes manageable, an iterative approximation is proposed. A parametric representation of the apex (interaction) point I is formulated with respect to the vertices P of the diffraction edge

$$I\left(s\right)=P_0+s\cdot \left(P_1-P_0\right).$$(2)

The parameter 0 ≤ s ≤ 1 defines the relative location of the apex point on the edge. By starting with an arbitrary initial value for each segment of the path, the distance between interaction points can be assumed constant while the angular constraint remains. This way, the matrix solution is approximated by solving linear equations only, and the locally correct apex points are determined again. With respect to distances between segments, the diffraction interaction points are evaluated and the current relative location values s is replaced, if a shorter path is found. Because a set of apex points represents the shortest path via a specific set of edges, the solution converges towards the global minimum of the overall path length. Bearing homogeneous propagation paths in mind, the situation from Figure 2 a is now modified slightly. Extending either the lower or the higher corner of the left building, the right building is no longer visible from either R or S. Thus, no valid reflection paths can be generated, regardless of the number of reflections. Only higher order diffraction paths can be found, which are nearly inaudible because the transmission loss of the diffractions into the shadow region can be assumed very large. By mirroring the source or receiver, reflected and diffracted paths can be constructed to a certain extent, as Torres et al. describe in [24]. In general, audible paths can be created by sequentially performing homogeneous path algorithms, e.g., first the ISM and then the ESEA. The mirrored image source S _i or image receiver R _i are constructed first and replace the source S or receiver R in the ESEA. Propagation paths with reflections that precede and succeed diffractions, i.e., fulfilling the sequence $S{I}_F^{\mathrm{*}}{I}_E^{\mathrm{*}}{I}_F^{\mathrm{*}}R$, can be determined this way. The sequence of interaction points is starting at source S and is ending at receiver R. I ^* describes zero or more interaction points, subscripts describe either images at faces F or edges E. Finally, a situation that cannot be covered by this combination of homogeneous path algorithms is easily constructed. If the right building is not directly visible from neither the source nor the receiver, as depicted in Figure 2b, an audible low-order path that starts or ends with reflections cannot be found. i.e. no path complies with sequences $S{I}_F^{+}{I}_E^{\mathrm{*}}R$ and $S{I}_E^{\mathrm{*}}{I}_F^{+}R$. In this case, the Kleene plus operator of I ⁺ represents at least one interaction point I. Only a single yet strongly diffracted path reaches the receiver. Regardless the order, images of source or receiver are not leading to paths with less diffracted components. However, with a practiced eye, a more significant propagation path can be found by a sound wave bouncing off the right building, that is subject to two diffractions, as depicted in Figure 2c.

This problem can only be solved by introducing diffraction into the image model, i.e., by using image edges to determine heterogeneous propagation paths. Arbitrary combinations of reflections and diffractions are constructed according to the visibility tree method described by Antani et al. [40]. Instead of applying extensive from-region and from-point visibility tests for sampled finite edges, we propose to find all theoretically possible shortest paths via scene primitives according to Keller’s law [4], and only perform validity tests between each interaction point only on those potentially audible candidates. For the application of the UTD and Maekawa diffraction model, no further visibility tests are required. However, for the integration of the BTMS model, the edge must be sampled in general and potential occlusion of edge parts during sound propagation must be retrofitted in order to determine the valid integration ranges for the calculation of the impulse response.

3.1 Diffraction for the image model

To integrate propagation paths with $S{I}_E^{+}{I}_F^{+}{I}_E^{+}R$, the concept of image edges is introduced. Together with classical image sources, all combinations of reflections and diffractions can be resolved. This leads to propagation paths with interaction points SI ^* R, whereby each I can be either a reflection I _F or an apex point I _E . In accordance with the generation of image sources in the ISM, an image edge can be created by mirroring an edge along a face as seen in Figure 2c.

In combination with an image source, the image edge can be used instead of the edge to calculate the apex points, as Tsingos et al. [42] already indicated. There, a direct path to the subsequent interaction point can be created and the influence of the reflecting face can be neglected as shown in Figure 2c. Accordingly, the reflection points can be determined by calculating the intersection points on the faces. Higher order reflections can be created by further mirroring the image edges as demonstrated in Figure 4a.

Figure 4 Determination of higher order propagation paths. In (a) the handling of multiple mirror operations of diffractions is depicted, in (b) higher order diffractions are applied and in (c) a reduced number of mirroring operations compared to (b) are applied due to the preserving of angles in Euclidean transforms. (d) depicts the determination of apex points and visibility aspects of image edges in three dimensional space. Analogous to the image source method, the mirrored apex point must lie behind the reflection plane.

When subsequent higher order diffractions are desired, further image edges must be created, as depicted in a straightforward manner in Figure 4b. Hence, by using images of both sources and edges, all combinations of reflections and diffractions can be constructed. In three-dimensional space, an edge is not a single point as it may appear in the two-dimensional view of buildings of Figures 2a–4b, but a line segment. Since an edge can be geometrically described by its start and end point, the image of the edge is fully described by mirroring these points, as shown in Figure 4d. To determine the apex points among reflections, the image edges are used in the ESEA instead. Furthermore, mirroring is an application of an Euclidean transformation. The incident angle as well as the diffracted angle of the original edge and the mirrored edge are preserved. This property can be exploited, because the tridiagonal structure of the ESEA only considers local segments with normalized direction vectors, which allows a reduction of the problem as indicated in Figure 4c. In contrast to the double mirroring as proposed initially by Tsingos et al. [42] and depicted in Figure 4b, only the prior edge of the regarded reflection must be determined. This leads to a significantly reduced number of image edges. The length of the mirrored edge is preserved, too, and the relative location of the apex point s is maintained. Thus, mirrored apex points directly lead to the corresponding apex points. The modified ESEA for the calculation of the apex points is

$$\begin{array}{ll}\overrightarrow{I_{1,\mathrm{in}}{S}_{\mathrm{out}}}\cdot \overrightarrow{E_{1,\mathrm{in}}}& =\overrightarrow{I_{1,\mathrm{out}}{I}_{2,\mathrm{in}}} \cdot (-\overrightarrow{E_{1,\mathrm{out}}})\\ \overrightarrow{I_{2,\mathrm{in}}{I}_{1,\mathrm{out}}}\cdot \overrightarrow{E_{2,\mathrm{in}}}& = \overrightarrow{I_{2,\mathrm{out}}{I}_{3,\mathrm{in}}}\cdot (-\overrightarrow{E_{2,\mathrm{out}}})\\ & ⋮\\ \overrightarrow{I_{n,\mathrm{in}}{I}_{n-1,\mathrm{out}}}\cdot \overrightarrow{E_{n,\mathrm{in}}}& =\overrightarrow{I_{n,\mathrm{out}}R} \cdot (-\overrightarrow{E_{n,\mathrm{in}}})\end{array},$$(3)where the incoming normalized edge direction E _k,in represents the actual direction of the edge E _k and the incoming interaction point I _k,in represents the interaction point I _k of the k-th edge of the propagation path. In case of two subsequent edges k and k + 1, the outgoing parameters E _k,out and I _k,out equal the corresponding incoming parameters. If the k-th edge is followed by a face, E _k,out and I _k,out are represented by the direction of the corresponding mirrored image edge and by the mirrored apex point on the image edge. S _out contains the image source of the last face that appeared prior to the first edge. Instead, if the first interaction point of the propagation path is a diffraction, S _out equals the source S. The modified ESEA can be resolved by the use of the non-linear equation system in combination with the parametric representation of the original and mirrored interaction points. These are

$$I_i\left(s_i\right)=P_{\mathrm{start}}+s_i(P_{\mathrm{end}}-P_{\mathrm{start}})$$(4)and

$$I_{i,k}\left(s_i\right)=P_{\mathrm{start},k}+s_i(P_{\mathrm{end},k}-P_{\mathrm{start},k})$$(5)with the start and end points of the ith edge P _start and P _end and the corresponding image edge P _start,k and P _end,k , respectively. The parameter s _i is the relative location on the original and image edge. The valid range of s _i is between 0 and 1 for a completely visible image edge, but is reduced if corresponding parts of the original edge do not lie in front of the wall plane, as depicted in Figure 4d. The final non-linear equation system can again be approximated efficiently using an iterative Newton scheme.

3.2 Acceleration

Applying the image edge model on a computer for all mesh items of a scenario is an approach that quickly reaches operational restrictions concerning memory and performance. The foremost reason is the exponential complexity of the image algorithm, that does not only grow with the number of polygons, but as well with the number of edges [31]. To accelerate the computation and decrease memory consumption, two measures of different class are proposed.

Firstly, it is desirable to dismiss invalid propagation paths during the search as soon as possible. This way, the final result remains unchanged, but computational load can be reduced significantly. A second acceleration approach has a progressive nature and aborts further processing of paths, if the expected contribution becomes irrelevant. The basis of the decision whether a path is relevant or irrelevant is context-dependent and must be justified carefully.

The pathfinder is constructing an item matrix that holds references among all faces and edges. In outdoor scenarios, many of these items may never construct a valid subsection of a propagation path for any given source and receiver location. This can happen, for example, if faces are representing opposite facades of a building, or two arbitrary edges are entirely occluded by other buildings. If a pre-processing of the geometry data is performed and an illumination matrix that contains visibility among geometrical items is generated, the construction of the propagation tree can be drastically thinned out. The purge of the propagation tree results in reduction of computation time as it allows to skip connections of illuminable mesh items during the generation of propagation path candidates at all image orders, which comes at the expense of a single pre-processing calculation. Pre-processing acceleration is not discussed in detail, but can also include geometrical hierarchy methods, like binary-space partitioning, during construction of the illumination matrix, which requires dynamic adaption in case of to non-stationary geometry. In the context of acoustics in urban environments, buildings (i.e. mesh items) are considered static here.

The second class of acceleration options, the progressive methods, are aiming at dismissing paths during geometrical construction based on acoustic considerations, like signal-to-noise ratio and perception thresholds. By applying a sound pressure level threshold and an angular threshold for diffraction, most paths can be culled with limited spectral error [31]. These threshold values are user-defined and must meet the goals of the application, i.e., one would select conservative thresholds for high quality simulations, but can allow more relaxed values for real-time rendering purposes.

A multitude of different aspects may be decisive for the implementation of such progressive algorithms, however the most apparent may be the consideration of an indiscriminate broadband or frequency-selective sound pressure level threshold, which has been implemented to demonstrate the scope of acceleration. In principle, it can be assumed that sound is always mitigated during transmission along a propagation path, hence the probability of relevance decreases with interaction at boundaries and distance, if non-planar wave forms are emitted.

For those, penalty levels can be introduced for each interaction, that conservatively approximates the transmission loss per section. If more than one source is present in the scenario, each corresponding sound power level becomes a factor that can be integrated straightforwardly and that creates an energetic balance among propagation paths of different origin. A simple formula can be generated, that follows ISO 9613-2 [19]. If the estimated sound pressure level L at the receiver falls below a threshold level, the path is marked perceptually irrelevant. Modifications have been included that assume spherical spreading loss to the first edge, and then apply the even stronger spreading attenuation after diffraction, as shown in equation (6)

$$\begin{array}{ll}L\le L_W& -11 \mathrm{dB}-20 {\log }_{10}\left(r_1\right)\\ & -10 {\log }_{10}\left(\frac{r_2\left(r_1+r_2\right)}{r_1}\right)-L_{\mathrm{Penalty}}\end{array},$$(6) r ₁ (in m) is the length of the propagation path from the source to the first edge and r ₂ (in m) is the length from the first edge all the way to the receiver for simplification reasons. L _Penalty, however, is the accumulated level loss caused by each diffraction and reflection of the propagation path, which is based on the assumption of a minimum expected attenuation. This value can be derived from a statistical evaluation of the absorption at surfaces and the average or best-case diffraction attenuation. Equation (6) is also applicable at the pre-processing stage. If the minimum distance between a pair of edges and/or faces in the illumination matrix exceeds the threshold value, the connection can be dismissed. To determine the value conservatively, we use sphere bounding boxes per mesh item.

Additionally, the accumulated angle of all occurrences of diffraction gives a sound idea about the expected energy transfer, i.e. paths that make a full circle around buildings can doubtfully carry enough energy to be relevant in the presence of other paths. Diffraction components that are subject to multiple deep diffractions into the shadow regions are likely to be attenuated significantly, e.g., a path traveling around a building may not contribute to a relevant scale, as shown in Figure 5a. We therefore allow an abortion criterium to be applied that compares the accumulated angle φ _A against a user-defined threshold value.

Figure 5 (a) Two diffracted propagation paths. On the left building, the path is diffracted around the entire building with an accumulated diffraction angle φ A  > π. The right building shows a diffracted path with φ A > 2π. Depending on the angular threshold value, these paths are dismissed. (b) Wedge with aligned faces F 1 and F 2. Propagation paths to edges that are completely located in the illuminated region of both aligned faces (red area) are diffracted neither from the shadow region of F 2 nor from the shadow region of F 2.

Calamia et al. proposed another angular culling procedure that compares the angular distance between diffraction angle and the nearest zone boundary, which is either the shadow or reflection boundary, with a user-defined threshold ϵ [31]. We depict the deviation angle to the shadow boundary by α and to the reflection boundary by β, respectively. If either one of the angle differences becomes small, i.e. φ _B  = min|α, β| < ϵ holds, the coherent superposition of the direct or reflected sound wave with the diffraction contribution becomes significant. A generalized formulation for higher order diffraction

$${\phi }_B=\sum _i\mathrm{}\min |{\alpha }_i,{\beta }_i|$$(7)is introduced, where paths are culled with φ _B  ≥ ϵ.

Furthermore, diffraction components into the illuminated region (back-scattering) have a potentially low impact on the overall sound field due to the presence of coherent direct (and reflected) energy. Apart from the issue that omitting back-scattering can cause discontinuities in the overall sound field at the reflection boundary [31], we optionally included this acceleration feature in the pre-processing step. The potential speed-up is tremendous, and it is important in order to achieve real-time rates for outdoor scenarios. In our approach, subsequent mesh items that lie entirely inside the illuminated region of the previous edge can be dismissed, as depicted in Figure 5b.

In general, the consequences of the described accelerations options on the auralization result and their audibility (as well as the impact on measures like plausibility and such) must be investigated in the future.

3.3 Pathfinder

The algorithm for the construction of propagation paths – the pathfinder – is able to determine audible paths with every possible combination of higher order reflection and diffraction components. It is a combination of the image source method and the modified equation system of equal angles, which handles the image edges. The implementation can be divided into four major steps, as depicted in Figure 6.

Figure 6 Block diagram of the construction of propagation paths (pathfinder).

The first step intends to construct an audibility matrix based on the geometry of the urban environment. For the discretization, the employment of a half-edge compatible mesh data structure [43] was chosen. The propagation simulation considers the mesh to be static and, consequently, allows for geometry pre-processing. At the expense of memory consumption, computation time is saved to accelerate the pathfinder algorithm, a necessary step to achieve real-time update rates of moving sound sources and receivers.

As a side note, it should be mentioned that dynamic mesh items cannot be handled in this pipeline due to the high complexity of re-computing the pre-processing step. To incorporate occlusion of dynamic geometrical objects like vehicles or variable sound barrier heights, one could run an a-posteriori segment-wise intersection test and replace a dynamically blocked segment with an occlusion algorithm. The rotation-based edge diffraction method would enable real-time rates for diffraction around the convex hull of an occluder, but is in fact not able to add specular reflections or any further interaction with other mesh items. Further methods are applicable in the same faction that incorporate acoustic diffraction based on a-priori knowledge of the geometric occluder. For example, the diffraction filtering neural network generated by a machine-learning algorithm that matches filter coefficients for known geometric objects such as screens, as described by Pulkki and Svensson [44], can be combined with our propagation paths.

Also, the approach by Rungta et al. using diffraction kernels for dynamic interactive situations is an option to add dynamic occlusion, which uses wave-based pre-calculation [45]. Both are reported to generate plausible auralization results and present a meaningful extension to our pathfinder algorithm.

The second step, performed after the geometry pre-processing, requires knowledge on the source’s location. An ordered propagation tree can be constructed, as depicted in Figure 7. The root node S is represented by the sound source itself, and all illuminable mesh items F _i and E _x from its location are building the child nodes. Iteratively, further audible items are added as leaf nodes until a corresponding abortion criteria is reached. The abortion criteria is not necessarily defined as the order of images, but can also be an Euclidean distance metric or the approximation of an accumulated diffraction angle.

Figure 7 Ordered propagation with the source being the root of the tree. The environment consists of two faces F 1 and F 2 as well as of two half-edges E a and E b . Each node contains one mesh item which can be either a face or a half-edge. Faces can additionally be linked to further image sources or image edges depending on the position of the node in the tree. The receiver is connected to colored items that represent a valid path, grey-scale items are not forming valid paths.

Figure 8 An urban environment consisting of buildings at a T junction with asymmetric locations of a source and a receiver.

The pathfinder algorithm, the third step of the procedure, begins with the interconnection of mesh items that interact acoustically with a wave front, i.e. an edge or a face, that finally reaches the receiver. In case of a face, interaction is considered a specular reflection and an image is generated from the parent item. In case of a half-edge, the interaction is considered a diffraction and geometrical wedge parameters are constructed from the linked face as well as the opposite face, which is connected via the opposite half-edge. Half-edges of the mesh data structure are used instead of simple edges because each half-edge provides a fixed orientation and a fixed assignment to a certain face. Although hereby each edge is included twice, the validation of propagation will remove one half-edge and retain the other, if a propagation via this diffracting edge can be found. As the receiver position is introduced, each node can be finalized by performing a validation test regarding audibility, i.e. item-to-item line-of-sight assertion. A valid propagation path represents a sequence of nodes beginning at source via given mesh items down to the receiver. The previously mentioned acoustic reflection and diffraction models are applied here to determine the interaction points and hence seal the propagation path’s properties in the process. The apex points are calculated by using image sources and image edges following the solution of the modified ESEA. Face-related reflection points are then determined by using the projection-space image apex points.

It is worth mentioning, that, technically, the image edges with respect to all face items of the mesh are independent of actual source and receiver locations. A certain set of corresponding parameters can be generated a-priori. However, because the apex points depend on source and receiver data and it is straight-forward to implement image instances during path finding, the edge’s images are not generated during pre-processing but are directly applied in this last step.

An acceleration realization that depends on the source position, in particular the streamlining of the tree by back-face culling and determination of image source locations is proposed in the second step as it drastically reduces memory consumption and computational load without affecting the resulting paths. For further tests, the depth of the tree has been limited to a combined interaction order separated into maximum number of considered reflections and diffractions, respectively. Although this represents an exponentially growing structure, the memory consumption has been found manageable, as each node only holds pointers to a fixed number of items.

Apart from trivial disqualifying properties like back-face location of an interaction point, some of the remaining paths may still intersect with geometries. An additional intersection test is required to be performed on each section between a pair of interaction points of the residual paths. This results in a deterministic subset of valid propagation paths which are ready to be handed over to a processing unit that applies acoustic models and performs the auralization.

Due to the reciprocity of the propagation paths, the receiver and the source can be exchanged in the algorithms, leading to the same propagation paths with inverted order of interaction points. Adaption to dynamic movement of emitting or immitting entities or integration of multiple entities can be incorporated by re-considering and re-performing the third step. Again, the reciprocity principle holds and the direction of propagation path finding can be applied as deemed appropriate.

3.4 Results

The algorithm has been implemented in the programming language C++ and is available as an open-source project named ITAGeometricalAcoustics¹.

An open-science package including the geometry file shown in Figure 9 and a precompiled Windows^® application that makes use of the library as well as the MATLAB^® scripts evaluating the perceptual culling results can be found in the supplementary data. The algorithm has been applied to a simple street canyon situation as depicted in Figure 8.

Figure 9 Visualization of all 2550 propagation paths up to 4th order from source S to receiver R. Yellow dyed segments depict stronger diffraction components, green indicates more reflections in the path.

This non-optimized implementation leads to very high memory consumption even for relatively low orders. Thus, this rather simple environment contains a manageable number of faces and edges that has been chosen over a more complex situation, because it allows non-accelerated calculations up to order 4 of combined reflections or diffractions.

The test environment consists of three buildings located at a T-junction. The total amount of faces is n _f  = 26, and the total amount of edges to be found is n _e  = 52. The center of the T-junction marks the origin of the coordinate system. The roads have a width of 12 m and are considered infinitely long. The buildings have a height of 20 m and a depth of 15 m. The sloping roof building has an additional height of 7.5 m. All results relate to a source position S = (−1 m; 20 m; 1 m) with the position in the side street. A slightly off-centered receiver position R = (50 m; 5 m; 1.7 m) is intentional to avoid symmetries in the propagation paths.

First, the number of propagation paths is analyzed in terms of the reflection and diffraction order. Path-reducing operations are included and compared with the non-reducing implementations as reference result. The algorithm and its accelerations are tested with respect to computational run-time. Finally, the implementation of the algorithm is statistically evaluated by comparing the computed incident angles with reflected and diffracted angles, which are subject to angular displacements due to the iterative approximation in the solution of the modified ESEA.

The resulting propagation paths of the scenario are displayed in Figure 9. In the absence of direct line-of-sight and invalid first-order image sources, it is proven that mainly diffracted paths dominate the scene. Additionally, paths crawl over the rooftop surface of the L-shaped buildings, which are mostly subject to back-scattering, i.e. the next propagation segment (or receiver) is not in the shadow region but in the illuminated area of a canonical diffraction situation (cf. Fig. 5b).

The total amount of various combined reflection and diffraction paths is shown in Table 1 for each scenario.

By considering one or two diffractions, occasional valid paths can be found. However, with increasing order paths are assembled that are of questionable importance to the overall sound field at the receiver (cf. Fig. 9). Optimizing the result in the sense that irrelevant paths should be removed, a progressive evaluation is performed as proposed in the previous chapter. Perceptual path length culling, an accumulated diffracted angle threshold φ _A = π and the exclusion of non-shadowed edges that always lead to back-scattered diffractions have been applied. A sound power level L _W  = 90 dB of a large Diesel vehicle is assumed. The default reflection factor of $R=\sqrt{0.8}$ is applied at all surfaces following ISO 9613-2 [19]. The level penalty for one reflection has been set to ΔL _refl = 0.97 dB. Additionally, a conservatively chosen penalty for each diffraction is applied at a level of ΔL _diffr = 2.5 dB. It is determined at half the attenuation level empirically found at the shadow boundary according to Maekawa [17]. The path is neglected after its sound pressure level drops by 50 dB.

The resulting 37 propagation paths after perceptual culling are depicted in Figure 10. The relative contribution of specular reflections (that are also subject to diffraction) is significantly higher than in Figure 9, as paths traveling over the roof top of the L-shaped buildings are dismissed.

Figure 10 Visualization of remaining 37 propagation paths after perceptual culling from source S to receiver R. Yellow dyed segments depict stronger diffraction components, green indicates more reflections in the path.

In general, the number of paths could be vastly reduced, as shown in Table 2. Figure 11a shows the original impulse train, Figure 11b the pruned impulses after dynamic range culling and Figure 11c the remaining result with all perception-based reductions applied.

Figure 11 Normalized impulse series of T-junction propagation paths under different reduction approaches. (a) No path reduction, (b) dynamic range culling, (c) perceptual culling.

Figure 12 Accumulated run-times for the scenario excluding pre-processing. The dark blue bars depict the inclusion of back-face culling, the yellow bars the additional inclusion of occlusion detection and light blue bars further include the path number reducing methods.

From the initial 2550 propagation paths, a remaining count of 37 are selected that appear practically relevant as to the visualization of Figure 10. While the dynamic range culling has a moderate influence due to the limited spatial extent of the geometrical mesh, the accumulated diffraction angle threshold as well as the removal of edge-edge combinations that lead to double scattering into illuminated diffraction areas (neglecting the non-shadowed edges) reveal a significant impact on the final result.

After applying all perceptual acceleration methods, the total sound path numbers are considered to be in a range where real-time processing is applicable. The omission of the second wave of impulses after 0.3 ms (cf. Fig. 11c) appears harsh when viewed individually. It should be noticed that this example is an extreme case where no direct sound is present. However, it is assumed that the consequence on the auralization result is negligible considering loudness masking effects in the context of an urban scenario with several sources present, some of which contribute with direct sound. At this point, perceptual studies must follow and individual selection strategies for parameters must be validated for various scenarios and purposes. A general recommendation cannot be drawn at this point.

3.5 Run-time of the algorithm

Run-times for the scenario for combined diffractions and reflections up to order four were measured. The ESEA is solved by using initial relative apex points at the mid of the corresponding edges and are refined with five iterations. The obtained mean and standard deviation of the angular displacement error was 1.7° and 6.5°, respectively.

The benchmark was computed on a PC with Windows^® 10 Enterprise operating system. The binary was compiled with Visual Studio 2013 and executed on an Intel^® CoreTM i7-7700 CPU with 3.60 GHz and 8 GB DDR3-RAM. No parallelization optimizations have been used in the code.

Figure 10 depicts the measured run-times with the brute-force method as well as the optimizations from Table 2. From the third order onwards, the computation time is reduced by adding the occlusion detection and the path number reducing methods. In case of the fourth order, the computation time is reduced from more than 120 s to less than one second.

4 Proof of concept

Due to lack of availability of reference data for a large urban scene, the pathfinder and filter synthesis was checked for consistency with basic performance features on a rather simple case. This does not replace a continued validation process with more complex cases, which will be subject to future work. However, the C++ source code, employed applications, configurations and scripts to reproduce the results have been made available under public licenses (see dataset in Supplementary Material). The image edge model has been checked with the help of the round robin scenarios of the Benchmark for Room Acoustical Simulation (BRAS) database² [46]. The setups of scene 5 and scene 6 encourage the consideration of diffraction, for which geometrical propagation paths have been determined both manually and with the image edge model approach. The paths then are processed by means of acoustic propagation models including source directivity, specular reflections, edge diffractions, spherical spreading loss and transmission. To attribute validity to the image edge model, attention must be drawn to the timing of incident wave fronts at the receiver as well as the magnitude of the propagation frequency response Furthermore, notches are expected to be present in the frequency domain that theoretically form a comb filter structure due to coherent superposition of the source’s signal via paths of different delays. However, this effect is deteriorated by the directivity spectra of different directions concerning the source, and are additionally faded as the impact of the diffraction is stronger for those paths reflected on the ground, because the image’s path over the aperture of the barrier is steeper compared to the shortest path.

The scenario of a thin medium-density fiberboard of 25 mm strength has been evaluated for the source position LS01 and receiver position MP01, which are placed at a height of 1.2 m in the front and the back of the board in a symmetrical manner (see [46] for further details).

As shown in Figure 13, the image edge model reliably determines the potentially dominant propagation paths of the scenario including combinations of reflections and diffractions. In total, six main propagation paths are considered relevant in this situation. Because the fiberboard can be considered as an infinitesimally thin plate in the acoustics simulation, double diffractions at the edges are merged into one and are modeled as a single reflection at the infinite plate. However, the image edge algorithm has been executed on the actual geometrical model and hence had to find propagation paths of second order diffraction. Table 3 shows a list of relevant propagation paths for this situation. Back-scattering paths from the corners formed by the ground floor and the fiberboard are dismissed as well as paths circulating around the plate at the vertical edges, which can only be found if diffraction orders greater than 3 are pursued. The path transmitted through the fiberboard has been added manually to the result using the mass law formula.

Figure 13 Automatic geometrical propagation path determination based on the image edge model for BRAS scene 5 (LS01, MP01) [46].

Reflections have been considered without losses, and two different diffraction models have been applied, namely the energetic Maekawa approximation based on empirical data [17] and the Uniform Theory of Diffraction (UTD) based on solutions from optics [47]. The latter is an approximation of the electro-magnetic diffraction solution that is applied for acoustic waves. It includes phase information and delivers a smooth transition from the shadow region into the illuminated region. The directivity of the loudspeaker has been integrated as well. The provided BRAS dataset has been queried for the outgoing directions to the first interaction point of the propagation paths, e.g. reflection or diffraction points. The impulse responses measured at a distance of 2 m have been transformed into the DFT domain and multiplied with the separated transmission filters of the paths. The result then have been time shifted to invert the additional delay introduced by the directivity dataset in order to equalize the additional delay.

To match the simulation results with the provided measurements, the consideration of the loudspeaker frequency response and directivity component is imperative [46], otherwise the fine structure modulating the comb filter structure of the different propagation paths is not comparable. From Figure 14 it becomes obvious, that even in the bandwidth of interest between 100 Hz and 4 kHz significant variation in magnitude is present. Also the texture of impulses in the time domain are evidently influenced by the response of the loudspeaker.

Figure 14 Selection of loudspeaker transfer functions required for the calculation of BRAS scene 5 (LS01, MP01) [46].

Figure 15 shows the simulation results together with the calibrated measurement in the frequency domain in Decibel relative to a free-field transmission with the same equipment. Values are representing the impulse response re. free-field caused by the setup, i.e. floor reflections and diffractions at fiberboard. The comb filter structure ranging from about 300 Hz up to 2 kHz can clearly be recognized and proves that the structure of the attenuation filter matches in the timing as well as magnitude. On the one hand, the increase of deviations in the frequency region below 100 Hz can be explained by room modes in the anechoic chamber that are not well-suppressed [46], by violation of the plane wave approximation of the incident wave, and by limitations of the UTD and Maekawa models. On the other hand, deviations in the higher frequency range are likely to be subject to the approximations in the used diffraction models, for which an infinitely thin barrier is assumed. They may also be explained with missing reflection impedances and the shifted interference pattern due to position uncertainties that can have a significant effect on the comb filter structure revealed in the frequency domain. The simplification of the barrier being an infinitely long and infinitesimally thin partition ignores effects of double diffraction and interaction with boundaries that becomes apparent with decreasing wavelength. If the wavelength is supposed to be greatly larger compared to the plate thickness, the presented frequency range up to 4 kHz can be assumed valid.

Figure 15 Measurement and simulation results of BRAS scene 5 (LS01, MP01) in the frequency domain [46].

Figure 16 depicts the simulation results together with the calibrated measurement in the time domain. The impulse response has been equalized with the magnitude of a free-field transmission with the same equipment in order to maintain the temporal structure. Again, the timing of the incident wave fronts at the receiver position are well met (cf. Tab. 3). The amplitudes differ minimally for the first wave fronts but deteriorate henceforth which can be explained by room modes and interaction with equipment that was not part of the simulation.

Figure 16 Measurement and simulation result of BRAS scene 5 (LS01, MP01) [46] in the time domain.

5 Conclusion and outlook

A complex diffraction pathfinder algorithm was presented which allows for a drastic increase of computational efficiency in auralization of large urban scenarios. The image edge model is introduced and applied to multiple combinations of reflections and diffractions. The energy loss associated with each path leads to exclusion due to irrelevant contribution to the total sound field. After this way of culling, only the relevant paths remain and are stored with their parameters of distances and reflection or diffraction factors.

The results of the computational performance as concerns real-time processing are promising up to second order, which in fact covers the largest contributions to the receiver sound pressure. Whether neglecting higher orders is appropriate or not, however, must be studied in listening experiments for a number of typical scenes.

This approach is well suited for calculation of large distance point-to-point reflection and diffraction effects. It is not applicable for multiple diffraction where the contribution edges are on close distance compared with the wavelengths. This limiting condition can be accepted in urban scenarios but not in indoor scenarios, where object sizes are in the order of magnitude of the wavelength.

The approach pathfinder and the implemented diffraction filters were put to test by comparing an open-data benchmark scenario with a finite-size screen. The diffraction filters based on Maekawa and UTD gave very good agreement with measurement results, A comparison with measurements of a more complex situation including separable propagation paths of higher reflection and diffraction orders would reveal the potential of the method, which is to date not available to the best of the authors knowledge.

Another research direction in the future is to further optimize the culling algorithms. In augmented reality, the existing background noise level well serves as a limiting factor in the culling procedure. In full Virtual Reality (VR), the background noise relevant for one contributing path is in fact the sum over all other paths. The convergence of the culling algorithm is a critical problem, and this will have to be studied in follow-up research. The real-time performance is subject to investigation as well but this depends on the speed of the sources and listeners in the scene. Our experience shows that with slow listener movements, update rates of 1 s with of up to third order may still lead to plausible and smooth auralizations.

Supplementary materials

Supplementary file supplied by the authors. Access here

Conflict of interest

The authors declare that they have no conflict of interest.

Data Availability Statement

The research data associated with this article are included within the article.

Acknowledgments

The authors would like to express their gratitude to the German Research Foundation (DFG, Deutsche Forschungsgesellschaft) for funding the project Auralization of Urban Environments Real time simulation of Diffraction under grant number VO 600/39-1, which made this contribution possible.

References

All Tables

All Figures

	Figure 1 Propagation path along three buildings determined by the image edge model.
In the text

Figure 2 Comparison of the construction of propagation paths without and with image edges. (a) By use of image source S 1, a first order reflection path can be found. One second order diffraction path is composed of two diffractions at edges E a and E b of the left building. (b) Neither source S nor receiver R are directly illuminating the right building. Irrespective of the order of reflection, neither image sources nor mirrored receivers directly lead to valid reflection paths. (c) By introducing image edges, a valid propagation paths consisisting of reflections and diffractions is found.

In the text

	Figure 3 Constraint of equal incoming and outgoing diffraction angles and the relationship between direction vectors shown as top-down and side view of a two-building scene with an elevated receiver.
In the text

Figure 4 Determination of higher order propagation paths. In (a) the handling of multiple mirror operations of diffractions is depicted, in (b) higher order diffractions are applied and in (c) a reduced number of mirroring operations compared to (b) are applied due to the preserving of angles in Euclidean transforms. (d) depicts the determination of apex points and visibility aspects of image edges in three dimensional space. Analogous to the image source method, the mirrored apex point must lie behind the reflection plane.

In the text

Figure 5 (a) Two diffracted propagation paths. On the left building, the path is diffracted around the entire building with an accumulated diffraction angle φ A  > π. The right building shows a diffracted path with φ A > 2π. Depending on the angular threshold value, these paths are dismissed. (b) Wedge with aligned faces F 1 and F 2. Propagation paths to edges that are completely located in the illuminated region of both aligned faces (red area) are diffracted neither from the shadow region of F 2 nor from the shadow region of F 2.

In the text

	Figure 6 Block diagram of the construction of propagation paths (pathfinder).
In the text

Figure 7 Ordered propagation with the source being the root of the tree. The environment consists of two faces F 1 and F 2 as well as of two half-edges E a and E b . Each node contains one mesh item which can be either a face or a half-edge. Faces can additionally be linked to further image sources or image edges depending on the position of the node in the tree. The receiver is connected to colored items that represent a valid path, grey-scale items are not forming valid paths.

In the text

	Figure 8 An urban environment consisting of buildings at a T junction with asymmetric locations of a source and a receiver.
In the text

	Figure 9 Visualization of all 2550 propagation paths up to 4th order from source S to receiver R. Yellow dyed segments depict stronger diffraction components, green indicates more reflections in the path.
In the text

	Figure 10 Visualization of remaining 37 propagation paths after perceptual culling from source S to receiver R. Yellow dyed segments depict stronger diffraction components, green indicates more reflections in the path.
In the text

	Figure 11 Normalized impulse series of T-junction propagation paths under different reduction approaches. (a) No path reduction, (b) dynamic range culling, (c) perceptual culling.
In the text

	Figure 12 Accumulated run-times for the scenario excluding pre-processing. The dark blue bars depict the inclusion of back-face culling, the yellow bars the additional inclusion of occlusion detection and light blue bars further include the path number reducing methods.
In the text

	Figure 13 Automatic geometrical propagation path determination based on the image edge model for BRAS scene 5 (LS01, MP01) [46].
In the text

	Figure 14 Selection of loudspeaker transfer functions required for the calculation of BRAS scene 5 (LS01, MP01) [46].
In the text

	Figure 15 Measurement and simulation results of BRAS scene 5 (LS01, MP01) in the frequency domain [46].
In the text

	Figure 16 Measurement and simulation result of BRAS scene 5 (LS01, MP01) [46] in the time domain.
In the text