A comparison between the high-frequency Boundary Element Method and Surface-Based Geometrical Acoustics

Jonathan A. Hargreaves 1 Acoustics Research Group Newton Building, University of Salford, Salford, M5 4WT, United Kingdom

ABSTRACT The audible frequency range covers many octaves in which the wavelength changes from being large with respect to dominant features of a space to being comparatively much smaller. This makes numerical prediction of a space’s acoustic response, e.g. for auralisation, extremely challenging if all frequencies are to be represented accurately. Different classes of algorithm give the best balance of accuracy to computational cost in different frequency bands – ‘wave solvers’ such as Boundary Element Method (BEM) at low frequencies and Geometrical Acoustics (GA) methods at high frequencies. But combining their output data can be an awkward process due to their very different formulations. This is particularly important for early reflections, which give crucial spatial perceptual cues. Hence there is a need for a unified full audible bandwidth algorithm for early reflections. This paper will describe ongoing research to develop such an algorithm by exploiting synergies between high-frequency BEM and GA. It will describe how appropriately chosen oscillatory basis functions in BEM can produce leading-order GA behaviour at high frequencies and explore how interactions between these compare to the same interactions arising in a surface-based Geometrical Acoustics scheme.

1. INTRODUCTION

Since its inception 25 years ago, auralisation [1] has become an important tool for acoustic engineers to communicate the sonic benefits of designs to stakeholders. This is especially true in architectural and automotive applications. In order for an auralisation to be accurate, however, the numerical simulation on which it is based, and the spatial audio encoding and rendering processes used to present it to the listener, must all be accurate too. This places challenging requirements on simulation algorithms. Presently it is commonplace to run different algorithms in different frequency bands, but this raises questions of how their output data should be combined, since the listener needs to hear all frequencies together. Postprocessing techniques have been proposed to combine impulse response data from different algorithms [2], but these are not perfect and do not address the fundamental dissimilarities between the algorithms. An energetic raytracing method has almost nothing in common with the grid-based pressure sampling used in Finite-Difference-Time- Domain (FDTD), for example. This is particularly important for early reflections, which give crucial spatial perceptual cues – for late time the wave field becomes chaotic at high frequencies and the benefits are less clear. There is, therefore, a need for a new framework with a unified full audible bandwidth algorithm for computing early reflections.

This paper describes ongoing research to develop such an algorithm by exploiting synergies between Boundary Element Method (BEM) and Geometric Acoustics (GA). It will describe how

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

appropriately chosen oscillatory basis functions in BEM can produce leading-order GA behaviour at high frequencies, and how these might be assembled into a wave solver. It will introduce elements of the ‘Wave-Matching’ BEM, a new High-Frequency BEM (HF-BEM) algorithm and it’s similarity to what Svensson & Savioja term “Surface-based” GA (SB-GA) [3]. Variants of SB-GA that consider wave direction are the most sophisticated, and this has been suggested independently by several groups [4–7]. Notably, a Galerkin BEM implementation of SB-GA is published in these proceedings [8] which emphasises this similarity. Both formulations can be solved by marching-on- in-reflection order, a property that makes them suitable for modelling early reflections.

Section 2 outlines a vision for how these algorithms could operate together. Section 3 addresses some initial issues of superficial dissimilarity and discusses how the total solution is represented as a sum of reflections. Sections 4 and 5 then look in more detail at the behaviour of the propagators in HF-BEM and SB-GA, and the matrices they produce. Finally, conclusions are drawn in section 6.

2. VISION

Room acoustic simulation is currently dominated by Geometrical Acoustic (GA) solvers. These models are efficient and have been widely used for approaching 30 years [3], but it is also widely known that they are not accurate in all scenarios. This is usually at lower frequencies or in smaller rooms, where the modal density remains low up to a higher frequency.

‘Wave-based’ methods, such as Finite Element Method (FEM), FDTD or BEM, are extremely accurate [9] but their computational cost scales badly with frequency, meaning there is an upper limit to their practical use. Moreover, it is well known that at late time, the sound field in rooms becomes diffuse and chaotic at high frequencies. There is, therefore, little to be gained from running a computationally expensive wave method in this region; an energy-based algorithm that gives an ensemble average of likely responses is more appropriate [7].

This leads to different algorithms being used for different frequency bands. This is depicted in Figure 1, with the status quo shown on the left. A distinction is made in the high-frequency GA regime between early reflections (yellow region), which are sparse in time and perceptually important, and late reflections (blue region), which are dense in time. Commonly, a deterministic algorithm such as Image Source Method (ISM) is employed for early reflections [3]. But this has a computational cost that rises rapidly with reflection order, so later time is modelled with a computationally cheaper method such as stochastic raytracing [10]. This is less accurate but is acceptable because only the energy decay envelope is important perceptually. A method is required to transition energy between these algorithms, represented by the curved arrow. Note that the raytracing is shown as if it occurs after the ISM, but it actually starts in parallel with it and energy is gradually transferred to it due to scattering. Then, at some pre-set transition order, the ISM method is terminated, and all energy is moved to the raytracing model. An algorithm to do this adaptively between ISM and Radiosity, a type of zeroth-order SB-GA, has also been proposed [11].

Figure 1: Illustrations of spectrograms of a Room Impulse Response, with regions of low and high reflection density (left to right) and modal density (bottom to top) indicated. a) How this is divided between simulation algorithms in the status quo; b) the vision explored herein.

The low modal density in the low frequency region (green) means it must be handled by a ‘wave’ solver, such as BEM. ‘Hybrid’ approaches [9,12–14] do what Figure 1a depicts, that is, run a wave solver for low frequencies and GA for higher frequencies. But this is while efficient, it is not as straightforward as it seems due to the requirement for a crossover between the impulse responses generated by two algorithms (red line). For late time, a conventional crossover filter is adequate, but for early time more care is required to merge the data correctly. This is because the reflection pattern is sparse in time and perceptually important. Aretz et al [2] suggested a non-linear method that avoids cancellations at the crossover frequency when the algorithms predict different phases for reflections, but in unpublished work by this author this has produced ‘pre-ringing’ artefacts.

The vision depicted in Figure 1b aims to addresses these issues. It features a unified full-audible bandwidth algorithm for early reflections, which feeds two late-time algorithms. Notably, this structure avoids the need to apply a crossover filter to early reflections. It is suggested that these should be BEM at low frequencies and SB-GA at high frequencies, which is essentially BEM with a ray propagator. For the unified early-time algorithm, HF-BEM is suggested. While ambitious, this is substantially more tractable than producing a single unified algorithm to cover the entire impulse response, because it is limited to sparse reflections. In particular, this demarcation would allow the efficiency gains seen in Hybrid Numerical Asymptotic BEM (HNA-BEM) [15] to be leveraged. For specific classes of scattering problem, HNA-BEM has been shown to reduce the number of degrees of freedom necessary so as to scale with log frequency [16], or even be frequency independent [17]. But such problems are characterized by a small number of dominant wave directions, and this only holds for rooms during the sparse early reflection regime. After that, energy must be transferred to other algorithms more suitable for dense wave fields, and SB-GA and standard BEM are hard to beat in this application. Most results in this paper are considered in the frequency domain, but all can be converted into time-domain equivalents using the Convolution Quadrature framework [18].

Two transfer processes (curved arrows) are now required. For HF-BEM to CQ-BEM this is straightforward. Both methods solve for pressure distributions satisfying the wave equation, so their results can simply be summed arithmetically. SB-GA differs because it is based on energy, so transfer will likely involve stochastic energy statements [19].

3. SOLUTIONS AS A SUM OF REFLECTIONS

BEM has been suggested for unification with GA because both are concerned with reflection from boundaries, which leads naturally to the synergies that are discussed herein. FEM and FDTD, in contrast, discretise the solution within the acoustic domain so are fundamentally different.

There are, however, some differences between BEM and GA, notably to do with how visibility and shadow zone concepts are handled. This section will first unpick these differences, arguing that they mainly arise due to differences in terminology, and then outlines how algorithms that solve reflection orders separately operate.

3.1 Visibility and Shadow Regions GA methods assume that sound propagates in straight lines like light. Indeed, the majority of GA algorithms originate from computer graphics. The notion of visibility is carried over from vision; if a sound path is occluded then it cannot be heard. This creates shadow zones . Related to this is the notion that reflections are only visible in certain regions.

Figure 2 presents a simple scattering problem to illustrate this. An incident wave arriving from the top left hits a polygonal obstacle, producing two reflected beams and a shadow zone behind. Pressure due to these phased beams can be written in the form 𝑃(𝐱) = 𝐴× 𝑉(𝐱) × e i𝑘𝐝̂∙𝐱 , where 𝐴 is amplitude, 𝑉(𝐱) is a visibility function that is either 0 or 1, 𝐝 ̂ is a direction unit vector and 𝑘= 2𝜋𝑓𝑐 0 ⁄ is wavenumber, where 𝑓 is frequency in Hz and 𝑐 0 is sound speed in metres per second. Figure 2b depicts the incident pressure field as it would be termed in a GA method. It’s direction 𝐝 ̂ is indicated by the arrow. The obstacle is deemed to occlude it; hence its visibility function is 1 in Ω + −Ω shad and 0 in Ω − ∪ Ω shad . Reflections arising from the illuminated faces produces the two reflections shown in (c). These again have their direction vectors depicted as arrows. Their visibility

Figure 2: Scattering from an obstacle decomposed into phased beams following the convention used

in GA (b & c) and the convention used in BEM (e & f). Arrows indicate propagation directions. a) geometry and denotations; b) incident pressure with shadow (GA); c) reflected pressure (GA); d) total pressure (same for both); e) incident pressure no shadow (BEM); f) scattered pressure (BEM).

functions are 1 in Ω top and Ω left respectively and 0 elsewhere. Total pressure is shown in (d) and is the sum of these. The checkerboard patterns there are due to interference.

BEM decomposes total pressure differently. It defines incident pressure to be ‘in the absence of the obstacle’, hence no shadowing is seen in (e). Instead, the scattered pressure it computes * , shown in (f), includes an antiphase beam that is visible in Ω − ∪ Ω shad , which cancels out the incident wave in the shadow region. Such behaviour could reasonably be termed a shadow wave . In the ‘Wave-Matching’ BEM [20], distinction between incoming and outgoing sound causes these to arise separately from reflections. A shadow wave cancels out the incoming wave, then a reflected wave is radiated with amplitude and phase dependent on the boundary condition. Total pressure (d) is again equal to the sum of (e) and (f), showing that these two decompositions are equivalent.

Figure 3: Diagram illustrating how a scattering problem is formulated in Galerkin BEM.

3.2 Solution by Marching on in Reflection Order GA algorithms operate by marching on in reflection order. ISM involves repeated re-reflection of sources in boundaries. Rays and beams reflect of boundaries, then propagate on and hit other boundaries, where they are reflected again and the process repeats.

BEM is not usually solved by marching on in reflection order, but it can be. Time-domain BEM [21] is often solved by marching on in time, a variant of the same process. Figure 3 illustrates the structure of this for Galerkin BEM. The core of the algorithm is the discretisation scheme, being the

* To be clear, subplots e and f of Figure 2 are not computed using BEM. They are geometric beams grouped in the manner that BEM decomposes a pressure field. In a BEM result, all three beams would see diffraction since they radiate from finite apertures, and it is diffraction of the shadow wave that smoothes the boundary of the shadow zone. Higher- order diffraction would occur too, which for this example has been shown to orbit around the obstacle [16].

boundary mesh and the interpolation (basis) functions defined upon it. This allows the pressure distribution, or energy for SB-GA, on the boundary to be represented by a finite set of coefficients. Mathematicians refer to this as the approximation space . These coefficients are found, either for the incident sound (left) or re-reflected sound (looped arrow at bottom) by the Testing Integral. This computes an inner product over the boundary between the incoming sound field and a set of testing functions (also called the dual-to-range in BEM literature [22]). These allow choice in how error in the solution is minimised. Choosing them to match the basis functions in the discretisation scheme 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 function space. We opt to do the former for reasons of simplicity, symmetry, and consistency with an energy interpretation presented in ref. [23]. Finally, the scattering integral computes how the pressure or energy on the boundary propagates to other locations, be them receivers in the acoustic domain Ω + or other parts of the boundary Γ . This integral is the well-known Kirchhoff-Helmholtz boundary integral equation that arises from application of Greens’ theorem to the acoustic wave equation [24,25].

A differentiating feature of SB-GA and HF-BEM, compared to standard BEM, is that wave propagation direction is expressly considered in the discretisation scheme. Incoming and outgoing directions 𝐤 ̂ − and 𝐤 ̂ + are parameterised by a unitless quantity 𝜒 according to:

𝐤 ̂ + (𝜒) = 𝐭 ̂ 𝑛 𝜒+ 𝐧̂ 𝑛 √1 −𝜒 2 , 𝐤 ̂ − (𝜒) = 𝐭 ̂ 𝑛 𝜒−𝐧̂ 𝑛 √1 −𝜒 2 .

(1)

Here 𝐭 ̂ 𝑛 and 𝐧̂ 𝑛 are unit vectors that are respectively tangential and normal to the n th element on the boundary. See fig. 1 in ref. [8] for a diagram. 𝜒 is related to the angle 𝜃 a ray makes with 𝐧̂ 𝑛 by 𝜒= sin 𝜃 . But while a geometrical view of wave propagation requires 𝜃 and 𝜒 to respectively be bounded by |𝜃| ≤ 𝜋 2 ⁄ and |𝜒| ≤1 , as is done in ref. [8], a full wave solution allows evanescent waves too. These are included by allowing |𝜒| > 1 , as is done in ref. [20], which makes the surface- normal components of 𝐤 ̂ − and 𝐤 ̂ + imaginary.

It is well known from beamforming theory that a finite aperture, be it for sensing or radiation, causes sidelobes, which equivalently are smudging in 𝜒 [26]. Geometrical methods ignore this and treat directionality as perfect (a Dirac delta function in angle) regardless of how small the radiating or receiving area is. A consequence of this is that position and direction can be treated separately in SB-GA whereas they are intertwined in HF-BEM.

The SB-GA discretisation equation from ref. [8] is:

𝑄

𝑁

𝑤 ± (𝐱(𝜇), 𝐤 ̂ ± (𝜒)) = ∑∑𝐰 ± [𝑛𝑄+ 𝑞]

𝑏 𝑛 (𝜇)𝑎 𝑞 (𝜒). (2)

𝑞=0

𝑛=0

Here 𝑁 is the number of spatial basis functions 𝑏 𝑛 , typically equal to the number of nodes or elements in the mesh, 𝑄 is the number of angular degrees of freedom, 𝑎 𝑞 (𝜒) are angular basis functions, and 𝐱 is position on the boundary, which is parameterised by a unitless coordinate 𝜇 . 𝑤 ± is the distribution of sound power, either incoming − or outgoing + , versus position 𝜇 and angle 𝜒 that is being solved for, and 𝐰 ± is the vector of coefficients used to represent it.

For comparison, the HF-BEM discretisation equation from ref. [20] is * :

𝑄

𝑁

𝑏 𝑛 (𝜇)e i𝑘𝐤̂ ± (𝜒 𝑞 )∙𝐱 (3)

𝑃 ± (𝒙(𝜇)) = ∑∑𝒘 ± [𝑛𝑄+ 𝑞]

𝑞=0

𝑛=0

* This equation is equivalent to the plane wave scheme in ref. [20]. The notation has been adapted to match ref. [8].

The main differences between eq. 2 and 3 are that:

i) Pressure (with phase) is being solved for (instead of power) and there is no explicit

dependence on direction. This will arise through beamforming. ii) A set of discrete directions 𝜒 𝑞 is used instead of angular distributions 𝑎 𝑞 (𝜒) . These will,

however, each lead to radiation over a distribution of angles due to the aforementioned smudging effect. The exact nature of this is dictated by the spatial Fourier transform of 𝑏 𝑛 . iii) The angular basis function 𝑎 𝑞 (𝜒) has been replaced by an oscillatory term that depends on

both 𝜒 and 𝐱 . It is this that (under integration) performs the beamforming. These differences will be explored further in later sections.

3.3 Algorithm Structure The algorithm for solving the SB-GA scheme in ref. [8] and the HF-BEM scheme in ref. [20] are extremely similar. All wave propagation operations are represented by matrix multiplications. A difference is that ref. [8] implements directional sources and receivers for SB-GA, while ref. [20] only implements omnidirectional sources and receivers, but the directional source framework demonstrated in ref. [9] could be readily added. For directional receivers in 2D, this would replace the orthogonality properties derived in ref. [24] with the equivalent 2D versions derived in section 4.2.2 of ref. [20]. If that were done, then both algorithms would take the following structure:

1) First the directional radiation from the source, represented by a vector of coefficients 𝐰 s , is

propagated to the boundary where it is mapped onto a distribution of incoming power (or pressure) for the zeroth reflection order, represented by 𝐰 0

− . 2) This is reflected in the boundary, including scattering and absorption as appropriate, to give

+ . 3) The outgoing distribution is propagated through the domain until it again reaches the

an outgoing distribution of power (or pressure) represented by 𝐰 0

boundary (looped arrow at the bottom of Figure 3 where it is mapped onto a distribution of incoming power (or pressure) for the first reflection order, represented by 𝐰 1

− .

± are summed with respect to reflection order 𝑖 to obtain the total boundary distributions, or a Neumann series is used to solve for the total steady state power in one step [5]. 5) Directional responses at receivers are computed, represented by 𝐰 r . This comprises:

4) Either steps 2 and 3 are repeated up the required number of reflection orders, and 𝐰 𝑖

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 in ref. [8] to validate the SB-GA algorithm for early reflections. The next section in this paper goes a step further and visualises pressure fields arising from individual coefficients in 𝐰 + , so their similarity with geometric beams can be discussed.

Figure 4: Diagram illustrating how a reflectance boundary condition operates.

The process above can be summarised by the following pair of matrix equations, which appear in both refs [8] and [20] (slightly differing notation in the latter):

− = { 𝐌 −1 𝐓 sb 𝐰 s for 𝑖= 0 𝐌 −1 𝐓 bb 𝐰 𝑖−1

+ = 𝐌 −1 𝐑𝐰 𝑖

− . (4)

𝐰 𝑖

+ for 𝑖> 0 , 𝐰 𝑖

Here 𝐓 sb is a matrix describing source-to-boundary propagation and 𝐓 bb describes boundary-to- boundary propagation. 𝐌 is called the mass matrix. It compensates for non-orthonormality in the testing integral inner product. It involves both the basis and the testing functions and is given this name due to its similarity with the mass matrix in FEM. 𝐑 is a reflectance matrix that describes boundary absorption and scattering. The purpose of this is illustrated in Figure 4.

The total directional distribution at a receiver is then given by:

+ 𝒊 . (5)

𝐰 r = 𝐓 sr 𝐰 s + 𝐓 br ∑𝐰 𝑖

Here 𝐓 sr is a matrix describing direct source-to-receiver propagation and 𝐓 br describes boundary- to-receiver propagation. 𝐰 r is the set of coefficients for a Fourier series for SB-GA in 2D, or a Fourier-Bessel series for 2D BEM.

Thus, it is seen that identical numerical frameworks arise from the two algorithms, a major step towards the vision in Figure 1b. The next sections will consider how the wave propagators compare.

4. DIRECTIONALITY OF RADIATION

The geometric phased beams shown in Figure 2 are plane waves times a visibility function. For the reflected and shadow beams, this visibility function is the projection of the extent of the boundary section that radiated them in their direction of propagation. The statement for radiation by a single SB-GA basis function reflects this. In the notation of ref. [8] this is:

𝑤 𝑛,𝑞 (𝐱, 𝐤 ̂ ) = ∫𝛿(cos −1 (𝐤 ̂ ∙𝐑 ̂ ))𝑉(𝐱, 𝐲)𝑏 𝑛 (𝜇 𝐲 )𝑎 𝑞 (𝜒 𝐲 ) 1

. (6)

𝑅 cos 𝜃 𝐲 𝑑Γ 𝐲 Γ

The left-hand side of this equation is a power density function, in position and angle, at a receiver point. Its first argument, 𝐱 , is the position of the receiver and the second argument, 𝐤 ̂ is the direction in which that the user wants to know the power flux. Subscripts 𝑛 and 𝑞 indicate that only the contribution from the basis function pair 𝑏 𝑛 × 𝑎 𝑞 is being considered here. This is one step removed from an elemental entry in 𝐓 br ; that additionally includes the receiver directivity model.

The integral on the right sums contributions to 𝐱 from all points 𝐲 on the boundary Γ . This can be thought of as a bundle of rays (orange) emanating from the boundary and converging on point 𝐱 , as depicted in Figure 5. Of course, if radiation from a single basis function pair is being considered, then the range of this integral will be limited to the support (non-zero range) of the spatial basis function 𝑏 𝑛 ( 𝜇 𝐲 is shorthand for the value of 𝜇 that gives 𝐲 ).

𝐑= 𝐱−𝐲 is the vector of the ray from 𝐲 to 𝐱 , which has length 𝑅= |𝐑| and direction 𝐑= 𝐑 ̂ 𝑅 ⁄ . The ray nature of the assumed propagation means than 𝐤 ̂ + = 𝐤 ̂ = 𝐑 ̂ . The former relation has

Figure 5: Geometry for radiation from a boundary segment to a receiver point.

already been substituted to get eq. 6. 𝐤 ̂ + was linked to 𝜒 by eq. 1, hence 𝜒 𝐲 = sin 𝜃 𝐲 , where 𝜃 𝐲 is as shown in the figure and sin 𝜃 𝐲 = 𝐭 ̂ 𝐲 . 𝐤 ̂ + . 𝐤 ̂ is an argument of the expression, however. This means it cannot be assumed that 𝐤 ̂ is aligned with 𝐑 ̂ , but the ray nature of the propagation model means that 𝑤 𝑛,𝑞 (𝐱, 𝐤 ̂ ) should be non-zero when 𝐤 ̂ = 𝐑 ̂ . The term 𝛿(cos −1 (𝐤 ̂ ∙𝐑 ̂ )) in the integrand enforces this. It is, therefore, associated with reception at the receiver. 𝑉(𝐱, 𝐲) is a visibility function that is 0 if the ray path is occluded and 1 if it is not.

The remaining two terms are cos 𝜃 𝐲 and 1 𝑅 ⁄ . The 1 𝑅 ⁄ arises due to spherical spreading from the boundary sources. It converts the ‘per radian’ of the boundary power density at 𝐲 into a ‘per metre’ tangential to the wavefront at 𝐱 . The cos 𝜃 𝐲 term is necessary because the power density includes a ‘per metre’ perpendicular to the ray, whereas the spatial integral is with respect to boundary length. Multiplication by cos 𝜃 𝐲 compensates for this. The same term appears in the so- called ‘form factor’ calculations in 3D radiosity [3]. It can be computed by cos 𝜃 𝐲 = 𝐧̂ 𝐲 . 𝐤 ̂ + .

The equivalent statement for the HF-BEM scheme in ref. [20] is * :

+ (𝐱) = ∫𝑏 𝑛 (𝜇 𝐲 )e i𝑘𝐤 ̂ + (𝜒 𝑞 )∙𝐲 [ 𝜕𝐺

(𝐱, 𝐲) −i𝑘𝐧̂ 𝐲 ∙𝐤 ̂ + (𝜒 𝑞 ) × 𝐺(𝐱, 𝐲)] 𝑑Γ 𝐲 Γ

. (7)

𝑃 𝑛,𝑞

𝜕𝑛 𝐲

+ (𝑘|𝐱−𝐲|) is the free space Green’s function for of the Helmholtz equation in 2D, with 𝐻 0

Here 𝐺(𝐱, 𝐲) = i 4 ⁄ 𝐻 0

+ being a zeroth order Hankel function of the first (outgoing) kind. This radiates as a monopole distribution. The notation 𝜕𝜕𝑛 𝐲 ⁄ is shorthand for 𝐧̂ 𝐲 ∙∇ 𝐲 . Application of this to 𝐺(𝐱, 𝐲) yields a dipole, hence 𝜕𝐺𝜕𝑛 𝐲 ⁄ radiates as a dipole distribution. The term i𝑘𝐧̂ 𝐲 ∙𝐤 ̂ + (𝜒 𝑞 ) has arisen †

from application of 𝜕𝜕𝑛 𝐲 ⁄ to e i𝑘𝐤̂ + (𝜒 𝑞 )∙𝐲 . Equations 6 and 7 have similarities. Terms in eq. 6 that do not have equivalents in eq. 7 are the delta function, which is omitted because pressure at the receiver is not directional, and the visibility function 𝑉(𝐱, 𝐲) , which is a purely geometrical concept. Terms that carry over exactly are the cos 𝜃 𝐲 , which is present in the dipole term 𝜕𝐺𝜕𝑛 𝐲 ⁄ and explicitly in the second term as 𝐧̂ 𝐲 ∙𝐤 ̂ + , and the 1 𝑅 ⁄ , which is embedded in the Hankel function ‡ . The remaining terms to be explained are 𝑎 𝑞 (𝜒 𝐲 ) in eq. 6 and e i𝑘𝐤̂ + (𝜒 𝑞 )∙𝐲 in eq. 7. These are both concerned with directionality – stated explicitly in the former and occurring through beamforming in the latter – meaning they play analogous roles. It is the spatial Fourier transform of 𝑏 𝑛 that dictates the nature of this smearing in angle, leading to the notion that it should be chosen based both on its Fourier transform as well as its direct interpolation properties in the spatial domain. But its relation to 𝑎 𝑞 is more complex, because the power arising from eq. 6 is the square of pressure in eq. 7. This leads to cross-terms involving spatial cross-correlations [19], meaning 𝑎 𝑞 is most likely related to 𝑏 𝑛 through a Wigner transform [27]. This is, however, an open research question.

The result of evaluating eq. 7 at a grid of receiver plots is plotted in Figure 6. This is shown for both discontinuous rectangular (top row) and continuous Hann (bottom row) spatial basis functions. To understand why the spreading of these waves, which is essentially aperture diffraction, is related to the Fourier transform of 𝑏 𝑛 , note that eq. 7 would simplify to 𝑃 𝑛,𝑞

+ (𝐱) = e i𝑘𝐤 ̂ + (𝜒 𝑞 )∙𝐱 if 𝑏 𝑛 weren’t present. This is the perfectly directional geometrical behaviour. The complex exponential term also gives eq. 7 the form of a Fourier transform, hence this is the classic signal processing windowing problem. Multiplication by a rectangular window smudges the directionality through convolution with a sinc function in 𝜒 . The high sidelobes of this cause the strong directional spreading in Figure 6a. The Hann window has comparatively lower sidelobes, hence less spreading is seen in Figure 6d.

* Notation adapted to match eq. 3 and eq. 6 † Better results are achieved by instead here using a Dirichlet to Neumann map [31]. But this is almost equivalent to

× i𝑘𝐧̂ 𝐲 ∙𝐤 ̂ + for the Hann spatial basis functions considered over the page, hence this simpler form is stated here.

‡ The large argument approximation of 𝐻 0 (𝑘𝑅) is √2 𝜋𝑘𝑅 ⁄ . Noting that 𝑤∝|𝑝| 2 , this leads to a 1/𝑅 term.

Figure 6: Decomposition of the result of eq. 7 (a&d) into a geometric beam (b&e) and a correction

term (c&f) for discontinuous rectangular (a-c) and continuous Hann (d-f) spatial basis functions.

Figure 6 also shows in (b) and (e) a phased beam, arising from the plane wave times a visibility function that is the geometric projection of the spatial basis function in the direction of propagation. Parts (c) and (f) are the different between these, with (c) being the classical Maggi-Rubinowicz aperture diffraction term. This has a reputation for being inaccurate for edge diffraction calculations [28], but that is associated with unrealistic boundary conditions that are used in its application to that problem. Here, as a building block for HF-BEM solutions, it is perfectly valid and provides a means to more efficiently evaluate these oscillatory boundary integral in 3D [29].

=o [PO0|/(v/ka ‘Normale secibver tates Rj

A peculiarity is that the geometric beams in parts (b) and (e) of Figure 6 do not decay with distance – they are akin to a bundle of rays all propagating in the same direction and not diverging – but the geometric propagator in eq. 6 includes a 1/𝑅 term due to spherical spreading. To investigate this inconsistency, eq. 7 was evaluated at rings of receivers of increasing radius with a few values of 𝑘𝑎 , where 𝑎 is the length of the element.

Results from this are presented in Figure 7 for the Hann envelope function. Directionality of radiation is shown in (a). This shows mainly that the angular smudging effect is reduced for large 𝑘𝑎 . This is as expected since it means the wavelength is small compared to the aperture [26]. It also confirms that the peak values are scaled by cos 𝜃= √1 −𝜒 2 , matching eq. 6. Trends with receiver ring radius are shown in (b). These are all seen to follow a 𝒪(1 √𝑅 ⁄ ) trend which, taking into account the squaring required to get from pressure to power, matches the 1/𝑅 term in eq. 6.

a)

b)

Figure 7: a) Directionality and b) on-axis distance attenua tion for pressure radiated by a HF-BEM

element. Pressure is normalised by √𝑘𝑎× √1 −𝜒 2 . Legend is common to both plots.

Hence the trends in eq. 6 are confirmed to exist in eq. 7, even if the geometric terms in Figure 6 b) & e) do not exhibit them. Increasing cancellation with diffraction must occur at greater distances. Ref. [29] also uses a decomposition based on a solid angle term and a correction [30]. This might provide more synergy with SB-GA since it matches the ‘form factor’ terms that occur in 3D [3].

5. SOURCE AND RECEIVER RECIPROCITY

The same directionality-smudging effects also affect the testing integral in Figure 3, but there is insufficient space to discuss this fully here. Instead, it is emphasised that both the Wave-Matching HF-BEM in ref. [20] and the SB-GA formulation in ref. [8] fully satisfy reciprocity, so everything that has been written about boundary to receiver mappings in section 4 equally applies to source to boundary mappings. The source of this reciprocity, which is not present in standard BEM, is a dual interpretation of the Kirchhoff-Helmholtz Boundary Integral Equation (KHBIE) derived in ref. [23] and summarized by Figure 8. The Wave-Matching BEM in [20] exploits this and essentially uses the KHBIE for both radiation and Galerkin testing, and this is the source of the reciprocity.

a)

b)

Figure 8: Complimentary interpretations of the Kirchhoff-Helmholtz boundary integral equation.

a): the standard ‘Wave Field Synthesis’ interpretation of monopole and dipole boundary

sources radiating waves which sum to form the wave at an observer point 𝐱 . b): a ‘microphone array’ interpretation, computing the cross-intensity between the measured wave

and a ‘testing wave’ (blue), being a contracting spherical wave which coalesces at 𝐱 .

To emphasise this equivalence, the source to boundary map in fig. 4 of ref. [8] is reproduced here – this time with an omnidirectional source – and mirrored in HF-BEM. For explanation of the meaning of these plots, and the geometry they pertains to, see ref. [8]. The HF-BEM results (b & c) are seen to follow the same trend as SA-GA (a) and the red ISM trendline, though the HF-BEM map extends into the evanescent region |𝜒| > 1 because it is not constrained by geometrical assumptions. Increasing 𝑘𝑎 to 100 (c) is seen to tighten the distribution compared to 𝑘𝑎= 10 (b), as is expected by the increased aperture to wavelength ratio [26].

Figure 9: Projection of radiation from a point source onto the boundary coordinate space

(𝑥, 𝜒) using: a) SB-GA, b) HF-BEM with 𝑘𝑎= 10 , c) HF-BEM with 𝑘𝑎= 100

6. CONCLUSIONS

This paper has explored synergies between the High-Frequency Boundary Element Method (HF- BEM) and Surface-Based Geometrical Acoustics (SB-GA). Interest in this is motivated by the desire to create a single unified full-audible bandwidth simulation algorithm for early time.

This has been possible due to prior efforts to align the form of these two algorithms, which are described in ref. [8] in these proceedings and ref. [20]. The latter reports the ‘The Wave-Matching Boundary Integral Equation’, which contains several novel features congruent with the espoused vision, notably the ability to be solved by marching on in reflection order. It also has an attractive energy meaning and fully satisfies reciprocity, which the majority of BEM formulations do not.

This progress is early steps towards the stated target, and many challenges to do with analysing canonical problems and implementing efficient numerical integration remain, but this is progress in that direction, nonetheless. It is hoped that by pointing out these synergies, this paper inspires more interest in this promising research area.

7. REFERENCES

[1] M. Kleiner, B.-I. Dalenbäck, P. Svensson, Auralization-An Overview, J. Audio Eng. Soc. 41 (1993) 861–875. [2] M. Aretz, R. Nöthen, M. Vorländer, D. Schröder, Combined Broadband Impulse Responses Using Fem and Hybrid Ray-Based Methods, in: EAA Symp. Auralization, Espoo, 2009: pp. 1–6. [3] L. Savioja, U.P. Svensson, Overview of geometrical room acoustic modeling techniques, J. Acoust. Soc. Am. 138 (2015) 708–730. https://doi.org/10.1121/1.4926438. [4] S. Siltanen, T. Lokki, S. Kiminki, L. Savioja, The room acoustic rendering equation, J. Acoust. Soc. Am. 122 (2007) 1624–1635. https://doi.org/10.1121/1.2766781. [5] 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. https://doi.org/10.1016/j.jcp.2012.05.028. [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. https://doi.org/10.1121/1.1416201. [7] R.S. Langley, A wave intensity technique for the analysis of high frequency vibrations, J. Sound Vib. 159 (1992) 483–502. https://doi.org/10.1016/0022-460X(92)90754-L. [8] A. Emthyas, J.A. Hargreaves, Ray tracing in Galerkin Boundary Integral form, in: Proc. Internoise, Glasgow, 2022. [9] 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. https://doi.org/10.1121/1.5096171. [10] M. Vorländer, Auralization : Fundamentals of Acoustics, Modelling, Simulation, Algorithms

and Acoustic Virtual Reality, 2nd ed., Springer, Cham, Switzerland, 2020. [11] G. Marbjerg, J. Brunskog, C.-H. Jeong, E. Nilsson, Development and validation of a

combined phased acoustical radiosity and image source model for predicting sound fields in roomsa), J. Acoust. Soc. Am. 138 (2015) 1457–1468. https://doi.org/10.1121/1.4928297. [12] M. Aretz, M. Vorländer, Combined wave and ray based room acoustic simulations of audio

systems in car passenger compartments, Part I: Boundary and source data, Appl. Acoust. 76 (2014) 82–99. https://doi.org/10.1016/j.apacoust.2013.07.021. [13] M. Aretz, M. Vorländer, Combined wave and ray based room acoustic simulations of audio

systems in car passenger compartments, Part II: Comparison of simulations and measurements, Appl. Acoust. 76 (2014) 52–65. https://doi.org/10.1016/j.apacoust.2013.07.020. [14] J.E. Summers, K. Takahashi, Y. Shimizu, T. Yamakawa, Assessing the accuracy of

auralizations computed using a hybrid geometrical-acoustics and wave-acoustics method, J.

Acoust. Soc. Am. 115 (2004) 2514. https://doi.org/10.1121/1.4809339. [15] I. Graham, E. Spence, S. Chandler-Wilde, S. Langdon, Numerical-asymptotic boundary

integral methods in high-frequency scattering, Acta Numer. 21 (2012) 89–305. https://doi.org/10.1017/S0962492912000037. [16] S. Langdon, S.N. Chandler-Wilde, A Galerkin boundary element method for high frequency

scattering by convex polygons, SIAM J. Numer. Anal. 45 (2007) 610–640. https://doi.org/10.1137/06065595X. [17] O.P. Bruno, C.A. Geuzaine, An O (1) integration scheme for three-dimensional surface

scattering problems, J. Comput. Appl. Math. 204 (2007) 463–476. https://doi.org/10.1016/j.cam.2006.02.050. [18] 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. https://doi.org/10.1002/NME.6844. [19] J.W. Rouse, A boundary element based method for accurate prediction of the surface

pressure cross-spectral density matrix, J. Acoust. Soc. Am. 141 (2017) 3697–3697. https://doi.org/10.1121/1.4988058. [20] 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. https://doi.org/10.1016/J.WAVEMOTI.2018.07.003. [21] J.A. Hargreaves, T.J. Cox, A transient boundary element method model of Schroeder diffuser

scattering using well mouth impedance., J. Acoust. Soc. Am. 124 (2008) 2942–2951. https://doi.org/10.1121/1.2982420. [22] T. Betcke, E. van ’t Wout, P. Gélat, Computationally Efficient Boundary Element Methods

for High-Frequency Helmholtz Problems in Unbounded Domains, in: Mod. Solvers Helmholtz Probl., Birkhäuser, 2017: pp. 215–243. [23] J.A. Hargreaves, Y.W. Lam, Acoustic Cross - Energy Measures and Their Applications, in:

22nd Int. Congr. Sound Vib., Florence, 2015: pp. 12–16. [24] J.A. Hargreaves, Y.W. Lam, An Energy Interpretation of the Kirchhoff-Helmholtz Boundary

Integral Equation and its Application to Sound Field Synthesis, Acta Acust. United with Acust. 100 (2014) 912–920. https://doi.org/10.3813/AAA.918770. [25] J.A. Hargreaves, Y.W. Lam, Corrigendum to An energy interpretation of the Kirchhoff-

Helmholtz boundary integral equation and its application to sound field synthesis, Acta Acust. United with Acust. 104 (2018) 1134–1134. https://doi.org/10.3813/AAA.919278. [26] J.A. Hargreaves, Acquisition of Bi-Directional Reflectance Functions by Nearfield

Acoustical Holography - a preliminary study, in: Proc. Int. Congr. Acoust., Aachen, 2019: p. 8. https://doi.org/10.18154/RWTH-CONV-239163. [27] S.C. Creagh, M. Sieber, G. Gradoni, G. Tanner, Diffraction of Wigner functions, J. Phys. A

Math. Theor. 54 (2021) 015701. https://doi.org/10.1088/1751-8121/abc72a. [28] G.M. Jebsen, H. Medwin, On the failure of the Kirchhoff assumption in backscatter, J.

Acoust. Soc. Am. 72 (1982) 1607. https://doi.org/10.1121/1.388496. [29] J.A. Hargreaves, Y.W. Lam, S. Langdon, A transformation approach for efficient evaluation

of oscillatory surface integrals arising in three-dimensional boundary element methods, Int. J. Numer. Methods Eng. 108 (2016) 93–115. https://doi.org/10.1002/nme.5204. [30] M. Albani, Boundary Diffracted Wave and Incremental Geometrical Optics: A Numerically

Efficient and Physically Appealing Line-Integral Representation of Radiation Integrals. Aperture Scalar Case, IEEE Trans. Antennas Propag. 59 (2011) 586–594. https://doi.org/10.1109/TAP.2010.2096404. [31] E. van ’t Wout, P. Gélat, T. Betcke, S. Arridge, A fast boundary element method for the

scattering analysis of high-intensity focused ultrasound., J. Acoust. Soc. Am. 138 (2015) 2726–37. https://doi.org/10.1121/1.4932166.