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

The ISM was first applied to acoustics by Carslaw [52] and later implemented using a digital computer by Allen and Berkley [22] 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 omidirectional receiver inside of the room at positions $\mathbf{x}_{\text{s}}=(x_{\text{s}},y_{\text{s}},z_{\text{s}})$ and $\mathbf{x}_{\text{r}}=(x_{\text{r}},y_{\text{r}},z_{\text{r}})$, given in Cartesian coordinates, respectively. The RTF can be expressed as a weighted sum of free-space Green’s functions attenuated by the walls, that is:

where $k$ is the wavenumber. The triples $\mathbf{p}$ and $\mathbf{q}$ determine the position of the images and the attenuation factor $\beta(\mathbf{p},\mathbf{q})$. The triple $\mathbf{p}=(p_{x},p_{y},p_{z})$ is used to determine whether the image source is mirrored w.r.t. walls at $x=0$, $y=0$ or $z=0$, and its elements take values of $0$ or $1$, resulting in a set $\mathcal{P}$ of eight combinations. Each element of the other triple $\mathbf{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 [47], resulting in a set $\mathcal{Q}$ of $(2N_{m}+1)^{3}$ combinations. $\beta(\mathbf{p},\mathbf{q})=\beta_{x_{1}}^{|q_{x}-p_{x}|}\beta_{x_{2}}^{|q_{x}|}\beta_{y_{1}}^{|q_{y}-p_{y}|}\beta_{y_{2}}^{|q_{y}|}\beta_{z_{1}}^{|q_{z}-p_{z}|}\beta_{z_{2}}^{|q_{z}|}$ and $\beta_{x_{1}},\beta_{x_{2}},\beta_{y_{1}},\beta_{y_{2}},\beta_{z_{1}},\beta_{z_{2}}$ are the angle-independent reflection coefficients of the walls, and the subscripts $a_{1},a_{2}$ ($a\in\{x,y,z\}$) of $\beta_{(\cdot)}$ correspond to the boundaries at $a=0$ and the boundaries at $a=L_{a}$, respectively. For an image source with position

the Green’s function can be written as

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

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.

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. [49] 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 $\mathbf{x}_{\text{s}}$ as shown in Fig. 1. The sound field observed at an arbitrary point $\mathbf{r}$ outside of the sphere enclosing the entire speaker can be computed as [55]:

where $(d,\theta,\phi)$ denote the spherical coordinates associated with $\mathbf{r}$, viz. the distance $d$, the inclination angle $\theta$ and the azimuth angle $\phi$. 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^{(\text{s})}_{n,m}(k)$ represent the spherical harmonic source directivity coefficients. In practice, the infinite summation in Eq.(10) is truncated to a maximum order $N=\left\lceil kr_{0}^{(\text{s})}\right\rceil$ where $r_{0}^{(\text{s})}$ denotes the radius of the sphere that is large enough to surround the loudspeaker [56, 57].

The room transfer function from the source $\mathbf{x}_{\text{s}}$ to any point $\mathbf{y}^{(\text{r})}$ with spherical coordinates $(d_{y}^{({\text{r}})},\theta_{y}^{({\text{r}})},\phi_{y}^{({\text{r}})})$ relative to $\mathbf{x}_{\text{r}}$ can be expressed using the GISM as [49]:

where $V=\left\lceil k\,d_{y}^{(r)}\right\rceil$ denotes the truncation limit that produces a small error [58] and $j_{v}(\cdot)$ denotes the $n$th-order spherical Bessel function of the first kind. The reverberant mode coupling coefficients $\gamma_{v,u}^{n,m}(k)$ are given by

where $\mathbf{R}^{\text{sI}\to\text{r}}_{\mathbf{p},\mathbf{q}}$ has spherical coordinates $(||\mathbf{R}^{\text{sI}\to\text{r}}_{\mathbf{p},\mathbf{q}}||,\Omega^{\text{sI}\to\text{r}}_{\mathbf{p},\mathbf{q}})$, and $\Omega^{\text{sI}\to\text{r}}_{\mathbf{p},\mathbf{q}}=(\theta^{\text{sI}\to\text{r}}_{\mathbf{p},\mathbf{q}},\phi^{\text{sI}\to\text{r}}_{\mathbf{p},\mathbf{q}})$. 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 [59] 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 [48]. Using $\mathbf{x}_{\mathbf{0}}$ with spherical coordinates $(d_{x_{0}},\theta_{x_{0}},\phi_{x_{0}})$ instead of $\mathbf{R}^{\text{sI}\to\text{r}}_{\mathbf{p},\mathbf{q}}$ for brevity, the single-path mode-coupling coefficients $\alpha_{v,u}^{n,m}(\cdot)$ in Eq.(12) can be derived using the translation of fields [59],

with $W_{1}$ and $W_{2}$ denoting Wigner 3 - j symbols,

and $\xi=\sqrt{(2n+1)(2v+1)(2l+1)/4\pi}$.

The Fast Source Room Receiver (FSRR) method, proposed by Luo and Kim [50], 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.4.

Using our notation, the FSRR method RTF is given by

where $B_{\mathbf{p},\mathbf{q}}$ randomly takes value between $1$ and $-1$, $M_{\mathbf{p},\mathbf{q}}$ is a function fitted to frequency-dependent absorption coefficients, $\tilde{C}$ is a directivity coefficient obtained at 1 m from either the source or the receiver. The modified triples $\tilde{\mathbf{p}},\tilde{\mathbf{q}}$ are defined in Sec. III.4 (Eq. (30)), $\Omega^{\text{sI}\to\text{r}}_{\tilde{\mathbf{p}},\tilde{\mathbf{q}}}$ is the direction of a vector pointing from a source image to the receiver, and $\Omega^{\text{r}\to\text{sI}}_{\tilde{\mathbf{p}},\tilde{\mathbf{q}}}$ 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.4.

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 $\mathbf{x}_{\text{s}}$, e.g., a vibrating piston on one surface, and a microphone located at $\mathbf{x}_{\text{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.[51], the terms $j_{v}(kd_{y}^{(\text{r})})Y_{v,u}(\theta_{y}^{(\text{r})},\phi_{y}^{(\text{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 [50] 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.

The directivity of a transducer mounted on a device can either be 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 $\mathbf{x}_{\text{s}}$ or a microphone at $\mathbf{x}_{\text{r}}$ as shown in Fig. 2, acoustic diffraction occurs due to the wave interaction with the speakers. As shown in previous works [55, 60, 51], the total sound radiation of a source device can be quantified using directivity coefficients defined on a transparent sphere with radius $r_{0}^{(\text{s})}$. Using reciprocity, the receiver directivity can be obtained analogously to the source directivity [51]. 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. (10) with truncated order $N$ as

where $\mathbf{r}_{j}$ denotes the the $j$th position on the transparent sphere, $(\theta_{j},\phi_{j})$ the corresponding direction, $P_{n,m}(r_{0}^{(\text{s})},k)=C^{(\text{s})}_{n,m}(k)h_{n}(kr_{0}^{(\text{s})})$ is the spherical wave spectrum [55] and $j\in\{1,2,\ldots,J\}$. Different choices of the distribution of $(\theta_{j},\phi_{j})$ can be found in the literature [59]. 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^{(\text{s})}_{n,m}(k)=P_{n,m}(r_{0}^{(\text{s})},k)/h_{n}(k\,r_{0}^{(\text{s})})$. The receiver directivity coefficients $C^{(\text{r})}_{n,m}(k)$ can be obtained following the same procedure by making use of reciprocity [51].

The following parameterization of the transfer function between two transducers with arbitrary directivities in the free field, $H^{(\text{F})}$, was proposed by Xu et al. [51]:

where $\alpha_{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^{(\text{r})}_{v,u}(k)$ and the directivity coefficients $C^{(\text{r})}_{v,u}(k)$ obtained via reciprocity for the receiver was derived in [51]

Since the summation over source images is a linear operation and the mirroring effects are already included in $\gamma_{v,u}^{n,m}$ in Eq. (12), the total RTF can be written by substituting $\gamma_{v,u}^{n,m}(k)$ for $\alpha_{v,u}^{n,m}(k)$ in Eq. (17), yielding:

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

The final proposed RTF is therefore expressed as

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. Similarly to Xu et al. [51], for a receiver region without any physical objects and an observation point $\mathbf{y}^{(\text{r})}$ with spherical coordinates $(d_{y}^{(\text{r})},\theta_{y}^{(\text{r})},\phi_{y}^{(\text{r})})$ relative to $\mathbf{x}_{\text{r}}$, as shown in Fig. 1, the receiver directivity coefficients can be written as

By substituting Eq. (22) into Eq. (21), and utilizing the property of spherical harmonics $Y_{v,u}^{*}(\theta_{y}^{(\text{r})},\phi_{y}^{(\text{r})})=(-1)^{u}Y_{v,-u}(\theta_{y}^{(\text{r})},\phi_{y}^{(\text{r})})$, one can obtain the same expression of the RTFs as the GISM shown in Eq. (11).

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 [61], 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 $\mathbf{R}^{\text{sI}\to\text{r}}_{\mathbf{p},\mathbf{q}}$, for the following analysis. The mode coupling coefficients from Eq. (20) for a single reflection path can be expressed as follows,

where $m^{\prime}=(-1)^{p_{x}+p_{y}}m$ and $\varLambda(\mathbf{p},m,n)=(-1)^{(p_{y}+p_{z})m+p_{z}n}$ are the additional effect of the mirroring on the spherical harmonics [49]. Assuming $||\mathbf{R}^{\text{sI}\to\text{r}}_{\mathbf{p},\mathbf{q}}||$ is large, we employ the large argument of the spherical Hankel function, i.e., $h^{(2)}_{l}(k||\mathbf{R}^{\text{sI}\to\text{r}}_{\mathbf{p},\mathbf{q}}||)\approx i^{l+1}e^{-ik||\mathbf{R}^{\text{sI}\to\text{r}}_{\mathbf{p},\mathbf{q}}||}/k||\mathbf{R}^{\text{sI}\to\text{r}}_{\mathbf{p},\mathbf{q}}||$. The modified mode coupling coefficients in far field can be expressed as

Note that the product of two spherical harmonics can be expressed using the formula [62],

By applying Eq. (26) to Eq. (25), we can express the mode coupling coefficients as follows,

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

where DEISM-LC is used to refer to the Low-Complexity (LC) solution.

The terms $\varLambda(\mathbf{p},m,n)Y_{n,m^{\prime}}(\Omega^{\text{sI}\to\text{r}}_{\mathbf{p},\mathbf{q}})$ are a result of mirroring of spherical harmonics. Since the source images are used in Eq. (28), the additional factors $\varLambda(\mathbf{p},m,n)$ and $m^{\prime}$ need to be associated with the source spherical harmonics. To further simplify the expressions in Eq. (28), we propose to replace $\varLambda(\mathbf{p},m,n)Y_{n,m^{\prime}}(\Omega^{\text{sI}\to\text{r}}_{\mathbf{p},\mathbf{q}})$ by 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 $\mathbf{R}^{\text{sI}\to\text{r}}_{\mathbf{p},\mathbf{q}}$.

First, we note that the vector pointing from the receiver to the source image $\mathbf{R}^{\text{r}\to\text{sI}}_{\mathbf{p},\mathbf{q}}$ is the inverse of the vector $\mathbf{R}^{\text{sI}\to\text{r}}_{\mathbf{p},\mathbf{q}}$ given in Eq. (4). In addition, we define the vector pointing from the source to the receiver images $\mathbf{R}^{\text{s}\to\text{rI}}_{\mathbf{p},\mathbf{q}}$, which can be expressed similarly to Eq. (4) as

In spherical coordinates, the vectors $\mathbf{R}^{\text{r}\to\text{sI}}_{\mathbf{p},\mathbf{q}}$ and $\mathbf{R}^{\text{s}\to\text{rI}}_{\mathbf{p},\mathbf{q}}$ have directions $\Omega^{\text{r}\to\text{sI}}_{\mathbf{p},\mathbf{q}}$ and $\Omega^{\text{s}\to\text{rI}}_{\mathbf{p},\mathbf{q}}$, respectively. These two vectors are reversed reflection paths only when there are 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 $\mathbf{R}^{\text{s}\to\text{rI}}_{\tilde{\mathbf{p}},\tilde{\mathbf{q}}}$ (with directions $\Omega^{\text{s}\to\text{rI}}_{\tilde{\mathbf{p}},\tilde{\mathbf{q}}}$) that represents the reversed path of $\mathbf{R}^{\text{r}\to\text{sI}}_{\mathbf{p},\mathbf{q}}$ with its subscript indices $\tilde{\mathbf{p}},\tilde{\mathbf{q}}=(\tilde{p}_{x},\tilde{p}_{y},\tilde{p}_{z}),(\tilde{q}_{x},\tilde{q}_{y},\tilde{q}_{z})$ taking the following expressions

where $a\in\{x,y,z\}$, and $p_{a}^{\prime},q_{a}^{\prime}$ represent the image indices that are symmetric to $p_{a},q_{a}$ w.r.t the image indices $p_{a},q_{a}=0,0$. The relation between $p_{a}^{\prime},q_{a}^{\prime}$ and $p_{a},q_{a}$ can be summarized as

where $\lfloor\cdot\rfloor$ denotes the floored division.

Secondly, 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 $\Omega^{\text{sI}\to\text{r}}_{\mathbf{p},\mathbf{q}}$ in Eq. (28) and $\Omega^{\text{s}\to\text{rI}}_{\tilde{\mathbf{p}},\tilde{\mathbf{q}}}$ introduced above

Since $\mathbf{R}^{\text{r}\to\text{sI}}_{\mathbf{p},\mathbf{q}}=-\mathbf{R}^{\text{sI}\to\text{r}}_{\mathbf{p},\mathbf{q}}$, i.e., $\Omega^{\text{sI}\to\text{r}}_{\mathbf{p},\mathbf{q}}$ and $\Omega^{\text{r}\to\text{sI}}_{\mathbf{p},\mathbf{q}}$ have opposite directions, the following relation is obtained

Finally, by substituting Eqs. (33) and (34) in the approximated solution in Eq. (28), we obtain the simplified expression

The RTF of the DEISM-LC is then given by

The DEISM-LC is more computationally efficient since it does not require the summation over $l$ in the mode coupling coefficient in Eq. (12). The DEISM in Eq. (21) has a time complexity of $\mathcal{O}(IN^{2}V^{2}(N+V))$ per wavenumber $k$, where $I$ is the number of images, $N$ and $V$ are the maximum SH order used to describe directivities of the source and receiver. However, the time complexity of DEISM-LC is reduced to $\mathcal{O}(IN^{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.1, the FSRR utilizes directivity coefficients $\tilde{C}^{(\text{s})}_{n,m}(k),\tilde{C}^{(\text{r})}_{v,u}(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_{\mathbf{p},\mathbf{q}},M_{\mathbf{p},\mathbf{q}}$ used in Eq. (15) and $i^{n},i^{v}$ used in Eq. (36). This difference is also illustrated by an example described in Sec. IV.2.

A summary of the algorithm proposed in Sec. III.3 is listed as pseudocode in 1, 2 and 3. A Python implementation of the algorithm, along with an example, is available online at https://github.com/audiolabs/DEISM. The pseudocode of the far-field approximation method is not provided here for brevity.

The calculation of the Wigner 3j 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^{(\text{s})}_{n,m}(k)$ and $C^{(\text{r})}_{n,m}(k)$ are precalculated and stored in a lookup table. The procedures for generating the lookup tables are listed in Alg. 1.

In the provided implementation, we separate the whole procedure of the ISM into two parts, viz. the image generation given in Alg. 2 and the single-path transfer function calculation using DEISM given in Alg. 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.