Simulating room transfer functions between transducers mounted on audio devices using a modified image source method

Xu, Zeyu; Herzog, Adrian; Lodermeyer, Alexander; Habets, Emanuël A. P.; Prinn, Albert G.

doi:10.1121/10.0023935

The image source method (ISM) is often used to simulate room acoustics due to its ease of use and computational efficiency. The standard ISM is limited to simulations of room impulse responses between point sources and omnidirectional receivers. In this work, the ISM is extended using spherical harmonic directivity coefficients to include acoustic diffraction effects. These effects occur in practice when transducers are mounted on audio devices of finite spatial extent, e.g., modern smart speakers with loudspeakers and microphones. The proposed method is verified using finite element simulations of various loudspeaker and microphone configurations in a shoebox-shaped room. It is shown that the accuracy of the proposed method is related to the sizes, shapes, number, and positions of the devices inside a room. A simplified version of the proposed method, which can significantly reduce computational effort, is also presented. The proposed method and its simplified version can simulate room transfer functions more accurately than currently available image source methods and can aid the development and evaluation of speech and acoustic signal processing algorithms, including speech enhancement, acoustic scene analysis, and acoustic parameter estimation.

I. INTRODUCTION

Room acoustic simulation has been an active research field since the 1960s.¹ It is of paramount importance for, e.g., auralization and virtual acoustics,^2–8 evaluation of acoustic signal processing algorithms,^9–13 and data generation for data-driven methods.^14–21 Room acoustics simulation methods can be broadly categorized into geometric methods,^22–26 wave-based methods,^27–32 hybrid methods,^33,34 and, more recently, deep-learning–based methods.^8,35,36 This paper presents an extension to the well-known geometrical method, the image source method (ISM), which enables the inclusion of acoustic diffraction effects from audio devices of finite extent, e.g., smart speakers on which the transducers are mounted, in shoebox-shaped rooms. The extension to arbitrary room shapes, which has been proposed by Borish³⁷ and realized in well-known toolboxes,^4,38 and the consideration of sound scattering by room boundary surfaces is outside the scope of this paper.

Modern smart speakers are a typical embodiment of an audio device that contains loudspeakers and microphones. Generating a reverberant room acoustic simulation that includes the directivities of loudspeakers and microphones due to acoustic diffraction effects is essential for developing and testing new acoustic signal processing algorithms. When using wave-based methods, the acoustic diffraction effects caused by the loudspeaker enclosure are implicitly modeled. However, because of the much higher computational cost of solving wave-based models, geometrical methods are often preferred. While assuming that sound waves behave like light rays does not apply to long wavelengths, Aretz et al.³⁹ have shown that the ISM can be used to approximate wave-based method solutions above the Schroeder frequency⁴⁰ in rectangular rooms by specifying angle-dependent reflection coefficients. Thus, the ISM can be used to obtain acceptable solutions across most of the audible frequency range above the Schroeder frequency.

The ISM has been implemented in many widely used toolboxes,^4,38,41–43 along with various approaches for reducing the computational complexity, for example, parallelization based on graphics processing units (GPUs).^43,44 Directional properties of the source and receiver were also incorporated into the ISM to produce more realistic room impulse responses,^{38,41,42,45–47} and efforts have been made to include acoustic diffraction effects due to the human head⁴⁵ and rigid spherical baffles.^48,49 Jarrett et al.⁴⁸ extended the ISM to include rigid spherical microphones using an analytical expression of the diffraction effects in the spherical harmonic domain. Samarasinghe et al.⁵⁰ further extended the ISM, viz., the generalized ISM (GISM), to facilitate parameterization of the reverberant transfer function between a directional source and a directional receiver. More recently, an extension of the ISM to incorporate loudspeaker directivities in the far field in the spherical harmonic domain has been reported.⁵¹ However, the existing approaches cannot simulate a transfer function between a directive source and a directive receiver that includes near-field diffraction effects.

In this work, the GISM⁵⁰ is extended to enable the incorporation of the local acoustic diffraction effects that result from the interaction of sound waves with audio devices of finite spatial extent, which are used to excite and capture the acoustic field in a room. First, we review the basic concepts of the ISM and its extension to the spherical harmonic domain in Sec. II. In Sec. III, we present our proposed method, which incorporates the source directivity coefficients computed from solutions of an exterior radiation problem and receiver directivity coefficients obtained using reciprocity.⁵² The pseudocode of the proposed method is also included at the end of the section. In addition, we derive a simplified version of the proposed method to reduce computational complexity. The finite element method (FEM) is used to verify the proposed method by comparing room transfer functions (RTFs) with the source and receiver mounted on either different devices or on the same device. The experimental setups and the evaluation of the proposed method are presented in Sec. IV. The effect of varying the parameters of the device, e.g., its position in the room, physical shape, and size, on the accuracy of the proposed method, is investigated. Finally, conclusions are drawn in Sec. V.

II. THE IMAGE SOURCE METHOD

A. Standard image source method

The ISM was first applied to acoustics by Carslaw⁵³ and later implemented using a digital computer by Allen and Berkley²² to model the impulse response between a point source and an omnidirectional receiver in a reverberant rectangular room.

To derive the standard ISM, consider a rectangular room with length, width, and height given by L_x, L_y, and L_z, as well as an acoustic point source and an omnidirectional receiver inside of the room at positions

x_{s} = (x_{s}, y_{s}, z_{s})

and

x_{r} = (x_{r}, y_{r}, z_{r})

⁠, given in Cartesian coordinates, respectively. The room transfer function can be expressed as a weighted sum of free-space Green's functions attenuated by the walls, that is

H_{ISM} (x_{r} | x_{s}, k) = \sum_{p \in P} \sum_{q \in Q} β (p, q) G (x_{r} | x_{p, q}^{sI}, k),

(1)

where k is the wavenumber. The triples p and q determine the position of the images and the attenuation factor

β (p, q)

⁠. The triple

p = (p_{x}, p_{y}, p_{z})

is used to determine whether the image source is mirrored with respect to walls at x = 0, y = 0, or z = 0, and its elements take values of 0 or 1, resulting in a set

P

of eight combinations. Each element of the other triple

q = (q_{x}, q_{y}, q_{z})

takes a value between

[- N_{m}, N_{m}]

⁠, and it includes higher order reflections, where N_m is used to limit computational complexity and circular convolution errors,⁴⁸ resulting in a set

Q

of

{(2 N_{m} + 1)}^{3}

combinations.

β (p, q) = β_{x_{1}}^{| q_{x} - p_{x} |} β_{x_{2}}^{| q_{x} |} β_{y_{1}}^{| q_{y} - p_{y} |} β_{y_{2}}^{| q_{y} |} β_{z_{1}}^{| q_{z} - p_{z} |} β_{z_{2}}^{| q_{z} |}

and

β_{x_{1}}, β_{x_{2}}, β_{y_{1}}, β_{y_{2}}, β_{z_{1}}, β_{z_{2}}

are the angle-independent reflection coefficients of the walls, and the subscripts a₁, a₂ (⁠

a \in {x, y, z}

⁠) of

β_{(\cdot)}

correspond to the boundaries at a = 0 and the boundaries at a = L_a, respectively. For an image source with position

x_{p, q}^{sI} = [\begin{matrix} x_{s} - 2 p_{x} x_{s} + 2 q_{x} L_{x} \\ y_{s} - 2 p_{y} y_{s} + 2 q_{y} L_{y} \\ z_{s} - 2 p_{z} z_{s} + 2 q_{z} L_{z} \end{matrix}],

(2)

the Green's function can be written as

G (x_{r} | x_{p, q}^{sI}, k) = \frac{e^{- i k | | R_{p, q}^{sI \to r} | |}}{4 π | | R_{p, q}^{sI \to r} | |},

(3)

where i denotes the imaginary unit. The vector pointing from every image of the source to the receiver is expressed as

R_{p, q}^{sI \to r} = x_{r} - x_{p, q}^{sI} = [\begin{matrix} x_{r} - x_{s} + 2 p_{x} x_{s} - 2 q_{x} L_{x} \\ y_{r} - y_{s} + 2 p_{y} y_{s} - 2 q_{y} L_{y} \\ z_{r} - z_{s} + 2 p_{z} z_{s} - 2 q_{z} L_{z} \end{matrix}] .

(4)

B. Extensions of the standard ISM

Many extensions to the ISM have been proposed to facilitate more realistic modeling of room acoustics. In this section, we discuss two important extensions of the ISM relevant to the present study, namely, using angle-dependent reflection coefficients and including source and receiver directivities.

1. Angle-dependent reflection coefficients

In Eq. (1), angle-independent reflection coefficients are used. Aretz et al.³⁹ have shown that including the angle of incidence when specifying the boundary conditions in the ISM can result in solutions that agree with FEM solutions. As such, angle-dependent reflection coefficients are used in this work.

For each reflection path represented via the indices

(p, q)

in the ISM, three incident angles in Cartesian coordinates are needed, i.e.,

θ_{a} (p, q) = arccos (\frac{a_{R_{p, q}^{sI \to r}}}{| | R_{p, q}^{sI \to r} | |}),

(5)

where

a \in {x, y, z}, a_{R_{p, q}^{sI \to r}}

is the corresponding coordinate of the vector

R_{p, q}^{sI \to r}

and

θ_{a} (p, q)

is the incident angle on the walls perpendicular to the

a -

axis. The incident angle

θ_{a} (p, q)

is the same for the reflections at the walls perpendicular to the same axis in a rectangular room.³⁹ For a uniform normalized impedance, ζ, the angle-dependent reflection coefficients are given by⁵⁴

β_{a_{2}} = β_{a_{1}} = \frac{ζ cos (θ_{a} (p, q)) - 1}{ζ cos (θ_{a} (p, q)) + 1} .

(6)

From this point onward, the attenuation term due to the reflections,

β (p, q)

⁠, contains the angle-dependent reflection coefficients defined in Eq. . Note that, while the reflection properties of surfaces are generally frequency dependent, we use a frequency-independent impedance ζ for all walls in this work without loss of generality.

2. Extension to the spherical harmonic domain

In practice, the sound source and receiver in a room usually exhibit directional properties that are far more complex than those of a point source or an omnidirectional receiver. Various methods for incorporating source and receiver directivities into the ISM have been presented in the literature.^{15,18,38,41,42,45–47,49} In this section, we present the spherical microphone array impulse response (SMIR) generator proposed by Jarrett et al.⁵⁵

The directivity of a source or receiver depends not only on the directional properties of the transducer itself but also on acoustic wave interaction with the physical object on which the transducer is mounted. The SMIR generator extends the ISM to the spherical harmonic domain, which facilitates the incorporation of the acoustic diffraction effect due to a rigid sphere and enables the simulation of RTFs for microphones mounted on such a sphere. When considering the receiver at

x_{r}

as the center of the spherical array with an omnidirectional microphone positioned at

y^{(r)}

⁠, the SMIR generator computes the RTF via

H_{SMIR} (y^{(r)} | x_{s}, k) = \sum_{p \in P} \sum_{q \in Q} β (p, q) G_{N} (y^{(r)} | {\tilde{R}}_{p, q}^{r \to sI}, k),

(7)

where

G_{N} (\cdot | \cdot)

denotes the Green's function fulfilling Neumann boundary conditions on the sphere and

{\tilde{R}}_{p, q}^{r \to sI}

is the vector pointing from the center of the array to the source image in spherical coordinates. Detailed expressions of

H_{SMIR} (y^{(r)} | x_{s}, k)

can be found in the literature.⁵⁵ It is worth noting that the Green's function given in Eq. and the Neumann Green's function in Eq. use different coordinate systems for the input arguments.

The acoustic diffraction on the sphere affects the directivity of the microphones on the sphere. Further directional properties of the source and receiver are usually implemented as a gain factor for each reflection path in the ISM. For example, Eq. is modified by considering the source directivity⁴⁹ yielding the expression

\begin{matrix} H_{SMIR +} (y^{(r)} | x_{s}, k) = \sum_{p \in P} \sum_{q \in Q} β (p, q) g (θ_{p, q}^{e, o}, ϕ_{p, q}^{e, o}, k) \\ \times G_{N} (y^{(r)} | {\tilde{R}}_{p, q}^{r \to sI}, k), \end{matrix}

(8)

where

θ_{p, q}^{e, o}, ϕ_{p, q}^{e, o}

are the emission inclination and azimuth angles of the image source with indices

(p, q)

with respect to the orientation of the source directivity pattern. The angles

θ_{p, q}^{e, o}, ϕ_{p, q}^{e, o}

can be obtained by calculating the angles between the source orientation, e.g., the main direction of the directivity pattern at the cell with indices

(p, q)

⁠, and the vector from the source image with indices

(p, q)

to the receiver. For first-order directivity patterns,

g (\cdot)

⁠, can be expressed as⁴⁹

g (θ_{p, q}^{e, o}, ϕ_{p, q}^{e, o}, k) = δ + (1 - δ) cos θ_{p, q}^{e, o} sin ϕ_{p, q}^{e, o},

(9)

where

0 \leq δ \leq 1

⁠. The relation between the emission and impinging angles, i.e., the mirrored emission angles after all reflections, depends on the number of reflections by walls perpendicular to different axes. A detailed formulation of this relation is provided by Hafezi et al.;⁴⁹ however, the phase information of the source directivity is missing in this approach.

C. Generalized image source method

The method introduced in Eq. (8) to model directional sources and receivers is restricted to directivity patterns with known analytical expressions. In practice, directivity patterns may exhibit complex shapes without analytical expressions. Samarasinghe et al.⁵⁰ further extended the ISM to facilitate the parameterization of the transfer function between a source with arbitrary directivity and a receiver region. The directional properties of the receiver region can be formed using beamforming techniques based on the received sound field. In this section, we describe the GISM, which is the foundation of our proposed method.

The GISM decomposes the transfer function into three parts: (i) an outgoing sound field from the directional source, (ii) an incident sound field captured by the receiver that can be directional, and (iii) the mode coupling coefficients that couple the source and receiver directivities in the spherical harmonic domain.

Consider an outgoing sound field from a directional source, for example, a loudspeaker, located at

x_{s}

⁠, as shown in Fig. 1. The sound field observed at an arbitrary point r outside of the sphere enclosing the entire speaker can be computed as⁵⁶

P (r, k) = \sum_{n = 0}^{\infty} \sum_{m = - n}^{n} C_{n, m}^{(s)} (k) h_{n} (k d) Y_{n, m} (θ, ϕ),

(10)

where

(d, θ, ϕ)

denote the spherical coordinates associated with r, viz., the distance d, the inclination angle θ, and the azimuth angle

ϕ

⁠. Moreover,

h_{n} (\cdot)

denotes the spherical Hankel function of order n,

Y_{n, m}

denotes the spherical harmonic function of order n and mode m, and

C_{n, m}^{(s)} (k)

represent the spherical harmonic source directivity coefficients. In practice, the infinite summation in Eq. is truncated to a maximum order

N = ⌈ k r_{0}^{(s)} ⌉

⁠, where

r_{0}^{(s)}

denotes the radius of the sphere that is large enough to surround the loudspeaker.^57,58

FIG. 1.

View large Download slide

Illustration of the scenario modeled by Samarasinghe et al. (Ref. 50). A directional source is positioned inside a rectangular room. The room transfer function from the source located at $x_{s}$ to an observation point $y^{(r)}$ around the receiver located at $x_{r}$ can be simulated.

The room transfer function from the source

x_{s}

to any point

y^{(r)}

with spherical coordinates

(d_{y}^{(r)}, θ_{y}^{(r)}, ϕ_{y}^{(r)})

relative to

x_{r}

can be expressed using the GISM as⁵⁰^,

\begin{matrix} H_{GISM} (y^{(r)} | x_{s}, k) = \sum_{v = 0}^{V} \sum_{u = - v}^{v} \sum_{n = 0}^{N} \sum_{m = - n}^{n} C_{n, m}^{(s)} (k) γ_{v, u}^{n, m} (k) \\ \times j_{v} (k d_{y}^{(r)}) Y_{v, u} (θ_{y}^{(r)}, ϕ_{y}^{(r)}), \end{matrix}

(11)

where

V = ⌈ k d_{y}^{(r)} ⌉

denotes the truncation limit that produces a small error⁵⁹ and

j_{v} (\cdot)

denotes the nth-order spherical Bessel function of the first kind. The reverberant mode coupling coefficients

γ_{v, u}^{n, m} (k)

are given by

\begin{matrix} γ_{v, u}^{n, m} (k) = \sum_{p \in P} \sum_{q \in Q} β (p, q) {(- 1)}^{(p_{y} + p_{z}) m + p_{z} n} \\ \times α_{v, u}^{n, {(- 1)}^{p_{x} + p_{y}} m} (R_{p, q}^{sI \to r}), \end{matrix}

(12)

where

R_{p, q}^{sI \to r}

has spherical coordinates

(| | R_{p, q}^{sI \to r} | |, Ω_{p, q}^{sI \to r})

⁠, and

Ω_{p, q}^{sI \to r} = (θ_{p, q}^{sI \to r}, ϕ_{p, q}^{sI \to r})

⁠. The additional factors

{(- 1)}^{(p_{y} + p_{z}) m + p_{z} n}

and

{(- 1)}^{p_{x} + p_{y}}

result from the mirror effect of the image sources⁶⁰ on the spherical harmonics, which is a more general description compared to the use of emission and impinging angles as applied in the SMIR+ method.⁴⁹ Using

x_{0}

with spherical coordinates

(d_{x_{0}}, θ_{x_{0}}, ϕ_{x_{0}})

⁠, instead of

R_{p, q}^{sI \to r}

for brevity, the single-path mode coupling coefficients

α_{v, u}^{n, m} (\cdot)

in Eq. can be derived using the translation of fields,⁶⁰

\begin{matrix} α_{v, u}^{n, m} (x_{o}) = 4 π i^{v - n} {(- 1)}^{m} \sum_{l = | n - v |}^{n + v} i^{l} h_{l}^{(2)} (k d_{x_{0}}) \\ \times Y_{l, m - u} (θ_{x_{0}}, ϕ_{x_{0}}) W_{1} W_{2} ξ, \end{matrix}

(13)

with W₁ and W₂ denoting Wigner 3-j symbols,

W_{1} = (\begin{matrix} n & v & l \\ 0 & 0 & 0 \end{matrix}), W_{2} = (\begin{matrix} n & v & l \\ - m & u & m - u \end{matrix}),

(14)

and

ξ = \sqrt{(2 n + 1) (2 v + 1) (2 l + 1) / 4 π}

⁠.

D. Fast source room receiver method

The fast source room receiver (FSRR) method, proposed by Luo and Kim,⁵¹ also simulates the RTF in the spherical harmonic domain. The RTF parameterization includes an approximation of the spherical harmonic description of the source and receiver directivities and utilizes a simplification of the cross-modal coupling, i.e., the mode coupling coefficients in Eq. (12). The method also uses a randomized sign parameter. Notably, the FSRR method uses reversed reflection paths—a technique explored in this work for deriving a simplified version of the proposed method, which is presented in Sec. III D.

Using our notation, the FSRR method RTF is given by

\begin{matrix} H_{FSRR} (x_{r} | x_{s}, k) = \sum_{p \in P} \sum_{q \in Q} B_{p, q} M_{p, q} \frac{e^{- i k | | R_{p, q}^{r \to sI} | |}}{| | R_{p, q}^{r \to sI} | |} \\ \times \sum_{n = 0}^{N} \sum_{m = - n}^{n} {\tilde{C}}_{n, m}^{(s)} (k) Y_{n, m} (Ω_{\tilde{p}, \tilde{q}}^{s \to rI}) \\ \times \sum_{v = 0}^{V} \sum_{u = - v}^{v} {\tilde{C}}_{v, u}^{(r)} (k) Y_{v, u} (Ω_{p, q}^{r \to sI}), \end{matrix}

(15)

where

B_{p, q}

randomly takes value between 1 and –1,

M_{p, q}

is a function fitted to frequency-dependent absorption coefficients, and

\tilde{C}

is a directivity coefficient obtained at 1 m from either the source or the receiver. The modified triples

\tilde{p}, \tilde{q}

are defined in Sec. III D [Eq. ],

Ω_{\tilde{p}, \tilde{q}}^{sI \to r}

is the direction of a vector pointing from a source image to the receiver, and

Ω_{\tilde{p}, \tilde{q}}^{r \to sI}

is the direction of the vector corresponding to the reversed path of the same reflection. Definitions of the two vectors are given in Sec. III D.

III. PROPOSED METHOD

We have seen that the SMIR generator relies on analytical descriptions of directivity, the GISM simulates either a directive source or a directive receiver, and the FSRR method approximates the RTF between directional sources and receivers. In this section, we propose a method for simulating RTFs that can include near-field diffraction effects of both directional sources and directional receivers.

A. Problem formulation

We consider the following scenario: one source device and one receiver device are positioned in a rectangular room, as depicted in Fig. 2. The source contains a sound-radiating transducer located at $x_{s}$ ⁠, e.g., a vibrating piston on one surface, and a microphone located at $x_{r}$ is mounted on the receiver. The two devices can have arbitrary shapes and materials, and the transducers can exhibit arbitrary directivity patterns. The acoustic field generated by the source transducer interacts with both devices and the walls and is captured by the microphone. We aim to incorporate the local diffraction effects around the devices into the RTFs based on the ISM. The existing methods are unable to resolve the problem under consideration due to the following reasons:

As shown by Xu et al.,⁵² the terms $j_{v} (k d_{y}^{(r)}) Y_{v, u} (θ_{y}^{(r)}, ϕ_{y}^{(r)})$ corresponding to the incident sound field in Eq. (11) cannot model the diffraction around the receiver. Therefore, even utilizing the acoustic reciprocity principle, the GISM can only model the RTFs when either the source or the receiver does not include local diffraction effects.
The FSRR method⁵¹ extends the ISM to include the directivity of both the source and receiver using a far-field approximation and, therefore, cannot model near-field directivities.

FIG. 2.

View large Download slide

Illustration of the scenario considered for the proposed method. Two speakers are positioned inside a rectangular room. The source speaker has a vibrating piston on one of its surfaces located at $x_{s}$ ⁠. The receiver speaker has a microphone on the top surface located at $x_{r}$ ⁠. The room transfer function from $x_{s}$ to $x_{r}$ is of interest. Not shown here is the room in which the speakers are placed. The two spherical regions marked in light gray are nonoverlapping.

B. Source and receiver directivities

The directivity of a transducer mounted on a device can be either numerically simulated or practically measured at discrete points on a sphere around a local origin. When sound waves are generated or captured by the mounted transducers, e.g., a circular piston at

x_{s}

or a microphone at

x_{r}

as shown in Fig. 2, acoustic diffraction occurs due to the wave interaction with the speakers. As shown in previous works,^52,56,61 the total sound radiation of a source device can be quantified using directivity coefficients defined on a transparent sphere with radius

r_{0}^{(s)}

⁠. Using reciprocity, the receiver directivity can be obtained analogously to the source directivity.⁵² As an example, consider the source device. If the sound field is sampled at J discrete points on the transparent sphere, one can reformulate Eq. with truncated order N as

P (r_{j}, k) = \sum_{n = 0}^{N} \sum_{m = - n}^{n} P_{n, m} (k, r_{0}^{(s)}) Y_{n, m} (θ_{j}, ϕ_{j}),

(16)

where

r_{j}

denotes the jth position on the transparent sphere,

(θ_{j}, ϕ_{j})

is the corresponding direction,

P_{n, m} (r_{0}^{(s)}, k) = C_{n, m}^{(s)} (k) h_{n} (k r_{0}^{(s)})

is the spherical wave spectrum,⁵⁶ and

j \in {1, 2, \dots, J}

⁠. Different choices of the distribution of

(θ_{j}, ϕ_{j})

can be found in the literature.⁶⁰ If the number of samples is larger than

{(N + 1)}^{2}

⁠, the spherical wave spectrum can be calculated from the sampled sound field by solving an over-determined linear system using a least squares approach. The directivity coefficients of the source can then be determined by

C_{n, m}^{(s)} (k) = P_{n, m} (r_{0}^{(s)}, k) / h_{n} (k r_{0}^{(s)})

⁠. The receiver directivity coefficients

C_{n, m}^{(r)} (k)

can be obtained following the same procedure by making use of reciprocity.⁵²

C. Room transfer function

The following parameterization of the transfer function between two transducers with arbitrary directivities in the free field,

H^{(F)}

⁠, was proposed by Xu et al.:⁵²

\begin{matrix} H^{(F)} (x_{r} | x_{s}, k) = \sum_{v = 0}^{V} \sum_{u = - v}^{v} \sum_{n = 0}^{N} \sum_{m = - n}^{n} C_{n, m}^{(s)} (k) \\ \times α_{v, u}^{n, m} (k) D_{v, u}^{(r)} (k), \end{matrix}

(17)

where

α_{v, u}^{n, m} (k)

represent the free-field mode coupling coefficients. Note that the two transparent spheres shown in Fig. 2 are nonoverlapping in this work. However, when both the source and receiver are on the same device, one can measure or simulate the direct path to avoid the overlapping of the two transparent spheres. The following relation between

D_{v, u}^{(r)} (k)

and the directivity coefficients

C_{v, u}^{(r)} (k)

obtained via reciprocity for the receiver was derived in⁵²

D_{v, u}^{(r)} (k) = \frac{i {(- 1)}^{u}}{k} C_{v, - u}^{(r)} (k) .

(18)

Since the summation over source images is a linear operation and the mirroring effects are already included in

γ_{v, u}^{n, m}

in Eq. , the total RTF can be written by substituting

γ_{v, u}^{n, m} (k)

for

α_{v, u}^{n, m} (k)

in Eq. , yielding

\begin{matrix} H_{DEISM} (x_{r} | x_{s}, k) = \sum_{v = 0}^{V} \sum_{u = - v}^{v} \sum_{n = 0}^{N} \sum_{m = - n}^{n} C_{n, m}^{(s)} (k) \\ \times γ_{v, u}^{n, m} (k) D_{v, u}^{(r)} (k) . \end{matrix}

(19)

The modified reverberant mode coupling coefficients can be defined, based on Eq. and the reciprocal relation in Eq. , as follows:

{\tilde{α}}_{v, u}^{n, m} (k) = \frac{i {(- 1)}^{u}}{k} γ_{v, u}^{n, m} (k) .

(20)

The final proposed RTF is therefore expressed as

\begin{matrix} H_{DEISM} (x_{r} | x_{s}, k) = \sum_{v = 0}^{V} \sum_{u = - v}^{v} \sum_{n = 0}^{N} \sum_{m = - n}^{n} C_{n, m}^{(s)} (k) \\ \times {\tilde{α}}_{v, u}^{n, m} (k) C_{v, - u}^{(r)} (k), \end{matrix}

(21)

where, for clarity, the proposed method is identified using the acronym diffraction enhanced image source method (DEISM).

The GISM RTF, i.e., the simulation of a directional source to a local receiving region, can be recovered from the proposed RTF as follows. Similar to Xu et al.,⁵² for a receiver region without any physical objects and an observation point

y^{(r)}

with spherical coordinates

(d_{y}^{(r)}, θ_{y}^{(r)}, ϕ_{y}^{(r)})

relative to

x_{r}

⁠, as shown in Fig. 1, the receiver directivity coefficients can be written as

C_{v, u}^{(r)} = - i k j_{v} (k d_{y}^{(r)}) Y_{v, u}^{*} (θ_{y}^{(r)}, ϕ_{y}^{(r)}) .

(22)

By substituting Eq. into Eq. , and utilizing the property of spherical harmonics

Y_{v, u}^{*} (θ_{y}^{(r)}, ϕ_{y}^{(r)}) = {(- 1)}^{u} Y_{v, - u} (θ_{y}^{(r)}, ϕ_{y}^{(r)})

⁠, one can obtain the same expression of the RTFs as the GISM shown in Eq. .

D. Far-field approximation

The computational cost of the DEISM RTF is expected to increase with frequency since the directivity patterns tend to be more complicated at higher frequencies,⁶² resulting in higher orders of the directivities, i.e., larger N and V, required in the DEISM. As either N or V increases, the number of summations over l in the mode coupling coefficients in Eq. (12) increases, thus consuming more computational resources. Therefore, it is necessary to investigate the possibility of simplifying the DEISM in Eq. (21) while maintaining sufficient accuracy. To this end, a far-field approximation is applied to each reflection path of the RTF in Eq. (21) in this section.

Without loss of generality, we consider a reflection path from the source image to the receiver, denoted by the vector

R_{p, q}^{sI \to r}

⁠, for the following analysis. The mode coupling coefficients from Eq. for a single reflection path can be expressed as follows:

{\tilde{α}}_{v, u}^{n, m} (k) = β (p, q) Λ (p, m, n) \frac{4 π i}{k} i^{v - n}

(23)

\begin{matrix} \times {(- 1)}^{m' + u} \sum_{l = | n - v |}^{n + v} i^{l} h_{l}^{(2)} (k | | R_{p, q}^{sI \to r} | |) \\ \times Y_{l, m' - u} (Ω_{p, q}^{sI \to r}) W_{1} W_{2} ξ, \end{matrix}

(24)

where

m' = {(- 1)}^{p_{x} + p_{y}} m

and

Λ (p, m, n) = {(- 1)}^{(p_{y} + p_{z}) m + p_{z} n}

are the additional effects of the mirroring on the spherical harmonics.⁵⁰ Assuming

| | R_{p, q}^{sI \to r} | |

is large, we employ the large argument of the spherical Hankel function, i.e.,

h_{l}^{(2)} (k | | R_{p, q}^{sI \to r} | |) \approx i^{l + 1} e^{- i k | | R_{p, q}^{sI \to r} | |} / k | | R_{p, q}^{sI \to r} | |

⁠. The modified mode coupling coefficients in far field can be expressed as

\begin{matrix} {\tilde{α}}_{v, u}^{n, m} (k) \approx β (p, q) Λ (p, m, n) \frac{4 π}{k} i^{v - n} \\ \times {(- 1)}^{m' + u + 1} \frac{e^{- i k | | R_{p, q}^{sI \to r} | |}}{k | | R_{p, q}^{sI \to r} | |} \sum_{l = | n - v |}^{n + v} {(- 1)}^{l} \\ \times Y_{l, m' - u} (Ω_{p, q}) W_{1} W_{2} ξ . \end{matrix}

(25)

Note that the product of two spherical harmonics can be expressed using the formula⁶³

\begin{matrix} Y_{n, m} (θ, ϕ) Y_{v, - u} (θ, ϕ) = {(- 1)}^{n + m + v - u} \\ \times \sum_{l = | n - v |}^{n + v} {(- 1)}^{l} Y_{l, m - u} (θ, ϕ) W_{1} W_{2} ξ . \end{matrix}

(26)

By applying Eq. to Eq. , we can express the mode coupling coefficients as follows:

\begin{matrix} {\tilde{α}}_{v, u}^{n, m} (k) \approx β (p, q) Λ (p, m, n) \frac{4 π}{k} i^{v - n} \\ \times {(- 1)}^{n + v - 1} \frac{e^{- i k | | R_{p, q}^{sI \to r} | |}}{k | | R_{p, q}^{sI \to r} | |} \\ \times Y_{n, m} (Ω_{p, q}^{sI \to r}) Y_{v, - u} (Ω_{p, q}^{sI \to r}) . \end{matrix}

(27)

The transfer function of one reflection path using the far-field approximation model can thus be written as

\begin{matrix} H_{DEISM - LC}^{p, q} (x_{r} | x_{s}, k) \approx - β (p, q) \frac{4 π}{k} \frac{e^{- i k | | R_{p, q}^{sI \to r} | |}}{k | | R_{p, q}^{sI \to r} | |} \\ \times \sum_{n = 0}^{N} \sum_{m = - n}^{n} Λ (p, m, n) i^{- n} {(- 1)}^{n} C_{n, m}^{(s)} (k) Y_{n, m'} (Ω_{p, q}^{sI \to r}) \\ \times \sum_{v = 0}^{V} \sum_{u = - v}^{v} i^{v} {(- 1)}^{v} C_{v, - u}^{(r)} (k) Y_{v, - u} (Ω_{p, q}^{sI \to r}), \end{matrix}

(28)

where DEISM-LC is used to refer to the low-complexity (LC) solution.

The terms $Λ (p, m, n) Y_{n, m'} (Ω_{p, q}^{sI \to r})$ are a result of mirroring of spherical harmonics. Since the source images are used in Eq. (28), the additional factors $Λ (p, m, n)$ and $m'$ need to be associated with the source spherical harmonics. To further simplify the expressions in Eq. (28), we propose to replace $Λ (p, m, n) Y_{n, m'} (Ω_{p, q}^{sI \to r})$ with a simplified expression, which can be achieved as shown in the following by using the direction that corresponds to the reversed traveling path of the vector $R_{p, q}^{sI \to r}$ ⁠.

First, we note that the vector pointing from the receiver to the source image

R_{p, q}^{r \to sI}

is the inverse of the vector

R_{p, q}^{sI \to r}

given in Eq. . In addition, we define the vector pointing from the source to the receiver images

R_{p, q}^{s \to rI}

⁠, which can be expressed similarly to Eq. as

R_{p, q}^{s \to rI} = [\begin{matrix} x_{r} - x_{s} - 2 p_{x} x_{r} + 2 q_{x} L_{x} \\ y_{r} - y_{s} - 2 p_{y} y_{r} + 2 q_{y} L_{y} \\ z_{r} - z_{s} - 2 p_{z} z_{r} + 2 q_{z} L_{z} \end{matrix}] .

(29)

In spherical coordinates, the vectors

R_{p, q}^{r \to sI}

and

R_{p, q}^{s \to rI}

have directions

Ω_{p, q}^{r \to sI}

and

Ω_{p, q}^{s \to rI}

⁠, respectively. These two vectors are reversed reflection paths, only when there is an odd number of reflections between parallel walls, as shown by the gray arrows in Fig. 3. To take the even number of reflections into account, as shown by the black arrows in Fig. 3, we use the modified vector

R_{\tilde{p}, \tilde{q}}^{s \to rI}

(with directions

Ω_{\tilde{p}, \tilde{q}}^{s \to rI}

⁠) that represents the reversed path of

R_{p, q}^{r \to sI}

with its subscript indices

\tilde{p}, \tilde{q} = ({\tilde{p}}_{x}, {\tilde{p}}_{y}, {\tilde{p}}_{z}), ({\tilde{q}}_{x}, {\tilde{q}}_{y}, {\tilde{q}}_{z})

taking the following expressions:

{\tilde{p}}_{a}, {\tilde{q}}_{a} = {\begin{matrix} p_{a}, q_{a} & if | 2 q_{a} - p_{a} | mod 2 = 1, \\ p_{a}^{'}, q_{a}^{'} & if | 2 q_{a} - p_{a} | mod 2 = 0, \end{matrix}

(30)

where

a \in {x, y, z}

⁠, and

p_{a}^{'}, q_{a}^{'}

represent the image indices that are symmetric to p_a, q_a with respect to the image indices

p_{a}, q_{a} = 0, 0

⁠. The relation between

p_{a}^{'}, q_{a}^{'}

and p_a, q_a can be summarized as

2 q_{a} - p_{a} = p_{a}^{'} - 2 q_{a}^{'},

(31)

q_{a}^{'} = ⌊ \frac{p_{a} - 2 q_{a}}{2} ⌋,

(32)

where

⌊ \cdot ⌋

denotes the floored division.

FIG. 3.

An example of reversed reflection paths along the x -direction used in the ISM. For an odd number of reflections, R 1 , 0 r → sI and R 1 , 0 s → rI share the same-length, but reversed reflection path. However, for an even number of reflections, e.g., R 0 , − 1 r → sI, the reversed path R 0 , 1 s → rI is symmetric with respect to the original room cell instead of R 0 , − 1 s → rI.

View large Download slide

An example of reversed reflection paths along the $x -$ direction used in the ISM. For an odd number of reflections, $R_{1, 0}^{r \to sI}$ and $R_{1, 0}^{s \to rI}$ share the same-length, but reversed reflection path. However, for an even number of reflections, e.g., $R_{0, - 1}^{r \to sI}$ ⁠, the reversed path $R_{0, 1}^{s \to rI}$ is symmetric with respect to the original room cell instead of $R_{0, - 1}^{s \to rI}$ ⁠.

Second, due to the properties of the mirroring effect on the spherical harmonics, one can write the following relation between the spherical harmonics associated with the spherical direction

Ω_{p, q}^{sI \to r}

in Eq. and

Ω_{\tilde{p}, \tilde{q}}^{s \to rI}

introduced above:

Λ (p, m, n) Y_{n, m'} (Ω_{p, q}^{sI \to r}) = Y_{n, m} (Ω_{\tilde{p}, \tilde{q}}^{s \to rI}) .

(33)

Since

R_{p, q}^{r \to sI} = - R_{p, q}^{sI \to r}

⁠, i.e.,

Ω_{p, q}^{sI \to r}

and

Ω_{p, q}^{r \to sI}

have opposite directions, the following relation is obtained:

Y_{v, - u} (Ω_{p, q}^{sI \to r}) = {(- 1)}^{v} Y_{v, - u} (Ω_{p, q}^{r \to sI}) .

(34)

Finally, by substituting Eqs. and in the approximated solution in Eq. , we obtain the simplified expression

\begin{matrix} H_{DEISM - LC}^{p, q} (x_{r} | x_{s}, k) \approx - β (p, q) \frac{4 π}{k} \frac{e^{- i k | | R_{p, q}^{r \to sI} | |}}{k | | R_{p, q}^{r \to sI} | |} \\ \times \sum_{n = 0}^{N} \sum_{m = - n}^{n} i^{n} C_{n, m}^{(s)} (k) Y_{n, m} (Ω_{\tilde{p}, \tilde{q}}^{s \to rI}) \\ \times \sum_{v = 0}^{V} \sum_{u = - v}^{v} i^{v} C_{v, u}^{(r)} (k) Y_{v, u} (Ω_{p, q}^{r \to sI}) . \end{matrix}

(35)

The RTF of the DEISM-LC is then given by

\begin{matrix} H_{DEISM - LC} (x_{r} | x_{s}, k) = \sum_{p \in P} \sum_{q \in Q} - β (p, q) \frac{4 π}{k} \frac{e^{- i k | | R_{p, q}^{r \to sI} | |}}{k | | R_{p, q}^{r \to sI} | |} \\ \times \sum_{n = 0}^{N} \sum_{m = - n}^{n} i^{n} C_{n, m}^{(s)} (k) Y_{n, m} (Ω_{\tilde{p}, \tilde{q}}^{s \to rI}) \\ \times \sum_{v = 0}^{V} \sum_{u = - v}^{v} i^{v} C_{v, u}^{(r)} (k) Y_{v, u} (Ω_{p, q}^{r \to sI}) . \end{matrix}

(36)

The DEISM-LC is more computationally efficient since it does not require the summation over l in the mode coupling coefficient in Eq. . The DEISM in Eq. has a time complexity of

O (I N^{2} V^{2} (N + V))

per wavenumber k, where I is the number of images, and N and V are the maximum spherical harmonic (SH) order used to describe directivities of the source and receiver. However, the time complexity of DEISM-LC is reduced to

O (I N^{2} V^{2})

per wavenumber k, which can lead to approximately N + V times faster computations compared to DEISM. As the frequency increases, larger N or V is usually necessary to accurately describe the directivity. Therefore, one might find DEISM-LC more practical for simulating high frequencies.

While the simplified method described here is similar to that of the FSRR method, we note that there are significant differences between Eq. (36) and Eq. (15). As stated in Sec. III A, the FSRR utilizes directivity coefficients ${\tilde{C}}_{n, m}^{(s)} (k), {\tilde{C}}_{v, u}^{(r)} (k)$ that are obtained at 1 m from the source or the receiver. However, the directivity coefficients used in DEISM-LC are not limited to this assumption as long as the transparent spheres enclose the devices. Furthermore, there are inherent dissimilarities between DEISM-LC and FSRR due to the different coefficients, e.g., $B_{p, q}, M_{p, q}$ used in Eq. (15) and iⁿ, i^v used in Eq. (36). This difference is also illustrated by an example described in Sec. IV B.

Algorithm 1.

Calculate lookup tables.

1: Calculate the Wigner 3-j symbols

W_{1} (n, v, l)

and

W_{2} (n, v, l, m, u)

from Eq. .

2: Calculate the spherical harmonic coefficients

P_{n, m} (k, r_{0})

shown in Eq. for both the source and receiver using the least squares approach.

3: Calculate the directivity coefficients of the source and receiver, viz.,

C_{n, m}^{(s)} (k)

and

C_{n, m}^{(r)} (k)

⁠, using the relation

P_{n, m} (r_{0}, k) = C_{n, m} (k) h_{n} (k r_{0})

⁠. Then, store the directivity coefficients as matrices

C^{(s)} (n, m, k)

and

C^{(r)} (n, m, k)

⁠.

E. Algorithm summary

A summary of the algorithm proposed in Sec. III C is listed as pseudocode in 1, 2, and 3. A Python implementation of the algorithm, along with an example, is available online.⁶⁴ The pseudocode of the far-field approximation method is not provided here for brevity.

The calculation of the Wigner 3-j symbols can be computationally expensive if they are computed for each reflection path. Since they are identical for each reflection path, to reduce the computational effort, we precalculate them for all the given spherical harmonic orders and modes $n, m, v, u, l$ ⁠, and store them in a lookup table. In addition, the directivity coefficients $C_{n, m}^{(s)} (k)$ and $C_{n, m}^{(r)} (k)$ are precalculated and stored in a lookup table. The procedures for generating the lookup tables are listed in Algorithm 1.

ALGORITHM 2.

Image generation.

1:

p = (p_{x}, p_{y}, p_{z}) \in P = {0, 1}^{3}

(A triple introduced in Table I).

2:

q = (q_{x}, q_{y}, q_{z}) \in Q = {- N_{m}, \dots, 0, \dots, N_{m}}^{3}

(A triple introduced in Table I).

3:

A = P \times Q

(A set containing all possible images.)

4: for

(p, q) \in A

do

5: if

| | 2 q_{x} - p_{x} | | + | | 2 q_{y} - p_{y} | | + | | 2 q_{z} - p_{z} | | \leq N_{o}

then {

N_{o} =

maximum reflection order}.

6: Calculate the vector pointing from the source image to the receiver

R_{p, q}^{sI \to r}

in Eq. and its spherical coordinates

(d_{x_{0}}, θ_{x_{0}}, ϕ_{x_{0}})

⁠.

7: Calculate the incident angles

θ_{a} (p, q)

with respect to different walls using Eq. .

8: Calculate the reflection coefficients

β_{x_{1}}, β_{x_{2}}, β_{y_{1}}, β_{y_{2}}, β_{z_{1}}, β_{z_{2}}

using Eq. .

9: Calculate the wall attenuation

β (p, q) = β_{x 1}^{| q_{x} - p_{x} |} β_{x 2}^{| q_{x} |} β_{y 1}^{| q_{y} - p_{y} |} β_{y 2}^{| q_{y} |} β_{z 1}^{| q_{z} - p_{z} |} β_{z 2}^{| q_{z} |}

⁠.

10: else

11:

A = A ∖ {(p, q)}

(Remove current image if it exceeds N_o.)

12: end if

13: end for

In the provided implementation, we separate the whole procedure of the ISM into two parts, viz., the image generation given in Algorithm 2 and the single-path transfer function calculation using DEISM given in Algorithm 3. The computation in the spherical harmonic domain is generally also expensive, and therefore, we use parallel computation for all reflection paths after all images are calculated.

ALGORITHM 3.

DEISM core.

1: P(k) = 0 (Initialize the sound field at

x_{r}

⁠.)

2: for

(p, q) \in A

do

3: for n = 0 to N do

4: for

m = - n

to n do

5:

Γ = {(- 1)}^{(p_{y} + p_{z}) m + p_{z} n}

[Part of mirror effect of spherical harmonics in Eq. .]

6:

m' = {(- 1)}^{p_{x} + p_{y}} m

[The modified mode m of

α_{v, u}^{n, m}

in Eq. .]

7: for v = 0 to V do

8: for

u = - v

to v do

9:

α_{v, u}^{n, m'} = 0

[Initialize the mode coupling coefficient in Eq. .]

10: for

l = | n - v |

to n + v do

11:

ξ = \sqrt{(2 n + 1) (2 v + 1) (2 l + 1) / 4 π}

⁠.

12:

α_{v, u}^{n, m'} + = 4 π i^{v - n} {(- 1)}^{m^{'}} i^{l} h_{l}^{(2)} (k d_{x_{0}}) Y_{l, m' - u} (θ_{x_{0}}, ϕ_{x_{0}}) W_{1} (n, v, l)

13:

\times W_{2} (n, v, l, m', u) ξ

(Update the mode coupling coefficient.)

14: end for

15:

P (k) + = Γ β (p, q) C^{(s)} (n, m, k) α_{v, u}^{n, m'} \frac{i {(- 1)}^{u}}{k} C^{(r)} (v, u, k)

[Summation over

A

constitutes the room transfer function in Eq. .]

16: end for

17: end for

18: end for

19: end for

20: end for

IV. EVALUATION OF THE PROPOSED METHODS

In the following, we describe the experimental setups used to evaluate the performance of DEISM and DEISM-LC [see Eqs. (21) and (36), respectively]. The simulated RTFs are plotted and compared for different scenarios. The differences between the DEISM and the FEM solutions are evaluated using error metrics, and the difference between DEISM and DEISM-LC is analyzed. In addition, the effects of incorporating directivities into the RTFs by comparing the DEISM-LC with the original ISM, and the differences between DEISM-LC and FSRR are also demonstrated by comparing them with FEM.

A. Test setup

The proposed methods require a description of the directivities of the devices, i.e., a loudspeaker (or microphone) in the spherical harmonic domain. These data can be obtained from measurements of the pressure field in the local vicinity of the transducer and its enclosure. Alternatively, directivity patterns obtained from simulations could be used. In this work, the FEM is used to simulate the free-field pressure fields surrounding the transducer and its enclosure. Aretz et al.³⁹ have demonstrated that solutions of the ISM with angle-dependent reflection coefficients are in good agreement with FEM solutions above the Schroeder frequency. Therefore, the FEM is used in this work to evaluate the performance of the proposed method.

This section presents the details of generating the simulated transducer directivities and the RTFs.

1. Transducer directivities

The three speaker shapes considered in this work are shown in Fig. 4. For each speaker, the loudspeaker driver is modeled as a vibrating piston with frequency-independent unit acceleration. For each microphone, a monopole source with frequency-independent unit volume velocity is placed at the microphone position, which is located on the top and close to the front of each enclosure.

FIG. 4.

View large Download slide

Speaker shapes considered in this work. From left to right, these speakers are referred to in the text as cuboidal (cub.), spherical (sph.), and cylindrical (cyl.), where the terms in brackets are abbreviations also used in this work. The three speakers share the same height and the same-sized vibrating piston. Note that the speakers are not to scale in this figure.

One must ensure that free-field conditions are used to capture the loudspeaker and microphone directivities. In a measurement setup, one might use an anechoic chamber. When using the FEM, a suitable non-reflecting boundary condition is required. In this work, a perfectly matched layer⁶⁵ is used to reduce unwanted reflections from the edge of the computational domain. Note that care must be taken to ensure that the perfectly matched layer sufficiently reduces the unwanted reflections—this is achieved in this work by performing a priori mesh studies. Each speaker is meshed using triangular elements, the surrounding air (with a sound speed of 343 m/s and density of 1.2 kg/m³) is meshed using tetrahedral elements, and the perfectly matched layer is meshed using hexahedral elements. An average element size that ensures that there are 10 degrees of freedom (or five quadratic elements) per wavelength is used. Quadratic shape functions are used. The free-field transfer function is simulated from 20 Hz to 1 kHz with steps of 2 Hz.

For each speaker, two sizes are considered, which are referred to as Size 1 and Size 2, where Size 2 is the smallest. The sound pressure field is simulated on a transparent sphere around the speaker for each size and shape. Note that, for a given speaker shape, the radius of the transparent sphere changes depending on the size of the speaker and the position of the transducer on the speaker. These changes occur because we place the transducer at the center of the coordinate system. For the source on the speakers of Size 1, the radius of the transparent sphere is 0.4 m, while for Size 2, it is 0.2 m. For the receivers on the speakers, a radius of 0.5 m is used for Size 1, and a radius of 0.25 m is used for Size 2.

The complex-valued pressure field simulated on the surface of the transparent sphere is used to compute the spherical harmonic directivity coefficients. For the loudspeaker directivity, this involves a straightforward implementation of spherical harmonic decomposition.⁵⁶ To obtain the receiver directivities, we use acoustic reciprocity⁵² (see Sec. III B).

An example of the directivity patterns is shown in Fig. 5, where the sound pressure level (SPL) and phase are plotted for the spherical and cuboidal speakers with Size 1. The vibrating piston faces the positive direction of the x axis, which is hereafter referred to as $+ x$ ⁠. The pressure field on the XY-plane is reconstructed using Eq. (16) using simulated directivity coefficients $C_{n, m}^{(s)}$ with maximum spherical harmonic order N = 5. It can be seen that both the SPL and phase plots show differences between the two selected speaker shapes, especially at angles that are in the shadow region of the speakers. One can also observe a significant phase jump in the shadow region when a cuboidal speaker is used, which is less prominent for a spherical speaker.

FIG. 5.

View large Download slide

Plots of sound pressure level and phase reconstructed using simulated speaker directivity coefficients of a large (Size 1) spherical and cuboidal speaker. The pressure field is reconstructed with maximum spherical harmonic order 5. The results are evaluated on the XY-plane at three frequencies.

2. Room transfer functions

A rectangular room with dimensions $L_{x} \times L_{y} \times L_{z} = 4 m \times 3 m \times 2.5 m$ is considered. The air in the room has a speed of sound of 343 m/s and a density of 1.2 kg/m³. A uniform frequency-independent normalized impedance of 18, which gives an absorption coefficient of approximately 0.2, is specified at the walls of the room. This value is chosen because it gives a reverberation time of approximately 0.29 s at the fundamental frequency of 43.88 Hz, which gives a Schroeder frequency of approximately 198 Hz. The reverberation time was estimated using the modal model described by Morse and Bolt.⁶⁶ Note that, in general, the surface impedance is a complex-valued frequency-dependent quantity. Our choice of a real-valued frequency-independent impedance does not limit the generality of the method described here.

Five source–receiver configurations are simulated. The center position of each transducer, for each configuration, is given in Table II. Four of the configurations contain two devices—one with a loudspeaker and the other with a microphone. The remaining Configuration 3 has a single device on which the loudspeaker and microphone are placed. The transfer functions of the direct path between the source and receiver are simulated using the FEM for Configuration 3. The microphone position $x_{r}^{C 3}$ in Table II depends on the size, shape, and rotation of the speaker and takes the values given in Table III. A two-dimensional (2D) illustration of the speaker position and facing configurations is shown in Fig. 6. The source and receiver are both at least 1 m away from the walls in Positions 1–3. Positions 4–5 allow a closer distance between the walls and at least one speaker to examine the performance of the DEISM in challenging configurations.

TABLE I.

Parameters and descriptions used in ISM.

Parameter	Descriptions
$L_{x}, L_{y}, L_{z}$	Room dimensions: length, width, and height
$β_{x 1}, β_{x 2}, β_{y 1}, β_{y 2}, β_{z 1}, β_{z 2}$	Reflection coefficients of the walls
$p = (p_{x}, p_{y}, p_{z}) \in P$	Triple determining whether the image source is mirrored
	With respect to walls at x = 0, y = 0, or z = 0, elements take values of 0 or 1
$q = (q_{x}, q_{y}, q_{z}) \in Q$	Triple determining higher order reflections
	Each element can take values between $[- N_{m}, N_{m}]$
$x_{p, q}^{sI}$	Source images
$x_{p, q}^{rI}$	Receiver images
$R_{p, q}^{sI \to r} = x_{r} - x_{p, q}^{sI}$	Vector pointing from source images to the receiver
$R_{p, q}^{r \to sI} = - R_{p, q}^{sI \to r}$	Vector pointing from receiver to source images
$R_{\tilde{p}, \tilde{q}}^{s \to rI}$	The reversed vector corresponding to the same reflection path of $R_{p, q}^{r \to sI}$

TABLE II.

Configurations of the transducers and the corresponding speaker orientations used in the room simulation. The position vector $x_{r}^{C 3}$ for Configuration 3 is listed for different shapes and sizes in TABLE III.

Configuration	Source, orientation	Receiver, orientation
1	[1.1, 1.1, 1.3], $+ x$	[2.9, 1.9, 1.3], – x
2	[1.1, 1.1, 1.3], $+ x$	[1.9, 1.6, 1.4], – x
3	[1.1, 1.1, 1.3], $+ x$	$x_{r}^{C 3}, + x$
4	[0.4, 1.1, 1.3], $+ x$	[2.1, 1.6, 1.3], – x
5	[0.4, 1.1, 1.3], $+ x$	[2.5, 2.6, 1.3], – y

Configuration	Source, orientation	Receiver, orientation
1	[1.1, 1.1, 1.3], $+ x$	[2.9, 1.9, 1.3], – x
2	[1.1, 1.1, 1.3], $+ x$	[1.9, 1.6, 1.4], – x
3	[1.1, 1.1, 1.3], $+ x$	$x_{r}^{C 3}, + x$
4	[0.4, 1.1, 1.3], $+ x$	[2.1, 1.6, 1.3], – x
5	[0.4, 1.1, 1.3], $+ x$	[2.5, 2.6, 1.3], – y

TABLE III.

Positions $x_{r}^{C 3}$ for simulations with the source and receiver on the same speaker (Configuration 3).

Speaker shape	Size 1 (large)	Size 2 (small)
Cuboid	[1.05, 1.1, 1.5]	[1.075, 1.1, 1.4]
Spherical	[1.05, 1.1, 1.473]	[1.075, 1.1, 1.387]
Cylindrical	[1.05, 1.1, 1.5]	[1.075, 1.1, 1.4]

Speaker shape	Size 1 (large)	Size 2 (small)
Cuboid	[1.05, 1.1, 1.5]	[1.075, 1.1, 1.4]
Spherical	[1.05, 1.1, 1.473]	[1.075, 1.1, 1.387]
Cylindrical	[1.05, 1.1, 1.5]	[1.075, 1.1, 1.4]

FIG. 6.

View large Download slide

A 2D view of the positions of the mounted transducers inside the room. The black dots and gray dots represent the approximate positions of the source and receiver transducers, respectively. The arrows of the dots indicate the directions in which the speakers face, and the numbers indicate the source and receiver positions.

Each model is meshed using tetrahedral elements with an average size that ensures that 10 degrees of freedom per wavelength are used. Quadratic shape functions are used. The RTFs for each configuration are simulated, from 20 Hz to 1 kHz, in steps of 2 Hz.

B. Preliminary comparisons with ISM and FSRR

In this section, we compare the DEISM-LC to existing RTF simulation methods. To this end, we choose the original ISM with omnidirectional transducers (referred to as ISM-Omni) and the FSRR method, as introduced in Eqs. and , respectively. We demonstrate the effects of incorporating transducer directivities into the RTFs by comparing the ISM-Omni with the DEISM-LC. We are also interested in the differences between the FSRR method and DEISM-LC as they are computationally more efficient compared to the DEISM. In Fig. 7, the magnitudes and phases of the RTFs using the ISM-Omni, DEISM-LC, FSRR, and FEM are plotted. As an example, the results are shown for large cuboidal speakers in Configuration 2 with maximum reflection order 25 and maximum spherical harmonic order 5 for both the source and receiver. The magnitudes of the RTFs are given in terms of SPL in decibels and the unwrapped phases in radians. Unless specified, the magnitude and phase responses use the same units in the subsequent sections. The frequency-dependent SPL is defined as⁵⁴

SPL (k) = 20 {log}_{10} (\frac{P_{rms} (k)}{P_{0}}),

(37)

where

P_{rms}

is the root mean square (RMS) pressure calculated from the RTFs and

P_{0} = 2 \times 10^{- 5} Pa

is the reference pressure. The FEM solutions are the reference RTFs in this section and throughout the following evaluations.

FIG. 7.

View large Download slide

A comparison of the (unshifted) magnitudes and phase responses between the ISM-Omni, DEISM-LC, FSRR, and FEM solutions. The results are shown for large cuboidal speakers with Configuration 2. Max. reflection order is 25 and max. SH orders are 5 for both the source and receiver.

The RTFs of the ISM-Omni were generated using an omnidirectional source and receiver at the same locations as those used in Configuration 2. The directivity coefficients of the source and receiver are the same and contain only a monopole component based on Eq. , i.e.,

C_{0, 0}^{(s)} (k) = C_{0, 0}^{(r)} (k) = - i k j_{0} (0) Y_{0, 0}^{*} (θ, ϕ) = \frac{- i k}{\sqrt{4 π}},

(38)

where

(θ, ϕ)

denotes the directions in spherical coordinates of either the point source or the omnidirectional receiver relative to their local origins. For the RTFs of the FSRR given in Eq. , the directivity coefficients

{\tilde{C}}_{n, m}^{(s)} (k), {\tilde{C}}_{n, m}^{(r)} (k)

are extrapolated from

C_{n, m}^{(s)} (k), C_{v, u}^{(r)} (k)

to the measurement sphere with a radius of 1 m, and we choose

M_{p, q} = β (p, q)

⁠.

Overall, the differences between the DEISM-LC and FEM are small, compared to those between the ISM-Omni and FEM and those between the FSRR and FEM. This indicates that the DEISM-LC can generate more realistic RTFs when compared to the other image source methods for the considered scenario. The large offset in magnitude responses between the DEISM-LC and ISM-Omni mainly come from the different source types, i.e., a point source used in ISM-Omni and a vibrating piston used in DEISM-LC. However, the locations and Q-factors of the resonant peaks and the phase responses also vary greatly. This shows that the directivities of the source and receiver, including sound diffraction around the speakers, can affect the RTFs to a large extent.

The differences between the DEISM-LC and FSRR are mainly attributed to the inherent mismatch in their RTF parameterizations, i.e., between Eqs. (36) and (15). The magnitude responses are different, not only in the amplitude, but also in the locations and Q-factors of the resonant peaks. The phase response of the FSRR even contains positive values at lower frequencies corresponding to non-causal behaviors, which results from the random sign $B_{p, q}$ used in Eq. (15).

Since the inherent differences between DEISM-LC and other methods are shown to be large, i.e., the ISM-Omni and FSRR, we limit the following evaluations to DEISM, DEISM-LC, and FEM.

C. Varying distance between source and receiver

In this section, we present a comparison between the FEM and DEISM solutions for the small spherical speaker when the distance between the source and receiver is changing, corresponding to Configurations 1–3 in Table II. The two distances in Configurations 1 and 2 are around 1.97 and 0.95 m, and the distances in Configuration 3 are around 0.1 and 0.2 m for the small and large speakers, respectively. The magnitude and phase are compared in Fig. 8. To avoid visual overlap of the curves for different configurations, the SPLs and phases are plotted by shifting the curves vertically, as denoted by the arrows.

FIG. 8.

View large Download slide

Magnitudes and phase responses for small spherical speakers with varying distance between speakers. Max. reflection order is 25 and max. SH orders are 5 for both the source and receiver. The magnitudes and phases are shifted by vertical arrows for better distinction of the curves. Config., configuration.

We see that good agreement is obtained between the FEM and DEISM solutions. In general, greater differences are found at the notches of the RTFs, than at the peaks. Additionally, the differences between DEISM and DEISM-LC solutions are much smaller than those between the FEM and DEISM solutions in the notches of the RTFs. Since Configuration 3 contains only one speaker in the room, the corresponding RTFs show fewer variations across the frequency range when compared to the other two configurations. Similar agreement between FEM and DEISM was found for the other speaker shapes. For brevity, the respective plots are not presented here.

D. Source and receiver on the same device

Next, we compare solutions for cases where the source and receiver are placed on the same device (Configuration 3). For this comparison, we analyze data for the three speaker shapes considered.

The respective magnitude and phase solutions are shown in Fig. 9. Once again, good agreement is obtained between the FEM and DEISM solutions. The magnitudes and phases of the RTFs of the cuboidal and cylindrical speakers show very similar behaviour. In contrast, those of the spherical speaker exhibit clear differences when compared to the RTFs of the other speakers. Note that the main differences appear at the notches of the RTF magnitudes, which correspond to frequencies around 444 and 650 Hz. In addition, the phase curves of the cuboidal and cylindrical speakers also exhibit significant jumps at some frequencies, which cannot be observed for the spherical speaker. The differences between the curves indicate that the edges of the cuboidal and cylindrical speakers can induce additional diffraction and lead to distinct characteristics in the RTFs.

FIG. 9.

View large Download slide

Magnitudes and phase responses for large speakers with Configuration 3 while varying shapes of the speakers. Max. reflection order is 25 and max. SH orders are 5 for both the source and receiver. The magnitudes and phases are shifted by vertical arrows for better distinction of the curves.

E. Quantitative analysis

In this section, we quantify the differences between the FEM and DEISM/DEISM-LC solutions and DEISM and DEISM-LC solutions. The following error metrics are used for the comparison: The root mean square log spectral distance (LSD)⁶⁷

e_{lsd} = \sqrt{\frac{1}{K} \sum_{k \in K} | 20 {log}_{10} (\frac{| H_{Test} (k) |}{| H_{FEM} (k) |}) |^{2}},

(39)

and the normalized root mean square (NRMS) phase error in radians⁴⁸

e_{p} = \sqrt{\frac{1}{K} \sum_{k \in K} (∠ H_{Test} (k) - ∠ H_{FEM} (k))^{2}},

(40)

where

K

is the set of the K wavenumbers used in simulations,

H_{Test} (k)

is the simulated RTF of either the DEISM or DEISM-LC at wavenumber k, and

H_{FEM} (k)

is the simulated RTF of the FEM at wavenumber k. The log spectral distance and the phase errors are presented for different positions, speaker shapes, and sizes in Tables IV and V, respectively. In both tables, lower values indicate better performance. General observations can be made for Positions 1–3 in the two tables:

The smaller the speaker, the smaller the difference between the DEISM and FEM solutions.
The fewer speakers, i.e., for Configuration 3, the smaller the difference.

TABLE IV.

Root mean square log spectral distance in decibel between DEISM and FEM, and between DEISM-LC and FEM. S., small. Cub., sph., and cyl. denote cuboidal, spherical, and cylindrical, respectively.

Configuration	Cub.	S. cub.	Sph.	S. sph.	Cyl.	S. cyl.
1 DEISM	1.918	1.214	2.268	1.143	2.062	0.950
1 DEISM-LC	1.669	1.070	2.206	1.125	2.007	0.916
2 DEISM	1.974	1.051	2.521	1.003	1.900	0.964
2 DEISM-LC	1.855	1.003	2.060	0.984	1.809	0.943
3 DEISM	0.423	0.204	0.841	0.170	0.920	0.154
3 DEISM-LC	0.427	0.206	0.883	0.174	0.952	0.155
4 DEISM	7.018	8.187	6.457	7.418	6.350	7.287
4 DEISM-LC	6.864	8.149	6.461	7.298	6.240	7.211
5 DEISM	6.887	6.810	7.077	6.309	7.011	6.425
5 DEISM-LC	6.792	6.807	6.986	6.277	6.968	6.397

Configuration	Cub.	S. cub.	Sph.	S. sph.	Cyl.	S. cyl.
1 DEISM	1.918	1.214	2.268	1.143	2.062	0.950
1 DEISM-LC	1.669	1.070	2.206	1.125	2.007	0.916
2 DEISM	1.974	1.051	2.521	1.003	1.900	0.964
2 DEISM-LC	1.855	1.003	2.060	0.984	1.809	0.943
3 DEISM	0.423	0.204	0.841	0.170	0.920	0.154
3 DEISM-LC	0.427	0.206	0.883	0.174	0.952	0.155
4 DEISM	7.018	8.187	6.457	7.418	6.350	7.287
4 DEISM-LC	6.864	8.149	6.461	7.298	6.240	7.211
5 DEISM	6.887	6.810	7.077	6.309	7.011	6.425
5 DEISM-LC	6.792	6.807	6.986	6.277	6.968	6.397

TABLE V.

Normalized root mean square phase error in radians between DEISM and FEM, and between DEISM-LC and FEM. S., denotes small. Cub., sph., and cyl. denote cuboidal, spherical, and cylindrical, respectively.

Configuration	Cub.	S. cub.	Sph.	S. sph.	Cyl.	S. cyl.
1 DEISM	0.258	0.179	0.307	0.174	0.278	0.162
1 DEISM-LC	0.265	0.153	0.300	0.176	0.267	0.181
2 DEISM	0.251	0.138	0.246	0.128	0.255	0.123
2 DEISM-LC	0.246	0.149	0.259	0.142	0.260	0.144
3 DEISM	0.051	0.026	0.104	0.032	0.093	0.024
3 DEISM-LC	0.052	0.026	0.109	0.032	0.097	0.024
4 DEISM	1.796	1.764	1.834	1.776	1.801	1.778
4 DEISM-LC	1.791	1.763	1.828	1.774	1.781	1.776
5 DEISM	1.795	1.926	1.818	1.906	1.822	1.921
5 DEISM-LC	1.791	1.928	1.816	1.907	1.821	1.923

Configuration	Cub.	S. cub.	Sph.	S. sph.	Cyl.	S. cyl.
1 DEISM	0.258	0.179	0.307	0.174	0.278	0.162
1 DEISM-LC	0.265	0.153	0.300	0.176	0.267	0.181
2 DEISM	0.251	0.138	0.246	0.128	0.255	0.123
2 DEISM-LC	0.246	0.149	0.259	0.142	0.260	0.144
3 DEISM	0.051	0.026	0.104	0.032	0.093	0.024
3 DEISM-LC	0.052	0.026	0.109	0.032	0.097	0.024
4 DEISM	1.796	1.764	1.834	1.776	1.801	1.778
4 DEISM-LC	1.791	1.763	1.828	1.774	1.781	1.776
5 DEISM	1.795	1.926	1.818	1.906	1.822	1.921
5 DEISM-LC	1.791	1.928	1.816	1.907	1.821	1.923

The errors of the DEISM solutions increase when the speakers are located close to one wall of the room, i.e., for Configurations 4 and 5. The large errors may be attributed to the mirror source approximation since the speakers can have more complicated wave interaction with the nearby walls, which is not considered in this model, and can lead to additional RTF fluctuations. For devices with a physical extent, the assumption of a source image might not hold and could require more complex modeling.

In Configuration 3, the DEISM outperforms DEISM-LC for all speaker shapes and sizes. However, this does not necessarily hold for Configurations 1 and 2. This is mainly caused by the notches of the RTF magnitudes, where the DEISM produces much lower values than DEISM-LC, as shown in Figs. 8 and 9.

To provide insight into the performance of the methods in larger rooms, we consider the case in which the room's length in the x-dimension is doubled. For this test, only the cylindrical driver is considered in Configurations 1–3. The error levels obtained are presented in Table VI. The error levels for the larger room are similar to the values given in Tables IV and V, and imply similar conclusions.

TABLE VI.

Root mean square log spectral distance in decibel and normalized root mean square phase error in radians between DEISM and FEM, and between DEISM-LC and FEM in a larger room. S., small.

	RMS LSD (decibel)		NRMS phase error (radians)
Configuration	Cyl.	S. cyl.	Cyl.	S. cyl.
1 DEISM	1.177	0.967	0.185	0.128
1 DEISM-LC	1.190	1.036	0.227	0.126
2 DEISM	0.987	0.928	0.123	0.131
2 DEISM-LC	1.188	0.972	0.153	0.119
3 DEISM	0.783	0.319	0.151	0.038
3 DEISM-LC	0.893	0.320	0.155	0.038

	RMS LSD (decibel)		NRMS phase error (radians)
Configuration	Cyl.	S. cyl.	Cyl.	S. cyl.
1 DEISM	1.177	0.967	0.185	0.128
1 DEISM-LC	1.190	1.036	0.227	0.126
2 DEISM	0.987	0.928	0.123	0.131
2 DEISM-LC	1.188	0.972	0.153	0.119
3 DEISM	0.783	0.319	0.151	0.038
3 DEISM-LC	0.893	0.320	0.155	0.038

F. Evaluation of DEISM-LC

In Sec. IV E, we showed that in some cases, DEISM-LC achieves a smaller difference to the FEM than DEISM, especially at the notches of the RTFs. It is worth noting that the frequencies of the notches might not be accurate since the maximum resolution of the frequencies is 2 Hz. Therefore, directly comparing the DEISMs with the FEM might not be enough to determine whether the simplified DEISM performs worse or better. In the following, we investigate the difference between the DEISMs by computing the relative error using the

ℓ^{2}

norm, defined as

e_{ℓ^{2}} = \frac{| | H_{DEISM} - H_{DEISM - LC} | |_{2}}{| | H_{DEISM} | |_{2}},

(41)

where

H_{(\cdot)}

is the complex-valued vector with the transfer function values

H_{(\cdot)} (k)

of all wavenumbers. In this comparison, we consider the error in free-field and reverberant conditions.

The errors in free-field conditions are given in Fig. 10, where the source–receiver distance and the error are plotted on logarithmic scales. One can observe that the difference decreases to around 1% at a distance of 10 m and to around 0.1% at a distance of 100 m. It can be concluded that DEISM-LC asymptotically converges to the DEISM solution as the distance between the source and receiver increases, which is expected from the derivations in Sec. III D.

FIG. 10.

View large Download slide

Relative error [Eq. (41)] between the DEISM and DEISM-LC solutions evaluated for the direct path only while changing the distance between the source and receiver. Both the distance and relative errors are plotted on a logarithmic scale.

In Fig. 11, the relative error is evaluated by changing the maximum reflection order in the DEISM for the cuboidal speakers. Note that for the curves of Configuration 3, the data of reflection order zero is missing since the direct path is simulated using the FEM. The difference between the DEISM solutions reduces for the first few reflection orders, but remains constant for higher reflection orders. Similar behaviors were also found for other speaker shapes, which are not presented here for brevity.

FIG. 11.

View large Download slide

Relative error [Eq. (41)] between the DEISM and DEISM-LC solutions as a function of reflection order for the cuboidal speakers at Positions 1, 2, and 3.

G. Summary

The presented comparisons and error analyses show that the proposed method and its simplified version can simulate RTFs that are in good agreement with FEM solutions in Configurations 1–3. Less agreement is found when a speaker is placed close to the walls, i.e., in Configurations 4 and 5. This might be attributed to the comb filtering caused by the direct path and reflection between the wall and speaker when the distance is less than 0.5 m.⁶⁸

The errors between the DEISM and FEM solutions decrease when smaller speakers are used in the simulation. Additionally, the lowest error is obtained when the source and receiver are placed on the same speaker, which reduces the number of speakers in the simulation, and hence, the number of reflecting surfaces. These observations support the hypothesis that the errors come mostly from the omission of speakers with physical extent in the ISM model, i.e., scattering effects between the walls and the speaker's surfaces are not included.

We have also shown that the simplified version (DEISM-LC) converges to its complete version (DEISM) in the free field and yields a small bias when used to simulate RTFs. The simplified version, DEISM-LC, yields a better agreement with the FEM simulations as opposed to the existing FSRR method.

In addition, the computation of a single RTF using FEM took approximately 2 days on a desktop computer with 256 GB of random access memory (RAM), while the computation of the same RTF using DEISM took approximately 5 min on a laptop with 16 GB of RAM. Moreover, DEISM-LC was around ten times faster than DEISM for the scenarios considered in this paper. This reduction of the time needed to compute RTFs can facilitate applications that require large amounts of data, e.g., the generation of training data for data-driven approaches.

V. CONCLUSION

A RTF simulation method, which makes use of spherical harmonic directivity coefficients, and which is based on the image source method, has been presented. The directivity coefficients contain information about the sound diffraction around transducers mounted on audio devices, and the method presented aims to capture this effect. The proposed method (DEISM) and its simplified version (DEISM-LC) were compared with FEM solutions and evaluated for different acoustic scenarios. The results show that the proposed approaches can simulate RTFs with good agreement to FEM simulations for shoebox-shaped rooms when devices are not placed close to the room boundaries. A rule of thumb for indicating the regions of the room in which the DEISM provides results that are significantly different from FEM results is a topic for future studies.

Future work could include extending the proposed method to incorporate the sound reflections between the enclosures of multiple devices, or between the enclosures and room boundaries. Such extensions might reduce the error incurred when multiple devices are used or when devices are placed close to the room boundaries. Another valuable extension of the proposed approach would be to apply the DEISM and DEISM-LC to rooms with more general geometries, with the aim of efficiently generating more realistic room acoustics simulations.

REFERENCES

1.

M. R.

Schroeder

, “

Novel uses of digital computers in room acoustics

,”

J. Acoust. Soc. Am.

33

(

11

),

1669

(

1961

).

https://doi.org/10.1121/1.1936681

Google Scholar

2.

M.

Kleiner

,

B.-I.

Dalenbäck

, and

P.

Svensson

, “

Auralization-an overview

,”

J. Audio Eng. Soc.

41

(

11

),

861

–

875

(

1993

), available at https://www.aes.org/e-lib/browse.cfm?elib=6976.

Google Scholar

3.

L.

Savioja

, “

Modeling techniques for virtual acoustics

,” Ph.D. thesis,

Helsinki University of Technology

,

Espoo, Finland

,

1999

.

Google Scholar

4.

D.

Schröder

, “

Physically based real-time auralization of interactive virtual environments

,” Ph.D. thesis,

RWTH Aachen

,

Aachen, Germany

,

2011

.

Google Scholar

5.

M.

Vorländer

, “

Computer simulations in room acoustics: Concepts and uncertainties

,”

J. Acoust. Soc. Am.

133

(

3

),

1203

–

1213

(

2013

).

https://doi.org/10.1121/1.4788978

Google Scholar

6.

C.

Schissler

and

D.

Manocha

, “

Interactive sound propagation and rendering for large multi-source scenes

,”

ACM Trans. Graph. (TOG)

36

(

4

),

2:1

–

2:12

(

2016

).

https://doi.org/10.1145/3072959.2943779

Google Scholar

7.

M.

Vorländer

,

Auralization: Fundamentals of Acoustics, Modelling, Simulation, Algorithms and Acoustic Virtual Reality

, 2nd ed. (

Springer

,

Berlin/Heidelberg, Germany

,

2020

).

Google Scholar

8.

Z.

Tang

,

H.-Y.

Meng

, and

D.

Manocha

, “

Learning acoustic scattering fields for dynamic interactive sound propagation

,” in

IEEE Virtual Reality and 3D User Interfaces (VR)

,

Lisboa

,

Portugal

(

2021

), pp.

835

–

844

.

Google Scholar

9.

E.

Mabande

,

K.

Kowalczyk

,

H.

Sun

, and

W.

Kellermann

, “

Room geometry inference based on spherical microphone array eigenbeam processing

,”

J. Acoust. Soc. Am.

134

(

4

),

2773

–

2789

(

2013

).

https://doi.org/10.1121/1.4820895

Google Scholar

10.

D.

Cherkassky

and

S.

Gannot

, “

Blind synchronization in wireless acoustic sensor networks

,”

IEEE/ACM Trans. Audio. Speech. Lang. Process.

25

(

3

),

651

–

661

(

2017

).

https://doi.org/10.1109/TASLP.2017.2655259

Google Scholar

11.

C.

Evers

and

P. A.

Naylor

, “

Acoustic slam

,”

IEEE/ACM Trans. Audio. Speech. Lang. Process.

26

(

9

),

1484

–

1498

(

2018

).

https://doi.org/10.1109/TASLP.2018.2828321

Google Scholar

12.

S.

Das

and

T.

Bäckström

, “

Enhancement by postfiltering for speech and audio coding in ad hoc sensor networks

,”

JASA Express Lett.

1

(

1

),

015206

(

2021

).

https://doi.org/10.1121/10.0003208

Google Scholar

13.

A.

Herzog

and

E. A. P.

Habets

, “

Distance estimation in the spherical harmonic domain using the spherical wave model

,”

Appl. Acoust.

193

,

108733

(

2022

).

https://doi.org/10.1016/j.apacoust.2022.108733

Google Scholar

14.

D.

Diaz-Guerra

,

A.

Miguel

, and

J. R.

Beltran

, “

Robust sound source tracking using SRP-PHAT and 3D convolutional neural networks

,”

IEEE/ACM Trans. Audio. Speech. Lang. Process.

29

,

300

–

311

(

2021

).

https://doi.org/10.1109/TASLP.2020.3040031

Google Scholar

15.

F. B.

Gelderblom

,

Y.

Liu

,

J.

Kvam

, and

T. A.

Myrvoll

, “

Synthetic data for dnn-based doa estimation of indoor speech

,” in

IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)

,

Toronto, Ontario, Canada

(

2021

), pp.

4390

–

4394

.

Google Scholar

16.

F.

Hübner

,

W.

Mack

, and

E. A. P.

Habets

, “

Efficient training data generation for phase-based DOA estimation

,” in

IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)

, Toronto, Canada (

2021

), pp.

456

–

460

.

Google Scholar

17.

J.

Wechsler

,

W.

Mack

, and

E. A. P.

Habets

, “

End-to-end signal-aware direction-of-arrival estimation using weighted steered-response power

,”

European Signal Processing Conference (EUSIPCO)

, Belgrade, Serbia (

2022

), pp.

41

–

45

.

Google Scholar

18.

P.

Srivastava

,

A.

Deleforge

, and

E.

Vincent

, “

Realistic sources, receivers and walls improve the generalisability of virtually-supervised blind acoustic parameter estimators

,” in

International Workshop on Acoustic Signal Enhancement (IWAENC)

,

Bamberg

,

Germany

(

2022

), pp.

1

–

5

.

Google Scholar

19.

P.

Srivastava

,

A.

Deleforge

,

A.

Politis

, and

E.

Vincent

, “

How to (virtually) train your speaker localizer

,” in Proceedings of Interspeech, Dublin, Ireland (

2023

), pp.

1204

–

1208

.

Google Scholar

20.

A.

Herzog

,

S. R.

Chetupalli

, and

E. A. P.

Habets

, “

AmbiSep: Ambisonic-to-ambisonic reverberant speech separation using transformer networks

,” in

International Workshop on Acoustic Signal Enhancement (IWAENC)

,

Bamberg

,

Germany

(

2022

), pp.

1

–

5

.

Google Scholar

21.

P.-A.

Grumiaux

,

S.

Kitić

,

L.

Girin

, and

A.

Guérin

, “

A survey of sound source localization with deep learning methods

,”

J. Acoust. Soc. Am.

152

(

1

),

107

–

151

(

2022

).

https://doi.org/10.1121/10.0011809

Google Scholar

22.

J. B.

Allen

and

D. A.

Berkley

, “

Image method for efficiently simulating small-room acoustics

,”

J. Acoust. Soc. Am.

65

(

4

),

943

–

950

(

1979

).

https://doi.org/10.1121/1.382599

Google Scholar

23.

P. M.

Peterson

, “

Simulating the response of multiple microphones to a single acoustic source in a reverberant room

,”

J. Acoust. Soc. Am.

80

(

5

),

1527

–

1529

(

1986

).

https://doi.org/10.1121/1.394357

Google Scholar

24.

J. C. B.

Torres

,

L.

Aspöck

, and

M.

Vorländer

, “

Comparative study of two geometrical acoustic simulation models

,”

J. Braz. Soc. Mech. Sci. Eng.

40

(

6

),

1

–

15

(

2018

).

https://doi.org/10.1007/s40430-018-1226-1

Google Scholar

25.

Z.

Tang

,

L.

Chen

,

B.

Wu

,

D.

Yu

, and

D.

Manocha

, “

Improving reverberant speech training using diffuse acoustic simulation

,” in

IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)

, Barcelona, Spain (

2020

), pp.

6969

–

6973

.

Google Scholar

26.

E.

Bezzam

,

R.

Scheibler

,

C.

Cadoux

, and

T.

Gisselbrecht

, “

A study on more realistic room simulation for far-field keyword spotting

,” in

Asia-Pacific Signal and Information Processing Association Annual Summit and Conference (APSIPA ASC)

, Auckland, New Zealand (

2020

), pp.

674

–

680

.

Google Scholar

27.

D. T.

Murphy

, “

Digital waveguide mesh topologies in room acoustics modelling

,” Ph.D. thesis,

University of York

,

York, UK

,

2000

.

Google Scholar

28.

K.

Kowalczyk

, “

Boundary and medium modelling using compact finite difference schemes in simulations of room acoustics for audio and architectural design applications

,” Ph.D. thesis,

Queen's University Belfast

,

Belfast, UK

,

2010

.

Google Scholar

29.

T.

Sakuma

,

S.

Sakamoto

, and

T.

Otsuru

,

Computational Simulation in Architectural and Environmental Acoustics

(

Springer, Tokyo,

2014

).

Google Scholar

30.

B.

Hamilton

and

S.

Bilbao

, “

FDTD methods for 3-D room acoustics simulation with high-order accuracy in space and time

,”

IEEE/ACM Trans. Audio. Speech. Lang. Process.

25

(

11

),

2112

–

2124

(

2017

).

https://doi.org/10.1109/TASLP.2017.2744799

Google Scholar

31.

J. A.

Hargreaves

,

L. R.

Rendell

, and

Y. W.

Lam

, “

A framework for auralization of boundary element method simulations including source and receiver directivity

,”

J. Acoust. Soc. Am.

145

(

4

),

2625

–

2637

(

2019

).

https://doi.org/10.1121/1.5096171

Google Scholar

32.

T.

Yoshida

,

T.

Okuzono

, and

K.

Sakagami

, “

A parallel dissipation-free and dispersion-optimized explicit time-domain FEM for large-scale room acoustics simulation

,”

Buildings

12

(

2

),

105

(

2022

).

https://doi.org/10.3390/buildings12020105

Google Scholar

33.

M.

Aretz

, “

Combined wave and ray based room acoustic simulations of small rooms

,” Ph.D. thesis,

RWTH Aachen

,

Aachen, Germany

,

2012

.

Google Scholar

34.

M. R.

Thomas

, “

Wayverb: A graphical tool for hybrid room acoustics simulation

,” Ph.D. thesis,

University of Huddersfield

,

Huddersfield, UK

,

2017

.

Google Scholar

35.

A.

Ratnarajah

,

Z.

Tang

, and

D.

Manocha

, “

TS-RIR: Translated synthetic room impulse responses for speech augmentation

,” in

IEEE Automatic Speech Recognition and Understanding Workshop (ASRU)

,

Cartagena

,

Columbia

(

2021

), pp.

259

–

266

.

Google Scholar

36.

A.

Luo

,

Y.

Du

,

M. J.

Tarr

,

J. B.

Tenenbaum

,

A.

Torralba

, and

C.

Gan

, “

Learning neural acoustic fields

,” in

Conference on Neural Information Processing Systems (NeurIPS)

,

New Orleans, LA

(

2022

), available at https://openreview.net/pdf?id=D21DRzkZbSB.

Google Scholar

37.

J.

Borish

, “

Extension of the image model to arbitrary polyhedra

,”

J. Acoust. Soc. Am.

75

(

6

),

1827

–

1836

(

1984

).

https://doi.org/10.1121/1.390983

Google Scholar

38.

R.

Scheibler

,

E.

Bezzam

, and

I.

Dokmanić

, “

Pyroomacoustics: A Python package for audio room simulation and array processing algorithms

,” in

IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)

,

Calgary

,

Alberta

,

Canada

,

351

–

355

(

2018

).

Google Scholar

39.

M.

Aretz

,

P.

Dietrich

, and

M.

Vorländer

, “

Application of the mirror source method for low frequency sound prediction in rectangular rooms

,”

Acta Acust. united Acust.

100

(

2

),

306

–

319

(

2014

).

https://doi.org/10.3813/AAA.918710

Google Scholar

40.

M. R.

Schroeder

and

K. H.

Kuttruff

, “

On frequency response curves in rooms. Comparison of experimental, theoretical, and Monte Carlo results for the average frequency spacing between maxima

,”

J. Acoust. Soc. Am.

34

(

1

),

76

–

80

(

1962

).

https://doi.org/10.1121/1.1909022

Google Scholar

41.

E. A. P.

Habets

, “

Room impulse response generator

,”

Technical report

(

2000

)

1

–

21

.

https://doi.org/10.5281/zenodo.4096348

Google Scholar

42.

A.

Wabnitz

,

N.

Epain

,

C.

Jin

, and

A.

Van Schaik

, “

Room acoustics simulation for multichannel microphone arrays

,” in

Proceedings of the International Symposium on Room Acoustics

, Melbourne, Australia,

Citeseer

(

2010

), pp.

1

–

6

.

Google Scholar

43.

D.

Diaz-Guerra

,

A.

Miguel

, and

J. R.

Beltran

, “

gpuRIR: A python library for room impulse response simulation with GPU acceleration

,”

Multimed. Tools Appl.

80

(

4

),

5653

–

5671

(

2021

).

https://doi.org/10.1007/s11042-020-09905-3

Google Scholar

44.

Z.-h.

Fu

and

J.-w.

Li

, “

GPU-based image method for room impulse response calculation

,”

Multimed. Tools Appl.

75

(

9

),

5205

–

5221

(

2016

).

https://doi.org/10.1007/s11042-015-2943-4

Google Scholar

45.

M.

Kompis

and

N.

Dillier

, “

Simulating transfer functions in a reverberant room including source directivity and head-shadow effects

,”

J. Acoust. Soc. Am.

93

(

5

),

2779

–

2787

(

1993

).

https://doi.org/10.1121/1.405825

Google Scholar

46.

T.

Betlehem

and

M.

Poletti

, “

Sound field of a directional source in a reverberant room

,”

J. Acoust. Soc. N. Z.

25

(

4

),

12

–

22

(

2012

).

Google Scholar

47.

F.

Brinkmann

,

V.

Erbes

, and

S.

Weinzierl

, “

Extending the closed form image source model for source directivity

,” in

DAGA, München

(

2018

).

Google Scholar

48.

D. P.

Jarrett

,

E. A. P.

Habets

,

M. R. P.

Thomas

, and

P. A.

Naylor

, “

Rigid sphere room impulse response simulation: Algorithm and applications

,”

J. Acoust. Soc. Am.

132

(

3

),

1462

–

1472

(

2012

).

https://doi.org/10.1121/1.4740497

Google Scholar

49.

S.

Hafezi

,

A. H.

Moore

, and

P. A.

Naylor

, “

Modelling source directivity in room impulse response simulation for spherical microphone arrays

,” in

European Signal Processing Conference (EUSIPCO)

,

Nice

,

France

(

2015

), pp.

574

–

578

.

Google Scholar

50.

P. N.

Samarasinghe

,

T. D.

Abhayapala

,

Y.

Lu

,

H.

Chen

, and

G.

Dickins

, “

Spherical harmonics based generalized image source method for simulating room acoustics

,”

J. Acoust. Soc. Am.

144

(

3

),

1381

–

1391

(

2018

).

https://doi.org/10.1121/1.5053579

Google Scholar

51.

Y.

Luo

and

W.

Kim

, “

Fast source-room-receiver acoustics modeling

,” in

European Signal Processing Conference (EUSIPCO)

(

2021

), pp.

51

–

55

.

Google Scholar

52.

Z.

Xu

,

A.

Herzog

,

A.

Lodermeyer

,

E. A. P.

Habets

, and

A. G.

Prinn

, “

Acoustic reciprocity in the spherical harmonic domain: A formulation for directional sources and receivers

,”

JASA Express Lett.

2

,

124801

(

2022

).

https://doi.org/10.1121/10.0016542

Google Scholar

53.

H.

Carslaw

, “

Some multiform solutions of the partial differential equations of physical mathematics and their applications

,”

Proc. London Math. Soc.

s1-30

(

1

),

121

–

165

(

1898

), https://academic.oup.com/plms/article/s1-30/1/121/2917781 (Last viewed December 2023).

https://doi.org/10.1112/plms/s1-30.1.121

Google Scholar

54.

H.

Kuttruff

,

Room Acoustics

(

CRC Press

,

Boca Raton, FL

,

2017

).

Google Scholar

55.

D. P.

Jarrett

and

E. A. P.

Habets

, “

On the noise reduction performance of a spherical harmonic domain tradeoff beamformer

,”

IEEE Signal Process. Lett.

19

(

11

),

773

–

776

(

2012

).

https://doi.org/10.1109/LSP.2012.2220351

Google Scholar

56.

E. G.

Williams

,

Fourier Acoustics: Sound Radiation and Nearfield Acoustical Holography

, 1st ed. (

Academic Press

,

London

,

1999

).

Google Scholar

57.

W.

Zhang

,

T. D.

Abhayapala

,

R. A.

Kennedy

, and

R.

Duraiswami

, “

Insights into head-related transfer function: Spatial dimensionality and continuous representation

,”

J. Acoust. Soc. Am.

127

(

4

),

2347

–

2357

(

2010

).

https://doi.org/10.1121/1.3336399

Google Scholar

58.

P. N.

Samarasinghe

and

T. D.

Abhayapala

, “

Room transfer function measurement from a directional loudspeaker

,” in

IEEE International Workshop on Acoustic Signal Enhancement (IWAENC)

,

Xi'an

,

China

(

2016

), pp.

1

–

5

.

Google Scholar

59.

D.

Ward

and

T.

Abhayapala

, “

Reproduction of a plane-wave sound field using an array of loudspeakers

,”

IEEE Trans. Speech Audio Process.

9

(

6

),

697

–

707

(

2001

).

https://doi.org/10.1109/89.943347

Google Scholar

60.

B.

Rafaely

,

Fundamentals of Spherical Array Processing

(

Springer

,

New York

,

2015

), Vol. 8.

Google Scholar

61.

J.

Ahrens

and

S.

Bilbao

, “

Computation of spherical harmonics based sound source directivity models from sparse measurement data

,”

Forum Acusticum

(

2020

), pp.

2019

–

2026

.

Google Scholar

62.

J. G.

Tylka

,

R.

Sridhar

, and

E.

Choueiri

, “

A database of loudspeaker polar radiation measurements

,” in 139th Audio Engineering Society Convention, New York (

2015

), pp.

1

–

4

.

Google Scholar

63.

A. R.

Edmonds

,

Angular Momentum in Quantum Mechanics

(

Princeton University Press

,

Princeton, NJ

,

2016

).

Google Scholar

64.

See https://github.com/audiolabs/DEISM for the implementation of the algorithm described in Sec. III E.

65.

J.-P.

Berenger

, “

A perfectly matched layer for the absorption of electromagnetic waves

,”

J. Comput. Phys

114

(

2

),

185

–

200

(

1994

).

https://doi.org/10.1006/jcph.1994.1159

Google Scholar

66.

P. M.

Morse

and

R. H.

Bolt

, “

Sound waves in rooms

,”

Rev. Mod. Phys.

16

(

2

),

69

–

150

(

1944

).

https://doi.org/10.1103/RevModPhys.16.69

Google Scholar

67.

A.

Gray

and

J.

Markel

, “

Distance measures for speech processing

,”

IEEE Trans. Acoust, Speech, Signal Process.

24

(

5

),

380

–

391

(

1976

).

https://doi.org/10.1109/TASSP.1976.1162849

Google Scholar

68.

J.

Paasonen

,

A.

Karapetyan

,

J.

Plogsties

, and

V.

Pulkki

, “

Proximity of surfaces – Acoustic and perceptual effects

,”

J. Audio Eng. Soc.

65

(

12

),

997

–

1004

(

2017

).

https://doi.org/10.17743/jaes.2017.0039

Google Scholar

Simulating room transfer functions between transducers mounted on audio devices using a modified image source method

I. INTRODUCTION

II. THE IMAGE SOURCE METHOD

A. Standard image source method

B. Extensions of the standard ISM

1. Angle-dependent reflection coefficients

2. Extension to the spherical harmonic domain

C. Generalized image source method

D. Fast source room receiver method

III. PROPOSED METHOD

A. Problem formulation

B. Source and receiver directivities

C. Room transfer function

D. Far-field approximation

E. Algorithm summary

IV. EVALUATION OF THE PROPOSED METHODS

A. Test setup

1. Transducer directivities

2. Room transfer functions

B. Preliminary comparisons with ISM and FSRR

C. Varying distance between source and receiver

D. Source and receiver on the same device

E. Quantitative analysis

F. Evaluation of DEISM-LC

G. Summary

V. CONCLUSION

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Simulating room transfer functions between transducers mounted on audio devices using a modified image source method

I. INTRODUCTION

II. THE IMAGE SOURCE METHOD

A. Standard image source method

B. Extensions of the standard ISM

1. Angle-dependent reflection coefficients

2. Extension to the spherical harmonic domain

C. Generalized image source method

D. Fast source room receiver method

III. PROPOSED METHOD

A. Problem formulation

B. Source and receiver directivities

C. Room transfer function

D. Far-field approximation

E. Algorithm summary

IV. EVALUATION OF THE PROPOSED METHODS

A. Test setup

1. Transducer directivities

2. Room transfer functions

B. Preliminary comparisons with ISM and FSRR

C. Varying distance between source and receiver

D. Source and receiver on the same device

E. Quantitative analysis

F. Evaluation of DEISM-LC

G. Summary

V. CONCLUSION

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only