# Directional reverberation time and the image source method for rectangular parallelepipedal rooms

Stefan Bilbao, Benoit Alary

► To cite this version:

Stefan Bilbao, Benoit Alary. Directional reverberation time and the image source method for rectangular parallelepipedal rooms. Journal of the Acoustical Society of America, 2024, 155 (2), pp.1343-1352. 10.1121/10.0024975. hal-04457732

HAL Id: hal-04457732

https://hal.science/hal-04457732v1

Submitted on 14 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# Directional reverberation time and the image source method for rectangular parallelepipedal rooms

Stefan Bilbao¹,ª) and Benoit Alary²)
simple means of prediction in the very simple case of the box-shaped room using IS methods.

In this article, we propose a closed-form expression to obtain directional reverberation times based on the IS method in the case of a parallelepipedal room. For the sake of generality, here, we work with angle-dependent wall reflection coefficients,³⁰ but the resulting expression may be simplified by using angle-independent reflection coefficients, and examples are provided in such cases. (Impedance-based absorption is also referred to as “directional/direction-dependent” absorption in the literature³⁰ and, here, we use the term angle-dependent absorption to make a clear distinction with directional reverberation, which is the main topic of this study.) In either case, energy decay is purely exponential in the limit of large times along a given direction, a fact that has been indicated previously.²⁰ This general formula for the directivity of late reverberation in a parallelepiped is potentially useful for applications applying IS that could benefit from approximating characteristics of late reverberation using statistically motivated methods, such as artificial reverberation, instead of the full simulation of an impulse response.¹⁹,²²⁻²⁴ In practical applications, the formula may be used to efficiently obtain directional decay characteristics by sampling around the sphere and inform the design of a reverberator capable of synthesizing rooms with directional reverberation, such as a long corridor. In addition, the general expression presented here may be employed to yield an easily computed approximation to the full omnidirectional EDC. This omnidirectional approximation is distinct from those presented previously¹⁹ and holds in the case of angle-dependent reflection coefficients.

In Sec. II, the general problem of the solution of the wave equation over a parallelepipedal domain is presented, including frequency-independent impedance boundary conditions over all surfaces. In Sec. III, the pressure-based IS method is reviewed, including the definition of reflection coefficients in the cases of angle-dependent and angle-independent reflections. In Sec. IV, a closed-form expression for the EDC and reverberation time as a function of direction is derived, which is valid in the late response. After integration over the sphere, a direct expression for omnidirectional reverberation time also results. In Sec. V, various numerical results are presented. First, plots of late reverberation times are displayed as a function of angle for rooms of various shapes and different choices of wall impedances. Comparisons between directional reverberation times in the cases of angle-dependent and angle-independent reflection are then shown, illustrating trough-like features in the case of angle-dependent reflection (corresponding to angles at which no energy is reflected in the IS method). Next, the formula for directional decay time is validated against the direct output of the IS method for a variety of room geometries and impedance conditions. Finally, comparisons between EDCs drawn from the IS method against curves drawn from a closed-form expression are compared in the omnidirectional case. A public code repository demonstrating how to implement the directional decay time formula is accessible online.⁴³

## II. THE WAVE EQUATION OVER A RECTANGULAR PARALLELEPIPED

Room acoustics simulation is generally concerned with numerical solutions to the three-dimensional wave equation over an enclosure:

$$ \frac{1}{c^2} \partial_t^2 p - \nabla^2 p = f \delta^{(3)}(\mathbf{r} - \mathbf{r}^{(S)}). \quad (1) $$

Here, $p(\mathbf{r}, t)$ is the acoustic pressure (in Pa) as a function of time $t \ge 0$, and spatial coordinate $\mathbf{r} = [r_1, r_2, r_3] \in \mathcal{D} \subset \mathbb{R}^3$. $c$ is the wave speed in m s⁻¹. $\partial_t$ represents partial differentiation with respect to time, $t$, and $\nabla^2$ is the three-dimensional spatial Laplacian operator. Finally, $f(t)$ is a driving term, assumed to be acting on the field pointwise at $\mathbf{r} = \mathbf{r}^{(S)} \in \mathcal{D}$, where $\delta^{(3)}$ represents a three-dimensional Dirac delta function. The function $f(t)$, defined for $t \ge 0$, is assumed to be impulsive in nature and square integrable such that

$$ \int_{0}^{\infty} f^{2}(t) dt = 1. \quad (2) $$

This value of the squared norm has been chosen to be one here to simplify subsequent calculations of reverberation time.³¹ The source term in Eq. (1) is omnidirectional but can easily be generalised to include source directivity.³²

In this simple study, the domain $\mathcal{D}$ is defined as a rectangular parallelepiped, centered at $\mathbf{r} = \mathbf{0}$, with side lengths $L_1, L_2$, and $L_3$,

$$ \mathcal{D} = [-L_1/2, L_1/2] \times [-L_2/2, L_2/2] \times [-L_3/2, L_3/2]. \quad (3) $$

See the left side of Fig. 1. Over the six bounding faces of the parallelepiped, which are defined as

$$ \partial \mathcal{D}_{\nu}^{\pm} = \{ \mathbf{r} \in \mathcal{D} | r_{\nu} = \pm L_{\nu} / 2 \}, \quad \nu = 1, \dots, 3, \quad (4) $$

boundary conditions of purely resistive type are defined such that

$$ p = \pm \rho c z_{\nu}^{\pm} v_{\nu} \quad \text{for} \quad \mathbf{r} \in \partial \mathcal{D}_{\nu}^{\pm} \quad \text{and} \quad \nu = 1, \dots, 3, \quad (5) $$

where $v_\nu$ is the component of the particle velocity in direction $\nu$. Here, the six nondimensional constants $z_\nu^\pm$, $\nu = 1, \dots, 3$, are assumed real and nonnegative. These are the wall impedances, normalized against the characteristic impedance of air, $\rho c$. For the purposes of this contribution, the impedances are assumed to be purely real. In general, wall impedances exhibit frequency dependence, in which case Eq. (5) is replaced by a frequency domain relationship between pressure and velocity at a given surface, and impedances are, in general, complex valued. Such conditions are locally reactive and do not model extended reaction conditions due, e.g., to the presence of porous media.³³

A response is assumed to be drawn from the acoustic field at coordinates $\mathbf{r} = \mathbf{r}^{(R)} \in \mathcal{D}$ as

$$ p^{(R)}(t) = p(\mathbf{r}^{(R)}, t). \quad (6) $$
<figure><img src="image_4.png" /><figcaption>FIG. 1. (Color online) (Left) Rectangular parallelepiped with source and receiver coordinates, **r**<sup>(S)</sup> and **r**<sup>(R)</sup> as indicated by red and yellow points, respectively; (center) an array of image volumes with IS locations indicated in blue; and (right) cross section in the (r<sub>1</sub>, r<sub>2</sub>) plane with the qth IS coordinates, **r**<sup>(q)</sup>, and displacement vector, **r**<sub>Δ</sub><sup>(q)</sup>, as indicated are shown.</figcaption></figure>

For general choices of wall impedances, $z_{\nu}^{\pm}$, a closed-form solution for $p^{(R)}(t)$ in terms of the driving function, $f(t)$, is not available, in general.

### III. THE IMAGE SOURCE METHOD

In the special cases where $z_{\nu}^{\pm} = +\infty$, corresponding to a perfectly reflecting wall condition on all sides of the parallelepiped, one approach to obtaining an exact solution at a single output location, **r**<sup>(R)</sup>, is through the use of the IS method.¹ The IS method is still employed, however, even in cases of finite values of the wall impedance. In this case, the solution generated is no longer exact, and the IS method is approximate. Careful comparisons against wave-based solutions have been presented.³⁴

The method consists of the redefinition of the problem domain to $\mathcal{D} = \mathbb{R}^3$ and the insertion of so-called image sources at locations **r**<sup>(q)</sup> reflected multiple times through the parallelepipedal boundaries and indexed by the integer-valued three-vector **q** = [q<sub>1</sub>, q<sub>2</sub>, q<sub>3</sub>] ∈ $\mathbb{Z}^3$,

$$ \mathbf{r}^{(\mathbf{q})} = [r_{1}^{(q_1)}, r_{2}^{(q_2)}, r_{3}^{(q_3)}], \qquad (7) $$

where

$$ r_{\nu}^{(q_{\nu})} = q_{\nu}L_{\nu} + (-1)^{q_{\nu}}r_{\nu}^{(S)} \quad \text{for } \nu = 1, \dots, 3. \qquad (8) $$

(See Fig. 1 center.) Assuming point sources, the resulting solution at coordinates **r**<sup>(R)</sup> is of the form

$$ p_{\text{IS}}^{(\text{R})}(t) = \sum_{\mathbf{q} \in \mathbb{Z}^3} \xi^{(\mathbf{q})}(t), \qquad (9) $$

where

$$ \xi^{(\mathbf{q})}(t) = a^{(\mathbf{q})} \frac{f(t - r_{\Delta}^{(\mathbf{q})}/c)}{4\pi r_{\Delta}^{(\mathbf{q})}}. \qquad (10) $$

Here, $r_{\Delta}^{(\mathbf{q})}$ is the distance from IS, **q**, to the receiver location, **r**<sup>(R)</sup>, such that

$$ r_{\Delta}^{(\mathbf{q})} = |\mathbf{r}^{(\mathbf{q})} - \mathbf{r}^{(\mathbf{R})}|. \qquad (11) $$

The signal $\xi^{(\mathbf{q})}(t)$ represents the contribution to the computed pressure waveform, $p_{\text{IS}}^{(\text{R})}(t)$, due to image source, **q**. In practical computations, the triply infinite sum in Eq. (9) is truncated to a finite set of terms (normally falling within a spherical region centered at **r** = **r**<sup>(R)</sup>). Note that $p_{\text{IS}}^{(\text{R})}(t)$ is, in general, not equal to $p^{(\text{R})}(t)$ from Eq. (6), which is the exact solution to the wave equation over a parallelepiped, except under perfectly reflecting conditions. In this article, we employ a pressure-based IS method rather than an energy-based IS method.⁶

The dimensionless constants, $a^{(\mathbf{q})}$, must be determined through an analysis of the number of “reflections” undergone by the wave emitted by the image source in its transit to the receiver location. To this end, it is useful to decompose the vector, $\mathbf{r}_{\Delta}^{(\mathbf{q})}$, in terms of its length, $r_{\Delta}^{(\mathbf{q})}$, and the unit-length direction vector, $\mathbf{u}^{(\mathbf{q})} = [u_1^{(\mathbf{q})}, u_2^{(\mathbf{q})}, u_3^{(\mathbf{q})}]$, as

$$ \mathbf{r}_{\Delta}^{(\mathbf{q})} = r_{\Delta}^{(\mathbf{q})} \mathbf{u}^{(\mathbf{q})}, \qquad (12) $$

where

$$ \mathbf{u}^{(\mathbf{q})} = [\cos(\theta^{(\mathbf{q})}) \sin(\phi^{(\mathbf{q})}), \sin(\theta^{(\mathbf{q})}) \sin(\phi^{(\mathbf{q})}), \cos(\phi^{(\mathbf{q})})]. \qquad (13) $$

Here, $\theta^{(\mathbf{q})} \in [0, 2\pi)$ and $\phi^{(\mathbf{q})} \in [0, \pi]$ are the azimuthal and polar angles of the qth image source relative to the receiver location, respectively. Next, the number of reflections undergone by the qth image source through each of the six faces of the original parallelepipedal domain must be determined. The six numbers will be denoted as $N_{\nu}^{(\mathbf{q}),\pm}$, $\nu = 1, \dots, 3$ and can be written in terms of the components of **q** as

$$ N_{\nu}^{(\mathbf{q}),+} = \left| \frac{q_{\nu}}{2} \right|, \quad N_{\nu}^{(\mathbf{q}),-} = \left| \frac{q_{\nu}}{2} \right| \quad \text{for } \nu = 1, \dots, 3. \qquad (14) $$

Here, $[\cdot]$ and $\lfloor \cdot \rfloor$ indicate ceiling and floor operations, respectively, taking sign into consideration. Note that $N_{\nu}^{(\mathbf{q}),+} + N_{\nu}^{(\mathbf{q}),-} = |q_{\nu}|$.

The constants, $a^{(\mathbf{q})}$, in which the reflectances undergone by image source **q**, are accumulated, may now be assembled as
$$a^{(\mathbf{q})} = \prod_{\nu=1}^{3} \prod_{\alpha=+,-} \left( \beta_{\nu}^{(\mathbf{q}),\alpha} \right)^{N_{\nu}^{(\mathbf{q}),\alpha}}, \quad (15)$$

where the reflectances, $\beta_{\nu}^{(\mathbf{q}),\alpha}(\mathbf{u}^{(\mathbf{q})})$, are defined by

$$\beta_{\nu}^{(\mathbf{q}),\alpha}(\mathbf{u}^{(\mathbf{q})}) = \frac{z_{\nu}^{\alpha}|u_{\nu}^{(\mathbf{q})}| - 1}{z_{\nu}^{\alpha}|u_{\nu}^{(\mathbf{q})}| + 1}. \quad (16)$$

This definition takes into account the specular nature of the reflection assumed here and, thus, the reflection coefficients, $\beta_{\nu}^{(\mathbf{q}),\alpha}$, are angle dependent.⁶

If, instead of angle-dependent reflection, angle-independent reflection is assumed, then the coefficients, $\beta_{\nu}^{(\mathbf{q}),\alpha}$, become constants that are independent of angle or $\mathbf{u}^{(\mathbf{q})}$. In the present context, these are often written¹,¹⁹ in terms of a constant absorption, $a_{\nu}^{\alpha}$, at a given wall surface,

$$\beta_{\nu}^{(\mathbf{q}),\alpha} = \sqrt{1 - a_{\nu}^{\alpha}}. \quad (17)$$

The absorption may be derived from measurement but also may be related to an impedance at normal incidence, $z_{\nu}^{\alpha}$. In this article, we have used the approach outlined by London.³⁵ Many other approaches to determining angle-independent reflection coefficients are available³⁶⁻³⁹—the choice above is made for reasons of simplicity and does not affect the computational procedure for determining direction-dependent reverberation time given in Sec. IV.

It follows that the reflectances, $\beta_{\nu}^{(\mathbf{q}),\alpha}$, are bounded by one in magnitude such that

$$|\beta_{\nu}^{(\mathbf{q}),\alpha}| \le 1. \quad (18)$$

This holds in either the case of angle-dependent coefficients, from Eq. (16), or angle-independent coefficients, from Eq. (17). In the angle-dependent case, the reflection coefficients can be negative, leading to reflection with inversion for nearly tangential reflections for which $|u_{\nu}^{(\mathbf{q})}| < 1/z_{\nu}^{\alpha}$. In the angle-independent case, the reflectances are constrained to be positive in Eq. (17), although some authors employ negative values.¹⁹

For the purposes of the determination of reverberation times, it is also useful to define the so-called “energy” associated with the individual contribution to the IS solution:

$$e^{(\mathbf{q})}(t) = \left(\xi^{(\mathbf{q})}(t)\right)^2 = (a^{(\mathbf{q})})^2 \frac{f^2(t - r_{\Delta}^{(\mathbf{q})}/c)}{16\pi^2 (r_{\Delta}^{(\mathbf{q})})^2}. \quad (19)$$

# IV. REVERBERATION TIME

## A. Directional reverberation time

Now, consider a single image source at coordinates $\mathbf{r}_{\Delta} = r_{\Delta}\mathbf{u}$ relative to the receiver location, where, as above, $r_{\Delta}$ is the distance between the image source and receiver, and $\mathbf{u} = [u_1, u_2, u_3]$ is a unit-length direction vector, which is defined as

$$\mathbf{u} = [\cos(\theta) \sin(\phi), \sin(\theta) \sin(\phi), \cos(\phi)]. \quad (20)$$

Next, consider the definition of the contribution of the signal $\xi^{(\mathbf{q})}(t)$ to the total IS solution, $p_{\text{IS}}^{(\text{R})}(t)$, from Eq. (9), and the associated energy from Eq. (19). These may be rewritten for a general image source location, $\mathbf{r}_{\Delta}$, as

$$\xi(\mathbf{r}_{\Delta}, t) = a(\mathbf{r}_{\Delta}) \frac{f(t - r_{\Delta}/c)}{4\pi r_{\Delta}}, \quad (21a)$$

$$e(\mathbf{r}_{\Delta}, t) = a^2(\mathbf{r}_{\Delta}) \frac{f^2(t - r_{\Delta}/c)}{16\pi^2 r_{\Delta}^2}. \quad (21b)$$

Returning now to Eq. (15) of the amplitudes of the individual contributions to the IS solution, note that in the limit of large image source locations, $\mathbf{r}_{\Delta}$, one may approximate the numbers of reflections as

$$N_{\nu}^{\pm}(\mathbf{r}_{\Delta}) \approx \frac{r_{\Delta}|u_{\nu}|}{2L_{\nu}}. \quad (22)$$

This limiting “large radius” assumption is also used in other work on approximations to the omnidirectional EDC curve.¹⁹ In contrast here, we do not use a simplified first-order approximation to the components $u_{\nu}$ in terms of elevation and azimuthal angles.

One arrives, finally, at expressions for the squared amplitudes,

$$a^2(\mathbf{r}_{\Delta}) = \prod_{\nu=1}^{3} \prod_{\alpha=+,-} |\beta_{\nu}^{\alpha}(\mathbf{u})|^{r_{\Delta}|u_{\nu}|/L_{\nu}}, \quad (23)$$

where in the case of angle-dependent reflection [Eq. (16)],

$$\beta_{\nu}^{\alpha}(\mathbf{u}) = \frac{z_{\nu}^{\alpha}|u_{\nu}| - 1}{z_{\nu}^{\alpha}|u_{\nu}| + 1}, \quad (24)$$

and for angle-independent reflection, as in Eq. (17),

$$\beta_{\nu}^{\alpha}(\mathbf{u}) = \beta_{\nu}^{\alpha} = \sqrt{1 - a_{\nu}^{\alpha}}. \quad (25)$$

At this point, it is useful to introduce the time parameter, $\tau_{\Delta} = r_{\Delta}/c$, corresponding to the time delay of the contribution to $p_{\text{IS}}^{(\text{R})}(t)$, which is received from an image source at coordinates $\mathbf{r}_{\Delta}$ and, thus,

$$e(\tau_{\Delta}, \mathbf{u}, t) = a^2(\tau_{\Delta}, \mathbf{u}) \frac{f^2(t - \tau_{\Delta})}{16\pi^2 c^2 \tau_{\Delta}^2}, \quad (26)$$

where

$$a^2(\tau_{\Delta}, \mathbf{u}) = \prod_{\nu=1}^{3} \prod_{\alpha=+,-} |\beta_{\nu}^{\alpha}|^{c\tau_{\Delta}|u_{\nu}|/L_{\nu}}. \quad (27)$$

In particular, note that for a fixed direction $\mathbf{u}$, the expression for $a^2(\tau_{\Delta}, \mathbf{u})$ is a pure exponential function,

$$a^2(\tau_{\Delta}, \mathbf{u}) = e^{-K(\mathbf{u})\tau_{\Delta}}, \quad (28)$$

where

$$K(\mathbf{u}) = -c \sum_{\nu=1}^{3} \sum_{\alpha=+,-} \frac{\ln(|\beta_{\nu}^{\alpha}|) |u_{\nu}|}{L_{\nu}} \ge 0. \quad (29)$$
The time-integrated energy received may be defined, recall-
ing Eq. (2), as

$$
E(\tau_{\Delta}, \mathbf{u}) = \int_{0}^{\infty} e(\tau_{\Delta}, \mathbf{u}, t) dt = \frac{a^2(\tau_{\Delta}, \mathbf{u})}{16\pi^2 c^2 \tau_{\Delta}^2} = \frac{e^{-K(\mathbf{u})\tau_{\Delta}}}{16\pi^2 c^2 \tau_{\Delta}^2}. \quad (30)
$$

The pure exponential decay of the energy in a given direc-
tion and in the limit of large times in the case of the IS
method, defined over a parallelepiped, has been indicated
previously.²⁰

Now, following Kuttruff,¹⁸ consider a thin cone with its
apex at r<sub>Δ</sub> = 0, oriented along direction **u**, and subtending a
solid angle Ω. For small enough Ω, all image sources con-
tained within this cone can be considered to be at direction
**u** relative to the receiver location. The total number of
image sources over a time interval, [τ<sub>Δ</sub>, τ<sub>Δ</sub> + dτ<sub>Δ</sub>], for a
cone of solid angle Ω, and taking into account the density,
1/V, of image sources per unit volume, where V = L<sub>1</sub>L<sub>2</sub>L<sub>3</sub>
is the room volume, is

$$
d\Xi = \frac{\Omega c^3 \tau_{\Delta}^2 d\tau_{\Delta}}{V}. \qquad (31)
$$

Thus, the total energy received at **r**<sup>(R)</sup> from image sources
from direction **u** over this range of delays is *dE*, which is
defined as

$$
d\mathcal{E} = E(\tau_{\Delta}, \mathbf{u})d\Xi = \frac{c\Omega}{16\pi^2 V} e^{-K(\mathbf{u})\tau_{\Delta}} d\tau_{\Delta}. \quad (32)
$$

A late-time approximation to the EDC, as a function of inci-
dence direction **u**, may then be obtained by integrating *dE*
from time *t* to ∞, yielding

$$
\mathrm{EDC}_{\mathrm{late}}(\mathbf{u}, t) = \frac{c\Omega}{16\pi^2 V K(\mathbf{u})} e^{-K(\mathbf{u})t}. \quad (33)
$$

The pure exponential dependence on time *t* is important
here. On a dB scale, this curve is in the form of a line with
slope −10 log₁₀(*e*)*K*(*u*).

An expression for the directional RT₆₀(**u**), which are
defined here from the decay rates in the late response, is
obtained from *K*(**u**) through

$$
\mathrm{RT}_{60}(\mathbf{u}) = \frac{6 \ln(10)}{K(\mathbf{u})}. \qquad (34)
$$

Unsurprisingly, this expression, which pertains to late reverber-
ation, is independent of the source and receiver locations. This
expression has been derived using angle-dependent reflection,
leading to a dependence of the factors β<sub>ν</sub><sup>α</sup> on **u** as in Eq. (23). If
angle-independent reflection is employed, then these factors
become constants, but the directional behavior of the RT<sub>60</sub> will
persist as per the definition of K(**u**) from Eq. (16).

B. Reduction to an omnidirectional EDC

Equation (32), which gives the energy received at the
receiver location at time τΔ from direction **u**, could be

integrated over the sphere to yield a global omnidirectional
measure of energy decay. In this case, non-exponential
decay is expected, and examples are given by Kanev²⁰ under
particular configurations (such as, e.g., a room with only one
absorptive wall, leading to a 1/τΔ dependence).

Beginning from Eq. (32), for infinitesimal Ω, integrat-
ing over the unit sphere and backward time integrating ana-
lytically leads directly to the following expression for the
omnidirectional EDC that is valid in the late response:

$$
\mathrm{EDC}_{\mathrm{late}}(t) = \frac{c}{16\pi^2 V} \iint_{S^2} \frac{e^{-K(\mathbf{u})t}}{K(\mathbf{u})} d\Omega. \quad (35)
$$

This is simple to evaluate numerically, requiring only an integration over the unit sphere S² (which, by symmetry, can be reduced to an integration over a single octant). It relies on fewer simplifying assumptions than the closed-form expression of Lehmann and Johansson.¹⁹ In particular, (a) this form yields the EDC directly without requiring a further numerical backward time integration; (b) there is not a further first-order approximation to the components of the direction vector **u** in terms of angles; and (c) it holds, generally, for the cases of angle-dependent and angle-independent reflections—and indeed, results differ considerably according to this choice— see Sec. V D. Although it is intended to model the energy decay in the late reverberant tail, like the expression of Lehmann and Johansson, it does an excellent job of approximating the early decay of the IS method. See Sec. V D.

The omnidirectional EDC defined by Eq. (35) is non-
exponential, but notice from the directional representation of
the RT₆₀ from Eq. (34) that the integrand in Eq. (35) is larg-
est for directions **u** with a long RT₆₀(**u**) [or small value of
*K*(**u**)], implying that such directions will dominate in the cal-
culation of the omnidirectional EDC—a known result.¹⁸,²⁰

V. RESULTS

In this section, we show various examples of the
closed-form directional reverberation time expression from
Eq. (34) for different choices of room geometry and wall
conditions and validate the expression against the output of
the IS method, which is modified to compute the directional
EDC over a cone of small angle. In addition, we compare
omnidirectional EDC curves obtained directly from the IS
method with Eq. (35).

Angle-dependent [Eq. (16)] and angle-independent [Eq.
(17)] reflection coefficients will be employed here. The
wave speed, c, is set uniformly to c = 344 m s⁻¹.

A. RT₆₀ maps

It is direct, given Eq. (34) for RT₆₀ as a function of
direction **u** through *K*(**u**), as defined in Eq. (29), to examine
the directional dependence for a given choice of room
dimensions, *L*<sub>ν</sub>, *ν* = 1, ..., 3, and associated impedances,
*z*<sub>ν</sub><sup>±</sup>, *ν* = 1, ..., 3. See Fig. 2, which shows directional RT₆₀
over the unit sphere as well as over the planar cross sections,
*φ* = π/2, *θ* = 0, and *θ* = π/2, for a variety of choices of
<figure><img src="image_5.png" /><figcaption>FIG. 2. (Color online) Variation in RT₆₀ under three different choices of room dimensions and wall impedances. In each case, the expression (34) is represented as a surface plot over a sphere, where maximal and minimal values are indicated in the accompanying color bar. In addition, planar cross sections of the RT₆₀ for $\phi = \pi/2$, $\theta = 0$, and $\theta = \pi/2$ are shown, where the polar angle, $\phi$, and azimuthal angle, $\theta$, are as indicated in the top left plot.</figcaption></figure>

room dimensions and impedances and uses angle-dependent reflection coefficients. There are large variations, depending on the impedances and wall dimensions—normally, the axial directions are associated with locally maximal RT₆₀ times.

An interesting feature, visible in the plots at the bottom in Fig. 2, is the presence of “troughs” in the RT₆₀ at certain angles. Such angles correspond to those for which the reflection vanishes—this is possible when at least one of the normalized impedances, $z_{\nu}^{\pm}$, is greater than unity. This feature results from the use of angle-dependent reflection coefficients and disappears if angle-independent reflection coefficients are employed. As a simple example, consider the case depicted in the middle in Fig. 2 with uniform dimensions of $L_{\nu} = 8$ m and varying wall impedances, $z_{\nu}^{\alpha}$. In Fig. 3, the RT₆₀ over the circle $\phi = \pi/2$ is shown with troughs indicated at the angles $\theta = \sin^{-1}(1/z_{2}^{+})$ and $\theta = \sin^{-1}(1/z_{2}^{-})$.

## B. Angle-dependent vs angle-independent reflection coefficients

The directional dependence of the late reverberant tail depends strongly on the type of reflection model. In Sec. VA,

plots of RT₆₀ were shown in the case of angle-dependent reflection. As can be observed from Eq. (29), however, the variation in RT₆₀ persists in the case where the reflection coefficients, $\beta_{\nu}^{\alpha}$, are constants independent of angle. See Fig. 4, which illustrates angle-dependent and angle-independent

<figure><img src="image_6.png" /><figcaption>FIG. 3. (Color online) RT₆₀ for $\phi = \pi/2$ and as a function of angle $\theta \in [0, \pi/4]$ for the case depicted at the middle in Fig. 2. The locations of troughs are indicated by red dotted lines at angles $\theta = \sin^{-1}(1/z_{2}^{+})$ and $\theta = \sin^{-1}(1/z_{2}^{-})$.</figcaption></figure>
$$ 
\begin{aligned}
L_1 &= 6, & L_2 &= 7, & L_3 &= 11 \\
z_1^+ &= 20, & z_2^+ &= 0.4, & z_3^+ &= 20 \\
z_1^- &= 20, & z_2^- &= 72, & z_3^- &= 0.4
\end{aligned}
$$
$$ \phi = \pi/2 $$

<figure><img src="image_7.png" /><figcaption>FIG. 4. (Color online) Cross section of RT$_{60}$ for room geometry and wall impedances. Curves in the cases of angle-dependent reflection (in blue) and angle-independent reflection, using the form in Eq. (17) for all angles (in dashed red), are shown.</figcaption></figure>

predicted RT$_{60}$ curves over a cross section of the sphere. Notably absent are the troughs in the case of angle-independent reflection.

Beyond this distinction, however, the predicted curves are quite different. Here, the angle-independent reflection coefficients have been chosen to correspond to London's formula³⁵ at all angles for simplicity, but this difference remains regardless of how the reflection coefficients are chosen (perhaps derived from measured absorption and the assumption of random incidence³⁸). This signals that the calculation of global (omnidirectional) measures of reverberation derived from the IS method,¹⁹ such as that given in Sec. IV B, will be influenced strongly by this choice.

### C. Directional EDC validation

The validation against EDCs obtained through the IS method is also direct. In this case, the simplest computational procedure is to define a cone subtending a small but finite angle, $\psi$, relative to the direction vector, $\mathbf{u}$, and describing a region $D(\mathbf{u}, \psi)$, which is defined by

$$ D(\mathbf{u}, \psi) = \{\mathbf{r} \in \mathbb{R}^3 | (\mathbf{r}/|\mathbf{r}|) \cdot \mathbf{u} \ge \cos(\psi)\}. \quad (36) $$

The summation (9) used in the computation of the response may, thus, be modified to obtain a directional response, $p_{\text{IS}}^{(\text{R})}(\mathbf{u}, t)$,

$$ p_{\text{IS}}^{(\text{R})}(\mathbf{u}, t) = \sum_{\substack{\mathbf{q} \in \mathbb{Z}^3, \mathbf{r}_{\text{A}}^{(\text{q})} \in D(\mathbf{u}, \psi) \\ \psi \in \mathcal{C}^3}} \xi^{(\mathbf{q})}(t), \quad (37) $$

such that the sum above is limited to those image sources clustered around the direction, $\mathbf{u}$, relative to the receiver location. All other aspects of the IS computational procedure remain the same. Then, a directional EDC may be obtained directly as

$$ \text{EDC}(\mathbf{u}, t) = \int_{t}^{\infty} \left( p_{\text{IS}}^{(\text{R})}(\mathbf{u}, t') \right)^2 dt'. \quad (38) $$

It is direct to transfer the calculation above to discrete time at a sample rate of $F_s$ Hz, and by specifying a unit impulse to approximate the excitation, $f(t)$. Note that, in contrast with other authors,¹⁹ we do not normalize the EDC.

Energy decay curves for a particular choice of room geometry, a variety of incidence directions, and wall impedances are displayed in Fig. 5. These are calculated using a high sample rate of $F_s = 1$ MHz so as to achieve a good resolution and separation of impulses received from ISs when truncated to sampling instances. In addition, it is necessary to use a relatively small cone angle, $\psi$ (here, set to $\psi = \pi/200$), to ensure that variations in RT$_{60}$ over the cone are minimized. For such thin cones, the number of ISs contained is small in the early part of the response, explaining the jagged character of the EDC curves in Fig. 5. However, at later times, the number of image sources within the cone increases, and the EDC settles down to pure exponential decay with a slope on a log scale predicted by $K(\mathbf{u})$ as defined in Eq. (16), and as indicated in Fig. 5.

As a further cross-checking exercise, we illustrate differences in the approximation over the entire sphere, as shown in Fig. 6, using the room dimensions and impedances as given at right in Fig. 2 and under angle-dependent reflection conditions. This is a worst case with very large variations in wall reflectance. The differences represent less than 0.48% error over the majority of the sphere (92%) with larger differences observed near the axes. Larger differences occur in directions where very high IS orders are necessary to resolve very localized variations in RT$_{60}$ as seen in the cross section plots at right in Fig. 2. In these directions, it becomes more challenging to estimate the RT$_{60}$ by analysing the EDC slope from the generated directional RIR. The problem is especially severe near the troughs in the RT$_{60}$, as illustrated in Fig. 3, where very little energy is received, requiring extreme IS order.

The practical challenges of obtaining reliable RT$_{60}$ estimates using the IS method demonstrate the benefits of using the proposed formula to obtain measures at arbitrary angles. For sound reproduction, the formula may be sampled around the sphere to provide design parameters for a directional reverberator. However, as the perceptual threshold for directional decay characteristics remains unclear at this point, the appropriate level of detail is more likely to be informed by the available computing resources for a given practical application.⁴⁰ This is similar to how the absorbent filters of
<figure><img src="image_8.png" /><figcaption>FIG. 5. (Color online) Directional EDCs in dB at angles A–G for the IS method with geometric and impedance parameters as indicated. All curves are normalized to 0 dB. Superimposed on the EDCs in blue are lines of slope -10 $log_{10}(e)K(\mathbf{u})$, where $K(\mathbf{u})$ is as defined in Eq. (16).</figcaption></figure>

a reverberator are designed to ensure frequency-dependent decay without clear guidelines on minimum requirements.⁴¹

## D. Omnidirectional EDCs

As a final example, we revert to the reduced omnidirectional estimate [Eq. (35)] of the EDC curve given in Sec. IV B. Returning to the three examples shown in Fig. 2, the omnidirectional EDC curve calculated directly from Eq. (35) is depicted, alongside EDC curves drawn directly from the IS

method. See Fig. 7. The cases of angle-dependent and angle-independent reflection are displayed. In all cases, the match is very good (and to within 1–2 dB over the interval shown). As expected, the particular approximation³⁵ used to derive angle-independent reflectance from the impedances exhibited has a large bearing on the resulting EDC but will generally lead to a faster energy decay. In both cases, however, the approximation matches well and includes the changes in slope in the very early part of the decay curve.

<figure><img src="image_9.png" /><figcaption>FIG. 6. (Color online) Differences between the estimated directional RT₆₀ and formula (34) for the case depicted at the right in Fig. 2. The RT₆₀ was estimated from IS simulation up to order 400. The differences are between 0 and 120 ms. The mean difference is 4.88 ms and less than 25 ms over 99.24% of the sphere.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIG. 7. (Color online) Energy decay curves for the three choices of room geometries and impedances displayed in Fig. 2 under angle-independent reflection (top row) and angle-dependent reflection (bottom row). The EDC curve calculated directly from IS output is shown as a black line, and results computed from formula (35) are shown as red circles (angle-independent case) or blue circles (angle-dependent case).</figcaption></figure>

The surface integral used in Eq. (35) is numerically integrated at intervals of 0.1 s, using a resolution of $\pi/800$ in the elevation angle, $\phi$, and $\pi/400$ in azimuth, $\theta$. As mentioned previously, no numerical time integration is necessary, and no further adjustments have been applied (such as shifts of the computed values from the formula up or down for matching purposes).

## VI. CONCLUDING REMARKS

The main result here is a closed-form expression for reverberation time in a parallelepipedal room as a function of direction as given in Eqs. (34) and (16). It is valid in the case of angle-dependent reflection, where the reflectances, $\beta_{\nu}^{\alpha}$, are as defined in Eq. (23), and in the case of angle-independent reflection, in which case, the reflectances, $\beta_{\nu}^{\alpha}$, are held constant [and perhaps set according to absorbances as per Eq. (17)]. Furthermore, energy decay in late reverberation is purely exponential for any given incidence angle. This result has been cross-checked against EDCs obtained directly from the IS method. There are major differences in the resulting directional profile of reverberation time, including the presence of so-called troughs at which the reverberation time becomes zero at particular angles in the case of angle-dependent reflectances. In addition, a simple direct formula for the EDC in the omnidirectional case has been presented and holds for either angle-dependent or angle-independent reflection.

While useful in framing a test problem, the box-shaped region exhibits an extreme symmetry that leads to preferred directions in terms of decay time. It is clear that the magnitude of variation in the late reverberation time with direction will decrease when more realistic geometries are considered, but an open question is the extent of this decrease and how sensitive it is to small variations in geometry (i.e., can one

expect residual preferred directions in a room that is nearly box-shaped?). Some case studies$^{42}$ indicate that the sensitivity may be high.

A further aspect to consider is the inaccuracy of the IS method—one could envisage employing a wave-based method to test the directivity of reverberant characteristics further. Nonetheless, this work has valuable applications in creating a statistical approximation of the late reverberation behavior of the IS method, notably through the use of a directional artificial reverberator,$^{29}$ which may be specified using the decay parameters obtained using the proposed formula.

## ACKNOWLEDGMENTS

Thanks go to Lauri Savioja of the Department of Computer Science at Aalto University for many interesting discussions on this topic.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

<footnote>$^{1}$J. Allen and D. Berkley, “Image method for efficiently simulating small-room acoustics,” J. Acoust. Soc. Am. 65(4), 943–950 (1979).</footnote>

<footnote>$^{2}$J. Borish, “Extension of the image model to arbitrary polyhedra,” J. Acoust. Soc. Am. 75(6), 1827–1836 (1984).</footnote>

<footnote>$^{3}$B. Champagne, “Simulation of the response of multiple microphones to a moving point source,” Appl. Acoust. 42(4), 313–332 (1994).</footnote>
<footnote n="4">⁴D. Jarrett, E. Habets, M. Thomas, and P. Naylor, "Rigid sphere room impulse response simulation: Algorithm and applications," *J. Acoust. Soc. Am.* **132**(3), 1462–1472 (2012).</footnote>

<footnote n="5">⁵M. Vorländer, *Auralization: Fundamentals of Acoustics, Modelling, Simulation, Algorithms and Acoustic Virtual Reality* (Springer, Berlin, 2008).</footnote>

<footnote n="6">⁶L. Savioja and U. Svensson, "Overview of geometrical room acoustic modeling techniques," *J. Acoust. Soc. Am.* **138**(2), 708–730 (2015).</footnote>

<footnote n="7">⁷A. Krokstad, S. Strom, and S. Sørsdal, "Calculating the acoustical room response by the use of a ray tracing technique," *J. Sound Vib.* **8**(1), 118–125 (1968).</footnote>

<footnote n="8">⁸T. Funkhouser, N. Tsingos, I. Carlbom, G. Elko, M. Sondhi, J. E. West, G. Pingali, P. Min, and A. Ngan, "A beam tracing method for interactive architectural acoustics," *J. Acoust. Soc. Am.* **115**(2), 739–756 (2004).</footnote>

<footnote n="9">⁹C.-H. Jeong, J.-G. Ih, and J. Rindel, "An approximate treatment of reflection coefficient in the phased beam tracing method for the simulation of enclosed sound fields at medium frequencies," *Appl. Acoust.* **69**(7), 601–613 (2008).</footnote>

<footnote n="10">¹⁰S. Laine, S. Siltanen, T. Lokki, and L. Savioja, "Accelerated beam tracing algorithm," *Appl. Acoust.* **70**(1), 172–181 (2009).</footnote>

<footnote n="11">¹¹R. R. Torres, U. P. Svensson, and M. Kleiner, "Computation of edge diffusion for more accurate room acoustics auralization," *J. Acoust. Soc. Am.* **109**(2), 600–610 (2001).</footnote>

<footnote n="12">¹²S. Bilbao, B. Hamilton, J. Botts, and L. Savioja, "Finite volume time domain room acoustics simulation under general impedance boundary conditions," *IEEE/ACM Trans. Audio. Speech. Lang. Process.* **24**(1), 161–173 (2016).</footnote>

<footnote n="13">¹³F. Pind, A. Engsig-Karup, C. Jeong, J. Hesthaven, M. Mejling, and J. Stromann-Andersen, "Time domain room acoustic simulations using the spectral element method," *J. Acoust. Soc. Am.* **145**(6), 3299–3310 (2019).</footnote>

<footnote n="14">¹⁴H. Wang, M. Cosnefroy, and M. Hornikx, "An arbitrary high-order discontinuous Galerkin method with local time-stepping for linear acoustic wave propagation," *J. Acoust. Soc. Am.* **149**(1), 569–580 (2021).</footnote>

<footnote n="15">¹⁵A. Neidhardt, C. Schneiderwind, and F. Klein, "Perceptual matching of room acoustics for auditory augmented reality in small rooms—Literature review and theoretical framework," *Trends Hear.* **26**, 233121652210929 (2022).</footnote>

<footnote n="16">¹⁶M. Vorländer, D. Schröder, S. Pelzer, and F. Wefers, "Virtual reality for architectural acoustics," *J. Build. Perf. Sim.* **8**(1), 15–25 (2014).</footnote>

<footnote n="17">¹⁷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).</footnote>

<footnote n="18">¹⁸H. Kuttruff, *Room Acoustics*, 5th ed. (Taylor and Francis, London, 2009).</footnote>

<footnote n="19">¹⁹E. Lehmann and A. Johansson, "Prediction of energy decay in room impulse responses simulated with an image-source model," *J. Acoust. Soc. Am.* **124**(1), 269–277 (2008).</footnote>

<footnote n="20">²⁰N. G. Kanev, "Sound decay in a rectangular room with impedance walls," *Acoust. Phys.* **58**(5), 603–609 (2012).</footnote>

<footnote n="21">²¹V. Välimäki, J. Parker, L. Savioja, J. O. Smith, and J. S. Abel, "Fifty years of artificial reverberation," *IEEE Trans. Audio. Speech. Lang. Process.* **20**(5), 1421–1448 (2012).</footnote>

<footnote n="22">²²E. Accolti and F. Miyara, "Fast and controllable box-shaped room impulse response algorithm," in *13th Meeting on Work in Information Processing and Control*, Rosario, Argentina (September 16–18, 2009).</footnote>

<footnote n="23">²³E. A. Lehmann and A. M. Johansson, "Diffuse reverberation model for efficient image-source simulation of room impulse responses," *IEEE Trans. Audio. Speech. Lang. Process.* **18**(6), 1429–1439 (2010).</footnote>

<footnote n="24">²⁴A. Nogueiras Rodríguez and J. Colom Olivares, "A statistical approach to reverberation in non-diffusive rectangular rooms based on the image source model," in *IEEE Int. Conf. Acoust., Speech Sig. Proces.*, Vancouver, Canada (May 26–31, 2013), pp. 448–452.</footnote>

<footnote n="25">²⁵J.-D. Polack, "Playing billiards in the concert hall: The mathematical foundations of geometrical room acoustics," *Appl. Acoust.* **38**(2), 235–244 (1993).</footnote>

<footnote n="26">²⁶M. Nolan, E. Fernandez-Grande, J. Brunskog, and C.-H. Jeong, "A wave-number approach to quantifying the isotropy of the sound field in reverberant spaces," *J. Acoust. Soc. Am.* **143**(4), 2514–2526 (2018).</footnote>

<footnote n="27">²⁷M. Nolan, M. Berzborn, and E. Fernandez-Grande, "Isotropy in decaying reverberant sound fields," *J. Acoust. Soc. Am.* **148**(2), 1077–1088 (2020).</footnote>

<footnote n="28">²⁸B. Alary, P. Massé, S. J. Schlecht, M. Noisternig, and V. Välimäki, "Perceptual analysis of directional late reverberation," *J. Acoust. Soc. Am.* **149**(5), 3189–3199 (2021).</footnote>

<footnote n="29">²⁹B. Alary, A. Politis, S. Schlecht, and V. Välimäki, "Directional feedback delay network," *J. Audio Eng. Soc.* **67**(10), 752–762 (2019).</footnote>

<footnote n="30">³⁰P. M. Morse and R. H. Bolt, "Sound waves in rooms," *Rev. Mod. Phys.* **16**, 69–150 (1944).</footnote>

<footnote n="31">³¹Note, also, that for consistency, the units of $f(t)$ are kg s⁻². Units and the normalized value in Eq. (2) are immaterial in the analysis of room reverberation time.</footnote>

<footnote n="32">³²A. Pierce, *Acoustics: An Introduction to Its Physical Principles and Applications*, 3rd ed. (Springer Nature, Switzerland, 2019).</footnote>

<footnote n="33">³³J. Allard and N. Atalla, *Propagation of Sound in Porous Media: Modelling Sound Absorbing Materials*, 2nd ed. (Wiley, New York, 2009).</footnote>

<footnote n="34">³⁴M. Aretz, P. Dietrich, and M. Vorländer, "Application of the mirror source method for low frequency sound prediction in rectangular rooms," *Acta Acust. Acust.* **100**(2), 306–319 (2014).</footnote>

<footnote n="35">³⁵A. London, "The determination of reverberant sound absorption coefficients from acoustic impedance measurements," *J. Acoust. Soc. Am.* **22**(2), 263–269 (1950).</footnote>

<footnote n="36">³⁶E. T. Paris, "Sound absorption coefficients measured by reverberation and stationary-wave methods," *Nature* **120**, 806 (1927).</footnote>

<footnote n="37">³⁷Y. Makita and T. Hidaka, "Comparison between reverberant and random incident sound absorption coefficients of a homogeneous and isotropic sound absorbing porous material - experimental examination of the validity of the revised $\cos \theta$ law," *Acta Acust. Acust.* **66**, 214–220 (1988).</footnote>

<footnote n="38">³⁸ISO 354:2003, "Measurement of sound absorption in a reverberation room" (International Organization for Standardization, Geneva, Switzerland, 2003).</footnote>

<footnote n="39">³⁹Y. Zhao, J. Xu, J. L. Davy, Z. Liu, and M. Fard, "Prediction of random incidence sound absorption coefficients of porous materials," *Appl. Acoust.* **189**, 108625 (2022).</footnote>

<footnote n="40">⁴⁰B. Alary and A. Politis, "Frequency-dependent directional feedback delay network," in *Proc. IEEE ICASSP-2020* Barcelona, Spain (May 4–8, 2020), pp. 176–180.</footnote>

<footnote n="41">⁴¹K. Rawda, V. Välimäki, and S. J. Schlecht, "Improved reverberation time control for Feedback Delay Networks," in *Proc. 22th Int. Conf. Digital Audio Effects (DAFx-19)*, Birmingham, UK (September 2–6, 2019) (Birmingham City University, Birmingham, UK, 2019).</footnote>

<footnote n="42">⁴²N. G. Kanev, "Reverberation in a trapezoidal room," *Acoust. Phys.* **59**(5), 559–564 (2013).</footnote>

<footnote n="43">⁴³See https://github.com/BenoitAlary/ImageSource2DRT60 for example code for directional reverberation time from a shoebox configuration (Last viewed February 7, 2024).</footnote>