Spatial characterization of the sound field in a room is a challenging task, as it usually requires a large number of measurement points. This paper presents a probabilistic approach for sound field reconstruction in the modal frequency range for small and medium-sized rooms based on Bayesian inference. A plane wave expansion model is used to decompose the sound field in the examined domain. The posterior distribution for the amplitude of each plane wave is inferred based on a uniform prior distribution with limits based on the maximum sound pressure observed in the measurements. Two different application cases are studied, namely a numerically computed sound field in a non-rectangular two-dimensional (2D) domain and a measured sound field in a horizontal evaluation area of a lightly damped room. The proposed reconstruction method provides an accurate reconstruction for both examined cases. Further, the results of Bayesian inference are compared to the reconstruction with a deterministic compressive sensing framework. The most significant advantage of the Bayesian method over deterministic reconstruction approaches is that it provides a probability distribution of the sound pressure at every reconstruction point, and thus, allows quantifying the uncertainty of the recovered sound field.

Capturing the spatial properties of a sound field in a room is important in various applications, e.g., in sound field analysis,1 source localization,2,3 and room compensation and equalization for sound field control systems.4,5 Typically, an extensive number of measurements is required in order to determine the frequency response experimentally with basic sampling techniques within an extended region of a room.6 In order to minimize the experimental effort, recent research has focused on reducing the number of measurement points.

Several approaches to reconstruct the sound field in a closed domain based on a limited number of measurement points have been published during the last decades. Generally, the reconstruction of the sound field is based on a physical wave propagation model.7 The principal idea is to project the sound field onto an elementary wave basis, often consisting of analytic propagating wave functions and use this basis to predict the sound field throughout the domain.8,9 A set of plane waves, spherical waves, or also spherical or cylindrical harmonics are primarily used as the reconstruction basis.10–14 This also allows for an approximate reconstruction of the sound field at regions where no direct measurements are performed. It should be noted that no information about the geometry, the boundary properties, or the sound source that generates the sound field, is necessary for the reconstruction, since the model just incorporates a propagating wave model.18 

The unknown amplitudes of the wave propagation model need to be determined based on measured sound pressure data at a limited number of observation points. This results in an inverse acoustic problem, which is usually ill posed. Small perturbations, due to measurement noise in the observed data, can lead to large errors and hence, an unstable solution of the inverse problem for the unknown plane wave amplitudes. Standard deterministic inference methods try to overcome this ill-posedness by adding regularization. By this, the inverse problem is reformulated such that it is uniquely solvable.15 

The number of required measurement points can be reduced in the low frequency range by assuming some sort of sparsity, either in the spatial, temporal or modal domain. The sparsity justifies the use of sparse reconstruction approaches such as compressive sensing (CS).16–19 Mignot et al.17 developed a CS algorithm to interpolate room impulse response functions based on sound pressure measurements in a rectangular room at low frequencies. Verburg and Fernandez-Grande16 presented a method based on CS for the reconstruction of the frequency response at arbitrary points based on a plane wave expansion model. CS has been typically applied to rectangular rooms due to the necessary sparsity requirements.16,17 Recently, a data driven deep-learning based approach for sound field reconstruction in rectangular rooms was proposed by Lluíz et al.20 Pham Vu and Lissek introduced a method based on a greedy algorithm and a global curve-fitting technique to recover the modal parameters of a non-rectangular room and to reconstruct the entire sound field at low frequencies.21 Beyond that, numerical reconstruction approaches exist based on the inverse boundary element method (IBEM)22–24 or the inverse finite element method (IFEM).25 These approaches solve the inverse reconstruction problem to provide a point estimate, i.e., one value of the sound pressure or other field quantities at the given reconstruction points. Accordingly, the uncertainty of the reconstruction is not provided. This is appealing in situations where the observed data are not abundant.

This work examines the use of Bayesian inference (BI) for the reconstruction of a sound field in an enclosure, since Bayesian probability theory is well suited to rigorously handle a lack of experimentally observable information.26 The use of BI results in a probability distribution of the sound pressure at each reconstruction point in space. Examination of this probabilistic distribution allows obtaining knowledge about the uncertainty of the estimation and insights into its robustness. In contrast to deterministic methods, the BI approach inherently considers measurement noise in the statistical formulation and thus provides advantages for noisy measurement data.27 

Over the last years, BI has proved to be a viable tool to generally solve ill-posed inverse problems and has been applied to different problems in the field of acoustics.28 Beaton and Xiang29 employed BI to model the sound field in a room as a superposition of decaying sinusoids based on a time-domain Prony model, which has a similar structure to the model used in this work. Xiang et al.30 used BI to estimate the sound energy decay rate for coupled enclosures, where multiple decay slopes exist. Antoni28 employed BI to reconstruct the field radiated by sound sources. Moreover, statistical approaches are frequently adopted in signal processing applications26 and sound source localization.31,32 Beyond that, BI is used for the estimation of the amplitudes of a spherical harmonics sound radiation model for sound field control purposes33 and the estimation of porous material model parameters.34,35

This work introduces the probabilistic BI method for the purpose of spatial sound field reconstruction in enclosures. To the authors' knowledge, Bayesian approaches to reconstruct the sound field in a room have mostly focused on estimating modal parameters in the time-frequency domain. The spatio-spectral reconstruction has not received similar attention. The authors specifically examine the uncertainty of the obtained reconstructions with a focus on the modal frequency range.

The present article is structured as follows: The plane wave expansion model and the BI method to determine the posterior distributions for the plane wave amplitudes are described in Secs. II and III, respectively. In Sec. IV, the reconstruction results for a synthesized non-rectangular 2D domain problem and an experimentally captured sound field in a horizontal measurement area of a cuboid room are presented and discussed. Furthermore, the results of the BI approach are compared to a deterministic method to approximate the inverse acoustic problem. Finally, the conclusions and the benefit of the BI method over conventional methods are summarized in Sec. V.

The sound field inside a source free acoustic domain Ω is described by the homogeneous Helmholtz equation, where the harmonic time dependence for the sound pressure p(rm,t)={p(rm)eiωt} is assumed. The position where the sound pressure is measured is represented by the vector rm=[xm,ym,zm]T and t stands for the time variable. The sound pressure field inside an acoustic domain at a certain frequency f=ω/2π can be approximated by a truncated plane wave expansion8 

p(rm)j=1NAjei(kj·rm)+n,
(1)

where each term in the sum corresponds to one of the N plane waves which compose the sound field and the time dependence eiωt has been omitted.36 The complex-valued coefficient Aj represents the amplitude of the jth plane wave, each propagating in the direction specified by the wave number vector

kj=[kcosϕjsinθj,ksinϕjsinθj,kcosθj]T
(2)

in spherical coordinates. In the three-dimensional case, the wave number vector lies on the surface of a sphere with radius k=2πf/c.8 The polar angles ϕj and the azimuthal angles θj are sampled from a Fibonacci sphere in order to guarantee evenly distributed sampling of the sphere. The speed of sound is represented by c and n stands for noise present in measurements.

Sampling the sound pressure field at a discrete number of M measurement positions, the inverse acoustic problem can be expressed in matrix notation as

p=Ha+n,
(3)

where p=[p1,p2,,pM]T represents the measured sound pressure vector and HM×N is the so called sensing matrix containing the plane wave functions, with Hmj=eikj·rm. The vector a=[A1,A2,,AN]T contains the unknown complex plane wave amplitudes and n=[n1,n2,,nM]T consists of the measurement noise at each observation point. Since it is desirable to reduce the number of measurement points, the problem is typically underdetermined (N > M) for practically relevant cases and the number of unknown amplitudes Aj exceeds the number of measurement positions.15 

Bayesian probability theory is well suited to rigorously handle a lack of information and thus, a Bayesian approach has shown to be advantageous, when experimentally obtainable information is not abundant.26 

Within the modal frequency range, the sound field inside a room can be assumed to be composed by a limited number of waves.36 Thus, only the corresponding coefficients of a will be non-zero. CS is a deterministic approach to approximate the ill-posed inverse problem formulated in Eq. (3) for given sound pressure data, assuming a to be sparse.16–19 The problem of finding the optimal set of plane wave amplitudes ã=[A1̃,A2̃,,AÑ]T can be expressed in terms of a convex optimization problem19 

ã=argminaN||pHa||22subjectto||a||1η.
(4)

The problem in Eq. (4) minimizes the cost function depending on the squared sum of the residual and a 1-norm penalty on the coefficients. The penalty term on the sum of absolute values in a promotes sparse solutions. The parameter η represents an estimation on the upper bound of the noise present in the sensing process.16 Thus, the choice of η will influence the number of non-zero elements of ã. According to Verburg and Fernandez-Grande,16 an appropriate estimate for η is ||n||2, which is employed in this work. Note that the 1-norm represents a relaxation of the 0-norm sparsity penalty, which would lead to a non-convex and non-deterministic polynomial-time (NP) hard optimization problem. By replacing the 0-norm with the convex 1-norm, the problem can be efficiently solved with convex optimization.19 The reconstruction results obtained by CS are compared to the performance of the statistical BI method.

The estimated plane wave amplitudes Aj̃ can be used to reconstruct the sound pressure at a set of positions rR=[xR,yR,zR]T, where no measurements are available. Recalling the plane wave expansion model in Eq. (1), the reconstructed sound pressure p̃ can be formulated as

p̃(rR)j=1NAj̃ei(kj·rR).
(5)

This formulation allows not only the interpolation of the sound pressure field inside the domain, where the measurements are taken, but also the extrapolation at positions outside the measurement domain.

In general, BI is the process of evaluating a probability model to a set of data and estimating the result by probability distributions on the model parameters.27 The centerpiece of the underlying theory is Bayes' theorem for parameter estimation. Applied to the case of sound field reconstruction, the theorem is governed by

π(a|p)=π(p|a)π(a)π(p),
(6)

where π(a) denotes the prior probability density function (PDF) of the plane wave amplitudes a. This quantity includes the physically meaningful information on the parameters which is known a priori. The function π(p) is the evidence and π(p|a) expresses the likelihood of the observed sound pressure data. By combining all these quantities, the posterior probability distribution π(a|p) of the amplitudes is obtained.28 The Bayesian framework can be interpreted as updating one's prior knowledge π(a) about the model parameters once experimental data are incorporated through the likelihood function π(p|a) and obtaining the posterior probability distribution π(a|p) by applying Bayes' theorem.26 Since the evidence in Eq. (6) only serves as a normalization constant in the case of parameter estimation, Bayes' theorem for parameter estimation simplifies to27 

π(a|p)π(p|a)π(a).
(7)

BI is used to determine the posterior distributions of the plane wave amplitudes, given the observed sound pressures at the measurement points and taking the prior knowledge into account.

For complex engineering problems, a closed-form solution for the posterior distribution π(a|p) in Eq. (7) is not available. In order to circumvent this issue, the posterior distribution is approximated by sampling based methods, namely using the Markov Chain Monte Carlo (MCMC) method introduced by Metropolis et al.37 The general idea behind MCMC sampling is to obtain random samples from the unknown posterior by drawing a sequence of samples [a1,a2,,aS] called Markov chain, which proportionally represents the posterior distribution. For a large number of samples S, the approximate distribution converges towards the target posterior distribution. Numerous modifications of MCMC methods are present in the literature, which differ in the transition rule of the Markov chain.27,38

The efficient state-of-the-art gradient based No-U-Turn sampler (NUTS)39 is utilized in this work, which belongs to the class of Hamiltonian Monte-Carlo (HMC) algorithms. It suppresses the random walk behavior, which is present for the Metropolis–Hastings algorithm37 by using a physical analogy to propose the move from the current sampling state to the next state.40 This allows an efficient exploration of the posterior distribution compared to basic random walk sampling algorithms [27]. In the scope of this work, the Python interface PyStan of the open-source probabilistic programming environment Stan41,42 is utilized to perform the sampling. Since Stan does not support complex numbers, the plane wave propagation model is divided into real and imaginary part.

A uniform prior distribution is assigned to the complex amplitude coefficients. It a priori prescribes an equal probability for each amplitude value. This corresponds to an assignment according to the principle of maximum entropy.43 The bounds of the uniform distribution for the real and imaginary parts of the amplitudes are assigned based on the physically measured sound pressure data at the measurement positions as

{Aj},{Aj}U(2pmax,+2pmax),
(8)

where pmax represents the maximum sound pressure magnitude observed at the considered frequency at the point with maximum sound pressure among all measurement points. Preliminary simulations showed that this is a reasonable choice for the number of measurement points and in the considered frequency range.

The likelihood function is analogous to the prior distribution rigorously assigned according to the principle of maximum entropy, which results in a Gaussian likelihood distribution.26,43 It can be formulated as

π(p|a)=(2π)m/2det(Σ)m/2exp(12(pp̃)TΣ1(pp̃)),
(9)

where m stands for the number of unknown parameters, Σ=diag(σ2) is the diagonal covariance matrix and the vector p̃ contains the reconstructed sound pressure values for a specific set of plane wave amplitudes. The assigned likelihood distribution ensures no bias for which there is no prior evidence.26 It is worth mentioning that the Gaussian likelihood assignment is fundamentally different from assuming the statistics of the residual errors being Gaussian.26 In the case of sound field reconstruction, the likelihood represents the degree of belief that the measured sound pressure data would have been generated from a given set of complex-valued amplitudes.27 

Sampling of the posterior is performed independently and frequency-wise using the NUTS sampling algorithm. Six independent Markov chains with different initial states and 1000 samples for each chain per frequency step are performed in both presented analyses. Additionally, 1000 tuning samples are initially drawn in order to tune the NUTS algorithm's hyperparameters. These tuning samples are discarded as burn-in. The Gelman–Rubin convergence criterion27 is evaluated to qualitatively assess convergence of the sampling algorithm. The criterion compares the variance between the six independent Markov chains with different initial states to the variance within each individual chain. Hence, it indicates, if all chains converge towards the same target distribution, regardless of the initial state. The posterior distributions for the complex-valued amplitudes are used to reconstruct the sound pressure in the domain based on the plane wave superposition model formulated in Eq. (5).

The reconstructed real and imaginary parts of the sound pressure are assumed to be independent and normally distributed with mean values μ{p} and μ{p} and identical standard deviations.15 According to Aja-Fernández and Vegas-Sánchez-Ferrero44 the transformed sound pressure magnitude follows a Rice distribution

|p|R(ν,σ),whereν=μ{p}2+μ{p}2.
(10)

The BI reconstruction method is applied to both a numerical and an experimental study. At first, a generic 2D problem concerning a non-rectangular shaped domain is examined, where the reference sound field is computed numerically. Second, the reconstruction method is applied to an experimentally measured sound field in a horizontal evaluation area of a room. The inference is carried out using the proposed BI method as well as the deterministic CS approach for both cases in the modal frequency range. The CS optimization problem formulated in Eq. (4) was solved making use of an interior point method, with the CVXPY https://www.cvxpy.org/ (Nov. 20, 2021) toolbox for convex optimization.47 

Two error measures are introduced to assess the performance of the BI method and to compare the results with the deterministic method, namely the normalized root mean square error (NRMSE) in percent

NRMSE=1Kk=1K||pkp̃k||22||pk||22×100%
(11)

and the modal assurance criterion (MAC)45 

MAC=||pHp̃||22(pHp)(p̃Hp̃).
(12)

The vectors p,p̃K represent the reference and the reconstructed sound pressure fields, K stands for the number of points in the whole domain, and the superscript H indicates the Hermitian. The NRMSE is an indicator for the average reconstruction error on the whole domain in percent for every frequency under consideration. The MAC is a measure for the consistent correspondence between two vectors. Thus, it is a reliable measure to assess the spatial similarity between the reference sound pressure field and the reconstruction. It takes values between 0 and 1, where 0 indicates no similarity and 1 perfect similarity.

The geometry of the studied non-rectangular 2D domain is depicted in Fig. 1. The upper boundary of the domain is tilted such that the height is decreasing from Ly(x=0)=2m to Ly(x=Lx)=1m over the total length of Lx=3.4m in x-direction. The reference sound field is computed numerically by using the commercial Finite Element Method (FEM) software COMSOL Multiphysics https://www.comsol.com/ (Nov. 20, 2021). A normal particle velocity of v=103ms1 is applied at the boundaries at x=0m and y=0m, corresponding to a wall-piston excitation of the sound field. Hence, all room modes in the domain are excited. The fluid density ρ=1.2kgm3 and the speed of sound c = 340 ms1 are used. All boundaries are assumed to be sound hard. The domain is discretized with 11 124 triangular elements of quadratic order. Ten data points are arbitrarily chosen, whose locations are depicted as black crosses in Fig. 2. Gaussian noise with a signal-to-noise ratio (SNR) of 30 dB is added to the numerically calculated sound pressure data in order to emulate measurement noise artificially. Since this numerical problem is two-dimensional, the azimuthal angles θj in the plane wave expansion model formulated in Eq. (1) are set to θj=π/2, which simplifies the 3D model to 2D. Thus, the plane waves are only incoming from the 2D reconstruction plane. The polar angle ϕj is subdivided into 16 uniformly spaced values, which represent the directions of the incoming plane waves with unknown amplitudes Aj. These incidence directions of the plane waves are depicted as arrows in Fig. 1. This leads to a problem with 32 parameters, 16 real and imaginary parts of the amplitudes, which need to be inferred, given the sound pressure at the 10 data points for each frequency step. The reconstruction is applied in a frequency range between 30 Hz and 150 Hz, where the study is successively executed frequency-wise with a step size of 10 Hz.

FIG. 1.

Geometry of the generic 2D reconstruction problem with arrows representing the angles of incidence of the plane waves.

FIG. 1.

Geometry of the generic 2D reconstruction problem with arrows representing the angles of incidence of the plane waves.

Close modal
FIG. 2.

(Color online) Comparison between the FEM reference sound pressure field (top), the BI reconstruction with a uniform prior (middle), and the CS reconstruction (bottom) for the generic 2D problem at 50 Hz (left) and 150 Hz (right). For the BI reconstruction, the mean of the sound pressure magnitude posterior distribution is depicted in color and the standard deviation is superimposed as a gray scale contour line plot. The black crosses indicate the chosen data points for the reconstruction.

FIG. 2.

(Color online) Comparison between the FEM reference sound pressure field (top), the BI reconstruction with a uniform prior (middle), and the CS reconstruction (bottom) for the generic 2D problem at 50 Hz (left) and 150 Hz (right). For the BI reconstruction, the mean of the sound pressure magnitude posterior distribution is depicted in color and the standard deviation is superimposed as a gray scale contour line plot. The black crosses indicate the chosen data points for the reconstruction.

Close modal

Figure 2 shows the comparison between the reference sound pressure field, the results of, and CS at the frequencies of 50 Hz and 150 Hz. The reference sound field is depicted in the top, the BI reconstruction in the middle and the CS reconstruction in the bottom row. The frequency of 50 Hz is shown on the left and 150 Hz on the right hand side, respectively. The BI reconstruction in the middle row shows the mean of the fitted Rice distribution for the reconstructed sound pressure magnitude in color scale and the standard deviation as a superimposed gray scale contour line plot. This allows analyzing the standard deviation of the determined posterior distribution in the reconstruction domain and identify regions of high uncertainty of the reconstruction. It is observed that the regions of high standard deviation of the reconstruction coincide with regions where no data points are available. The standard deviation of the posterior distribution generally takes higher values for 150 Hz than for 50 Hz. This originates in the fact that the complexity of the sound field increases with the modal density and wavelengths are shorter compared to the microphone spacing. By comparing the results of the Bayesian approach and CS, it becomes evident that both methods are capable of accurately reconstructing the sound field in the investigated frequency range. However, the deterministic method just leads to a single sound pressure at every reconstruction point. In contrast, the BI method yields a posterior distribution for the sound pressure at every point and thus, additional information about the uncertainty of the reconstruction is obtained.

Since the reconstruction depends on the superimposed noise, the reconstruction is performed Navg = 50 times in order to average out the stochastic nature of the noise and to quantitatively assess the accuracy of the reconstruction methods. Figure 3(a) shows the mean and the standard error (SE) as the shaded area for the NRMSE and Fig. 3(b) shows the MAC for the 2D reconstruction problem. The SE is evaluated as

SE=σavgNavg,
(13)

where σavg represents the standard deviation for the NRMSE or the MAC resulting from the 50 evaluations.

FIG. 3.

(Color online) Mean and SE of (a) the NRMSE and (b) the MAC for the BI and CS reconstruction for the 2D generic problem in the studied frequency range between 30 Hz and 150 Hz for an ensemble of 50 realizations. Note that the MAC is only depicted over a small range of [0.98 1].

FIG. 3.

(Color online) Mean and SE of (a) the NRMSE and (b) the MAC for the BI and CS reconstruction for the 2D generic problem in the studied frequency range between 30 Hz and 150 Hz for an ensemble of 50 realizations. Note that the MAC is only depicted over a small range of [0.98 1].

Close modal

In general, both methods perform well for this specific problem, with mean NRMSE values below 6%. Also, the MAC mean values remain above 0.99 for both methods in the considered frequency range. Yet, with a maximum NRMSE mean of 5.4% at 150 Hz, the BI approach generally results in slightly smaller reconstruction errors compared to CS, where the maximum NRMSE mean is 5.8% at 85 Hz. Further, the SE of the NRMSE and the MAC is lower for BI, which indicates a more consistent reconstruction for noisy observation data. Both methods perform well for this investigated problem, with a marginal advantage of the Bayesian approach.

Experimental testing was performed at the Acoustic Technology Labs at the Technical University of Denmark and reference is made to Hahmann et al.46 regarding the dataset. An approximately cuboid room of dimensions Lx=4.41m,Ly=3.31m,Lz=2.97m, a reverberation time of around T60=2.7s, and a Schroeder frequency of fS=500Hz is investigated. A sketch of the room under investigation is depicted in Fig. 4. The walls, the floor, and the ceiling are made out of concrete with minor irregularities at the boundary, such as the door. No scattering elements are placed inside the room. A Universal Robots (Odense, Denmark) UR5 robotic arm is used to precisely position a microphone at 289 measurement positions in a horizontal measurement area parallel to the xy-plane, which is highlighted in blue in Fig. 4. By this, the frequency response on the defined plane and the corresponding reference sound field can be measured with a single microphone and high accuracy in position. The measurement area is located at x=[2.6m,3.4m],y=[0.85m,1.65m], and a height of z=1.58m. Since the wavelengths in the examined frequency range are considerably larger than the cross section of the robotic arm, it just has a negligible effect on the sound pressure measurements. The room is excited using 10 s logarithmic sweeps by a Dynaudio (Skanderborg, Denmark) BM6 loudspeaker positioned close to a room corner, in order to effectively excite all modes in the considered frequency range. A Brüel & Kjær (Nærum, Denmark) 0.5 in. free field microphone is used for measuring the sound pressure. Figure 5 shows the experimental setup in the room under investigation. The proposed reconstruction method presented in Sec. II is applied in order to reconstruct the sound field in the measurement plane of the room. Fifty plane waves are assumed for the wave propagation model, cf. Eq. (1). Twenty randomly chosen observation positions are considered for the reconstruction. They are marked with black crosses in Fig. 6. The reconstruction is executed frequency-wise in a frequency range between 50 Hz and 500 Hz.

FIG. 4.

(Color online) Geometry and dimensions of the cuboid room under investigation with the highlighted horizontal measurement area.

FIG. 4.

(Color online) Geometry and dimensions of the cuboid room under investigation with the highlighted horizontal measurement area.

Close modal
FIG. 5.

(Color online) Experimental setup to measure the reference sound pressure field inside the room with a loudspeaker placed close to a room corner to excite the room and a robotic arm, which precisely positions the microphone in the defined measurement area.

FIG. 5.

(Color online) Experimental setup to measure the reference sound pressure field inside the room with a loudspeaker placed close to a room corner to excite the room and a robotic arm, which precisely positions the microphone in the defined measurement area.

Close modal
FIG. 6.

(Color online) Comparison between the measured sound pressure field in the measurement area of the room (top), the BI reconstruction with a uniform prior (middle), and the CS reconstruction (bottom) at 150 Hz (left) and 400 Hz (right). For the BI reconstruction, the mean of the sound pressure magnitude posterior distribution is depicted in color and the standard deviation is superimposed as a gray scale contour line plot. The black crosses indicate the chosen data points, and the black dot highlights the random evaluation point for the posterior distributions shown in Fig. 7.

FIG. 6.

(Color online) Comparison between the measured sound pressure field in the measurement area of the room (top), the BI reconstruction with a uniform prior (middle), and the CS reconstruction (bottom) at 150 Hz (left) and 400 Hz (right). For the BI reconstruction, the mean of the sound pressure magnitude posterior distribution is depicted in color and the standard deviation is superimposed as a gray scale contour line plot. The black crosses indicate the chosen data points, and the black dot highlights the random evaluation point for the posterior distributions shown in Fig. 7.

Close modal

Figure 6 shows the comparison between the measured reference sound pressure field in the measurement area of the room and the reconstructed sound field using BI and CS for the exemplary frequencies of 150 Hz (left) and 400 Hz (right). The measured reference sound field is depicted in the top, the results of the BI reconstruction in the middle, and the CS reconstruction in the bottom row. The mean of the sound pressure magnitude is plotted in color and the superimposed gray scale contour line plot shows the standard deviation in the reconstruction area. For 150 Hz, both methods lead to an accurate reconstruction of the measured sound field. In the higher frequency range at 400 Hz, the BI approach results in a slightly more precise reconstruction than CS. This stems from the restriction of the uniform prior distribution of the plane wave amplitudes based on the maximum sound pressures observed at the measurement positions. The deterministic method especially leads to slight reconstruction inaccuracies at the boundaries of the domain, where no data points are located. This indicates a little benefit of the presented statistical method for the extrapolation of the sound field outside the measurement domain. At the boundary regions, where CS regularization is slightly imprecise, the standard deviation of the posterior distribution of BI is the highest. The Bayesian approach provides information about the uncertainty of the reconstruction through the posterior distribution and thus, enables to quantify the uncertainty of the reconstruction.

Figure 7 shows the histograms of the sound pressure magnitude posterior distribution and the fitted Rice distributions (black lines) evaluated at the random evaluation point marked with a black dot in Fig. 6 for 150 Hz on the left and for 400 Hz on the right side. The posterior distribution of the reconstructed sound pressure magnitude is wider for the higher frequency of 400 Hz, which is reflected in a higher standard deviation of the distribution. This indicates an increasing uncertainty of the reconstruction for higher frequencies. The slightly more accurate reconstruction performance of the BI approach is also observed when comparing the NRMSE and the MAC for the measured sound pressure field in the evaluation area.

FIG. 7.

(Color online) Histograms of the reconstructed posterior distributions of the sound pressure magnitude at the random evaluation point in the measurement area of the room, which is marked with a black dot in Fig. 6 at a frequency of (a) 150 Hz, and (b) 400 Hz. The solid black lines indicate the fitted Rice distributions.

FIG. 7.

(Color online) Histograms of the reconstructed posterior distributions of the sound pressure magnitude at the random evaluation point in the measurement area of the room, which is marked with a black dot in Fig. 6 at a frequency of (a) 150 Hz, and (b) 400 Hz. The solid black lines indicate the fitted Rice distributions.

Close modal

Figure 8(a) shows the NRMSE and Fig. 8(b) shows the MAC in the frequency range from 50 Hz up to 500 Hz. In general, the reconstruction accuracy decreases in the higher frequency range, which is reflected in higher NRMSE and lower MAC values. The modal density and the complexity of the sound field are increasing with frequency and the considered number of microphone positions are no longer sufficient for the reconstruction. In the frequency range below approximately 300 Hz, both methods lead to an accurate reconstruction of the measured sound pressure field. Yet, in the higher frequency range, the maximum NRMSE of BI amounts to 34.9% at 445 Hz whereas the highest NRMSE value of CS is 58.7% at 445 Hz. Further, the evaluation of the MAC confirms a slightly more accurate reconstruction and a lower variability of the Bayesian approach. At the power frequency of 50 Hz, the statistical method leads to a slightly more precise reconstruction than deterministic CS.

FIG. 8.

(Color online) (a) NRMSE, and (b) MAC for the BI and CS reconstruction for the measured sound pressure field in the measurement area of the room between 50 Hz and 500 Hz.

FIG. 8.

(Color online) (a) NRMSE, and (b) MAC for the BI and CS reconstruction for the measured sound pressure field in the measurement area of the room between 50 Hz and 500 Hz.

Close modal

Since the Bayesian approach allows incorporating prior information and it restricts the assigned uniform prior distribution based on the measured sound pressure values, it enables a more accurate reconstruction of the investigated sound field. The improvement is especially significant for the reconstruction of the sound field outside the measurement domain (extrapolation), and in the higher frequency range. BI primarily advances conventional reconstruction methods by providing a posterior distribution of the sound field and hence, allows quantifying the uncertainty of the reconstruction.

Yet, this additional uncertainty information has its price. The present implementations require tBI = 184 s (BI via NUTS) and tCS = 6.5 s (CS via interior point method) on average to evaluate each frequency step. Since the computational effort of BI is significantly higher than deterministic methods, it appears to be particularly suitable for cases where uncertainty quantification of the reconstruction is desired.

A Bayesian inference (BI) approach to reconstruct the sound field in an enclosure has been presented. The study considers the reconstruction in the modal frequency range for two different application cases. Namely, a numerically computed sound field in a non-rectangular 2D domain and a measured sound field in a horizontal evaluation area of a lightly damped room are analyzed. A plane wave expansion model is used to expand the sound pressure field in the reconstruction domain and the posterior distribution of the plane wave amplitudes are determined by Markov chain Monte Carlo sampling. Yet, sampling of the posterior distribution is computationally demanding, which leads to a higher computation time than deterministic reconstruction methods. In both examined cases, the Bayesian approach leads to slightly more accurate and robust reconstructions than the compressive sensing reference, especially in extrapolation regions. The most significant advantage of BI is that it allows quantifying the uncertainty at each reconstruction point via the posterior distribution. This uncertainty information can provide important insights. For example, future experiments can be designed in such a way that additional measurement points can be placed in regions of high uncertainty of the reconstruction.

The authors would like to thank Samuel A. Verburg Riezu for the valuable discussions and for the support regarding the measurements at the facilities of the Acoustic Technology Group at the Technical University of Denmark (DTU). The research is partly supported by Villum Foundation through Villum Young Investigator Grant No. 19179 for the project “Large-scale Acoustic Holography.”

1.
I. B.
Witew
,
M.
Vorländer
, and
N.
Xiang
, “
Sampling the sound field in auditoria using large natural-scale array measurements
,”
J. Acoust. Soc. Am.
141
(
3
),
EL300
306
(
2017
).
2.
T.
Nowakowski
,
J.
de Rosny
, and
L.
Daudet
, “
Robust source localization from wavefield separation including prior information
,”
J. Acoust. Soc. Am.
141
(
4
),
2375
2386
(
2017
).
3.
A.
Xenaki
,
E.
Fernandez-Grande
, and
P.
Gerstoft
, “
Block-sparse beamforming for spatially extended sources in a bayesian formulation
,”
J. Acoust. Soc. Am.
140
(
3
),
1828
1838
(
2016
).
4.
F. M.
Heuchel
,
E.
Fernandez-Grande
,
F. T.
Agerkvist
, and
E.
Shabalina
, “
Active room compensation for sound reinforcement using sound field separation techniques
,”
J. Acoust. Soc. Am.
143
(
3
),
1346
1354
(
2018
).
5.
S.
Spors
,
H.
Buchner
,
R.
Rabenstein
, and
W.
Herbordt
, “
Active listening room compensation for massive multichannel sound reproduction systems using wave-domain adaptive filtering
,”
J. Acoust. Soc. Am.
122
(
1
),
354
369
(
2007
).
6.
T.
Ajdler
,
L.
Sbaiz
, and
M.
Vetterli
, “
The plenacoustic function and its sampling
,”
IEEE Trans. Signal Process.
54
(
10
),
3790
3804
(
2006
).
7.
Y.
Haneda
,
Y.
Kaneda
, and
N.
Kitawaki
, “
Common-acoustical-pole and residue model and its application to spatial interpolation and extrapolation of a room transfer function
,”
IEEE Trans. Speech Audio Process.
7
(
6
),
709
717
(
1999
).
8.
E. G.
Williams
,
Fourier Acoustics: Sound Radiation and Nearfield Acoustical Holography
(
Academic Press
,
London, UK, San Diego, California
,
1999
).
9.
E. G.
Williams
, “
Regularization methods for near-field acoustical holography
,”
J. Acoust. Soc. Am.
110
(
4
),
1976
1988
(
2001
).
10.
B.
Rafaely
, “
Plane-wave decomposition of the sound field on a sphere by spherical convolution
,”
J. Acoust. Soc. Am.
116
(
4
),
2149
2157
(
2004
).
11.
J.
Hald
, “
Basic theory and properties of statistically optimized near-field acoustical holography
,”
J. Acoust. Soc. Am.
125
(
4
),
2105
2120
(
2009
).
12.
M.
Nolan
and
E.
Fernandez-Grande
, “
Volumetric reconstruction of acoustic energy flows in a reverberation room
,”
J. Acoust. Soc. Am.
145
(
3
),
EL203
EL208
(
2019
).
13.
E.
Fernandez-Grande
, “
Sound field reconstruction using a spherical microphone array
,”
J. Acoust. Soc. Am.
139
(
3
),
1168
1178
(
2016
).
14.
E.
Fernandez-Grande
, “
Sound field reconstruction in a room from spatially distributed measurements
,” in
Proceedings of 23rd International Congress on Acoustics
(
2019
), pp.
4961
4968
.
15.
J.
Kaipio
and
E.
Somersalo
, “
Statistical and Computational Inverse Problems
,” in
Applied Mathematical Sciences
, Vol.
160
(
Springer Science+Business Media Inc
,
New York, NY
,
2005
).
16.
S. A.
Verburg
and
E.
Fernandez-Grande
, “
Reconstruction of the sound field in a room using compressive sensing
,”
J. Acoust. Soc. Am.
143
(
6
),
3770
3779
(
2018
).
17.
R.
Mignot
,
G.
Chardon
, and
L.
Daudet
, “
Low frequency interpolation of room impulse responses using compressed sensing
,”
IEEE/ACM Trans. Audio Speech, Lang. Process.
22
(
1
),
205
216
(
2014
).
18.
N.
Antonello
,
E.
de Sena
,
M.
Moonen
,
P. A.
Naylor
, and
T.
van Waterschoot
, “
Room impulse response interpolation using a sparse spatio-temporal representation of the sound field
,”
IEEE/ACM Trans. Audio Speech Lang. Process.
25
(
10
),
1929
1941
(
2017
).
19.
P.
Gerstoft
,
C. F.
Mecklenbräuker
,
W.
Seong
, and
M.
Bianco
, “
Introduction to compressive sensing in acoustics
,”
J. Acoust. Soc. Am.
143
,
3731
3736
(
2018
).
20.
F.
Lluís
,
P.
Martínez-Nuevo
,
M. B.
Møller
, and
S. E.
Shepstone
, “
Sound field reconstruction in rooms: Inpainting meets super-resolution
,”
J. Acoust. Soc. Am.
148
(
2
),
649
659
(
2020
).
21.
T.
Pham Vu
and
H.
Lissek
, “
Low frequency sound field reconstruction in a non-rectangular room using a small number of microphones
,”
Acta Acustica
4
(
2
),
5
(
2020
).
22.
B.-K.
Kim
and
J.-G.
Ih
, “
On the reconstruction of the vibro–acoustic field over the surface enclosing an interior space using the boundary element method
,”
J. Acoust. Soc. Am.
100
(
5
),
3003
3016
(
1996
).
23.
C.-X.
Bi
,
Y.
Liu
,
Y.-B.
Zhang
, and
L.
Xu
, “
Sound field reconstruction using inverse boundary element method and sparse regularization
,”
J. Acoust. Soc. Am.
145
(
5
),
3154
3162
(
2019
).
24.
S. F.
Wu
, “
Methods for reconstructing acoustic quantities based on acoustic pressure measurements
,”
J. Acoust. Soc. Am.
124
(
5
),
2680
2697
(
2008
).
25.
R.
Anderssohn
,
S.
Marburg
, and
C.
Großmann
, “
Fem-based reconstruction of sound pressure field damped by partially absorbing boundary conditions
,”
Mech. Res. Commun.
33
(
6
),
851
859
(
2006
).
26.
N.
Xiang
, “
Model-based Bayesian analysis in acoustics - A tutorial
,”
J. Acoust. Soc. Am.
148
(
2
),
1101
1120
(
2020
).
27.
A.
Gelman
,
J. B.
Carlin
,
H. S.
Stern
,
D. B.
Dunson
,
A.
Vehtari
, and
D. B.
Rubin
,
Bayesian Data Analysis
in
Chapman & Hall/CRC Texts in Statistical Science
, 3rd ed. (
CRC Press
,
Hoboken
,
2013
).
28.
J.
Antoni
, “
A Bayesian approach to sound source reconstruction: Optimal basis, regularization, and focusing
,”
J. Acoust. Soc. Am.
131
(
4
),
2873
2890
(
2012
).
29.
D.
Beaton
and
N.
Xiang
, “
Room acoustic modal analysis using Bayesian inference
,”
J. Acoust. Soc. Am.
141
(
6
),
4480
4493
(
2017
).
30.
N.
Xiang
,
P.
Goggans
,
T.
Jasa
, and
P.
Robinson
, “
Bayesian characterization of multiple-slope sound energy decays in coupled-volume systems
,”
J. Acoust. Soc. Am.
129
(
2
),
741
752
(
2011
).
31.
D.
Bush
and
N.
Xiang
, “
A model-based Bayesian framework for sound source enumeration and direction of arrival estimation using a coprime microphone array
,”
J. Acoust. Soc. Am.
143
(
6
),
3934
3945
(
2018
).
32.
C. R.
Landschoot
and
N.
Xiang
, “
Model-based Bayesian direction of arrival analysis for sound sources using a spherical microphone array
,”
J. Acoust. Soc. Am.
146
(
6
),
4936
4946
(
2019
).
33.
D.
Caviedes-Nozal
,
F. M.
Heuchel
,
J.
Brunskog
,
N. A. B.
Riis
, and
E.
Fernandez-Grande
, “
A Bayesian spherical harmonics source radiation model for sound field control
,”
J. Acoust. Soc. Am.
146
(
5
),
3425
3435
(
2019
).
34.
C. J.
Fackler
,
N.
Xiang
, and
K. V.
Horoshenkov
, “
Bayesian acoustic analysis of multilayer porous media
,”
J. Acoust. Soc. Am.
144
(
6
),
3582
3592
(
2018
).
35.
J.-D.
Chazot
,
E.
Zhang
, and
J.
Antoni
, “
Acoustical and mechanical characterization of poroelastic materials using a bayesian approach
,”
J. Acoust. Soc. Am.
131
(
6
),
4584
4595
(
2012
).
36.
F.
Jacobsen
and
P. M.
Juhl
,
Fundamentals of General Linear Acoustics
(
Wiley
,
Chichester
,
2013
).
37.
N.
Metropolis
,
A. W.
Rosenbluth
,
M. N.
Rosenbluth
,
A. H.
Teller
, and
E.
Teller
, “
Equation of state calculations by fast computing machines
,”
J. Chem. Phys.
21
(
6
),
1087
1092
(
1953
).
38.
S.
Brooks
,
A.
Gelman
,
G. L.
Jones
, and
Meng
Xiao-Li
, eds.,
Handbook of Markov Chain Monte Carlo (Chapman & Hall/CRC Handbooks of Modern Statistical Methods)
(
CRC Press Taylor & Francis
,
Boca Raton, Fla
.,
2011
).
39.
M. D.
Hoffman
and
A.
Gelman
, “
The no-u-turn sampler: Adaptively setting path lengths in hamiltonian monte carlo
,”
J. Mach. Learn. Res.
15
,
1593
1623
(
2014
).
40.
M.
Betancourt
, “
A conceptual introduction to hamiltonian monte carlo
,” arXiv: Methodology (1701.02434) (
2017
).
41.
B.
Carpenter
,
A.
Gelman
,
M. D.
Hoffman
,
D.
Lee
,
B.
Goodrich
,
M.
Betancourt
,
M.
Brubaker
,
J.
Guo
,
P.
Li
, and
A.
Riddell
, “
Stan: A probabilistic programming language
,”
J. Stat. Softw.
76
(
1
),
1
32
(
2017
).
42.
Stan Development Team, Stan reference manual
, version 2.19, https://mc-stan.org (
2019
) (Last viewed Nov. 20, 2021).
43.
E.
Jaynes
, “
Prior probabilities
,”
IEEE Trans. Syst. Sci. Cybern.
4
(
3
),
227
241
(
1968
).
44.
S.
Aja-Fernández
and
G.
Vegas-Sánchez-Ferrero
,
Statistical Analysis of Noise in MRI: Appendix A
(
Springer International Publishing
,
Cham
,
2016
).
45.
M.
Pastor
,
M.
Binda
, and
T.
Harčarik
, “
Modal assurance criterion
,”
Procedia Eng.
48
,
543
548
(
2012
).
46.
M.
Hahmann
,
S. A.
Verburg
, and
E.
Fernandez-Grande
, “
Analysis of the sound field in a room using dictionary learning
,” in
Proceedings of the 23rd International Congress on Acoustics
(
2019
).
47.
S.
Diamond
and
S.
Boyd
, “
CVXPY: A python-embedded modeling language for convex optimization
,”
J. Mach. Learn. Res.
17
(
1
),
2909
2913
(
2016
).