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

Acquisition of Bi-Directional Reflectance Functions by
Nearfield Acoustical Holography – a preliminary study
Jonathan Andrew HARGREAVES1
1

University of Salford, United Kingdom

ABSTRACT
It is well known that material absorption and scattering is dependent on incidence and observation angle.
Despite this, the corresponding standardised coefficients, which are used to represent these mechanisms
within computational acoustic models, aggregate all such dependency into single random-incidence
parameters. This limits the accuracy that can be achieved with computational acoustic models – even if thesealgorithms were to capture the wave physics perfectly, which they often do not, the results would not match
physical reality because the input data is too low resolution. Bi-Directional Reflectance Functions are an
established way of describing boundary absorption and scattering in computer graphics that have beensuggested for use in acoustics. To date, several algorithms have been published that do or could use these insimulation, but no measurement methods are available to acquire them. There is also ambiguity over some
aspects of their definition e.g. whether finite panel size is included as a scattering mechanism. This paper
adopts a definition suitable for high-frequency Boundary Element Method algorithms that use oscillatory
basis functions to capture wave directions. It then proposes an acquisition method based on double-layer
Near-Field Acoustical Holography and assesses it accuracy using 2D simulations.
Keywords: Room Acoustics, Microphone Array Methods, Scattering

1. INTRODUCTION AND BACKGROUND
Methods to measure the acoustic absorption of materials have long been of interest and various
approaches are standardised e.g. ISO 354:2003 (random-incidence in a reverberation room), ISO
10534-2:2001 (normal-incidence in an impedance tube) or ISO 13472-1:2002 (free-field). Accuratematerial absorption data is known to be a prerequisite if computational room acoustic models are
required to be accurate (1–3). This is particular clear when algorithms such as Boundary Element
Method (BEM) are applied; these are capable of producing extremely accurate results if the material
data is precisely known, but will still produce inaccurate results when it is not (4). ISO 10534-2:2001
is presently the most accurate standardised technique for absorption measurement, and can capture
surface impedance as well as absorption coefficient, but applying it for existing spaces is destructive(5). There has therefore been a recent surge in interest in ‘in-situ’ methods (6) and, of these, Nearfield
Acoustical Holography (NAH) techniques show particular promise. For example the recent method byHald et al (7) can be applied in an untreated room and requires only compact samples and hardware. It
is also capable of measuring the angle-dependence of the absorption of homogeneous samples.
Scattering from boundaries or obstacles is also known to be an important mechanism with in Room
Acoustics, with specialist ‘diffusers’ often being installed to subjectively improve the acoustics formusic performance or critical listening spaces (8). Its inclusion in room acoustic modelling software
has also been shown to be necessary to achieve realistic results (1), and predicted values for room
acoustic metrics depend strongly on the strength and type of scattering mechanism chosen (2). It istherefore now included in all commercial room acoustic modelling software, usually following therandom-incidence scattering coefficient that is measured according to ISO 17497-1:2004. The current
situation remains far from ideal though; there is relatively little measured scattering data available andwhat does exist typically has high uncertainty attached to it, leaving practitioners to make ‘best
guesses’ based on what seems to have worked satisfactorily in the past. There is also the issue that no
detail is captured about the nature of the scattering; devices that redirect sound or that scatter in only
one plane, such as the extruded diffusers common in concert halls, are not well characterised.
1

j.a.hargreaves@salford.ac.uk

4885

1.1 Directional Scattering Measurement
It is perhaps no surprise therefore that the group who have previously been most interested inmeasuring higher-resolution angle-dependent scattering data are diffuser designers. D’Antonio
discussed the possibility of standardising this as early as 1992 (9), but a scalar ‘diffusion coefficient’ toquantity uniformity of scattering was ultimately standardised in ISO 174 97-2:2012 instead; this also
includes an alternative scattering coefficient definition. The standard does however include detail on
how to measure high-resolution angle-dependent scattering because that is the raw data from which thecoefficients are extracted. This is achieved using a so-called ‘goniometer’ setup, where both the sourceand sensors are in the quasi-far field of the object to be measured.
There are several problems with this approach in the author’s opinion. Firstly, there are issues ofpracticality. They rely on subtraction approaches to separate incident and scattered sound, so requirevery low background noise, time invariant conditions (e.g. temperature), an d either an anechoic
environment or a very large room with time-gating to eliminate contamination from other reflections.
The rigs are also large and unwieldly to apply in situ (10). Secondly, their design is fundamentally
incompatible with the way that the directional data they measure could be used in Geometrical
Acoustics (GA) simulation algorithms. These assume an unrealistic model of acoustic wave
propagation as geometric rays or beams, the reflection of which from a boundary is most closely
approximated by far-field (i.e. plane wave) excitation and measurement with an infinite-sized sample.
The latter restriction occurs because GA models reflection and scattering as a local phenomenon; theedge effects that occur in measurement due to finite sample size should not factor and require removalfrom the dataset. We therefore distinguish between local scattering, that could occur to surface
undulations and roughness, and the leading-order edge diffraction that occurs at the edges of panels.
The preferred approach is that the latter should be handled separately in the simulation algorithm (11),hence the task of measuring panel scattering is to quantify only the former. This is a useful distinction
but one that has been muddied by the common practice of adjusting scattering coefficient to mask theworst errors occurring due to absence of diffraction in GA simulation algorithms (12).
Goniometers are, by nature of their design, restricted to finite sample size and measurement andexcitation distances. The requirement from GA for far-field measurements of infinite-sized samples
are both contradictory and unachievable with the goniometer ‘one physical source / sensor per
measurement / observation angle’ paradigm. The coefficients defined in ISO 17497-2:2012 include
normalisation for finite sample size, but this is not the same as correcting the raw directional measured
data. Measurements of directional reflection and/or scattering with phase, as would be ideal for earlyreflections that provide critical spatial cues, is certainly out of the question.
1.2 Scattering Measurement by Nearfield Acoustical Holography
The issue of finite sample size is circumvented by measuring in the nearfield; edge effects are
relatively less significant due to proximity to the sample and can be corrected for if required (13).
Far-field data can also be recovered mathematically by applying the far -field approximation to the
Kirchhoff-Helmholtz boundary integral equation (14), often termed “nearfield to far-field transform”.
Usually this utilises a plane-wave decomposition, though Müller-Trapet (15) investigated representing
scattering using spherical harmonic functions; there it was assumed that this means the microphonearray also has to be spherical, though work by this author shows that this is not the case (16,17).
D’Antonio (9) was aware of the usefulness of NAH in 1992 through its pioneering application tosurface impedance measurement by Tamura (18,19) in 1990, but appears to have concluded that it
required sample isotropy; that was true for the method of Tamura, which assumed specular reflection,
but is not a limitation of NAH in general. Kleiner et al (20) actually used NAH to measure directional
scattering in 1995, though here were attempting to include finite sample size effects to as to mimic the
results acquired by goniometers, and viewed the difficulty in doing this as a limitation rather than anadvantage. Such early studies were limited by the availability and cost of multi-channel acquisition
hardware, which today is ubiquitous. Separation of incident and reflected waves was also an issue,
with researchers often resorting to subtraction techniques. Today, multichannel integrated
double-layer pressure-sensing (PP) microphone arrays are available, as are integrated pressure andparticle-velocity (PU) sensors, that enable this issue to be overcome.
1.3 Representation of Scattering in Acoustic Simulation Algorithms
Commercial room acoustic simulation algorithms are almost exclusively GA based, and thatparadigm will be assumed here. We therefore wish to represent how an incoming wave with amplitude

4886

𝐴i arriving from an incident angle (𝜃i , 𝜙i ) might give rise to a distribution of outgoing waves withamplitudes 𝐴o (𝜃o , 𝜙o ); here 𝜃 is azimuth and 𝜙 is polar angle. The function 𝑅 that relates thesemust be a function of both pairs of angles, hence 𝑅(𝜃o , 𝜙o , 𝜃i , 𝜙i ). Siltanen et al (21) refer to this aBi-Directional Reflectance Function (BDRF), a term transferred from computer graphics. It can also
predict the amplitude at a given outgoing angle, given the incoming distribution 𝐴i (𝜃i , 𝜙i ), by𝐴o (𝜃o , 𝜙o ) = 𝑅(𝜃o , 𝜙o , 𝜃i , 𝜙i ) ⨂ 𝐴i (𝜃i , 𝜙i ); here ⨂ represents a two-dimensional convolution over 𝜃iand 𝜙i . The BDRF is about the most general representation of geometric scattering one could have,especially if it is also deemed to be position dependent, and Siltanen et al argued how it encompasses
previous approaches and showed how it can be applied to a variety of GA algorithms.
BDRFs do however have limitations when applied to acoustics. When used with raytracing incomputer graphics, it is reasonable for a BDRF to vary spatially in an arbitrary way; this is becauselight’s very short wavelength means optical scattering structures are orders of magnitude smaller than
the geometry being modelled. For acoustics in contrast, wavelength is comparable to the size ofgeometric features in many cases, meaning angular and spatial variation of BDRFs are inextricablyinterlinked. Subject to these restrictions, BDRFs can best be understood as representing acoustic
scattering in some spatially averaged sense that uses a particular choice of spatial window function
defined on the boundary. Understanding the effect of this for both measurement and simulation
requires consideration beyond the GA paradigm, and some initial analysis will be presented herein. In
particular, the transition from diffraction due to panel size to scatterin g due to surface roughness isessentially a ‘mid-frequency’ problem, so techniques developed for high-frequency BEM (22) can be
useful. These algorithms interpolate pressure fields using oscillatory basis functions that are designed
to capture leading-order propagation directions (23), which may even be found using a geometricmethod (24). The radiation from these basis functions rather resembles geometric beams as frequency
increases; some can be stated as a geometric beam plus a correction term (25). One recent algorithm(26), like GA, uses reflectance boundary conditions and is solved by marching on in reflection order.
Following the acknowledgement that acoustic BDRFs must be understood in a spatially-averaged
sense over some surface ‘patch’, the category of algorithms to which they are most suited is what
Svensson and Savioja (27) class as “surface-based” GA. Those that discretise wave arrival andreflection angles, in addition to spatial discretization with a mesh, have been given various names
including ‘Acoustic Radiance Transfer’ (ART) (21) and ‘Dynamical Energy Analysis’ (DEA) (28).
BDRFs will also require discretisation over space and angle; notably, if this matches what is used in
the simulation algorithm then the convolution above becomes a straightforward matrix multiplication.
Siltanen et al proposed piecewise-constant angular discretisation, but Chappell and co-workers
developed superior hierarchical schemes based on Legendre polynomials in 2D (28) and Zernike
polynomials in 3D (29). Both groups used discontinuous spatial discretisation, however continuous
spatial discretisation will be employed here because it produces beams that are more geometric innature (25), aids convergence of high-frequency BEM (26) and reduces microphone array sidelobes.

2. MODELLING FRAMEWORK AND TESTBED
When using BDRFs in a numerical model, the primaryobjective is to replace a complicated boundary Γ, which is either
costly to simulated or can only be measured, with a simplified
boundary Γs with BDRFs that give equivalent scattering
behaviour and absorption. A numerical testbed is useful sincereflections computed from BDRFs on the simplified boundary
can be validated against reflections computed using a directmodel of the complicated boundary e.g. using BEM.
Figure 1a show the configuration. A sample with complicatedboundary Γ (black) lies below a fictional simplified boundaryΓs (dashed blue); this is shown as being planar but needn’t
necessarily be due to the findings in (16) and (26). An incoming
plane wave (green) excites the sample causing an outgoing wave
(purple) to be scattered. Both are sensed by a microphone array(blue dots) located on Γs . This arrangement of having the output
data (BDRFs) defined on the same surface that is formeasurement is unusual for NAH; it is more common to
back-propagate the measured data to the physical boundary that

a)

b)

c)

Figure 1 Procedure to produce a
local model for BDRF evaluation.

4887

produced it e.g. so that boundary vibration can be found. This step is unnecessary for the intended
application, however, and is usually ill-posed so avoiding it is advantageous.
The Kirchhoff-Helmholtz boundary integral equation, being the mathematical foundation of BEM,
is useful for analysing this problem. It says that if pressure and particle velocity are known on theentire simplified boundary, then the reflections from the complicated boundary can be perfectly
reproduced from that data. This is a useful benchmark but is impractical for entire scenarios since
generating this data requires solution with an accurate algorithm such as BEM over the entire
complicated boundary. If this can be done reasonable effort then it renders the entire encoding as
BDRFs process pointless.
Instead it is desirable to exploit the assumption that scattering is a local phenomenon and try andacquire BDRFs via local simulations that are much computationally cheaper. There is precedent forthis in diffusion coefficient measurement and simulation; no one expects to have to simulate an entire
concert hall to understand how a diffuser performs, and BEM simulations of small samples have beenshown to be very effective (8). The complicated boundary could of course be simply truncated but thiswould likely produce strong edge effects, hence Figure 1b&c propose a more sophisticated approach.
In Figure 1b the complicated boundary is truncated but embedded in an infinite planar baffle; thisreduces truncation effects but means that the reflected wave (orange) from these sections can becomputed analytically. To avoid meshing the planar baffle, symmetry is exploited in Figure 1c,reflecting both the obstacle and the incident wave, which then becomes the analytical reflected wave
(orange). This approach is widely exploited in BEM and means that only the truncated section of thecomplicated boundary needs to be meshed, yielding much-reduced computational cost.
2.1 Mathematical Formulation
The mathematical formulation used for the testbed will now be presented. The model is formulated
in two dimensions for simplicity; all quantities and geometry are assumed to be invariant of the otherthird dimension. All quantities are time harmonic with time dependence e−i𝜔𝑡 , where 𝜔 is frequency
expressed in radians per second. Acoustic waves in the air domain have pressure 𝑝 that satisfies
Helmholtz’ equation ∇2 𝑝(𝐱) + 𝑘 2𝑝(𝐱) = 0. Here 𝐱 is a point in 2D Cartesian space and 𝑘 = 𝜔⁄𝑐0 is
the wavenumber, where 𝑐0 is the speed of sound in air. The simplified boundary Γs is assumed to bê are respectively its tangential and normal unit vectors and 𝐧̂ isplanar for ease of analysis; 𝐭̂ and 𝐧deemed to point away from the sample being characterised into the air volume. It is offset a distancê ∙ 𝐱 = 𝑑 for all 𝐱 ∈ Γs .
𝑑 from the coordinate origin, which is taken to align with the sample, hence 𝐧
A key concept in this formulation is that total pressure 𝑝t can be decomposed into an incoming
wave with pressure 𝑝i and an outgoing wave with pressure 𝑝o, so 𝑝t (𝐱) = 𝑝i (𝐱) + 𝑝o (𝐱). The word“outgoing” has been used here because the more common words “reflected” and “scattered” have been
reserved for more specific behaviours. Incidence and reflected wave directions (green and orange in
Figure 1) will be parameterised by 𝑘t; the component of wavenumber tangential to the measurement
plane. This definition is chosen over one based on angle because it supports inhomogeneous(evanescent) plane waves, for which |𝑘t | > 𝑘. In the non-evanescent region |𝑘t | ≤ 𝑘, polar angle 𝜙can be found by 𝜙 = sin−1(𝑘t⁄𝑘). The incoming and outgoing waves have complex amplitudedensities 𝐴i (𝑘t ) and 𝐴o (𝑘t ) respectively. The pressures 𝑝i and 𝑝o can be found from these by:1 ∞𝑝i|o (𝐱) =∫ 𝐴 (𝑘 )ei[𝑘t𝐭̂∓𝑘n 𝐧̂]∙𝐱 𝑑𝑘t .
(1)2𝜋 −∞ i|o t
The term ∓ in the exponent is negative for the incoming wave and positive for the outgoing wave and𝑘n = √𝑘 2 − 𝑘t2 is the component of wavenumber normal to Γs . The positive imaginary branch of thesquare root is chosen when |𝑘t | > 𝑘 because this means the outgoing evanescent waves decay in thê, which is consistent with the notion that they arise from the sample below. It may bedirection 𝐧
noticed that eq. (1) takes the form of a spatial Fourier transform. Williams (30) calls this system ofparameterising wave direction a ‘wavenumber spectrum’ and it is the basis of classical NAH.

3. Spatial Windowing and BDRFs
A system of spatial discretisation will now be defined. For the reasons mentioned in section 1.3, it
is chosen to implement this using smooth spatial windowing functions. The intention is that thisscheme could also be used as an approximation space for the reflection-based high-frequency BEM
scheme in (26). This means the window functions must form a partition-of-unity, i.e. sum to one.

4888

A system of overlapping raised Hanning windows fulfils these criteria. These have been shown togive improved accuracy compared to non-smooth alternatives in high-frequency BEM (31) and atechnique to accelerate evaluation of integrals involving them was presented in (25). They are defined:
1 |𝜇| < 1⁄2
𝑤(𝐱) = [1⁄2 + 1⁄2 cos(2𝜋𝜇(𝐱))] × Π(𝜇(𝐱)), where Π(𝜇) = {0 |𝜇| > 1⁄2(2)

Here 𝜇 is a normalised local coordinate that runs −1⁄2 ≤ 𝜇 ≤ 1⁄2 over the support of the window.
̂ , where 𝑛 is an integer
The windows are length 𝐿 and are centered on vertices 𝐯𝑛 = 1⁄2𝑛𝐿𝐭̂ + 𝑑𝐧
window index, and 𝜇(𝐱) = 𝐭̂ ∙ [𝐱 − 𝐯𝑛 ]⁄𝐿. For each window complex amplitude densities 𝐶i,𝑛 (𝑘t ) and𝐶o,𝑛 (𝑘t ) are defined, allowing 𝑝i and 𝑝o to be found for 𝐱 ∈ Γs :∞1
𝑝i|o (𝐱) =∑ ∫ 𝐶i|o,𝑛 (𝑘t )𝑏(𝐱, 𝑘t )𝑑𝑘t , where 𝑏(𝐱, 𝑘t ) = 𝑤(𝐱)ei[𝑘t 𝐭̂∓𝑘n𝐧̂]∙[𝐱−𝐯𝑛 ] .
(3)
2𝜋
−∞
𝑛

As with eq. (1), the term ∓ is negative for the incoming wave and positive for outgoing. Comparingeq. (3) and (1), and noting that the windows form a partition-of-unity, it is apparent that the ‘correct’values are 𝐶i|o,𝑛 (𝑘t ) = 𝐴i|o (𝑘t )ei[𝑘t𝐭̂∓𝑘n 𝐧̂]∙𝐯𝑛 = 𝐴i|o (𝑘t )ei[𝑘t 𝐿𝑛⁄2∓𝑘n𝑑] . Allowing each patch to have its
own distribution is however more flexible and can result in individual 𝐶i|o,𝑛 being sparser than 𝐴i|o.
To understand this, imagine a spherical wave impinging on a planar surface. 𝐴i|o must capture thecurvature of the entire wavefront, so will contain significant amplitude over a wide range of 𝑘t. 𝐶i|o,𝑛on the other hand, need only represent a small section of the wavefront which, being quasi planar if thesource is some distance away, will result in significant amplitude over a fairly narrow range of 𝑘t.
In a GA algorithm, 𝐶i,𝑛 (𝑘t ) would be the angular distribution of rays or beams arriving at the 𝑛th
surface patch, and 𝐶o,𝑛 (𝑘t ) would be what leaves in response. It follows that the BDRF for the 𝑛th
surface patch, 𝑅𝑛 (𝑘t,o , 𝑘t,i ), must relate these so, with ⊗ representing convolution over 𝑘t,i:(4)

𝐶o,𝑛 (𝑘t,o ) = 𝑅𝑛 (𝑘t,o , 𝑘t,i ) ⊗ 𝐶i,𝑛 (𝑘t,i ).

This is our mathematically formalised definition of a patch based BDRF. It has been defined forcomplex pressure amplitudes but could readily be applied to energy, as proposed in (21), by applying
it to squared pressure amplitudes. If 𝐶i,𝑛 (𝑘t ) and 𝐶o,𝑛 (𝑘t ) are discretised somehow using a finite setof coefficients in vectors 𝐜i,𝑛 and 𝐜o,𝑛 , then it can be written in matrix form as 𝐜o,𝑛 = 𝐑 𝑛 𝐜i,𝑛 .

4. Measurement of BDRFs
The task of the microphone array is to capture instances of 𝐶i,𝑛 (𝑘t ) and 𝐶o,𝑛 (𝑘t ) such that, given
sufficient different incident wave conditions, 𝑅𝑛 (𝑘t,o , 𝑘t,i ) may be derived. It is important however torealise that the array does not measure these directly. Adopting the lexicon of variational BEM (32),𝐶i,𝑛 and 𝐶o,𝑛 are termed ‘Coefficients’; their values correctly represent 𝑝i and 𝑝o. The distributions
that the array will measure are termed ‘Projections’. Evaluating these in BEM would typically involvean inner-product integral along the measurement surface; a simple windowed beamformer would be:
𝑃i|o,𝑛 (𝑘t ) = ∫ 𝑝i|o (𝐱)𝑏∗ (𝐱, 𝑘t )𝑑𝑙𝐱 .

(5)

Γs

This also takes the form of a spatial Fourier transform, but 𝑝i|o (𝐱) has been multiplied by 𝑤(𝐱).
This becomes convolution in the 𝑘t domain, leading to 𝑃i|o,𝑛 (𝑘t ) = 𝐶i|o,𝑛 (𝑘t ) ⊗ 𝑊(𝑘t ), where 𝑊(𝑘t )is the wavenumber spectrum of the window; for the choice of window used herein it is:
𝐿

𝑘𝐿

𝐿

𝑘𝐿

𝐿

𝑘𝐿

𝑊(𝑘t ) = 2 sinc ( 2t ) + 4 sinc ( 2t + 𝜋) + 4 sinc ( 2t − 𝜋).

(6)This convolution with 𝑊(𝑘t ) has the effect of spectrally ‘smudging’ the desired spectrum
𝐶i|o,𝑛 (𝑘t ) giving the spectrum 𝑃i|o,𝑛 (𝑘t ) measured by the array. Tamura (18) was aware of this issue
but seemed primarily concerned with how it interacted with the high 𝑘t limit due to finite spatial
resolution; he mitigated this by using a dipole source that had a more favourable excitation spectrum.
Wang-Lin et al (33) studied the performance of Hanning and Tukey windows in this application in2017, showing that both outperformed the rectangular windowing of Tamura; the spectrum of this is
𝐿 × sinc(𝑘t 𝐿⁄2), which decays slower with 𝑘t than the Hanning window spectrum in eq. (6).

4889

a)

b)

c)

Figure 2: Effect of spatial window. Amplitude of spectral ‘smudging’ function for a) 𝐿 = 0.2m andb) 𝐿 = 1.0m. c) Effect on absorption coefficient measurement under dipole excitation at 2kHz.
Another reason why these authors did not run into problems is that their measure apertures werevery large, being 1.8m diameter in (18) and 2.48m in (33). At these sizes the ‘smudging’ effect isnegligible over the majority of the frequency range of interest. Portable NAH arrays are however muchsmaller. Hald et al (7) used an array that was only 0.2m across; the amplitude of 𝑊(𝑘t ) for 𝐿 =0.2m is shown in Figure 2a and can be seen to be significant over the full frequency range of interest.
In this study a compromise value of 𝐿 = 1.0m has been chosen; 𝑊(𝑘t ) for that is shown in Figure 2b.
In terms of measuring BDRFs using projections, spectral ‘smudging’ would occur for both 𝐶i,𝑛 and𝐶o,𝑛 on the left and right of eq. (4). An approximate BDRF found from 𝑃i,𝑛 and 𝑃o,𝑛 would therefore
experience both pre-convolution and post-deconvolution with 𝑊(𝑘t ) compared to one computed
exactly. This actually has zero effect if the incident wave is wavenumber pure and the sample is
specularly reflecting, but in all other cases it will cause approximation.
To illustrate this, a uniform and infinite porous material sample has been simulated. The material is
10cm thick, has a flow resitivity of 50,000Nm-4s and is modelled using the approach in (7) to givenormalised surface admittacne Ys (𝑘t ). The sample reflects specularly and the exact BDRF is given by𝑅𝑛 (𝑘t,o , 𝑘t,i) = 𝛿(𝑘t,o − 𝑘t,i ) [𝑘n,i⁄𝑘 − Ys (𝑘t,i )]⁄[𝑘n,i ⁄𝑘 + Ys (𝑘t,i )].
In the first test, the incident wave was a single plane wave of unit amplitude and tangential
wavenumber 𝑘t,inc , hence 𝐶i,𝑛 (𝑘t,i ) = 𝛿(𝑘t,i − 𝑘t,inc ). The outgoing spectrum 𝐶o,𝑛 was found by
convolution of 𝐶i,𝑛 with 𝑅𝑛 , then 𝑃i,𝑛 and 𝑃o,𝑛 were found by convolving 𝑊 with 𝐶i,𝑛 and 𝐶o,𝑛respectively. 𝑘t,inc was varied and the BDRF was estimated as 𝑃o,𝑛 (𝑘t,inc )⁄𝑃i,𝑛 (𝑘t,inc ) , which isreasonable since the sample was known to reflect specularly, and this was compared to the exact value.
As expected, the BDRF, and the absorption coefficient that can be calculated from it, matched exactly.
In the second test, the sample was subject to excitation from a dipole source and the BDRF for all
angles was measured simultaneously on the assumption of specular reflection; this mimics theapproach used in (18) and (33). The incoming spectrum was set as 𝐶I,𝑛 (𝑘t,i ) = ei𝑘𝑛 𝑧0 , following thetrend given in (18) where 𝑧0 = 1m is the height of the source above the measurement plane, and thesame processing was performed. The results for this are shown in Figure 2c as absorption coefficient
extracted from the BDRFs; the exact analytical value is compared to values computed using 𝐿 = 0.2m
and 𝐿 = 1.0m. Both array sizes perform well around normal incidence ≡ 𝑘t = 0. However, at 𝑘t = 𝑘
the absorption coefficient varies rapidly as the incident wave transitions into the evanescent region,
and at this point the smudging due to the window functions has a noticeable effect; for the smaller
array, especially, there is significant oscillation and error. The true reflection coefficient for a finite
sample in real life is unlikely to vary as extremely as this – it is a quirk of the mathematical model ofinfinite plane waves and planar boundaries – but it is true that oblique incident angles are challenging.
It is interesting to note that the troublesome 75° measurement in (7) falls at sin 75∘ = 0.96 on a 𝑘t⁄𝑘scale, indicating how compressed oblique angles are in 𝑘t and therefore vulnerable to ‘smudging’.

5. SEPARATING INCOMING AND OUTGOING WAVES
The above approach required that 𝑝i and 𝑝o were known separately; the projection operator in eq.
(5) cannot differentiate between them. This is troublesome in reality since it can only be achieved bysubtraction, as applied in (20), which is error-prone and impractical in-situ. A more sophisticatedoperator that can discriminate between incoming and outgoing waves was defined in (26). It is a
generalisation and extension of the array designs of Hulsebos et al (34) and has the physicalinterpretation of sensing common energy flux between waves (35). For this problem it can be stated as:

4890

𝜕𝑝t
𝜕𝑏∗
(𝐱)𝑏∗(𝐱, 𝑘t ) − 𝑝t (𝐱)(𝐱, 𝑘t )] 𝑑𝑙𝐱 .
𝑃i|o,𝑛 (𝑘t ) = ∫ [
𝜕𝑛
𝜕𝑛

(7)

Γ𝐬

Note that this definition operates on total pressure 𝑝𝐭 rather than 𝑝i or 𝑝o. 𝜕⁄𝜕𝑛 is shorthand for̂ i.e. perpendicular to Γ𝐬 . 𝜕𝑝t ⁄𝜕𝑛 may be measured by a particle
spatial derivative in the direction 𝐧velocity sensor or a spaced-pair of microphones. 𝜕𝑏⁄𝜕𝑛 can be estimated as ∓𝑘n 𝑏 – note it takes
opposite sign for incoming and outgoing waves – but better differentiation between 𝑝i and 𝑝o can beachieved by defining it via a ‘Dirichlet to Neumann map’ (32); this also takes into account the effect of
the window function on the wave. The net result of this is thatthe measurements of pressure and particle velocity will bewindowed differently; particle velocity is multiplied by 𝑏∗ so is
windowed by the original window 𝑤(𝜇), whereas measured
pressure, which is multiplied by 𝜕𝑏∗ ⁄𝜕𝑛 , has a different
window 𝑤 ′ (𝜇) applied to it. Because 𝑤 ′(𝜇) is derived fromthe Dirichlet to Neumann map, it will depend on both 𝑘 and 𝑘t.
Figure 3 shows these for 1kHz with 𝑘t⁄𝑘 = 1⁄2. While 𝑤(𝜇)was purely real it turns out the optimal 𝑤 ′ (𝜇) is complex.
Equation (7) has been applied in a numerical test bed thatmodels an absorbing patch within an infinite rigid baffle using
BEM in the configuration depicted in Figure 1; details of andFigure 3: Pressure and particle
results from this will be included in the conference presentation.
velocity windows at 1kHz with
𝐿 = 1 and 𝑘t⁄𝑘 = 1⁄2

6. CONCLUSIONSThis paper has proposed a mathematical formalisation for Bi-Directional Reflectance Functions(BDRFs) for acoustics, including consideration of the effect of the finite patches they must be definedand measured over. It was seen that the common trade-off between spatial windowing and spectral‘smearing’ occurs in this application too, and its effect on BDRF measurement was analysed andexplored with a numerical test case. Knowledge and formulations from high-frequency variational
Boundary Element Method (BEM) have been adopted to assist in the design and analysis. Further work
could include attempting to apply other elements of this machinery to mitigate the spectral smudging
observed in this study, and to provide a validation approach for partial or whole geometries.

REFERENCES
1.

Vorländer M. International Round Robin on Room Acoustical Computer Simulations. In: Newman M,
editor. Proceedings of the 15th International Congress on Acoustics. Trondheim; 1995; 689–92.
2. Bork I. A Comparison of Room Simulation Software - The 2nd Round Robin on Room AcousticalComputer Simulation. Acta Acust united with Acust. 2000. 86 (6): 943–56.
3. Vorländer M. Computer simulations in room acoustics: Concepts and uncertainties. J Acoust Soc Am.
2013. 133 (3): 1203–13.
4. Hargreaves JA, Rendell LR, Lam YW. A framework for auralization of boundary element method
simulations including source and receiver directivity. J Acoust Soc Am. 2019. 145 (4): 2625–37.
5. Aretz M, Vorländer M. Combined wave and ray based room acoustic simulations of audio systems in
car passenger compartments, Part I: Boundary and source data. Appl Acoust. 2014. 76: 82–99.
6. Brandão E, Lenzi A, Paul S. A Review of the In Situ Impedance and Sound Absorption Measurement
Techniques. Acta Acust united with Acust. 2015. 101 (3): 443–63.
7. Hald J, Song W, Haddad K, Jeong C-H, Richard A. In-situ impedance and absorption coefficientmeasurements using a double-layer microphone array. Appl Acoust. 2019. 143: 74–83.
8. Cox TJ, D’Antonio P. Acoustic Absorbers and Diffusers: Theory, Design and Application. 2nd ed.
Taylor & Francis; 2009. 496 p.
9. D’Antonio P, Konnert J. The Directional Scattering Coefficient: Experimental Determination. J Audio
Eng Soc. 1992. 40( 12): 997–1017.
10. Ducourneau J, Faiz A, Chatillon J. New device for measuring mapping of sound scattering coefficientsof vertical uneven surfaces in a reverberant workplace. Applied Acoustics. 2015. 90: 21-30
11. Svensson UP, Fred RI, Vanderkooy J. An analytic secondary source model of edge diffraction impulse
responses. J Acoust Soc Am. 1999. 106 (5): 2331–44.
12. Christensen CL, Rindel JH. A new scattering method that combines roughness and diffraction effects.

4891

In: Forum Acusticum 2005. Budapest, Hungary; 2005;
13. Ottink M, Brunskog J, Jeong C-H, Fernandez-Grande E, Trojgaard P, Tiana-Roig E. In situ
measurements of the oblique incidence sound absorption coefficient for finite sized absorbers. J
Acoust Soc Am. 2016. 139 (1): 41–52.
14. Richard APA, Fernandez Grande E, Brunskog J, Jeong C-H. Characterization of acoustic scatteringfrom objects via near-field measurements. In: Euronoise 2018. Hersonissos, Crete. p. 2195–202.
15. Müller-Trapet M. Measurement of surface reflection properties : concepts and uncertainties. Phd
Thesis. RWTH Aachen; 2015.
16. Hargreaves JA, Lam YW. An Energy Interpretation of the Kirchhoff-Helmholtz Boundary Integral
Equation and its Application to Sound Field Synthesis. Acta Acust united with Acust. 2014. 100 (5):912–20.
17. Hargreaves JA, Lam YW. Corrigendum to An energy interpretation of the Kirchhoff-Helmholtz
boundary integral equation and its application to sound field synthesis. Acta Acust united with Acust.
2018. 104 (6): 1134–1134.
18. Tamura M. Spatial Fourier transform method of measuring reflection coefficients at oblique incidence.
I: Theory and numerical examples. J Acoust Soc Am. 1990. 88 (5): 2259.
19. Tamura M. Spatial Fourier-transform method for measuring reflection coefficients at oblique
incidence. II. Experimental results. J Acoust Soc Am. 1995. 97 (4): 2255–62.
20. Kleiner M, Gustafsson H, Backman J. Measurement of Directional Scattering Coefficients Using
Near-Field Acoustic Holography and Spatial Transformation of Sound Fields. In: Audio Engineering
Society Convention 99. 1995 Oct;
21. Siltanen S, Lokki T, Kiminki S, Savioja L. The room acoustic rendering equation. J Acoust Soc Am.
2007. 122 (3): 1624.
22. Graham I, Spence E, Chandler-Wilde S, Langdon S. Numerical-asymptotic boundary integral methods
in high-frequency scattering. Acta Numer. 2012. 21: 89–305.
23. Hargreaves JA, Lam YW, Langdon S, Hewett DP. A high-frequency BEM for 3D acoustic scattering.
In: 22nd International Congress on Sound and Vibration, ICSV 2015. 2015;
24. Groth SP, Hewett DP, Langdon S. A hybrid numerical–asymptotic boundary element method for highfrequency scattering by penetrable convex polygons. Wave Motion. 2018. 78: 32–53.
25. Hargreaves JA, Lam YW, Langdon S. A transformation approach for efficient evaluation of oscillatory
surface integrals arising in three-dimensional boundary element methods. Int J Numer Methods Eng.
2016. 108 (2): 93–115.
26. Hargreaves JA, Lam YW. The Wave-Matching Boundary Integral Equation — An energy approach toGalerkin BEM for acoustic wave propagation problems. Wave Motion. 2018. 87: 4–36.
27. Savioja L, Svensson UP. Overview of geometrical room acoustic modeling techniques. J Acoust Soc
Am. 2015. 138 (2): 708–30.
28. Chappell DJ, Tanner G, Giani S. Boundary element dynamical energy analysis: A versatile method forsolving two or three dimensional wave problems in the high frequency limit. J Comput Phys. 2012.
231 (18): 6181–91.
29. Bajars J, Chappell DJ, Søndergaard N, Tanner G. Transport of phase space densities throughtetrahedral meshes using discrete flow mapping. J Comput Phys. 2017. 328: 95–108.
30. Williams EG. Fourier acoustics: Sound radiation and nearfield acoustical holography. San Diego:
Academic Press; 1999. 306 p.
31. Peake MJ, Trevelyan J, Coates G. Novel basis functions for the partition of unity boundary element
method for Helmholtz problems. Int J Numer Methods Eng. 2013. 93 (9): 905–18.
32. Betcke T, van ’t Wout E, Gélat P. Computationally Efficient Boundary Element Methods forHigh-Frequency Helmholtz Problems in Unbounded Domains. In: Modern Solvers for Helmholtz
Problems. Birkhäuser; 2017. p. 215–43.
33. Lin W-L, Bi C-X, Opdam R, Zhang Y-B, Vorländer M. Performance of Spatial Windows in the Spatial
Fourier Transform Technique for the Angle-Dependent Reflection Factor Measurement. Acta Acust
united with Acust. 2017.103 (2): 349–53.
34. Hulsebos E, de Vries D, Bourdillat E. Improved microphone array configurations for auralization ofsound fields by Wave Field Synthesis. In: Audio Engineering Society Convention 110. Amsterdam;
2001 May;
35. Hargreaves JA, Lam YW. Acoustic Cross - Energy Measures and Their Applications. In: The 22ndInternational Congress on Sound and Vibration. Florence; 2015; p. 12–6.

4892