© The Author(s), Published by EDP Sciences, 2025

1 Introduction

Room acoustic design, auralizations, and virtual reconstructions of heritage spaces all make use of numerical models to predict room acoustics [1], [2]]. In any computational model, whether it be for auralizations of historic buildings or for predicting the modal frequencies of a guitar, proper material assignment is essential to obtaining meaningful results [3–5]. Model calibration, the processof modifying material properties to achieve agreement between simulated and measured results, is a central aspect of numerical techniques across disciplines. While more recent works have considered methods for model calibration in room acoustics simulations, simple yet numerically robust methods remain elusive.

Developing and calibrating a numerical model consists of three primary steps. First, selection of a numerical method to represent the underlying physics, such as finite-element methods (FEM) or ray-tracing techniques. The numerical method must reflect the desired simulated output data, e.g., FEM modeling for an output frequency response function. The choice of numerical method then determines the material properties which must be identified, such as the Young’s, bulk, and shear moduli in FEM structural simulations [4, 5]], surface impedances in acoustic boundary element method (BEM) simulations [6], or absorption and scattering coefficients for ray-tracing techniques [3], [7]]. Second, selection of target output values to compare numerical and measured results. Often, these are a subset of the entire simulation output, such as modal frequencies and their amplitudes or derived metrics such as reverberation time (RT) or clarity. Third, design of an optimization algorithm which adjusts the material properties to attain agreement between measurements and simulation to within some specified tolerance. Algorithms may include least-squares fits, gradient-descent solvers, or machine-learning based approaches. Importantly, the first two steps, namely the choice of numerical method and the choice of output values, strongly influence the practicality of these different optimization algorithms. Often, simpler formulations lead to straightforward solutions, including closed-form formulas.

Calibration algorithms vary widely across disciplines, often being well-adapted to the numerical method and the target outcomes. In room acoustics, the most common numerical methods presently are geometrical acoustics (GA) approaches (see [8] for a review), which require random incidence absorption and scattering coefficients. The most common output values for GA model calibration have included RT metrics such as early-decay time (EDT), T ₂₀ or T ₃₀, and clarity metrics such as C ₅₀ and C ₈₀ [3], [7], [9], [10]]. However, compared to other areas of computational acoustics, GA calibration algorithms remain extremely limited, either falling into guess-and-check manual adjustments or complicated non-linear optimization approaches with sometimes dubious results.

Although several references discuss calibrating absorption coefficient values “manually in an iterative process" [9], few have provided guidelines on best practices. Martellotta et al. [7] considered “case by case" modification of materials to achieve less than 1 just noticeable difference (JND) between predicted and measured RT. In their calibration procedure for the sound absorption of audiences in churches, adjustments began “from surfaces with the most uncertain behavior and covering the largest areas."

Postma and Katz’s [3] GA calibration framework likewise applied an empirical iterative method to adjust absorption and scattering coefficient values. First, they manually adjusted absorption coefficients of the material with the largest surface area until obtaining 1 JND of T ₂₀ and the EDT. Next, they calibrated C ₅₀ by adjusting scattering coefficients, requiring one or two more iterations. Final fine-tuning applied observations of local variations due to source-receiver locations in order to adjust other material properties. Calibration of a lecture hall and church required each about six iterations in all, approximately four for the initial calibration of T ₂₀ and EDT, and another one or two to calibrate C ₅₀ and to make other minor adjustments.

Rather than manual adjustments to absorption coefficients, Christensen et al. [9] applied a genetic algorithm to optimize the absorption coefficients in a GA model. When applied to an auditorium and optimized over seven different acoustic parameters including EDT, T ₁₅, T ₂₀, T ₃₀, C ₅₀, and C ₈₀, they found that the non-linear iterative approach stabilized after eight generations. The converged result decreased deviations between simulations and measurement data to less than 1 JND for several bands and for several parameters from 500 Hz to 8 kHz. However, they also found that while the algorithm could reduce deviations, without careful initialization and parameter constraints, the calibrated absorption coefficients “did not lead to realistic materials in most cases" [9]. In addition, the computational expense for such algorithms is quite large, to the point that the authors concluded that “we do not consider the calibration utility … fast enough to be a useful tool" [9]. This is especially true considering that an expert manual calibration may only require around six iterations instead of the eight used in their work. Nonetheless, the benefits of calibrating GA models using numerical algorithms cannot be understated.

A later work by Pilch [10] compared the effectiveness of three different non-linear optimization algorithms in calibrating GA models of five rooms with at most nine different materials. The work also considered uncertainties due to ray-tracing parameters, source-receiver location, and scattering coefficients. Simulation results suggested that absorption coefficient uncertainties had the greatest impact on C ₈₀ and EDT values. The optimization algorithms consequently used an objective function based on these output parameters. For all rooms, the work found that the automated method applied to a single source position achieved less than 1 JND between 250 Hz and 4 kHz for most rooms. However, cross-validated results for a second source position were less favorable, exceeding 1 JND for four different rooms in the range of 250 Hz–4 kHz. Additionally, the work did not consider whether the optimized absorption coefficients remained realistic or not.

Pilch’s and Christensen et al.’s works illustrate that numerical algorithms have potential to assist in calibrating material properties in room acoustic simulations, just as they are used in other areas of model calibration. However, there are several unique challenges to room acoustic model calibration that remain to be addressed. One of the most significant problems is the non-uniqueness of the solution. Other forms of model calibration typically make use of many measurable output parameters to inform a few material assignments, leading to an over-determined system. This ensures the solution space is sufficiently constrained so that the fitted material properties remain physically meaningful. In contrast, room acoustics has only a few measurable output parameters to inform many material assignments. The problem is exacerbated because the orthogonality of these few output parameters, such as RT or clarity, remains ambiguous. For example, a well-calibrated T ₂₀ will likely produce a reasonably tuned T ₃₀ or EDT, so that these measures cannot independently inform calibration. Consequently, with few independent output parameters to compare against, constraining an algorithm to produce a physically meaningful solution is challenging. This issue has been noted but not adequately addressed in previous works, whether it be for algorithms based on manual adjustments or numerical optimization.

Manual adjustment has the benefit that an experienced practitioner can ensure physically meaningful calibrated absorption coefficient values. Often, understanding the underlying physical problem can inform which materials to adjust and how much to adjust them [3]. However, it is also tedious and prone to human error and bias in the selected values. Calibration performed by two different people will likely result in two different solutions. On the other hand, while complicated numerical optimization approaches provide a programmatic approach, their precise behavior may be unclear, necessitating careful initialization and parameter ranges to ensure reliable results. As a result, there is a need for calibration approaches which are simple and intuitive to control while remaining numerically robust to avoid human errors and bias.

While many works in room acoustics have considered GA model calibration, Prawda et al. [11] more recently explored the reliability of the Sabine and Eyring equations for predicting the RT. Their work showed that the inaccuracy of these equations is largely due to incorrect absorption values and neglecting air absorption rather than the inherent limitations of the underlying physicalassumptions. In their study, they achieved agreement to within ±10% (2 JND) between predicted and measured RT values once the material absorption coefficients are known. While GA models are commonly employed for the most realistic numerical simulations, these results nonetheless suggest that calibration according to simpler RT formulations may be an effective first step. The primary advantage of RT equations is that they derive from linear relationships between the absorptive materials and the room reverberation properties. Consequently, they allow simplified, closed-form calibration solutions rather than requiring complicated non-linear optimization techniques.

This work develops two methods for calibrating the absorption coefficient values in RT equations. The approaches are based on constrained least-squares and maximum-likelihood estimation. In addition to exactly matching measured and simulated RT, stochastic simulations demonstrate that the proposed methods decrease deviations between the calibrated absorption coefficients and their true values compared to the initial estimate. Finally, case applications of the RT equation calibration methods to a scattering-delay-network reverberator and to three GA models highlight the portability of the proposed method to more complicated physical models.

2 Reverberation time formulas

Over the decades, researchers have proposed several RT formulas and their variants, although most may be generally expressed as [7]

$$\begin{array}{c}\hfill T_{60}=\frac{0.164V}{Sa+4mV}\end{array}$$(1)

where S is the total surface area, a is the absorption exponent, m is sound attenuation constant of air, and V is the volume of the room. The two most common variations are the Sabine equation, with [7, 12]]

$$\begin{array}{c}\hfill a_{\mathrm{Sab}}=\frac{1}{S}\sum _{n=1}^N{S}_n{\alpha }_n=\frac{1}{S}{\mathbf{S}}^T\mathit{\alpha }\end{array}$$(2)

where S = [S ₁, S ₂, ⋯, S _N ] ^T is a column vector containing the N surface areas of each wall in the room while α  = [α ₁, α ₂, ⋯, α _N ] ^T is a column vector containing the associated N absorption coefficient values for each surface. Note that using these definitions, the total surface area S = ∥S∥₁, the L ₁ or the “Manhattan" norm, i.e., S ≠ ∥S∥₂, the “magnitude" of the vector.

The Eyring formulation uses [7], [13]]

$$\begin{array}{c}\hfill a_{\text{Eyr}}=-ln(1-a_{\text{Sab}})\end{array}$$(3)

which converges to the Sabine value for small values of a and implies a RT of 0 for a _Sab = 1.

Both formulations relate the RT through a mapping of the total absorption S ^T α . For the case of more than one absorptive material (N >  1), there are potentially infinitely many solutions producing the same RT even for fixed surface areas S. This is because increasing the absorption coefficient of one surface may be counterbalanced by decreasing the absorption coefficient of another to produce the same RT.

An underdetermined system of equations is excellent from a room acoustic design standpoint; it implies that different solutions of absorptive treatments may achieve the same desired outcome. However, this property makes calibrating room acoustic models challenging, as matching measured and estimated RT is not sufficient for ensuring proper material assignment. In fact, considering the absorption coefficients α as the free variable, all possible combinations of α which yield the same value of T ₆₀ lie on the hyperplane defined by

$$\begin{array}{c}\hfill {\mathbf{S}}^T\mathit{\alpha }=b,\end{array}$$(4)

where b is a constant determined by the selected RT model. For example,

$$\begin{array}{c}\hfill b_{\text{Sab}}=\frac{0.164V}{T_{60}}-4mV\end{array}$$(5)

for the Sabine equation while

$$\begin{array}{c}\hfill b_{\text{Eyr}}=S(1-exp(-\frac{b_{\text{Sab}}}{S}))\end{array}$$(6)

for the Eyring equation.

3 Absorption coefficient calibration methods

Ideally, the absorption coefficient values used in RT equations and in GA models would follow from careful measurements. However, accurate in-situ measurements can be time-consuming, expensive and/or require specialized equipment [14], [15]]. Consequently, researchers and practitioners typically rely on previously published absorption coefficient values when modeling rooms.

While tabulated values may provide a reasonable estimate, they often differ from actual constructed surface properties. Absorption coefficient uncertainty has long been known to be a significant challenge in room acoustic models [16], [17]]. Calibrating room acoustic models requires (1) adjusting absorption coefficient values to achieve sufficient agreement between measured and simulated RT and (2) ensuring that the calibrated absorption coefficient values are closer to their true values than the original estimate, assuming that the original estimate is incorrect. Without this latter step, adjustments made to an acoustic model may not reflect real-life changes, requiring further calibration, modifications, and measurements.

Most works and methods on calibrating room acoustic models focus primarily on the first point. However, as discussed in the preceding section, RT equations and numerical simulations allow an infinite number of combinations which yield the same RT. This means that achieving agreement between measurements and simulations is a necessary, but not sufficient condition. Indeed, while the optimized calibration explored in Christensen et al. [9] could match measurement and simulation acoustic parameters, the converged values of absorptive materials were sometimes found to be not physically realistic if allowed to freely vary. Care must be taken to ensure that the calibrated values are closer to the ground-truth values than the initial estimate.

This section proposes two approaches to calibrate a set of estimated absorption coefficients $\tilde{\mathit{\alpha }}$ to yield the desired RT. The first applies a constrained least-squares estimate which can be shown to always improve upon the initial estimate with respect to the L ₂ norm. The second applies a constrained maximum-likelihood estimate which allows incorporation of varying material uncertainty. Both solutions assume measured or known values of T ₆₀, S, V, and m and an initial estimate $\tilde{\mathit{\alpha }}$ so that the only free parameters are the absorption coefficients α .

3.1 Least-squares calibration

Suppose one has an initial estimate of the absorption coefficient values $\tilde{\mathit{\alpha }}$, perhaps from previous measurements or from tabulated values in an engineering handbook such as [18–20]. When applied to the RT equations in Section 2, these values provide an initial estimate of the reverberation time ${\tilde{T}}_{60}$. Model calibration seeks to adjust the initial estimate $\tilde{\mathit{\alpha }}$ to a set of calibrated coefficients $\widehat{\mathit{\alpha }}$ that achieve agreement with the measured value T ₆₀.

Geometrically, the point $\tilde{\mathit{\alpha }}$ lies in an N-dimensional space, while the set of absorption coefficient combinations which yield the measured T ₆₀ lie in an (N − 1)-dimensional hyperplane. The least-squares estimate (LSE) seeks to find the point $\widehat{\mathit{\alpha }}$ which lies in this plane and is closest to the initial estimate. Formally,

$$\begin{array}{cc}\hfill \text{minimize}& ||\mathit{\alpha }-\tilde{\mathit{\alpha }}{||}_2^2\hfill \\ \hfill \text{subject to}& {\mathbf{S}}^T\mathit{\alpha }=b\hfill \\ \hfill & 0\le {\alpha }_n\le 1,\hfill \end{array}$$(7)

where b is the constant associated with the measured T ₆₀ (see Eqs. (5) and (6)). In other words, among all combinations of absorption coefficient values which produce the measured T ₆₀, the LSE finds the values $\widehat{\mathit{\alpha }}$ with the least amount of modification from the original estimate $\tilde{\mathit{\alpha }}$.

Momentarily setting aside the inequality constraints (see Sect. 3.3), Lagrange multipliers provide the solution to this constrained optimization problem. The Lagrangian may be expressed as

$$\begin{array}{c}\hfill \mathcal{L}(\mathit{\alpha },\lambda )=\frac{1}{2}{(\mathit{\alpha }-\tilde{\mathit{\alpha }})}^T(\mathit{\alpha }-\tilde{\mathit{\alpha }})+\lambda ({\mathbf{S}}^T\mathit{\alpha }-b)\end{array}$$(8)

with λ the Lagrange multiplier. Differentiating with respect to α and setting the result to 0 leads to the result

$$\hat{\mathit{\alpha }}=\tilde{\mathit{\alpha }}-\left(\frac{\tilde{b}-b}{{\mathit{S}}^T\mathit{S}}\right)\mathit{S}$$(9)

which is recognizable as the projection of $\tilde{\mathit{\alpha }}$ onto the plane S ^T α  = b. Thus, the LSE calibrated absorption coefficients are

$$\begin{array}{c}\hfill {\widehat{\alpha }}_n={\tilde{\alpha }}_n-C{S}_n\end{array}$$(10)

where C is a constant for the Sabine equation equal to

$$\begin{array}{c}\hfill C_{\text{Sab}}=\frac{0.164V}{{\sum }_{m=1}^N{S}_m^2}(\frac{1}{{\tilde{T}}_{60}}-\frac{1}{T_{60}})\end{array}$$(11)

and for the Eyring equation equal to

$$\begin{array}{c}\hfill C_{\text{Eyr}}=\frac{S}{{\sum }_{m=1}^N{S}_m^2}(exp(-\frac{b_{\text{Sab}}}{S})-exp(-\frac{{\tilde{b}}_{\text{Sab}}}{S}))·\end{array}$$(12)

The LSE formula contained in equation (10) highlights several important features of this calibration method. First, the LSE modifies each absorption coefficient proportional to its surface area S_n. Because the value of C is the same for all coefficients, absorption coefficients associated with larger surface areas receive greater adjustments than those with smaller surface areas. Interestingly, this concept is similar to the manual calibration approaches used by Martellotta et al. [7] and Postma and Katz [3], which suggested adjusting larger surfaces areas first as they have more significant impact on the RT result.

Second, the constant C, with units of absorption per unit area, represents the amount of absorption correction per unit area that must be applied to achieve the measured RT. Using this interpretation, the LSE first calculates the total amount of absorption correction necessary to obtain the measured T₆₀ and then applies this correction proportional to the surface areas associated with each absorption coefficient. Thus, unlike previously proposed numerical optimization calibration procedures, the LSE has a clear physical interpretation.

Figure 1 theoretically illustrates this calibration procedure for the case of a room covered with two different materials (N = 2). In the model, the floor of the 5 × 4 × 3 m room (V = 60 m³) is covered with a material with α₁ = 0.2. All other surfaces are covered with a material with α₂ = 0.4, leading to an T₆₀ = 0.164(60)/(0.2(20)+0.4(74)) = 0.29 s, neglecting air absorption. The blue circle marked by α^* in the α₁–α₂ plane represents this value. However, there are infinitely many combinations of α₁ and α₂ which could produce this same value of T₆₀; the solid black line represents all of these combinations.

Figure 1 Two-dimensional illustration of the least-squares calibration.

The red circle marked by $\tilde{\mathit{\alpha }}$ represents the estimated absorption coefficients used before measurement. The values of α ₁ = 0.4 and α ₂ = 0.6 would lead to a ${\tilde{T}}_{60}=0.19$ s. However, after a theoretical measurement of the actual T ₆₀ = 0.29 s, equation (10) allows calibration of the original coefficients $\tilde{\mathit{\alpha }}=(0.4,0.6)$ to the LSE represented by $\widehat{\mathit{\alpha }}=(0.34,0.36)$. The red “x" represents this solution; it falls on the space of all combinations of α yielding the measured T ₆₀ time. The dashed black line illustrates that $\widehat{\mathit{\alpha }}$ is the orthogonal projection of $\tilde{\mathit{\alpha }}$ onto this space.

Several important trends may be understood from this simplified two-dimensional case. First, since geometrically ${\Vert \tilde{\mathit{\alpha }}-{\mathit{\alpha }}^{\ast }\Vert }_2>{\Vert \widehat{\mathit{\alpha }}-{\mathit{\alpha }}^{\ast }\Vert }_2$ assuming ${\tilde{T}}_{60}\ne T_{60}$, the LSE-calibrated absorption values always decrease the overall deviation between the estimated and true values. For example, in the present case, ${\Vert \tilde{\mathit{\alpha }}-{\mathit{\alpha }}^{\ast }\Vert }_2=0.28$ while ${\Vert \widehat{\mathit{\alpha }}-{\mathit{\alpha }}^{\ast }\Vert }_2=0.14$. This feature is an essential requirement of a calibration approach because it guarantees improved absorption coefficient values.

Second, the values of the calibrated coefficients largely depend on the initial estimate $\tilde{\mathit{\alpha }}$. The deviation between the calibrated and actual values ${\Vert \widehat{\mathit{\alpha }}-{\mathit{\alpha }}^{\ast }\Vert }_2$ depends on the distance between $\tilde{\mathit{\alpha }}$ and α ^* in the direction tangent to the line S ^T α  = b. Consequently, care should always be taken to obtain the most reliable initial estimate of the absorption coefficient parameters $\tilde{\mathit{\alpha }}$ to achieve the best results. This result highlights that the RT calibration process improves with reliable initial estimates of material properties; it does not replace the need for careful measurements.

Finally, the projection always translates the vector $\tilde{\mathit{\alpha }}$ in the direction of −S as expressed in equations (9) and (10). This detail can be problematic because all coefficients are adjusted relative to their respective surface areas by a global correction represented by the constant C. Thus, the approach cannot achieve individual variations among calibrated values in ${\widehat{\alpha }}_n$. For example, consider the case given in Figure 1 but using ${\tilde{\alpha }}_1=0.2$ as the known true value. As a result, only the value ${\tilde{\alpha }}_2$ needs to be adjusted. However, the LSE will adjust both values in obtaining a solution. While the overall deviation will decrease, the error in α ₁ will increase while the error in α ₂ will decrease. This feature is undesirable in situations where some coefficients are known with great certainty and do not need further calibration or modifications. The following section explores this concept further and presents a maximum-likelihood calibration approach to overcome this limitation.

3.2 Maximum-likelihood calibration

The LSE calibration procedure follows from an initial estimate of absorption materials $\tilde{\mathit{\alpha }}$. However, often some values of absorption are known with greater certainty than others. Materials with higher uncertainty should see a greater relative change after calibration than those with little uncertainty, as discussed by Martellotta et al. [7] in their manual calibration approach. Incorporating the individual degrees of uncertainty into the estimate could provide better results than the LSE. This section develops a maximum-likelihood calibration procedure to meet this need.

The first step is to replace a deterministic description of the absorption coefficients with a probabilistic one. Under this model, rather than providing an initial fixed estimate for $\tilde{\mathit{\alpha }}$, one provides a probability density function. Because absorption is limited to 0 ≤ α ≤ 1, the truncated normal distribution, whose probability density function (PDF) is defined as [21]

$$\begin{array}{cc}\hfill f& (\alpha ;\mu ,\sigma )=\{\begin{array}{cc}\frac{1}{\mathrm{\Phi }(\mu ,\sigma )\sigma \sqrt{2\pi }}exp(-\frac{1}{2}{\left(\frac{\alpha -\mu }{\sigma }\right)}^2)\hfill & 0\le \alpha \le 1\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}\hfill \end{array}$$(13)

where

$$\begin{array}{c}\hfill \mathrm{\Phi }(\mu ,\sigma )=\frac{1}{2}[\text{erf}\left(\frac{1-\mu }{\sqrt{2}\sigma }\right)-\text{erf}\left(\frac{-\mu }{\sqrt{2}\sigma }\right)]\end{array}$$(14)

is a normalization factor, serves as a convenient parametrized model. Using the parameters μ and σ, the mean and standard deviation of the untruncated normal distribution, one may express the relative certainty of different absorption coefficient values. In this way, one speaks of the probability of the absorption coefficient falling between certain ranges.

To illustrate how the truncated normal distribution may represent material uncertainty, Figure 2 plots the PDFs for several parameter sets. The red curve marks the PDF for the case μ = 0.2 and σ = 0.02. In this case, the low standard deviation indicates a low absorption coefficient uncertainty. These parameter values indicate that the material has a 99% chance of having an absorption coefficient value between 0.2 ± 0.05. The green curve marks the PDF for μ = 0.5 and σ = 0.05. The higher σ value corresponds to a broader distribution width and higher uncertainty in the absorptive properties. For example, there is only a 68% chance that α has a value between 0.5 ± 0.05. For the blue curve, the even larger value of σ = 0.1 would translate to a 38% chance of the absorption coefficient falling between 0.7 ± 0.05.

Figure 2 Three truncated normal distributions with different parameters. Higher values of σ indicate lower certainty in the absorption coefficient value.

With a model of the probability of each individual absorption coefficient, one may describe the probability of any combination of α through a likelihood function

$$\begin{array}{c}\hfill L_N(\mathit{\alpha })=\prod _{n=1}^N{f}_n({\alpha }_n;{\mu }_n,{\sigma }_n)\end{array}$$(15)

where f_n are the individual PDFs of each absorption coefficient. The maximum probability occurs when α = μ, i.e., when assigning absorption coefficients to their untruncated mean value. However, calibration requires that the absorption coefficients produce the measured T₆₀. Consequently, one seeks to find the maximum likelihood on the hyperplane defined by S^Tα = b. In practice, it is easier to maximize the log-likelihood [22]

$$\begin{array}{cc}{l}_N(\mathit{\alpha })& =\log [L_N(\mathit{\alpha })]\end{array}$$(16)

$$\begin{array}{c}=\sum _{n=1}^N-\frac{1}{2}{\left(\frac{{\alpha }_n-{\mu }_n}{{\sigma }_n}\right)}^2+C_n\end{array}$$(17)

where

$$\begin{array}{c}\hfill C_n=-log\left(\mathrm{\Phi }({\mu }_n,{\sigma }_n){\sigma }_n\sqrt{2\pi }\right)\end{array}$$(18)

is a constant. The optimization is consequently

$$\begin{array}{cc}\hfill \text{maximize}& l_N(\mathit{\alpha })\hfill \\ \hfill \text{subject to}& {\mathbf{S}}^T\mathit{\alpha }=b\hfill \\ \hfill & 0\le {\alpha }_n\le 1.\hfill \end{array}$$(19)

Loosely speaking, among the set of absorption coefficients which yield the desired T ₆₀, the maximum-likelihood estimate (MLE) seeks to find the combination $\widehat{\mathit{\alpha }}$ which is most likely to occur given the probabilistic model.

Momentarily setting aside the inequality constraints (see Sect. 3.3), the solution follows from Lagrange multipliers. The Lagrangian is

$$\begin{array}{c}\hfill \mathcal{L}(\mathit{\alpha },\lambda )=(\sum _{n=1}^N-\frac{1}{2}{\left(\frac{{\alpha }_n-{\mu }_n}{{\sigma }_n}\right)}^2+C_n)+\lambda ({\mathbf{S}}^T\mathit{\alpha }-b).\end{array}$$(20)

Differentiating with respect to α and setting the result equal to zero yields

$$\begin{array}{c}\hfill {\widehat{\alpha }}_n-{\mu }_n=\lambda S_n{\sigma }_n^2.\end{array}$$(21)

Defining a new vector w by

$$\begin{array}{c}\hfill w_n=S_n{\sigma }_n^2\end{array}$$(22)

and multiplying the vector form of equation (21) by S ^T leads to the result

$$\begin{array}{cc}\hat{\mathit{\alpha }}& =\mathit{\mu }-\frac{\tilde{b}-b}{S^{T}w}w\\ & =\mathit{\mu }-\left(\frac{\tilde{b}-b}{S^{T}KS}\right)S_n{\sigma }_n^2,\end{array}$$(23)

where

$$\begin{array}{c}\hfill \mathbf{K}=\left[\begin{array}{cccc}{\sigma }_1^2& 0& \cdots & 0\\ [5pt]0& {\sigma }_2^2& \cdots & 0\\ [5pt]⋮& ⋮& \ddots & ⋮\\ [5pt]0& 0& \cdots & {\sigma }_N^2\end{array}\right]\end{array}$$(24)

is the covariance matrix and $\tilde{b}$ is the constant associated with using ${\tilde{\alpha }}_n={\mu }_n$, i.e., estimating T ₆₀ by using the mean values of the untruncated normal distributions. Thus, the MLE-calibrated absorption coefficients are

$$\begin{array}{c}\hfill {\widehat{\alpha }}_n={\mu }_n-D{S}_n{\sigma }_n^2\end{array}$$(25)

where D is a constant equal to

$$\begin{array}{c}\hfill D_{\text{Sab}}=\frac{0.164V}{{\sum }_{n=1}^N{\sigma }_n^2{S}_n^2}(\frac{1}{{\tilde{T}}_{60}}-\frac{1}{T_{60}})\end{array}$$(26)

for the Sabine equation and

$$\begin{array}{c}\hfill D_{\mathrm{Eyr}}=\frac{S}{{\sum }_{n=1}^N{\sigma }_n^2{S}_n^2}(exp(-\frac{b_{\text{Sab}}}{S})-exp(-\frac{{\tilde{b}}_{\text{Sab}}}{S}))\end{array}$$(27)

for the Eyring equation.

The solution to this MLE is the LSE but weighted by the variances σ_n² of each absorption coefficient. As a result, it can be interpreted as a weighted least-squares estimate (WLSE) [22]. This has several practical implications. First, when the variances or uncertainties of each material property are equal (σ_n² = σ²), the equation reduces to the LSE of equation (10) so that the LSE is the MLE with all absorption coefficients having equal uncertainties. Second, if a material is known with perfect certainty so that σ_n → 0, then that material property will be unchanged by the MLE so that ${\widehat{\alpha }}_n={\mu }_n$. Lastly, the MLE modifies each absorption coefficient proportional to its surface area and its uncertainty represented by σ _n . Thus, absorption coefficients with higher uncertainty will be modified more per unit area than those with lower uncertainty. This contrasts with the LSE, where the modification per unit area is the same for all coefficients. Similar to the LSE calibration, the closed-form MLE solution contrasts with numerical optimization procedures as it has clear intuitive meaning to its result.

To illustrate these concepts, Figure 3 shows the MLE for the same theoretical room as in Figure 1 except that half of the total surface area is covered with material with α₁^* = 0.6 while the other half is covered with material with α₂^* = 0.2. This point in the α₁–α₂ plane appears as a cyan dot. The orange line indicates all other combinations of α which produce the same T₆₀. In Figure 3, the underlain colormap plots the normalized likelihood function L_N for the case that μ = (0.6, 0.5) and σ = (0.05, 0.05). A red circle indicates this mean value. Because the variance or uncertainty for both materials is the same, contours of equal-likelihood appear as concentric circles about the mean value, indicated by thin white lines. Because σ₁ = σ₂, the MLE is the same as the LSE. Consequently, the projection is orthogonal and represents the nearest combination of α that satisfies the desired T₆₀. However, the LSE modifies the correct initial estimate of the first material’s absorption coefficient (μ₁ = α₁^* = 0.6) to an incorrect value of ${\widehat{\alpha }}_1=0.45$. Of course, this loss in accuracy for the first material is compensated by the fact that the second material’s calibrated absorption coefficient $\widehat{{\alpha }_2}=0.35$ is now much closer to its true value of α ₂ ^* = 0.2; the overall relative error is improved from ∥ μ  −  α ^*∥₂/∥ α ^*∥₂ = 47% to ${\Vert \widehat{\mathit{\alpha }}-{\mathit{\alpha }}^{\ast }\Vert }_2/{\Vert {\mathit{\alpha }}^{\ast }\Vert }_2=34\%$.

Figure 3 Maximum-likelihood estimation for calibrating absorption coefficient values for a two-dimensional problem.

However, by modifying the variance σ_n² of each absorption coefficient, the MLE can improve upon the LSE. Figure 3b now shows the case that σ₂ is increased to 0.1 while σ₁ = 0.05 remains the same. The underlying likelihood function stretches further in the α₂ direction to indicate this increased uncertainty. Contours of constant likelihood now appear as concentric ellipses. This case indicates a situation where the first material property is known with greater certainty than the second. Consequently, the MLE projection is no longer orthogonal;there is a larger adjustment to the uncertain α₂ value compared to the more certain α₁ value. The MLE for this case becomes $\widehat{\mathit{\alpha }}=(0.54,0.26)$, much closer to the true values of α ^* = (0.6, 0.2) than the LSE. The deviation between calibrated and true coefficients ${\Vert \widehat{\mathit{\alpha }}-{\mathit{\alpha }}^{\ast }\Vert }_2/{\Vert {\mathit{\alpha }}^{\ast }\Vert }_2$ decreases from the LSE value of 34% to 13%. Further increasing the uncertainty in α ₂ by increasing σ ₂ to 0.2 improves again upon the MLE estimate as shown by Figure 3c. The calibrated values $\widehat{\mathit{\alpha }}=(0.58,0.22)$ led to a deviation ${\Vert \widehat{\mathit{\alpha }}-{\mathit{\alpha }}^{\ast }\Vert }_2/{\Vert {\mathit{\alpha }}^{\ast }\Vert }_2$ of only 4%. Consequently, the MLE can significantly improve absorption coefficient calibration when uncertainty is properly incorporated.

This simple example illustrates the importance of understanding absorption coefficient uncertainty in calibrating RT equations. Although the initial estimate had a 47% relative error in the values of the absorption coefficients, judiciously applying material uncertainty placed the calibrated result to within 4%. This also highlights the difficulty in calibrating room acoustic models; all three solutions exactly reproduce the desired T ₆₀ time, yet the estimated coefficients varied significantly due to the undetermined system of equations.

3.3 Handling absorption coefficient constraints

The preceding developments neglected the additional optimization constraint that 0 ≤ α _n  ≤ 1. Handling these inequality constraints adds some additional complexity to determining a proper solution. Because the LSE and MLE calibrations lead to simple closed-form solutions, one should first apply these formulas. If all calibrated coefficients lie within the proper bounds, then the inequality constraints are inactive and may be ignored. With good initial estimates so that $\tilde{\mathit{\alpha }}$ lies close to α ^*, this case is assumed to be the most common in application.

However, in some situations, such as when large surface areas are covered in material with absorption values approaching either extremity (e.g., churches with primarily stone walls), the likelihood of the inequality constraints becoming active becomes higher. In the case that the calibration produced by the LSE or MLE without the inequality constraints violates 0 ≤ α _n  ≤ 1, then one or more of these constraints will be active and the solution requires further refinement. One possibility to avoid turning to a numerical solution is to increase the absorption coefficient value certainty of materials in violation by decreasing the variance σ _n . This in effect fixes the associated value of ${\widehat{\alpha }}_n$, although extreme care must be taken so that other calibrated values do not suffer from overcompensating for the fixed values. This approach is also problematic when the offending materials are the most uncertain and do require significantadjustments.

If this heuristic approach does not succeed, then using numerical optimization approaches such as gradient descent and its variants [22] becomes necessary. Nonetheless, a concave log-likelihood function, a closed-form solution to its gradient, and well-defined equality and inequality constraints allow straightforward implementation into most numerical software.

3.4 Stochastic simulations

The constrained optimization employed in the proposed calibration procedures always ensures that the calibrated set of coefficients $\widehat{\mathit{\alpha }}$ produces the measured T ₆₀. The illustrative examples given in the preceding sections demonstrate that both the LSE and MLE calibrations have the potential to not only match the target RT but also to yield improved agreement between calibrated and actual absorption coefficient values. However, it remains to evaluate their robustness over varying room conditions and material uncertainties. This section proposes stochastic (Monte-Carlo) simulations to assess each calibration scheme.

The stochastic simulations provide calibration results over many (∼ <  S P S D O U B L E D O L L A R >  10⁸) randomly generated room conditions. If, on average, the calibrated absorption coefficients are closer to the ground-truth values than the initial estimate, then the calibration method may be considered reliable. In contrast, a method which provides calibrated absorption coefficient values further from the ground-truth values may be considered unreliable. To quantify whether a calibration method is reliable or not, this section develops a calibration score (CS) based on results from the stochastic simulations. The CS will be positive for methods which bring calibrated coefficients closer to the ground-truth and negative for methods which lead to further deviations. A negative CS consequently warns of methods which can match target room acoustic parameters but may provide unrealistic or unphysical absorption coefficientvalues.

The numerical simulations model rooms in the following manner. Consider a room bounded by N surfaces each with surface area S _n . A room realization consists of randomly assigning a mean absorption coefficient value μ _n and a standard deviation value σ _n to each surface in the room. For the present work, μ _n are drawn from a uniform distribution over (2σ _n , 1 − 2σ _n ) while σ _n are drawn from a Rayleigh distribution whose scale parameter σ _R will vary over the simulations to analyze robustness under different uncertainty conditions. The standard deviations σ _n follow a Rayleigh distribution because this choice allows a simplified expression for the CS and because it provides a more realistic representation of model variance than a uniform distribution. These selected values of μ and σ represent the a-priori knowledge available to practitioners in reference compilations such as [18–20]. For further information, see Deetz and Boren [23] for a discussion on estimating material uncertainty σ from tabulated values.

Next, within a given room realization, a trial run consists of randomly drawing absorption coefficient values α ^* from the distributions defined by the now-fixed values of μ and σ . This in essence generates a room whose material properties will be similar to μ but will have individual variations according to the variance expressed in σ . These values represent the ground-truth absorption coefficients and produce the trial run’s “measured" T ₆₀. The calibration then employs this T ₆₀ value to produce calibrated absorption coefficients using the mean values as the initial estimate ($\tilde{\mathit{\alpha }}=\mathit{\mu }$).

Because the constrained optimization guarantees that the LSE and MLE methods provide absorption coefficients which exactly match the measured T ₆₀, of primary interest is whether the calibrated values $\widehat{\mathit{\alpha }}$ are closer to the ground-truth values α ^* compared to the initial estimate μ . To assess these differences, the mean-squared error (MSE) in the calibrated absorption coefficient values is given as

$$\begin{array}{c}\hfill \mathrm{MSE}=\frac{1}{N}||\widehat{\mathit{\alpha }}-{\mathit{\alpha }}^{\ast }{||}_2^2.\end{array}$$(28)

This value quantifies the calibration errors for a single trial run of a single room realization. The MSE is a random variable which differs between each trial run because α ^*, and consequently T ₆₀ and $\widehat{\mathit{\alpha }}$, change with each run. Averaging the MSE over many trial runs for a fixed room realization allows an estimate of its statistical properties. Over M trial runs, the expected value of the MSE (EMSE) may be approximated by a sample mean as

$$\begin{array}{c}\hfill \mathrm{EMSE}\approx \frac{1}{M}\sum _{m=1}^M{\mathrm{MSE}}_m\end{array}$$(29)

where MSE _m is the MSE of a single trial run. A smaller value of EMSE signifies an on-average higher accuracy of a given calibration method.

Next, the EMSE of a calibration scheme should be compared against that produced from no calibration. This baseline EMSE value follows by using $\widehat{\mathit{\alpha }}=\mathit{\mu }$. A calibration method with a smaller EMSE than this baseline value indicates a calibration which, on average, improves upon the initial estimate of $\tilde{\mathit{\alpha }}=\mathit{\mu }$. On the other hand, a calibration method with a larger EMSE than the baseline indicates that the calibrated absorption coefficients are, on average, worse than the initial estimate of $\tilde{\mathit{\alpha }}=\mathit{\mu }$.

While the EMSE values for the calibration methods will follow from the stochastic simulations, the baseline EMSE may be determined analytically as

$$\begin{array}{cc}\hfill {\mathrm{EMSE}}_{\mathrm{baseline}}& =E[\frac{1}{N}||\mathit{\mu }-{\mathit{\alpha }}^{\ast }{||}_2^2]\hfill \\ \hfill & =\frac{1}{N}\sum _{n=1}^{\mathit{NE}}[{({\mu }_n-{\alpha }_n^{\ast })}^2]\hfill \\ \hfill & =\frac{1}{N}\sum _{n=1}^N{\tilde{\sigma }}_n^2\hfill \end{array}$$(30)

where E[⋅] is the expectation operator and ${\tilde{\sigma }}_n^2$ is the variance of the nth truncated normal distribution. The variance of the truncated normal distributions relate to the variance of the untruncated normal distributions σ _n ² as [21]

$$\begin{array}{c}\hfill {\tilde{\sigma }}_n^2={\sigma }_n^2[1-\frac{\beta \varphi (\beta )-\gamma \varphi (\gamma )}{\mathrm{\Phi }({\mu }_n,{\sigma }_n)}-{\left(\frac{\varphi (\gamma )-\varphi (\beta )}{\mathrm{\Phi }({\mu }_n,{\sigma }_n)}\right)}^2],\end{array}$$(31)

where $\varphi (x)=1/\sqrt{2\pi }exp(-1/2x^2)$, γ = −μ _n /σ _n , and β = (1 − μ _n )/σ _n . In the case that σ _n is small, ${\tilde{\sigma }}_n\approx {\sigma }_n$, i.e., the truncated normal distribution has the same variance as a (untruncated) normal distribution.

After completing a realization by running the M trial runs and computing EMSE for each calibration, a new room realization follows by randomly selecting new values of μ and σ for each surface. Similar to how the MSE value varies for each trial run, the EMSE value varies for each room realization. The average absorption coefficient error (AACE) follows from an estimate of the expected value of EMSE by averaging EMSE values over R realizations so that

$$\begin{array}{c}\hfill \text{AACE}=\sqrt{\frac{1}{R}\sum _{r=1}^R{\text{EMSE}}_r}.\end{array}$$(32)

Whereas the EMSE represents the calibration’s robustness for a single choice of μ and σ , AACE represents its robustness across a variety of absorption coefficient combinations and uncertainty levels for a fixed room geometry.

Because the values of ${\tilde{\sigma }}_n$ depend on the choice of μ which is itself a random variable, analytically determining the baseline value for AACE is more involved. However, using ${\tilde{\sigma }}_n\approx {\sigma }_n$ yields an approximate result. Since the values of σ _n followed from a Rayleigh distribution with shape parameter σ _R , the sum of the variances

$$\begin{array}{c}\hfill Z=\sum _{n=1}^N{\sigma }_n^2\end{array}$$(33)

follows a Gamma distribution with shape parameters N and 1/2σ _R ². The expected value of Z is thus 2N σ _R ², so that

$$\begin{array}{cc}\hfill {\text{AACE}}_{\text{baseline}}& =\sqrt{E\left[{\text{EMSE}}_{\text{baseline}}\right]}\hfill \\ \hfill & =\sqrt{\frac{1}{N}E\left[\sum _{n=1}^N{\tilde{\sigma }}_n^2\right]}\hfill \\ \hfill & \approx \sqrt{\frac{1}{N}E[Z]}\hfill \\ \hfill & \approx \sqrt{2}{\sigma }_R.\hfill \end{array}$$(34)

Thus, it is clear that for an entire simulation for a fixed room, the parameter σ _R influences the overall uncertainty. For example, setting σ _R  = 0.1 will produce, on average, PDFs with σ ≈ 0.125, with resultant shapes similar to the blue dash-dot curve in Figure 2. The baseline AACE will then be ≈0.14. On the other hand, choosing a σ_R = 0.2 with PDF shapes more similar to the red solid curve in Figure 2 will lead to a baseline AACE of ≈0.03.

Lastly, with a suitable baseline value, the CS

$$\begin{array}{cc}\hfill {\text{CS}}_{{\sigma }_R}& =(1-\frac{\text{AACE}}{{\text{AACE}}_{\text{baseline}}})\hfill \\ \hfill & \approx (1-\frac{\text{AACE}}{\sqrt{2}{\sigma }_R}),\hfill \end{array}$$(35)

where the last line may be used for small values of σ _R , represents the improvement of a given calibration method relative to no calibration. The CS may be either positive or negative. A positive percentage indicates that a calibration improves upon the initial estimate of the absorption coefficient values, with a higher value indicating a more accurate approach. A negative percentage indicates that the calibration method yields a worse estimate of the absorption coefficient values relative to the initial estimate. Consequently, the CS provides a metric to assess whether a calibration method provides realistic absorption coefficient values or not.

As a reference calibration algorithm from the literature to compare against the two proposed methods, one simple method [24] scales absorption coefficient values so that their total absorption equals the total measured absorption (${\mathbf{S}}^T{\mathit{\alpha }}^{\ast }={\mathbf{S}}^T\widehat{\mathit{\alpha }}$). This yields the formula

$$\begin{array}{c}\hfill {\widehat{\alpha }}_n^{\mathrm{REF}}={\tilde{\alpha }}_{\mathit{nG}}\end{array}$$(36)

where G is a constant equal to Geometrically, this form of calibration finds the intersection of the vector in the direction as $\tilde{\mathit{\alpha }}$ with the hypersurface S ^Tα = b. For example, in the two-dimensional example in Figure 1, the calibrated values would fall at the point of the intersection between the line S^Tα = b and a line drawn between the origin (0, 0) and the estimated values $({\tilde{\alpha }}_1,\widetilde{{\alpha }_2})$. Unlike the LSE and MLE, this form of calibration does not weight the calibrated coefficients based on surface areas but on the total absorption deviation. Consequently, one anticipates higher errors for rooms with irregular surface area distributions.

$$\begin{array}{c}G=\frac{S^T{\mathit{\alpha }}^{\ast }}{S^T\tilde{\mathit{\alpha }}}=\frac{b}{\tilde{b}}.\end{array}$$(37)

4 Results

4.1 Stochastic simulations

Deriving the CS for a given room geometry and material variability requires several steps. The following paragraphs demonstrate the calculation process of the CS before presenting general results for a shoebox room and a lecture hall. First, consider the case of a simple shoebox room of dimensions 5 × 7 × 3 8 m with σ_R = 0.05. Table 1 lists the assigned absorption coefficient mean values μ_n and standard deviations σ_n for a single room realization. The associated PDFs for this room realization appear in Figure 4. From these PDFs, the randomly drawn values α_n^* for a single trial run appear in Table 1. These values produce three sets of calibrated absorption coefficients, using the reference calibration approach (Eq. (36)), LSE, and MLE methods.

Figure 4 (a) PDFs for the six material properties for a single room realization. (b) Statistical properties of the EMSE over the 104 room realizations incorporating 104 trial runs each. The boxes indicate the 25th and 75th percentiles while the whiskers indicate the extreme values for each case.

The MSE for this trial run appears at the foot of Table 1. Without effort at calibration, the MSE was 34.1×10⁻⁴, whereas the reference, LSE, and MLE methods decreased the MSE to 26.9×10⁻⁴, 23.6×10⁻⁴, and 20.4×10⁻⁴, respectively. Thus, for this specific room realization and trial run, the MLE performed best, followed by the LSE and then the reference calibration procedure.

Next, the EMSE follows by running 10⁴ trial runs (draws of α^*) for this room configuration (fixed values of μ and σ). This amount of trial runs was sufficient to consistently attain convergence between the numerically estimated value of EMSE_baseline and the analytic value given from equation (30) (26.8×10⁻⁴ for the present room configuration). Averaging over many trial runs gives a clearer picture of the relative performance of each calibration method.

Finally, EMSE values must be calculated for many different realizations of the same room geometry. After 10⁴ room realizations (10⁸ different simulated rooms in all), the AACE converges to values of 6.8×10⁻², 6.5×10⁻², 6.2×10⁻², and 5.8×10⁻². Note that the numerically estimated value of AACE_baseline falls close to the analytic approximation of $\sqrt{2}{\sigma }_R=7.1\times {10}^{-2}$. The aggregated statistics of the square root of the EMSE over the 10⁴ room realizations appear in Figure 4b. In each box-whisker plot, the AACE (roughly the mean value) appears as a red line; the box and whiskers mark the quartile ranges. The overlaid dashed black line indicates the analytical approximation to AACE${}_{\mathrm{baseline}}\approx \sqrt{2}{\sigma }_R=0.07$.

Among the 10⁴ room realizations, there were significant variations in calibration performances. However, averaged across realizations, no calibration tended to place the calibrated coefficients within 6.8 × 10⁻² of their true value for this specific choice of σ _R . For the calibration methods, the reference, LSE, and MLE improve upon this to achieve average deviations of only 6.5×10⁻², 6.2×10⁻² and 5.8×10⁻². From these numerically estimated AACE values, the CS become 5.5%, 8.8%, and 14.3% for the reference, LSE, and MLE methods, respectively. For this simple rectangular room geometry, the positive CS scores would indicate that each calibration method is effective at improving the estimate of absorption coefficient values from the initial guess. The relative magnitudes of the CS scores indicate that the MLE is the most accurate approach, followed the the LSE and the reference calibration method.

This same process is repeated for both the shoebox model with dimensions 5 × 7 × 3 m (N = 6, V = 105 m³) and a lecture hall (N = 23, V = 334 m³). These rooms are shown in Figure 5. Although in practice many surfaces could share similar absorption coefficient properties, the simulations allowed each surface to take on different values for the sake of evaluating calibration robustness. The simulations employed 10⁴ room realizations each with 10⁴ trial runs for three different levels of uncertainty: σ_R = 0.01, σ_R = 0.05, and σ_R = 0.1.

Figure 5 (a) Shoebox room and (b) lecture hall used in the stochastic simulations.

Table 2 records the calibration scores CS for the shoebox room. Several notable trends appear. First, the CS for all three methods remains positive for all three variance levels. This indicates that all three methods are effective calibration approaches for the room. Interestingly, the CS of the LSE remains fairly constant across variance uncertainty. The MLE’s CS decreases with increasing uncertainty, whereas the reference method’s CS increases with increasing uncertainty. Overall, the best methods for the room were the MLE, followed by the LSE, and lastly the reference method.

Table 3 reports the calibration scores CS for the lecture hall. Interestingly, the CS for the reference method is slightly negative. This result is unsurprising when considering that the lecture hall has large variances in room surface areas. For example, the ratio of the shoebox room’s largest and smallest surface is 2.3, whereas for the lecture hall it is 66.6. Because the reference method calibrates absorption coefficients without considering the relative surface area sizes, it will be prone to large errors when there are many surfaces with different sizes. Consequently, while the reference method can match measured and simulated and RT, for this room it is a poor calibration choice because the calibrated absorption coefficients, on average, are further from the true values than the initial estimate.

The MLE and LSE methods have positive CS scores like they did for the shoebox room. In fact, geometric arguments guarantee that the LSE will always have a positive CS. However, the CS scores are significantly lower than for the shoebox room. This result likely derives from the large number of material properties to predict based on a single output value.

4.2 Case study in scattering-delay-network reverberator calibration

Scattering Delay Networks (SDNs) are a class of efficient digital artificial reverberators (DARs) based on scattering junction nodes placed on the walls of a simulated room. The placement and delays between nodes, based on first-order image sources, ensure exactly rendered first-order reflections, while random-incidence absorption coefficients for each wall determine temporal energy losses. When De Sena et al. [25] introduced SDN reverberators, they showed that the RT of an SDN nearly exactly agrees with that produced from the image source method (ISM) and are in close agreement with the Sabine and Eyring RT equations. Previous works on calibrating DARs employed complicated algorithms, including machine learning, neural networks, or Bayesian optimization (see [26] for a review and further discussion). However, the close relationship between the RT of the SDN reverberator and that predicted by RT equations suggests that the proposed LSE and MLE methods may be useful tools in calibrating either SDN reverberators or an ISM.

In order to adapt the LSE and MLE methods to calibrating an SDN or an ISM, first consider that the SDN’s RT tends to undershoot the Sabine equation prediction or overshoot the Eyring equation prediction [25]. As a result, directly calibrating the absorption coefficients to either equation will result in deviations between the desired RT and the SDN’s output. However, adjusting the target RT used by the LSE or MLE method to compensate for this mismatch provides a straightforward SDN calibration procedure.

To begin with, assume that one has a measured RT as well as a set of initial absorption coefficient values $\tilde{\mathit{\alpha }}$ and their uncertainties $\tilde{\mathit{\sigma }}$ if using the MLE method. First, calibrate $\tilde{\mathit{\alpha }}$ through equation (10) or (25) by setting the initial target RT, $T_{60,\text{T}}^{(1)}$, to the measured RT. This step produces an initial set of calibrated coefficients ${\widehat{\mathit{\alpha }}}^{(1)}$. Applying these calibrated absorption coefficients to the SDN will produce an output RT which is close to, but not exactly equal to, the measured RT. Let $T_{60,\text{SDN}}^{(1)}$ be this RT derived from the first SDN run. Second, using the percent mismatch between the SDN’s output and measured RT, adjust the target RT to $T_{60,\text{T}}^{(2)}$ to compensate for the deviations. For example, if the SDN overshoots the target RT, then $T_{60,\text{T}}^{(2)}$ will be lower than $T_{60,\text{T}}^{(1)}$ to counter this effect. This new target RT produces the next iteration of calibrated absorption coefficients ${\widehat{\mathit{\alpha }}}^{(2)}$ and their associated SDN RT output $T_{60,\text{SDN}}^{(2)}$. Finally, using the secant method (finite-difference approximation to Newton’s method), continue to iteratively adjust the target RT and calibrated coefficients until obtaining convergence.

Equations (38)–(40) describe the numerical implementation of this procedure while Figure 6 illustrates the procedure. These equations use the reverberation rate R₆₀ = 1/T₆₀ because absorption is proportional to the reverberation rate but inversely proportional to the RT (see Eqs. (5) and (11)). Because the secant method extrapolates values based on the slope of the underlying function, using the reverberation rate led to a faster convergence compared to the RT when applying this method to calibrating GA models (see Sect. 4.3):

Fig. 6 SDN calibration flowchart using the LSE and MLE calibration methods.

$$\begin{array}{cc}\hfill R_{60,\text{T}}^{(1)}& =R_{60}\hfill \end{array}$$(38)

$$\begin{array}{cc}{R}_{60,\text{ T}}^{(2)}& =R_{60}(1-\frac{R_{60,\text{ SDN}}^{(1)}-R_{60}}{R_{60}})\end{array}$$(39)

$$\begin{array}{cc}{R}_{60,\text{ T}}^{(2)}& =R_{60}(1-\frac{R_{60,\text{ SDN}}^{(1)}-R_{60}}{R_{60}})\\ & \times \left(\frac{R_{60,\text{ T}}^{(n-1)}-R_{60,\text{ T}}^{(n-2)}}{R_{60,\text{ SDN}}^{(n-1)}-R_{60,\text{ SDN}}^{(n-2)}}\right),n>3.\end{array}$$(40)

It is interesting to consider the algorithm’s convergence for varying levels of average absorption S^Tα/S. Figure 7 compares the desired T₆₀ time, the initial $T_{60,\text{SDN}}^{(1)}$ derived from the SDN output, and the SDN output after using $T_{60}^{(1)}$ as the target RT. The gray region indicates a 1 JND value of ±5%. Both calibration steps employed the LSE optimization applied to the Sabine equation. The room geometry was a 5 × 7 × 3 m shoebox room with constant wall absorption coefficient and source and receiver locations at x_s = [4.5, 3.5, 2] m and x_r = [2, 2, 1.5] m, the same as used in Section 4 of [27]. Similar to results shown in [25] [27]], the SDN’s RT is too high for low absorption coefficient values while too low for high absorption coefficient values, only achieving values within 1 JND for values ranging from roughly α = 0.2 to α = 0.3. However, after a single correction step (Eq. (39)), the RT of the SDN’s output fell to within almost 1 JND for all averaged absorption coefficient values below α = 0.6. After the second iteration, the RT of the SDN’s output fell to within 1 JND for all absoprtion values. Because the SDN output is nearly identical to those of ISM for a rectangular room [25], one may likewise assume the approach will generalize well to the ISM. Consequently, although the proposed calibration methods directly apply to the Sabine and Eyring equations, they may be exploited to efficiently calibrate other numerical models and methods.

Figure 7 Target and simulated T 20 after the first and second calibration step using a LSE calibration of the Sabine equation.

4.3 Case studies in geometrical acoustics model calibration

While RT equations provide a simplistic yet reasonably realistic representation of room acoustic characteristics [11], most modern methods in room acoustic analysis employ either wave-based or GA models. Because these models are significantly more complicated than RT equations, previous works have employed complicated optimization algorithms to their calibration [9, 10]. However, the recent work of Prawda et al. [11] suggests that RT equations are often more reliable than assumed, able to achieve within 10% (2 JND) of measured RT when absorption coefficients are properly known. Additionally, both Christensen et al. and Pilch noted that a better initial guess led to better calibration results when using their optimization algorithms [9, 10]. Consequently, it is of interest to see if the proposed RT calibration methods may be beneficial in GA model calibrations. Although a rigorous consideration of this topic falls outside the scope of the present work, the following section presents a few cases studies highlighting the potential benefit of this approach.

4.3.1 Lecture hall

The fist step in calibrating the lecture hall under consideration (Fig. 8a) is to gather absorption coefficient data from available databases and handbooks to create an initial RT estimate based on the Sabine equation, including air absorption. This estimated RT appears as the blue dotted line in Figure 8b. Using best estimates of material properties from [18–20] gave an initial RT estimate within roughly 2–3 JND except for the 125 Hz band, whose value exceeded 5 JND. These initial coefficients serve as the μ values for the MLE estimation.

Second, comparing absorption coefficient ranges across databases and also visually assessing the similarity of materials present in the room compared to the available data led to material uncertainty σ assignment. The present work applied a three-tiered scale. First, materials with lower uncertainty, such as the linoleum floor or painted gypsum board walls, received a value of σ_n = 0.01. Second, materials with moderate uncertainty, such as the wood paneling and metal door assembly, received a value of σ_n = 0.05. Finally, materials with the highest uncertainty, such as the unique wooden desks and chairs and perforated ceiling panels, received a value of σ_n = 0.1.

Figure 8 (a) GA model of the lecture hall with 2 source and 3 receiver positions. Measured and first iteration of the GA-simulated (b) T20, (c) EDT, (d) C50, and (e) C80 averaged over the 6 source receiver pairs.

Third, the RT equation calibrated coefficients followed by using the measured T₂₀ averaged across two source and three receiver positions (six total measurements, see Fig. 8a) as the target RT for the MLE calibration. These RT-calibrated absorption coefficient values are then applied to a GA acoustics model in CATT-Acoustic v9 [28]. While the RT equation calibration provides absorption coefficient values, it does not include scattering data. Initial estimates for the scattering coefficients used the characteristic depth of each surface; for further details see [3].

The simulated results averaged across the six source-receiver pairs are shown in Figures 8b–8e. The results include a ±5% 1 JND value for the EDT and T₂₀ and a ±1 dB value for C₅₀ and C₈₀. The agreement is promising considering that the results represent the first iteration of a GA calibration. The percent deviation between the measured and simulated T₂₀ was within 1 JND for all bands, with a frequency averaged deviation of just 1.1%. The RT equation calibrated result slightly underestimated the EDT in all bands, although the frequency averaged deviation was just 4.7% <  1 JND. Both C₅₀ and C₈₀ also showed reasonable agreement, with only a few bands requiring further adjustment. These adjustments could be carried out, for example, by manually altering scattering coefficients and considering local variations in measured outputs as in [3] or using genetic algorithms as an optimizer as in [9, 10]].

4.3.2 Palais des Papes

The lecture hall used in the previous section represents a simple room with moderate absorption and whose sound field may be reasonably expected to behave diffusely due to the presence of many scattering objects. As a result, it is with little surprise that the Sabine-equation calibrated results agreed well with the GA software output without any further adaptations or adjustments. However, some rooms do not behave diffusely or have highly absorptive surfaces, which violate the Sabine equation’s underlying assumptions. For example, the Palais des Papes (Fig. 9a, see [29] for more details on the model and measurements) is a room with low absorption and which has rigid, parallel walls that lead to strong specular reflections. When applying the same calibration approach as in the lecture hall to this space, the results, which appear as the dash-dot red curves in Figure 9b–9e, only achieved good agreement for the 500 Hz octave band and above.

Figure 9 (a) GA model of the Palais des Papes with 4 source and 16 receiver positions. Measured and first iteration of the GA-simulated (b) T20, (c) EDT, (d) C50, and (e) C80 averaged over the 64 source-receiver pairs.

For lower bands, the Sabine equation underestimated the RT produced by the GA model.

Nonetheless, similar to the SDN calibration procedure, this mismatch may be exploited to iteratively adjust the absorption coefficients. The solid red curves in Figures 9b–9e show the results after the first iteration (Eq. (38)). In this case, all bands fell within 1 JND for T ₂₀. The frequency-averaged deviation was less than 0.1 JND for this metric. Additionally, the first correction step brought C ₅₀ and C ₈₀ to less than 1 JND for the 250 Hz band and above; the 125 Hz band fell to within less than 2 JND for these two metrics.

While T ₂₀, C ₅₀, and C ₈₀ fell within less than 1 JND for 250 Hz and above, the EDT was beyond 1 JND for 2 kHz and 4 kHz bands. Before the correction step, the EDT was only 1.3 JND and 2.1 JND too high for these bands, respectively. After the correction step, which decreased the total absorption in these bands to better match the measured T ₂₀, the simulated EDT rose to 1.6 JND and 3.7 JND. However, before the correction step, thesimulated T ₂₀ was too low by 0.3 JND and 1.3 JND. Consequently, one could comprise between these two metrics by using the original calibrated absorption coefficients for 2 kHz and 4 kHz bands and the corrected calibrated coefficients for all other bands.

4.3.3 Athenee Theater

The Athenee Theater [30, 31]] appearing in Figure 10 a represents an even more challenging room. The absorption and scattering properties of the stage and audience areas vary significantly, while narrow barriers between stalls lead to further deviations when using simple RT equations. After calibrating the initial coefficients of the thirty materials to the measured T₂₀, the simulation output is over 5 JND in all bands, with a frequency-averaged error of 8.2 JND (the 125 Hz band is excluded due to low measurement SNR). Even the first correction step (Eq. (39)) only brings in the simulated T₂₀ to an average 3.7 JND.

However, after applying the second correction step (Eq. (40)) the simulated T₂₀ fell to within 1 JND for all bands, with a frequency-averaged deviation of only 0.5 JND. Many of the other metrics had reasonable agreement between measured and simulated results by this correction step. For example, the simulated EDT fell within 1 JND for the 500 Hz, 1 kHz, and 2 kHz bands and within 2 JND for the 250 Hz and 4 kHz bands. Both C₅₀ and C₈₀ were within 1 JND for all bands except for the 1 kHz band. This result suggests that for the sake of model calibration, RT and clarity room acoustic metrics are non-orthogonal: obtaining agreement for one metric may result in reasonable agreement for the others.

Figure 10 (a) GA model of the Athenee Theater with 3 source and 43 receiver positions. Measured and first iteration of the GA-simulated (b) T20, (c) EDT, (d) C50, and (e) C80 averaged over the 129 source-receiver pairs.

5 Discussion

The results of this work highlight several important aspects of room acoustic model calibration. First, calibration requires not only achieving agreement between measured and simulated output but also ensuring that thecalibrated absorption coefficients lie closer to their ground-truth values than the initial estimate. This work proposed the calibration score to quantify this aim. Both the LSE and MLE calibration methods gave positive CS values from the statistical simulations. Consequently, these methods not only can match measured and estimated RT but also improve upon initial estimates of absorption coefficient values. This result is a distinct advantage compared to previous methods, which may produce incorrect absorption coefficient values despite achieving agreement with measured RT.

Another interesting aspect of calibration highlighted by the two methods is that the success of the calibration largely depends on the amount of initial data available. While previous works applying genetic-algorithm-based optimization methods noted that better initial estimates gave better results, the closed-form solutions given by the LSE and MLE illustrate how a better initial estimate and knowledge of uncertainty levels for the MLE case give improved results. Consequently, the calibration methods do not replace the need for careful measurement of material absorptive properties but rather emphasize the need for more reliable data.

Lastly, the presented case study suggests that simpler RT calibration methods may be beneficial tools in room acoustic model calibration, including for the ISM, SDN reverberators, or for GA models. For example, both [9] and [10] reported that their calibration methods converged better with a good initial guess. The proposed RT calibration could provide those values. Calibration of GA models would then follow as (1) collection of material properties and quantifying their uncertainty, (2) GA model pre-calibration using the proposed LSE or MLE methods, and (3) fine-tuning the GA model using a alternative optimization methods or manual tuning with regards to additional parameters beyond T ₆₀ [32].

6 Conclusions

This work proposed two methods for calibrating the absorption coefficients in RT equations. The first uses a constrained least-squares estimate, whereas the second employs a constrained maximum-likelihood estimate. The maximum-likelihood estimate allows users to include material uncertainty into the calibration process to more accurately predict property values. Stochastic simulations demonstrated that the proposed methods not only achieve the desired RT but also improve upon the initial estimates of the absorption coefficient values. For a simple shoebox room with six materials and high uncertainty conditions, the LSE method and MLE method improved upon the baseline method’s calibration score of 6.5% to 8.7% and 13.2%, respectively. For the lecture hall condition with 23 materials, the LSE method and MLE method improved upon the baseline method’s calibration score of −0.3% to 2.2% and 3.3%, respectively. Consequently, the proposed methods decrease the likelihood of providing unrealistic or unphysical values (evidenced by negative calibration scores) compared to works which only consider minimizing the differences between measured and simulated room acoustic parameters. In addition, the proposed methods result in simple closed-form expressions that significantly reduce the time needed for model calibration.

Application of the methods to iterative calibrations of an SDN reverberator and GA models highlighted the portability of the approach to more complicated but realistic room acoustic modeling applications. For the SDN calibration, the LSE method brought the simulated RT to within 1 JND of the target RT after only two iterations. For GA models, the MLE method applied to a lecture hall achieved 1 JND between the measured and simulated T ₂₀ across six octave bands with no need for further iterations. More complicated geometries such as the Palais des Papes and the Athenee Theater required one and two iterations, respectively, to achieve 1 JND between measured and simulated T ₂₀ values averaged across source-receiver pairs. Consequently, the proposed methods brought the simulated T ₂₀ to within 1 JND of the measured values across varying frequencies and room geometries. This result contrasts with machine-learning-based methods, which, for certain frequencies and room shapes, had deviations exceeding 1 JND. In addition, the proposed analytic methods used fewer iterations than manual tuning or machine-learning-based approaches, which required between four to eight iterations, respectively. Future work could include studying the application of these methods to improved calibration of GA and wave-based architectural acoustics models.

Cite this article as: Bellows S.D. & Katz B.F.G. 2025. Constrained least-squares and maximum-likelihood calibration of absorption coefficients in reverberation time equations. Acta Acustica, 9, 25. https://doi.org/10.1051/aacus/2025011.

Funding

This work was carried out in part in the context of the SONICOM project (https://www.sonicom.eu) that has received funding from the European Union’s Horizon 2020 research and innovation program under grant agreement No. 101017743.

Conflicts of interest

The authors declare no conflict of interests.

Data availability statement

An online repository at [33] contains freely available software implementations and demonstrations of the proposed methods. Additional data are available from the authors upon reasonable request.

References

All Tables

All Figures

	Figure 1 Two-dimensional illustration of the least-squares calibration.
In the text

	Figure 2 Three truncated normal distributions with different parameters. Higher values of σ indicate lower certainty in the absorption coefficient value.
In the text

	Figure 3 Maximum-likelihood estimation for calibrating absorption coefficient values for a two-dimensional problem.
In the text

	Figure 4 (a) PDFs for the six material properties for a single room realization. (b) Statistical properties of the EMSE over the 104 room realizations incorporating 104 trial runs each. The boxes indicate the 25th and 75th percentiles while the whiskers indicate the extreme values for each case.
In the text

	Figure 5 (a) Shoebox room and (b) lecture hall used in the stochastic simulations.
In the text

	Fig. 6 SDN calibration flowchart using the LSE and MLE calibration methods.
In the text

	Figure 7 Target and simulated T 20 after the first and second calibration step using a LSE calibration of the Sabine equation.
In the text

	Figure 8 (a) GA model of the lecture hall with 2 source and 3 receiver positions. Measured and first iteration of the GA-simulated (b) T20, (c) EDT, (d) C50, and (e) C80 averaged over the 6 source receiver pairs.
In the text

	Figure 9 (a) GA model of the Palais des Papes with 4 source and 16 receiver positions. Measured and first iteration of the GA-simulated (b) T20, (c) EDT, (d) C50, and (e) C80 averaged over the 64 source-receiver pairs.
In the text

	Figure 10 (a) GA model of the Athenee Theater with 3 source and 43 receiver positions. Measured and first iteration of the GA-simulated (b) T20, (c) EDT, (d) C50, and (e) C80 averaged over the 129 source-receiver pairs.
In the text