Proceedings of the Institute of Acoustics

 

 

Raytracing in Galerkin boundary integral form

 

Amal Emthyas¹, Acoustics Research Group, Newton Building, University of Salford, Salford, M5 4WT, United Kingdom
Jonathan A². Hargreaves, Acoustics Research Group, Newton Building, University of Salford, Salford, M5 4WT, United Kingdom

 

ABSTRACT

 

Raytracing is an established method for computing the late-time part of room impulse responses. But it has the drawback that only very crude Monte Carlo models of boundary scattering and diffraction are possible to include without losing its attractive computational cost scaling. This happens because higher-resolution models of these processes output multiple child rays for every parent ray received, making the number of rays grow with reflection order. An emerging solution is so-called ‘Surface-Based’ Geometrical Acoustics. Here the distribution of rays arriving at a boundary is mapped onto a predefined set of spatial elements and angular interpolation functions, producing a vector of coefficients. Re-radiation of subsequent reflections is then a matrix multiplication, with the steady state solution being solvable via a Neumann series. Since rays only ever propagate one reflection order before being collected, diffraction and scattering process that cause them to multiply can be included without issue. Here we present a new formulation based on a Galerkin Boundary Element Method (BEM). A unique feature is its ability to readily change interpolation functions, so their effect on accuracy and convergence can be assessed. In this preliminary work, it is verified against an Image Source Model for a rectangular room.

 

1. INTRODUCTION

 

Prediction models are at the heart of modern acoustic engineering. They are used in a diverse range of applications, from refining the acoustic design of classrooms and concert halls to predicting how noise exposure varies through an urban environment. They also allow Auralisation to be performed for buildings and spaces before they are built or long after they are lost [1] .

 

Geometrical Acoustics (GA) algorithms currently dominate in room acoustics applications [2]. This is an umbrella category for several algorithms, including Image Source Method (ISM) and Raytracing. These are common in their assumption that sound propagates geometrically, i.e., in straight lines in a homogeneous medium such as air in a room. Notably, GA methods do not automatically include diffraction, meaning it must be added in as an additional process where required. This contrasts with ‘low frequency’ solvers, such as the Finite Element Method (FEM) or Boundary Element Method (BEM), which directly include diffraction as part of how they solve the underlying Partial Differential Equation (PDE) that represents acoustic wave propagation within the domain under study. ISM constructs a solution from simple sources that individually satisfy the PDE, but the way it assembles them is only valid for rectangular (2D) or cuboid (3D) domains. Raytracing is related to the acoustic PDE via the eikonal approximation.

 

ISM is effective for reflections from smooth planar boundaries. Factors that constrain the accuracy include the inability of the standard algorithm to include boundary scattering or diffraction (the latter being zero in a rectangular or cuboid domain). These must be added on as additional processes [3,4]. But it’s main constraining factor is its exponential computational cost scaling with reflection order. This limits it to computing low-order early reflections only in most cases.

 

Raytracing is popular for late-time impulse response generation because its computational cost is linear with reflection order. However, all but the lowest-fidelity models of boundary scattering and diffraction lead to ray-splitting, and this cripples it with exponential cost scaling with reflection order. Recent work has used precomputed edge-to-edge visibility graphs to accelerate this, but this was only applied to ‘direct’ diffraction and the algorithm retains exponential complexity [5].

 

A final constraint is that the above methods bear little similarity to established methods for low frequencies (FEM and/or BEM), hindering prospects for a unified full-audible-bandwidth approach.

 

1.1 Surface-Based GA

 

An emerging category of models is what Savioja and Svensson termed “surface-based” GA [2]. These discretise the distribution of acoustic power flow versus position using a boundary mesh, just like BEM uses for low frequency problems. Indeed, it has been called “Energy-Intensity” BEM by some authors [6,7]. The oldest algorithm is this category is Acoustic Radiosity, developed by Kuttruff in the 1970s [8,9]. This does not consider the directional nature of sound arriving at the boundary and instead considers these to all scatter diffusely according to Lambert’s law, which is not a realistic model of real rooms. Results tend not to be too unreasonable for late reflections, where the sound field is diffuse anyway, but it is inaccurate for early reflections, which are sparse in direction. Usefully, however, it can be used to find T60 in rooms with non-uniform absorption [8].

 

More sophisticated algorithms also parameterise reflection direction [6,7,10–13]. This approach appears to have been developed independently by three different groups, but they can be shown to all be equivalent once choices in dimensionality, interpolation and notation are taken into account. One of the contributions of this paper is a framework that makes this equivalence clearer. A version of the algorithm that operates on volumetric meshes has also been demonstrated [14] but, while this opens up exciting applications in inhomogeneous and structure-borne sound problems, its benefits for homogeneous room acoustics problems are less clear.

 

A characteristic of the surface-based GA method is that the way the solution is interpolated – mathematicians would call this the ‘approximation space’ – is defined in advance. Spatially this is based on elements on a mesh, for which all the standard choices of interpolation (a.k.a. ‘shape’ or ‘basis’) function are possible. Common choices are piecewise-constant (PWC), piecewise-linear or piecewise-quadratic, where the latter two can be continuous or discontinuous but PWC is always discontinuous. To capture the distribution with respect to (w.r.t.) reflection angle, a set of PWC functions on non-overlapping angle ranges can be used – ref. [13] did this in 3D – or orthogonal polynomials such as Legendre or Chebyshev [10,11]. For the latter, it is more convenient to use an edge-tangential coordinate rather than angle, since this lies in the range -1 to +1. Franzoni and co workers used a half Fourier series in angle in 2001 and a full version in 2005, leading to problems of non-orthogonality and periodicity respectively. A cosine series in angle would have eradicated both these problems and is equivalent to Chebyshev polynomials in the tangential coordinate.

 

A key benefit of these algorithms is that they maintain a fixed number of degrees of freedom (DOF) irrespective of reflection order. Reflections are advanced by a matrix multiplication; hence computational cost is linear with reflection order. Crucially, this holds even if the raytracing process that occurs within a reflection order involves a lot of ray-splitting, as is required to model scattering or diffraction accurately [2,15]. The extra rays are absorbed into the predefined approximation space and the number of degrees of freedom is kept constant and computational cost linear. Other benefits are that steady-state can be easily solved for [10,11] and decay rates could be found too [8].

 

The discussion above of “approximation” makes it clear that this is present in the method. Repeated discretisation of high-order reflections will not give results identical to those from tracing the rays individually. This was viewed as a limitation by Pohl and Stephenson [15], who proposed a similar idea with their ‘ray recombination’ algorithm, but we do not see it that way. All raytracing algorithms acknowledge that individual rays do not describe sound intensity. Instead, it is their spatial distribution that captures this, hence receivers in classical raytracing must have volume. Nonetheless, the tracing of rays to high order makes it tempting to believe that each ray’s path over many reflection orders matters but it really does not – it is just one sample within a distribution. This is especially true when boundary scattering is present. In this case reflection density increases making individual reflections no longer important [2]. Surface-based GA exploits this to simplify the algorithm and provide computational benefits, plus the approximation is controllable. This work shows that early reflections can be captured accurately if the number of DOF is sufficiently high.

 

Some articles have suggested that surface-based GA is viewed as computationally expensive [9], but this view appears misplaced based on the discussion above and is not backed up by the all important computational cost scaling. What is true is that much of the computational effort is in the problem setup, so must be performed beforehand. We see such precomputation as an advantage. As Nosal et al [9] point out, listener position can be quickly recomputed once the boundary solution has been computed. Additionally, source position can be changed quite readily since the boundary to-boundary interaction matrices – computation of which are the main bottleneck – do not depend on source position. Identical arguments apply to standard BEM, and the algorithm presented herein matches exactly the framework in fig. 1 of ref. [16]. The similarity of surface-based GA to BEM is discussed further in this papers sister paper, #1005 in these proceedings [17].

 

1.2 Overview of this Paper

 

The main contribution of this paper is to cast surface-based GA into a Galerkin BEM framework. This expediates comparison of different basis and test functions since it moves factors required to compensate for these out of the main code body and into a matrix operation. It is also motivated by a desire to emphasise the synergy between surface-based GA and high-frequency BEM. The algorithm and matrix structure used herein matches exactly that of the high-frequency BEM in ref. [18]. This similarity is discussed further in paper #1005 in these proceedings [17].

 

Source and receiver directivity are also implemented using Fourier series of arbitrary order. Stacking of their coefficients into vectors means that all acoustic wave propagation operations are represented by a matrix multiplication, be it source to boundary, boundary to boundary, boundary to receiver, or direct source to receiver. This mimics the structure in fig. 1 of ref. [16].

 

Results are presented for early reflections and validated against an ISM that has been modified to use the same source and receiver directivity model. Agreement is seen to be very good, though differences exist between different discretisation schemes.

 

2. THEORY AND FORMULATION

 

This section presents the background theory behind the method. There is insufficient space for a full exposition, but the key concepts and features will be presented and dicsussed.

 

2.1 The Assumed Model of Acoustic Propagation

 

As with all GA, it is assumed that we are working in a high-frequency regime. The first starting assumption is that acoustic propagation can be represented by particles travelling in straight lines termed rays. An individual particle would have position 𝐱 and propagation direction 𝐤̂. The power distribution 𝑤(𝐱, 𝐤̂) of these particles is therefore a function of both position and direction.

 

The second starting assumption is that we will not consider phase – we only consider energy and power, albeit including its directional characteristics. This is also true for classical raytracing. It is justified by arguments of the chaotic nature of late-time sound fields and ensemble averaging [19].

 

The model will be formulated in two dimensions (2D) for simplicity. Within the domain Ω, position 𝐱 can vary in 2D and the propagation direction vector 𝐤̂ can point toward any angle within a 360° range. It will therefore be convenient to define 𝐤̂ = [cos 𝜃  sin 𝜃], where 𝜃 is a global angle relative to the x-axis of the global Cartesian coordinate system. But when 𝐱 lies on the boundary Γ the situation is slightly different. Now 𝐱 is only a function of boundary position, a 1D quantity (the boundary of a 2D domain is a 1D line) such as length 𝑙 in meters along the boundary from some

 

 

Figure 1: Geometry for source and boundary definitions: a) Incoming and outgoing directions: b) Power radiating from a point source; c) Geometry for point source to boundary mapping reference point. The orientation of the boundary Γ at 𝐱 is given by its normal vector 𝐧̂ and tangent vector 𝐭̂, which are perpendicular. 𝐧̂ is taken to point into the acoustic domain Ω.

 

It will be convenient to distinguish between incoming directions, for which 𝐧̂ ∙ 𝐤̂ < 0, and outgoing directions, for which 𝐧̂ ∙ 𝐤̂ > 0. These will be termed 𝐤̂⁻ and 𝐤̂⁺ respectively, and each occupies a 180° range. 𝐤̂^± are defined in the local coordinate system defined by 𝐧̂ and 𝐭̂ and can be parameterised by the angle 𝜃 made with the normal vector 𝐧̂. This is shown in Figure 1a.

 

Another option, which is closer to what is done is refs. [10,11], is to use a boundary tangential coordinate 𝜒=sin𝜃 to represent direction. This is bounded by |𝜒|  ≤ 1. 𝐤̂^± are then given by:

 

 

In this paper, 𝐤̂, 𝜃 and 𝜒 will be used interchangeably based on these definitions. Typically, 𝜃 is used when considering if propagation angles match, 𝜒 is used when evaluating angular integrals numerically, and 𝐤̂ is used to ease definitions by avoiding the need to align local coordinates.

 

2.2 Power Incident from a Point Source

 

Consider power radiating from a point source at position 𝐲₀, as depicted in Figure 1b. The source injects power 𝑊 into the domain. The units of 𝑊 is Watts/metre because the model is in 2D, and a 2D point source is equivalent to a 3D line source.

 

Consider now a receiver at position 𝐱. It receives a total power density of 𝑊 / 2𝜋𝑅 due to cylindrical spreading. Here 𝐑 = 𝐱 − 𝐲₀ and 𝑅 = |𝐑|. The unit of this density is Watts per metre², where this metre² is area tangential to the wavefront. The intensity at 𝐱 is this power density times the propagation direction = 𝐑 / 𝑅, so 𝐈 = 𝑊 / 2𝜋𝑅.

 

Now additionally consider that the source may be directional, so the power injected varies with direction 𝜃_s (which is defined in the global sense, i.e., relative to the Cartesian axes, because 𝐲₀ is in the domain Ω). The power injected is replaced by a directivity function 𝑊_s(𝜃_s). This is defined so that it’s integral w.r.t. 𝜃_s from 0 to 2𝜋 equals 𝑊. It follows that 𝑊_s(𝜃_s), therefore, has units Watts per metre per radian. For example, an omnidirectional source has 𝑊_s(𝜃_s) = 𝑊 / 2𝜋.

 

The power density at 𝐱 is now 𝑊_s(𝜃_s) / 𝑅. The 2𝜋 has been removed from this statement because of the units of 𝑊_s(𝜃_s). Intensity is 𝐈 = 𝑊_s(𝜃_s) / 𝑅. Note that a real line source would have nearfield effects, but these have been ignored because one of the starting assumptions is that the algorithm is being derived in the high frequency regime.

 

The directional power distribution 𝑤(𝐱, 𝐤̂) ≡ 𝑤(𝐱, 𝜃_r) is the acoustic quantity being propagated, so it must be defined for power incident from a point source. 𝜃_r is the evaluation angle at the receiver 𝐱 and 𝑤 has units Watts per metre² per radian; here the “per radian” is w.r.t. receiving angle 𝜃_r. The rays emanating from this source in the high-frequency regime are perfectly directional, so are represented by a Dirac delta function in angle 𝛿(𝜃_s − 𝜃_r). Hence, 𝑤 is given by:

 

The total power arriving at the receiver irrespective of direction is:

 

 

Here the ‘sifting’ property of the delta function has been exploited to remove the integral and set 𝜃_r = 𝜃_s. 𝑊^r is spatial power density with units Watts/metre², like intensity has – the ‘per radian’ has been removed by angular integration. Note that for an omnidirectional source, with 𝑊_s(𝜃_s) = 𝑊 / 2𝜋, this simplifies to 𝑊 / 2𝜋𝑅 as expected. A similar integral can be used to compute intensity.

 

2.3 Point Source to Boundary Mapping

 

The objective here is to map the directional power distribution 𝑤(𝐱, 𝐤̂) emanating from the source to an equivalent incoming distribution 𝑤⁻(𝐱, 𝐤̂⁻) based in the boundary coordinate system.

 

The source is radiating rays in all directions, and the boundary can accept rays from any direction within the domain, but only ones that align in both receiving position and angle are geometrically valid. This creates an equivalence between 𝜃_s and 𝜃_𝐱 (see Figure 1c), but this is not immediately straightforward to state because they are defined relative to different reference directions. But it can be expressed easily by requiring that cos⁻¹(𝐤̂⁻ ∙ )=0. This is satisfied by 𝐤̂⁻ = 𝐑̂, hence a boundary distribution including an angular term 𝛿( cos⁻¹(𝐤̂⁻∙𝐑̂)) will sift out a ray direction (blue arrow) matching the arrival direction from the source (orange arrow).

 

The other aspect that must accounted for is whether the incident ray arrives from the correct side of the boundary and whether the ray path is occluded. Such so-called “visibility” and “validity” checks [2] are encapsulated into a function 𝑉(𝐱,𝐲₀), which takes the value 0 or 1. However, this is always 1 in a convex domain, such as the rectangular room considered herein, so the detail of this is considered here no further. The incoming boundary power distribution is therefore:
 

 

2.4 Re-radiation from the Boundary

 

The objective here is to determine how the energy distribution 𝑤⁺ (𝐲, 𝐤̂⁺) leaving the boundary at 𝐲 propagates through the domain Ω and then maps onto a distribution 𝑤(𝐱,𝐤̂) at a point receiver 𝐱 or an incoming distribution 𝑤⁻(𝐱, 𝐤̂⁻) on another boundary Γ_𝐱. Figure 2 shows the geometry. The ray (shown in orange) connects 𝐲 and 𝐱 and the direction of propagation aligns with it, hence 𝐤̂ += 𝐤̂ = 𝐤̂ −= . But how this is parameterised as an angle varies:

• For at 𝐲 ∈ Γ_y the angle 𝜃_𝐲 is relative to the normal vector 𝐧̂_𝐲 at 𝐲;
• For a point receiver 𝐱 ∈ Ω the angle 𝜃_r is relative to the global Cartesian x-axis;
• For 𝐱 ∈ Γ_𝐱 the angle 𝜃_𝐱 is relative to the inverted normal vector −𝐧̂_𝐱 at 𝐱.

 

In all three cases, the units of 𝑤 are Watts per metre squared per radian, where the “per metre squared” is taken in the plane perpendicular to the direction of propagation. The formulaiton will be developed for a domain receiver 𝑤(𝐱, 𝐤̂) but this also applies to a boundary point 𝑤⁻(𝐱, 𝐤̂⁻).

 

 

Figure 2: Geometry for evaluating boundary to boundary interactions

 

All points on Γ_𝐲 can radiate to 𝐱 and the power density received will be the summation of their contributions. Because of this, there is an integral over boundary position. The boundary can radiate in any angle too, hence we also have an integral over 𝜃_𝐲. 𝑤(𝐱, 𝐤̂) therefore is given by:

 

 

In this, 𝐱 and 𝐤̂ are arguments of 𝑤. This means it cannot be assumed that 𝐤̂ is aligned with , but the statement needs to ‘sift’ out values of 𝐤̂ that are. The 𝛿(cos⁻¹(𝐤̂ ∙ )) term acheives this. Similarly, 𝐤̂⁺ is a function of the integration variable 𝜃_𝐲, but the geometric model of acoustic propagation means that only radiation from the boundary in the direction 𝐤̂ += contributes. The term 𝛿(cos⁻¹(𝐤̂ +∙ )) sifts this out. The 1 / 𝑅 arises due to spherical spreading from the boundary sources. It converts the ‘per radian’ of 𝑤⁺(𝐲, 𝐤̂⁺) at 𝐲 into a ‘per metre’ tangential to the wavefront at 𝐱. One of the ‘per metre’ of 𝑤⁺(𝐲, 𝐤̂⁺) at 𝐲 will be integrated out by the spatial integral over Γ_𝐲. But this ‘per metre’ of 𝑤⁺(𝐲, 𝐤̂⁺) is in the plane perpendicular to the direction of propagation, whereas the integral is w.r.t. boundary length. The cos 𝜃_𝐲 term compensates for this. The second delta function can be exploited to analytically evaluate the integral w.r.t 𝜃_𝐲 via the sifting property:

 

 

This is the directional energy distribution arriving at a point receiver in the domain or at a point on a boundary. cos 𝜃_𝐲 = 𝐧̂_𝐲 ∙ has be substituted into it for ease of numerical evaluation.

 

Eq. (5) is equivalent to the Frobenius-Perron operator that the Nottingham group [10,11] use as their starting point. The only significant different is that they integrate (and later sift) in momentum (equivalent to position and χ together), whereas we have integrated w.r.t θ_𝐲, instead of χ_𝐲, and sift in θ_𝐲 and θ_𝐱. This choice means they omit the cos θ_𝐲 and 1 / R, but a later coordinate transform re-introduces these and a cos θ_𝐱 term. This makes our formulation equivalent, as is the one from Duke.

 

2.5 Discretisation

 

We need to discretise 𝑤(𝐱, 𝐤̂) on the boundary. Other papers have usually discretised only the incoming or the outgoing sound energy, and then included somehow the fact that 𝑤⁻ is reflected to give 𝑤⁺. We will discretise 𝑤⁻ and 𝑤⁺ separately, and then account for boundary scattering via multiplication with a ‘scattering’ matrix.

 

The boundary Γ is meshed using a closed loop of 𝑁 boundary elements. Position on this is indexed by a unitless position parameter 0 ≤ 𝜇 ≤ 𝑁 where 𝑁 is the total number of mesh elements. The boundary quantities 𝑤⁻ and 𝑤⁺ are discretised by a weighted sum of basis functions:

 

 

Here 𝐰⁻ and 𝐰⁺ are vectors of coefficients for the incoming and outgoing waves 𝑤⁻ and 𝑤⁺ respectively. The two sets of indices 𝑛 and 𝑞 have been collapsed into a single index 𝑛𝑄 + 𝑞 . 𝑏_𝑛(𝜇) is the spatial basis function. This is chosen to be piecewise-constant (PWC), so 1 on element Γ_𝑛 and zero otherwise. 𝑎_𝑞(𝜒) is the angular interpolation function. It might be an orthogonal polynomial, such as Legendre or Chebyshev, or PWC discretisation of −1 ≤ 𝜒 ≤1 into 𝑄 regions.

 

The source directivity 𝑊_s(𝜃_s) will be represented by a Fourier series in angle 𝜃_s:
 

 

Here 𝑄_s is the source order and the set of coefficients 𝑎_𝑞 and 𝑏_𝑞 are parameters supplied by the user. For examples, An omnidirectional source has 𝑄_s = 0 and 𝑎₀ = 𝑊 / 2𝜋, where 𝑊 is the total source power in Watts / metre. A Fourier series in angle 𝜃r will also be used to represent directional energy arriving at receivers. This is akin to 2D Ambisonics. 𝑤(𝐱, 𝐤̂) will be represented by:

 

 

Here 𝑃_r is the receiver order. The coefficients 𝑎_𝑝 and 𝑏_𝑝 are output data computed by the algorithm by finding 𝑤 at 𝐱 and then mapping it onto the representation in eq. (9). This mirrors the framework used in ref. [16] for auralisation of BEM.

 

2.6 Galerkin BEM

 

There is insufficient space here to fully explain the details of the Galerkin BEM formulation. It involves a second pair of integrals over position on the boundary and receiving angle that must be evaluated, though another of these can be sifted out leaving just two boundary integrals to evaluate. These integrals include another pair of basis functions, sometimes called the ‘testing functions’ [18]. These allow choice in how error in the solution is minimised. Choosing them to match the basis functions in eq. (7) forms a ‘Bubnov–Galerkin method’, which minimises mean-squared error in the same approximation space as the solution. Choosing them to be different forms a ‘Petrov– Galerkin method’, which minimises mean-squared error in some other space. Rowbottom and Chappell recently opted to do the latter [20], using delta functions as basis in angle (so a finite set of angles instead of a distribution) and PWC testing functions. This produces an effect very similar to the ‘ray direction quantisation’ discussed by Pohl and Stephenson [15]. We opt to do the former for reasons of simplicity and symmetry, and synergy with ref. [18].

 

A feature of Galerkin BEM is that it does not compute the coefficients of eq. (7) directly. Instead, it computes a vector of ‘projections’, from which the coefficients can be found by solving a matrix equation. The matrix involved in this equation is an identity matrix of the approximation / testing space, and is commonly called the ‘mass matrix’, due to its similarity with the mass matrix in FEM. For an orthogonal discretisation and testing scheme, such as arises from pairs of PWC or Legendre functions, it is diagonal so has 𝑁𝑄 non-zeroes. For PWC discretisation in space and angle its entries are 𝑙_𝑚 × 2 / 𝑄, where 𝑙_𝑚 is the length of the m^th element. For PWC discretisation in space and Legendre in angle, its entries are 𝑙_𝑚 × 2 / (2𝑝+1). Chebyshev functions are not orthogonal with the chosen weighting in 𝜒, making their mass matrix block diagonal with 𝑁𝑄² non-zeroes. All mass matrices were therefore seen to be highly sparse and readily invertible. The ones encountered during this study all had condition numbers of 40 or below.

 

2.7 Algorithm Structure

 

The process for solving the discretised version of the scheme above is as follows:

1. First the directional power from the source, represented by 𝐰_s, is mapped onto the boundary as an incoming power distribution 𝑤⁻(𝐱, 𝐤̂⁻), represented by 𝐰₀⁻
2. This is reflected in the boundary, including scattering and absorption as appropriate, to give an outgoing power distribution 𝑤⁺(𝐱, 𝐤̂⁺), represented by 𝐰₀⁺
3. The outgoing power distribution 𝑤⁺(𝐱, 𝐤̂⁺) is propagated across the domain to the boundary, where it is mapped onto an incoming power distribution 𝑤⁻(𝐱, 𝐤̂⁻), represented by 𝐰₁⁻
4. Either steps 2 and 3 are repeated up the required number of reflection orders, and 𝐰𝑖± are summed with respect to reflection order 𝑖 to obtain the total boundary power distributions, or a Neumann series is used to solve for total steady state power distribution in one step [10].
5. Power distributions at receivers are computed, represented by 𝐰_r. This comprises:

   a) Computation of power direct from sources (represented by their 𝐰^s)

   b) Computation of power from the boundary (represented by Σ_𝑖 + 𝐰_𝑖⁺)

 

Step 5 can also be applied to individual reflection orders 𝐰_𝑖⁺, allowing the power density received from those to be visualised if that is of interest (e.g. for diagnostic purposes). This is done below.

 

The process above can be summarise by the following pair of matrix equations:
 

 

Here 𝐓_sb is a matrix describing source-to-boundary propagation and 𝐓_bb describes boundary-to-boundary propagation. 𝐌 is the mass matrix. 𝐑 is a reflectance matrix that describes boundary absorption and scattering. Boundaries are rigid in the examples herein, so 𝐑 = 𝐌.

 

Total power distribution at a receiver is then given by:

 

Here 𝐓_sr is a matrix describing direct source-to-receiver propagation and 𝐓_bb describes boundary to-receiver propagation. This form directly matches that used for BEM in eq. 32 and 41 of ref. [18].

 

3. RESULTS

 

A rectangular room was simulated. It had dimensions 4m x 3m and was centred on the origin of the coordinate system. The source was located at [-1, -0.5] and the receiver at [1, 1]. The setup is shown in Figure 3b. A 3rd order source directivity was defined with 𝑎₀ = 𝑎₁ = 𝑎₂ = 1, 𝑎₃ = 1/2, and all 𝑏_𝑞 = 0. This gives 𝑊_s(𝜃_s) (via eq. (8)) as shown in Figure 3a. The receiver is order 5. Direct propagation to the receiver gives the power distribution shown in Figure 3c. Power received directly from the source was 0.762 W/m and main lobe points in the direction of the source, as expected.

 

 

Figure 3: a) Source directivity polar plot, b) Room geometry used with source and receiver locations, c) Received power directivity for direct sound

 

For the initial results, a maximum element size of 0.25m was used, giving 𝑁 = 56 elements, and 𝑄 = 17 was used to discretise angle. This gave the system 952 degrees of freedom. Boundary integrals were evaluated by Gauss-Legendre quadrature and the 1 / 𝑅 singularity in eq. (6), which affects adjacent but non-coplanar elements, was regularised using the Sato transform [21].

 

Validation was against a cuboid room ISM code modified from the standard approach to include source and receiver directivity. This uses the same matrix mapping validated in Figure 3c but with image source directivity ‘flipping’ to respect reflection of sources in walls. A rigid boundary condition was assumed for simplicity. This, however, forces validation to be performed over finite reflection orders since neither algorithm converges to a finite steady state when energy is trapped.

 

3.1 Validation via Power Distribution on the Boundary

 

To validate the codebase and provide insight into the boundary power distribution 𝑤(𝜇, 𝜒) and choice of basis functions, plots of it are presented below. Here 𝑤(𝜇, 𝜒) is shorthand for the left-hand side of eq. (7) and no discrimination is necessary between incoming and outgoing waves because the boundary is planar and rigid, meaning no absorption or scattering occur.

 

Figure 4b gives an example of such a plot for PWC basis in both position and angle and Figure 4a illustrates its meaning. Figure 4b pertains to the lower boundary depicted in Figure 4a. This has 16 elements, hence the horizontal coordinate 𝜇 in Figure 4b runs from 0 to 16. 𝜒 parameterises reflection angle 𝜃 via 𝜒=sinθ and 𝜃 is depicted in Figure 4a. This has the range −^𝜋/2 ≤ 𝜃 ≤ 𝜋/2  so −1 ≤ 𝜒 ≤ 1. In ISM it is common to think only of the ray direct from the source to the receiver. This is shown towards the right of Figure 4a. It’s position and angle in boundary reflection space ( 𝜇, 𝜒 ) are indicated by the yellow cross on Figure 4b. But the source doesn’t radiate only to the receiver – it radiates to all angles, hitting all positions on the boundary. This leads to cases where 𝜃 may be negative (green cross) or zero (blue cross) – these are also overlaid Figure 4b. Evaluating arrival angle for all positions on the edge produces the red line, onto which all three crosses fall. This is the loci of the source arrival direction. It can be seen that the PWC discretisation scheme in position and angle causes a 2D ‘faceting’ of the distribution into small rectangles in (𝜇, 𝜒) space. The fading out around 𝜒 ≈ 0 is due to the choice of source directivity in Figure 3a.

 

This clearly shows approximation. The fine loci of the ISM method has been coarsely pixelized and its energy smudged across these. But this is not inappropriate for late-time responses, where the sound field becomes diffuse and individual reflections can no longer be distinguished. Here it illustrates that smaller element and higher 𝑄 cause less approximation since the ‘pixels’ are smaller.

 

Figure 4: Reflection of a source from a boundary: a) layout of a source, receiver and connecting rays in room; b) projection onto boundary (𝜇, 𝜒) space space with PWC basis funs.

 

Figure 5 shows the same plot as Figure 4b but for all four edges of the room and for PWC and Legendre angular basis functions. Using PWC functions makes the solution on the boundary discontinuous, with data quantised into blocks. But Legendre polynomials are continuous over the whole range −1 ≤ 𝜒 ≤ 1, yielding a series of smooth stripes. These exhibit some side-lobing in 𝜒, as is typical of orthogonal polynomials, but the peak in the distribution is clearly where it is expected to be. Peak values are different due to the differing ‘smudging’ in the two schemes.

 

Figure 5: Mapping of reflection angles onto the four edges of the rectangular room using (a-d) Piecewise-Constant functions, and (e-h) Legendre polynomials, as angular basis functions

 

3.2 Validation via Power Directivity at Receiver

 

Figure 6 shows the receiver directivity patterns for first and second order reflections. Figure 6a only involves the source-to-boundary and boundary-to-receiver propagation processes 𝐓_sb and 𝐓_br, but Figure 6b involves the boundary-to-boundary propagation process 𝐓_bb too, hence all three are validated. All three discretisation schemes match the ISM response well, Legendre and Chebyshev especially so. For these the error in total power (zeroth order directivity) was 0.1% for first-order reflections and 4.4% for second order. For PWC it was higher at 4.5% and 18.8% respectively.

 

That the Legendre and Chebyshev schemes gave identical results (to reported precision) was unexpected since ref. [11] reported that Chebyshev basis performed less well. To validate that this was not a coding error, the structure of the Chebyshev mass matrix was checked. Its version of the plots in Figure 5 (not shown) also showed slightly different sidelobes, proving it is different.

 

 

Figure 6: Received power directivity for: a) first order reflection, b) second order reflection

 

3.3 Convergence Study

 

In this section, element size and 𝑄 were varied to examine how the various schemes converge to the ISM result. The error metrics were normalised error versus ISM in total receiver power, as reported above, and a similar normalised RMS error over all the receiver coefficients. Thus, error is reported for receiver order 0 and receiver order 5.

 

Figure 7 shows results with 𝑄 fixed at 20 and element size varying from 0.1m to 1m. The Legendre and Chebyshev schemes again exhibit identical error trends. Error for PWC is again higher. Element size appears to have a fairly weak effect on accuracy, though convergence is visible in the RMS (directional to order 5) results for the Legendre and Chebyshev schemes. The PWC scheme shows little evidence of convergence in any case.

 

Figure 7: Normalised error (left) and RMS error (right) for vari ed element size and fixed 𝑄 = 20

 

Figure 8 kept the maximum element size constant at 0.2m and varied 𝑄 between 1 and 20. The Legendre and Chebyshev schemes show a clear convergence trend for the RMS (directional) results, including at second order. A weaker trend is also seen for total power. The PWC scheme, however, again shows little evidence on convergence in any case.

 

Figure 8: Normalised error (left) and RMS error (right) for var ied 𝑄 and fixed elem size 0.2m
 

One reason for poor convergence of the PWC scheme might be its discontinuous nature. This affects its convergence in standard BEM. But accuracy could also be compromised by application of Gauss-Legendre integration to evaluate matrix entries for these discontinuous functions. The boundary basis functions are smooth over each element, but the angular kernel is not, which could compromise integration accuracy. Refs. [9,19] addressed this in 3D using spherical triangles.

 

4. CONCLUSIONS

 

A surface-based geometrical acoustics raytracing algorithm has been presented in Galerkin boundary integral form. The inclusion of the ‘mass’ matrix was seen to simplify comparison between different interpolation schemes. Moreover, the framework and notation emphasise synergy with a high-frequency BEM, which we see as being a route towards a unified algorithm for simulation of the full audible frequency band. The scheme was validated against the image source method for early reflections in a rectangular room with rigid boundary conditions.

 

Future work is to add support for absorbing and scattering boundaries and study convergence w.r.t. reflection order. Other basis functions, e.g., continuous piecewise linear, could be studied too. Later, it would be valuable to add corner diffraction to the scheme and investigate T60 decay rates.

 

5. REFERENCES

 

1. M. Vorländer, Auralization : Fundamentals of Acoustics, Modelling, Simulation, Algorithms and Acoustic Virtual Reality, 2nd ed., Springer, Cham, Switzerland, 2020.
2. L. Savioja, U.P. Svensson, Overview of geometrical room acoustic modeling techniques, J. Acoust. Soc. Am. 138 (2015) 708–730.
3. P.T. Calamia, B.E. Markham, U.P. Svensson, Diffraction Culling for Virtual-Acoustic Simulations, Acta Acust. United With Acust. 94 (2008) 907–920.
4. A. Erraji, J. Stienen, M. Vorländer, The image edge model, Acta Acust. 5 (2021) 17.
5. C. Schissler, G. Mückl, P. Calamia, Fast diffraction pathfinding for dynamic sound propagation, ACM Trans. Graph. 40 (2021) 1–13.
6. L.P. Franzoni, D.B. Bliss, J.W. Rouse, An acoustic boundary element method based on energy and intensity variables for prediction of high-frequency broadband sound fields, J. Acoust. Soc. Am. 110 (2001) 3071.
7. J.W. Rouse, L.P. Franzoni, Improvement of a High-Frequency Broadband Energy-Intensity Boundary Element Method to Include High Resolution Specular Reflection, J. Comput. Acoust. 13 (2005) 99–125.
8. H. Kuttruff, A simple iteration scheme for the computation of decay constants in enclosures with diffusely reflecting boundaries, J. Acoust. Soc. Am. 98 (1995) 288–293.
9. E.-M. Nosal, M. Hodgson, I. Ashdown, Improved algorithms and methods for room sound field prediction by acoustical radiosity in arbitrary polyhedral rooms, J. Acoust. Soc. Am. 116 (2004) 970–980.
10. D.J. Chappell, G. Tanner, S. Giani, Boundary element dynamical energy analysis: A versatile method for solving two or three dimensional wave problems in the high frequency limit, J. Comput. Phys. 231 (2012) 6181–6191.
11. D.J. Chappell, G. Tanner, Solving the stationary Liouville equation via a boundary element method, J. Comput. Phys. 234 (2013) 487–498.
12. S. Siltanen, T. Lokki, S. Kiminki, L. Savioja, The room acoustic rendering equation, J. Acoust. Soc. Am. 122 (2007) 1624–1635.
13. S. Siltanen, T. Lokki, L. Savioja, Room Acoustics Modeling with Acoustic Radiance Transfer, Proc. Int. Symp. Room Acoust. 2010. (2010) 1–6.
14. J. Bajars, D.J. Chappell, N. Søndergaard, G. Tanner, Transport of phase space densities through tetrahedral meshes using discrete flow mapping, J. Comp. Phys. 328 (2017) 95–108.
15. A. Pohl, U.M. Stephenson, A Combination of the Sound Particle Simulation Method and the Radiosity Method, Build. Acoust. 18 (2011) 97–122.
16. J.A. Hargreaves, L.R. Rendell, Y.W. Lam, A framework for auralization of boundary element method simulations including source and receiver directivity, J. Acoust. Soc. Am. 145 (2019) 2625–2637.
17. J.A. Hargreaves, A comparison between the high-frequency Boundary Element Method and Surface-Based Geometrical Acoustics, in: Proc. Internoise, Glasgow, 2022.
18. J.A. Hargreaves, Y.W. Lam, The Wave-Matching Boundary Integral Equation — An energy approach to Galerkin BEM for acoustic wave propagation problems, Wave Motion. 87 (2018) 4–36.
19. J. Bajars, D. Chappell, N. Søndergaard, G. Tanner, Computing High-Frequency Wave Energy Distributions in two and three dimensions using Discrete Flow Mapping, in: Proc. 22nd Int. Congr. Sound Vib., Florence, 2015.
20. J. Rowbottom, D.J. Chappell, On hybrid convolution quadrature approaches for modeling time-domain wave problems with broadband frequency content, Int. J. Numer. Methods Eng. 122 (2021) 7581–7608.
21. M. Sato, S. Yoshiyoka, K. Tsukui, R. Yuuki, Accurate numerical integration of singular kernels in the two- dimensional boundary element method, in: C. Brebbia (Ed.), Bound. Elem. X, Vol. 1, Springer, Berlin, 1988: pp. 279–296.

 

---

¹ a.emthyas1@edu.salford.ac.uk

² j.a.hargreaves@salford.ac.uk