Sound field reconstruction in reverberant environments with rigid spherical microphone arrays

Abstract

To achieve accurate sound field reconstruction in reverberant environments, a method based on the principle of equivalent sources and a mixed-wave model is proposed with rigid spherical microphone arrays. The method assumes that the near-field sound sources are sparsely distributed and decomposes the sound field into the direct sound components modelled by monopole sources and the reflected sound components represented by plane waves, thereby constructing a mixed-wave model. The alternating direction method of multipliers is then used to alternately estimate the weights of monopole sources and plane waves, enabling the sound field decoupling. Finally, the target sound field in the region of interest is reconstructed using the obtained weights. Both numerical simulations and finite element method-based virtual validation experiment demonstrate that the proposed method generally expands the effective reconstruction region to twice that of Ambisonics, while improving reconstruction accuracy by 32% and 12.5% within a region spanning two to three times the array radius, compared with the plane wave decomposition method.

1. Introduction

Sound field reconstruction^1,2 estimates the acoustic information at arbitrary positions by post-processing the signals captured by microphone arrays, and this technique has been widely applied in areas such as personal sound zones,^3,4 virtual reality,^5,6 and noise control.^7,8 The effective reconstruction region and application scenarios are closely related to the geometry of microphone arrays. Common array configurations include linear arrays,^9–11 planar arrays,^12–14 and spherical arrays.^15–17 Compared to linear and planar arrays, spherical arrays, especially rigid spherical arrays, are more compact, flexible, and capable of capturing more spatial acoustic information.

Ambisonics,^18–20 an early proposed sound field reconstruction method with spherical microphone arrays, represents sound fields using spherical harmonics and reconstructs the sound field at unmeasured locations by exploiting the estimated expansion coefficients. However, due to the limited number of microphones and discrete spatial sampling, Ambisonics performs poorly in reconstructing the sound field far from the array. To overcome this limitation, Samarasinghe et al. ²¹ employed a distributed higher-order microphone array to expand the effective reconstruction region. The equivalent source method (ESM), another sound field reconstruction method, assumes that the target sound field is generated by a set of monopole sources or plane waves surrounding the region of interest, and estimates source strengths or plane wave weights via sparse optimization to achieve sound field reconstruction. Fernandez-Grande ²² used the equivalent sources emitting spherical waves to model the sound field recorded by rigid spherical microphone arrays and adopted the regularization technique to estimate source strengths, achieving accurate three-dimensional sound field reconstruction. Verburg et al. ²³ represented sound fields using a limited number of plane waves and reconstructed the room’s spatial frequency responses over the extended regions via compressive sensing. Compared to Ambisonics, ESM is not constrained by the spherical harmonic order truncation and thus provides a larger effective reconstruction region. However, the above ESM is based on the free-field assumption, which limits its applicability to reverberant environments, where sound reflections violate this assumption and degrade reconstruction accuracy.

To achieve accurate sound field reconstruction in reverberant environments, Koyama et al.^24,25 established a mixed wave model with monopole sources and plane waves. Damiano et al. ²⁶ modeled early reflections using the image source method, represented the late-reverberant components with a set of plane waves, and simulated the directivity using multipole sources. These mixed wave models are constructed based on linear and planar microphone arrays for one-dimensional and two-dimensional sound field reconstruction. In this paper, we extend the mixed wave modeling framework to three-dimensional space. The signals captured by rigid spherical microphone arrays are decoupled into direct and reflected components using monopole sources and plane waves. Then, sparse constraints are imposed, and an alternating direction method of multipliers (ADMM) is used to estimate the weights of monopole sources and plane waves, enabling accurate three-dimensional sound field reconstruction in reverberant environments.

The remainder of this paper is organized as follows. Sections 2 and 3 introduce the mixed wave model with rigid spherical microphone arrays and ADMM solver. Sections 4 and 5 present numerical simulations and finite element method-based virtual validation experiment, respectively. Section 6 draws the conclusions and outlines the perspectives for future research.

2. The mixed wave model with rigid spherical microphone arrays

As shown in Figure 1, a near-field source is located in a bounded space $\mathrm{Z}$ . The sound field induced by the source in this space satisfies the inhomogeneous Helmholtz equation

$\left({\nabla }^2+k^2\right)p\left(\mathbf{r},\omega \right)=-Q\left(\mathbf{r},\omega \right),$

(1)

where $p\left(\mathbf{r},\omega \right)$ represents the sound pressure at the position $\mathbf{r}=\left(r,\theta ,\varphi \right)$ , $r$ is the distance from the origin, $\theta \in \left[0,\mathrm{\pi }\right]$ is the elevation angle, $\varphi \in [0,2\mathrm{\pi })$ is the azimuth angle, $\omega$ is the circular frequency, $k$ is the wave number, $Q\left(\mathbf{r},\omega \right)$ is the source distribution in the space $\mathrm{Z}$ , and ${\nabla }^{2}f={\partial }^{2}f/\partial x^2+{\partial }^{2}f/\partial y^2+{\partial }^{2}f/\partial z^2$ is the Laplace operator. Therefore, $p\left(\mathbf{r},\omega \right)$ can be expressed as the sum of the particular solution $p_{\mathrm{p}}\left(\mathbf{r},\omega \right)$ and the homogeneous solution $p_{\mathrm{h}}\left(\mathbf{r},\omega \right)$ ,

$p\left(\mathbf{r},\omega \right)=p_{\mathrm{p}}\left(\mathbf{r},\omega \right)+p_{\mathrm{h}}\left(\mathbf{r},\omega \right).$

(2)

Figure 1.

The mixed wave model.

In Eq. (2), $p_{\mathrm{p}}\left(\mathbf{r},\omega \right)$ denotes the direct sound components, which can be modelled by monopole sources

$p_{\mathrm{p}}\left(\mathbf{r},\omega \right)={\sum }_{i=1}^{I}G\left(\mathbf{r}|{\mathbf{r}}_{\mathrm{S}i}\right)x\left({\mathbf{r}}_{\mathrm{S}i}\right),$

(3)

where I is the total number of monopole sources, $G\left(\mathbf{r}|{\mathbf{r}}_{\mathrm{S}i}\right)$ is the transfer function from the $i\hspace{0.5em}\text{th}$ monopole source at ${\mathbf{r}}_{\mathrm{S}i}=\left(r_{\mathrm{S}i},{\theta }_{\mathrm{S}i},{\varphi }_{\mathrm{S}i}\right)$ to the position $\mathbf{r}=\left(r,\theta ,\varphi \right)$ , and $x\left({\mathbf{r}}_{\mathrm{S}i}\right)$ is the strength (weight) of the $i\hspace{0.5em}\text{th}$ monopole source. Due to the rigid boundary of the spherical microphone array, acoustic scattering occurs on its surface. ²⁷ Therefore, $G\left(\mathbf{r}|{\mathbf{r}}_{\mathrm{S}i}\right)$ is modeled using the Neumann Green’s function ²⁸ to account for the scattering effect,

$G\left(\mathbf{r}|{\mathbf{r}}_{\mathrm{S}i}\right)={\displaystyle \sum }_{n=0}^{\infty }{\displaystyle \sum }_{m=-n}^n-4\mathrm{\pi }\mathrm{j}{h}_n^{\left(2\right)}\left(k{r}_{\mathrm{S}i}\right)\left[j_n\left(kr\right)-\frac{j_n^{\prime }\left(ka\right)}{{h_n^{\left(2\right)}}^{\prime }\left(ka\right)}{h}_n^{\left(2\right)}\left(kr\right)\right]Y_n^m\left(\theta ,\varphi \right)Y_n^{m*}\left({\theta }_{\mathrm{S}i},{\varphi }_{\mathrm{S}i}\right),$

(4)

where $h_n^{\left(2\right)}(·)$ is the spherical Hankel function of the second kind, $j_n(·)$ is the spherical Bessel function of the first kind, ${h_n^{\left(2\right)}}^{\prime }(·)$ and $j_n^{\prime }(·)$ are the derivatives of $h_n^{\left(2\right)}(·)$ and $j_n(·)$ , ${(·)}^{*}$ denotes the conjugate, $a$ is the radius of the spherical microphone array, $Y_n^m\left(\theta ,\varphi \right)$ is the spherical harmonic function in the direction of $\left(\theta ,\varphi \right)$ , and $\mathrm{j}=\sqrt{-1}$ is the imaginary unit. $p_{\mathrm{h}}\left(\mathbf{r},\omega \right)$ denotes the reflected sound components, which is represented by plane waves,

$p_{\mathrm{h}}\left(\mathbf{r},\omega \right)={\sum }_{l=1}^{L}W\left(\mathbf{r}|{\mathrm{\Omega }}_{\mathrm{S}l}\right)u\left({\mathrm{\Omega }}_{\mathrm{S}l}\right),$

(5)

where L is the number of plane waves, $W\left(\mathbf{r}|{\mathrm{\Omega }}_{\mathrm{S}l}\right)$ is the transfer function from the l th propagation direction ${\mathrm{\Omega }}_{\mathrm{S}l}$ to the position $\mathbf{r}=\left(r,\theta ,\varphi \right)$ , and $u\left({\mathrm{\Omega }}_{\mathrm{S}l}\right)$ is the weight of the l th plane wave. Similarly, taking scattering into consideration, $W\left(\mathbf{r}|{\mathrm{\Omega }}_{\mathrm{S}l}\right)$ is expressed as

$W\left(\mathbf{r}|{\mathrm{\Omega }}_{\mathrm{S}l}\right)={\displaystyle \sum }_{n=0}^{\infty }4\mathrm{\pi }{\mathrm{j}}^n\left[j_n\left(kr\right)-\frac{j_n^{\prime }\left(ka\right)}{{h_n^{\left(2\right)}}^{\prime }\left(ka\right)}{h}_n^{\left(2\right)}\left(kr\right)\right]{\displaystyle \sum }_{m=-n}^n{Y}_n^m\left(\theta ,\varphi \right)Y_n^{m*}\left({\theta }_{\mathrm{S}l},{\varphi }_{\mathrm{S}l}\right).$

(6)

Based on Eq. (2), Eq. (3), and Eq. (5), we construct a mixed wave model to express the pressures measured by rigid spherical microphone arrays,

$\mathbf{p}=\mathbf{Gx}+\mathbf{Wu},$

(7)

where $\mathbf{p}={[p({\mathbf{r}}_{\mathrm{M}1},\omega ),p({\mathbf{r}}_{\mathrm{M}2},\omega ),\cdots ,p({\mathbf{r}}_{\mathrm{M}Q},\omega )]}^{\mathrm{T}}\in {\mathrm{ℂ}}^{Q\times 1}$ is the vector of microphone pressures, ${\mathbf{r}}_{\mathrm{M}q}=\left(a,{\theta }_{\mathrm{M}q},{\varphi }_{\mathrm{M}q}\right)$ is the position of the $q\hspace{0.5em}\text{th}$ microphone, Q is the number of microphones, the superscript ^“T” denotes the transpose, $\mathrm{ℂ}$ denotes the set of complex numbers, $\mathbf{x}={[x\left({\mathbf{r}}_{\mathrm{S}1}\right),x\left({\mathbf{r}}_{\mathrm{S}2}\right),\dots ,x\left({\mathbf{r}}_{\mathrm{S}I}\right)]}^{\mathrm{T}}\in {\mathrm{ℂ}}^{I\times 1}$ and $\mathbf{u}={[u\left({\mathrm{\Omega }}_{\mathrm{S}1}\right),u\left({\mathrm{\Omega }}_{\mathrm{S}2}\right),\dots ,u\left({\mathrm{\Omega }}_{\mathrm{S}L}\right)]}^{\mathrm{T}}\in {\mathrm{ℂ}}^{L\times 1}$ are the vectors of the weights of monopole sources and plane waves, respectively, $\mathbf{G}=[{\mathbf{G}}_1,{\mathbf{G}}_2,\cdots ,{\mathbf{G}}_I]\in {\mathrm{ℂ}}^{Q\times I}$ with ${\mathbf{G}}_i={[G\left({\mathbf{r}}_{\mathrm{M}1}|{\mathbf{r}}_{\mathrm{S}i}\right),G\left({\mathbf{r}}_{\mathrm{M}2}|{\mathbf{r}}_{\mathrm{S}i}\right),\cdots ,G\left({\mathbf{r}}_{\mathrm{M}Q}|{\mathbf{r}}_{\mathrm{S}i}\right)]}^{\mathrm{T}}\in {\mathrm{ℂ}}^{Q\times 1}$ , and $\mathbf{W}=[{\mathbf{W}}_1,{\mathbf{W}}_2,\cdots ,{\mathbf{W}}_L]\in {\mathrm{ℂ}}^{Q\times L}$ with ${\mathbf{W}}_l={[W\left({\mathbf{r}}_{\mathrm{M}1}|{\mathrm{\Omega }}_{\mathrm{S}l}\right),W\left({\mathbf{r}}_{\mathrm{M}2}|{\mathrm{\Omega }}_{\mathrm{S}l}\right),\dots ,W\left({\mathbf{r}}_{\mathrm{M}Q}|{\mathrm{\Omega }}_{\mathrm{S}l}\right)]}^{\mathrm{T}}\in {\mathrm{ℂ}}^{Q\times 1}$ .

3. ADMM solver

The source distribution is usually spatially sparse, so $\mathbf{x}$ will have few non-zero elements. According to Vekua’s theory, ²⁹ the reflected sound component in a bounded region can be well approximated by a limited number of plane waves. When L in Eq. (5) is very large, the most elements of $\mathbf{u}$ can be approximated as 0. Therefore, to estimate $\mathbf{x}$ and $\mathbf{u}$ , we construct the optimization model,

$\begin{array}{l}\min \hspace{1em}{\lambda }_1{\Vert \mathbf{x}\Vert }_0+{\lambda }_2{\Vert \mathbf{u}\Vert }_0\\ \mathrm{s}.\mathrm{t}.\hspace{1em}\mathbf{p}=\mathbf{Gx}+\mathbf{Wu}\end{array},$

(8)

where ${\Vert ·\Vert }_0$ denotes the ${\mathcal{l}}_0$ norm, and ${\lambda }_1>0$ and ${\lambda }_2>0$ are the regularization parameters. ${\mathcal{l}}_0$ norm minimization problem is non-convex, which is challenging to solve, thus Eq. (8) is relaxed into a ${\mathcal{l}}_1$ norm minimization problem

$\begin{array}{l}\min \hspace{1em}{\lambda }_1{\Vert \mathbf{x}\Vert }_1+{\lambda }_2{\Vert \mathbf{u}\Vert }_1\\ \mathrm{s}.\mathrm{t}.\hspace{1em}\mathbf{p}=\mathbf{Gx}+\mathbf{Wu}\end{array}.$

(9)

ADMM ³⁰ is adopted to solve Eq. (9) due to its computational efficiency and flexibility, where each variable is updated alternately while keeping the other variables fixed. The steps are as follows. We introduce two additional variables, ${\mathbf{x}}_1$ and ${\mathbf{u}}_1$ , and convert Eq. (9) into

$\begin{array}{l}\min \hspace{1em}{\lambda }_1{\Vert {\mathbf{x}}_1\Vert }_1+{\lambda }_2{\Vert {\mathbf{u}}_1\Vert }_1+\frac{1}{2}{\Vert \mathbf{p}-\mathbf{Gx}-\mathbf{Wu}\Vert }_2^2\\ \mathrm{s}.\mathrm{t}.\hspace{1em}\mathbf{x}={\mathbf{x}}_1,\mathbf{u}={\mathbf{u}}_1\end{array}.$

(10)

Eq. (10) can be further transformed into a non-constrained optimization problem,

$\begin{array}{l}\mathcal{L}\left(\mathbf{x},{\mathbf{x}}_1,{\mathbf{v}}_x,\mathbf{u},{\mathbf{u}}_1,{\mathbf{v}}_u\right)\\ =\frac{1}{2}{\Vert \mathbf{p}-\mathbf{Gx}-\mathbf{Wu}\Vert }_2^2+{\lambda }_1{\Vert {\mathbf{x}}_1\Vert }_1+\text{Re}\left\{{\mathbf{v}}_x^{\mathrm{H}}\left(\mathbf{x}-{\mathbf{x}}_1\right)\right\}+\frac{{\rho }_1}{2}{\Vert \mathbf{x}-{\mathbf{x}}_1\Vert }_2^2\\ +{\lambda }_2{\Vert {\mathbf{u}}_1\Vert }_1+\text{Re}\left\{{\mathbf{v}}_u^{\mathrm{H}}\left(\mathbf{u}-{\mathbf{u}}_1\right)\right\}+\frac{{\rho }_2}{2}{\Vert \mathbf{u}-{\mathbf{u}}_1\Vert }_2^2\end{array},$

(11)

where ${\mathbf{v}}_x$ and ${\mathbf{v}}_u$ are the Lagrangian multipliers, $\text{Re}\{·\}$ represents taking the real part, the superscript ^“H” denotes the Hermitian conjugate, and ${\rho }_1>0$ and ${\rho }_2>0$ are penalty parameters. Let ${\mathbf{\mu }}_x={\mathbf{v}}_x/{\rho }_1$ and ${\mathbf{\mu }}_u={\mathbf{v}}_u/{\rho }_2$ . Eq. (11) becomes

$\begin{array}{l}\mathcal{L}\left(\mathbf{x},{\mathbf{x}}_1,{\mathbf{\mu }}_x,\mathbf{u},{\mathbf{u}}_1,{\mathbf{\mu }}_u\right)\\ =\frac{1}{2}{\Vert \mathbf{p}-\mathbf{Gx}-\mathbf{Wu}\Vert }_2^2+{\lambda }_1{\Vert {\mathbf{x}}_1\Vert }_1+\frac{{\rho }_1}{2}\left({\Vert \mathbf{x}-{\mathbf{x}}_1+{\mathbf{\mu }}_x\Vert }_2^2-{\Vert {\mathbf{\mu }}_x\Vert }_2^2\right)\\ +{\lambda }_2{\Vert {\mathbf{u}}_1\Vert }_1+\frac{{\rho }_2}{2}\left({\Vert \mathbf{u}-{\mathbf{u}}_1+{\mathbf{\mu }}_u\Vert }_2^2-{\Vert {\mathbf{\mu }}_u\Vert }_2^2\right)\end{array}.$

(12)

Eq. (12) can be solved iteratively. In the $\gamma \hspace{0.5em}\text{th}$ iteration, by taking the partial derivatives of $\mathcal{L}\left(\mathbf{x},{\mathbf{x}}_1,{\mathbf{\mu }}_x,\mathbf{u},{\mathbf{u}}_1,{\mathbf{\mu }}_u\right)$ with respect to each variable and setting these derivatives to 0, we can obtain

$\{\begin{array}{l}\begin{array}{l}{\mathbf{x}}^{\left(\gamma +1\right)}={\left({\mathbf{G}}^{\mathrm{H}}\mathbf{G}+{\rho }_1\mathbf{I}\right)}^{-1}\left({\mathbf{G}}^{\mathrm{H}}\mathbf{p}-{\mathbf{G}}^{\mathrm{H}}{\mathbf{Wu}}^{\left(\gamma \right)}+{\rho }_1{\mathbf{x}}^{\left(\gamma \right)}-{\rho }_1{\mathbf{\mu }}_x^{\left(\gamma \right)}\right)\\ {\mathbf{x}}_1^{\left(\gamma +1\right)}={\mathcal{T}}_{{\lambda }_1/{\rho }_1}\left({\mathbf{x}}^{\left(\gamma +1\right)}+{\mathbf{\mu }}_x^{\left(\gamma \right)}\right)\\ {\mathbf{\mu }}_x^{\left(\gamma +1\right)}={\mathbf{\mu }}_x^{\left(\gamma \right)}+{\mathbf{x}}^{\left(\gamma +1\right)}-{\mathbf{x}}_1^{\left(\gamma +1\right)}\end{array}\\ \begin{array}{l}{\mathbf{u}}^{\left(\gamma +1\right)}={\left({\mathbf{W}}^{\mathrm{H}}\mathbf{W}+{\rho }_2\mathbf{I}\right)}^{-1}\left({\mathbf{W}}^{\mathrm{H}}\mathbf{p}-{\mathbf{W}}^{\mathrm{H}}{\mathbf{Gx}}^{\left(\gamma \right)}+{\rho }_2{\mathbf{u}}_1^{\left(\gamma \right)}-{\rho }_2{\mathbf{\mu }}_u^{\left(\gamma \right)}\right)\\ {\mathbf{u}}_1^{\left(\gamma +1\right)}={\mathcal{T}}_{{\lambda }_2/{\rho }_2}\left({\mathbf{u}}^{\left(\gamma +1\right)}+{\mathbf{\mu }}_u^{\left(\gamma \right)}\right)\\ {\mathbf{\mu }}_u^{\left(\gamma +1\right)}={\mathbf{\mu }}_u^{\left(\gamma \right)}+{\mathbf{u}}^{\left(\gamma +1\right)}-{\mathbf{u}}_1^{\left(\gamma +1\right)}\end{array}\end{array},$

(13)

where ${\mathcal{T}}_{\alpha }(·)=\text{sign}(·)\max \left\{|·|-\alpha ,0\right\}$ is the soft-thresholding. Each variable is initialized as a zero vector. The iterations stop when $\left({{\Vert {\mathbf{x}}^{\left(\gamma +1\right)}-{\mathbf{x}}^{\left(\gamma \right)}\Vert }_2/\Vert {\mathbf{x}}^{\left(\gamma \right)}\Vert }_2\right)<{10}^{-3}$ and $\left({{\Vert {\mathbf{p}}^{\left(\gamma +1\right)}-{\mathbf{p}}^{\left(\gamma \right)}\Vert }_2/\Vert {\mathbf{p}}^{\left(\gamma \right)}\Vert }_2\right)<{10}^{-3}$ are satisfied, or when the maximum number of iterations ${\gamma }^{(\hspace{0.17em}\max \hspace{0.17em})}=1000$ is reached.

The estimated strengths $\widehat{\mathbf{x}}$ and weights $\widehat{\mathbf{u}}$ are used to reconstruct sound field,

${\mathbf{p}}^{\widehat{\mathrm{R}}\mathrm{e}}={\mathbf{G}}^{\widehat{\mathrm{R}}\mathrm{e}}\widehat{\mathbf{x}}+{\mathbf{W}}^{\widehat{\mathrm{R}}\mathrm{e}}\widehat{\mathbf{u}},$

(14)

where ${\mathbf{p}}^{\widehat{\mathrm{R}}\mathrm{e}}\in {\mathrm{ℂ}}^{V\times I}$ is the vector of reconstructed sound pressures at V positions ${\mathbf{r}}_{\text{Re}v}=\left(r_{\text{Re}v},{\theta }_{\text{Re}v},{\varphi }_{\text{Re}v}\right)$ ( $v=1,2,\dots ,V$ ), and ${\mathbf{G}}^{\widehat{\mathrm{R}}\mathrm{e}}\in {\mathrm{ℂ}}^{V\times I}$ and ${\mathbf{W}}^{\widehat{\mathrm{R}}\mathrm{e}}\in {\mathrm{ℂ}}^{V\times L}$ can be obtained by replacing ${\mathbf{r}}_{\mathrm{M}q}$ in $\mathbf{G}$ and $\mathbf{W}$ with ${\mathbf{r}}_{\text{Re}v}$ .

4. Numerical simulations

In this section, the performance of the proposed method is evaluated and compared with the Ambisonics method ¹⁸ and the plane wave decomposition method ³¹ through numerical simulations. Brüel & Kjær Type 8608 rigid spherical microphone array with 36 microphones is used. The radius is 0.0975 m. Figure 2 shows the microphone array and the mixed wave model. The center of the array is placed at the origin of the coordinate system, and a sound source is located on the y–z plane at $x=-1 \mathrm{m}$ . A spherical equivalent source surface is set for plane waves, which is 1 m from the center of the array, with angular intervals of $\mathrm{\Delta }\varphi =10°$ and $\mathrm{\Delta }\theta =10°$ . A rectangular equivalent source surface is used for spherical waves, which is at $x=-1 \mathrm{m}$ , with boundaries defined by $y_{\min }=-0.5 \mathrm{m}$ , $y_{\max }=0.5 \mathrm{m}$ , $z_{\min }=-0.5 \mathrm{m}$ and $z_{\max }=0.5 \mathrm{m}$ . The grid spacing is set to 0.05 m. ${\rho }_1={\rho }_2=1$ ³² , and ${\lambda }_1={\lambda }_2=0.1$ .To clearly demonstrate and compare performance, the reconstruction region ${\mathrm{Z}}^{\widehat{\mathrm{R}}\mathrm{e}}$ is set as a 1 m×1 m square in the x-y plane, centered at the origin. The theoretical sound pressures are given by the image source method. ³³ To evaluate the performance of each method, the reconstruction error is defined as

$\text{Error}=\frac{{\Vert {\mathbf{p}}^{\text{The}}-{\mathbf{p}}^{\text{Re}}\Vert }_2}{{\Vert {\mathbf{p}}^{\text{The}}\Vert }_2}\times 100\%,$

(15)

where ${\mathbf{p}}^{\text{The}}$ is the vector of theoretical pressures at the reconstruction positions, and ${\mathbf{p}}^{\text{Re}}$ is the vector of the reconstructed pressures.

Figure 2.

Spherical microphone array and the equivalent source model.

4.1. Sound field reconstruction in a semi-enclosed rectangular space

To simulate a semi-enclosed rectangular space, an infinitely long rectangular duct is modeled. Free boundary conditions are applied along the x-axis, indicating the space is infinitely extended in both the positive and negative x-directions. Along the y- and z-axes, rigid boundary conditions are set at y=0.75 m, y=-0.75 m, z=1.2 m and z=-1.2 m, with an absorption coefficient of 0.3. According to the Eyring formula, ³⁴ T₆₀=190 ms. A source is placed at (-1 m, 0, 0), with a strength of 94 dB (referenced to 20 μPa).

The proposed method, the Ambisonics method, and the plane wave decomposition method are used to reconstruct the sound field within the region ${\mathrm{Z}}^{\widehat{\mathrm{R}}\mathrm{e}}$ . Figure 3 shows the theoretical values of the real parts of the sound pressures at 1000 Hz in the reconstruction region, obtained using the 10th-order image source method. Figure 4 presents the reconstruction results and corresponding reconstruction error. The black dotted circle in the figures represents the outline of the used rigid spherical microphone array. As shown in Figure 4(a1) and (a2), the reconstruction result of Ambisonics is close to the theoretical sound field within the array radius, with reconstruction error below 10%. However, large distortions occur when far away the array due to the order truncation of spherical harmonics. The result of the plane wave decomposition method, as shown in Figure 4(b1), shows distinct plane wavefronts. Comparing Figure 4(a2) and (b2), it can be observed that, the area where reconstruction error is below 10% (i.e., the bright region) is larger in the map of the plane wave decomposition method. For the proposed method, as shown in Figure 4(c1) and (c2), the reconstructed real parts of the sound pressures coincide best with the theoretical values. Correspondingly, the area where reconstruction error is less than 10% is significantly larger than the other two methods. This demonstrates that the proposed method achieves more accurate reconstruction in a larger region, which benefits from the mixed wave model.

Figure 3.

The theoretical values of the real parts of the sound pressures at 1000 Hz in a semi-enclosed rectangular space.

Figure 4.

Reconstruction results and error in the semi-enclosed rectangular space.

As for computational cost, the proposed method requires longer running time for sound field reconstruction than Ambisonics and the plane wave decomposition method. Specifically, for the above case, the running time of the proposed method, Ambisonics and the plane wave decomposition method is 4.29 s, 1.09 s, and 2.07 s on the 3.6GHz AMD Ryzen7 7745HX CPU, respectively.

4.2. Sound field reconstruction in a closed rectangular space

The proposed method is further applied to reconstruct the sound field within a closed space, and compared with the Ambisonics method and the plane wave decomposition method. To simulate a closed space, rigid planar boundaries are applied at x=1.5 m, x=-1.5 m, y=0.75 m, y=-0.75 m, z=1.2 m and z=-1.2 m with an absorption coefficient of 0.3. The sound source position and strength remain the same as those used in Section 4.1. Figure 5 shows the closed space and the theoretical values of the real parts of the sound pressures at 1000 Hz, which are computed through the 10th-order image source method.

Figure 5.

The closed space and the theoretical values of the real parts of the sound pressures at 1000 Hz.

Figure 6 presents the sound field reconstruction results and reconstruction error at 1000 Hz in the closed rectangular space. By comparing the results and the error of the three methods, similar conclusions can be drawn. Among the three methods, the reconstructed field of the proposed method is the most consistent with the theoretical field. Correspondingly, the area with reconstruction error less than 10% (i.e., the bright region) is the largest for the proposed method, followed by the plane wave decomposition method. The Ambisonics method maintains reconstruction error below 10% only within the region close to array. These results demonstrate that the proposed method outperforms the plane wave decomposition method and the Ambisonics method.

Figure 6.

Reconstruction results and error in the closed rectangular space.

Monte Carlo trials are conducted to analyze the reconstruction performance under different frequencies and absorption coefficients. The frequency ranges from 500 Hz to 2000 Hz in a step of 500 Hz. The absorption coefficient $\alpha$ varies from 0.1 to 0.7 in a step of 0.2, corresponding to reverberation time T₆₀ of 568 ms, 159 ms, 82 ms, and 47 ms. For each frequency and absorption coefficient, 50 trials are considered. In each trial, a monopole source in the mixed wave model is randomly selected as the sound source, with a source strength of 94 dB (referenced to 20 μPa).

Figure 7 shows average reconstruction error of the three methods within regions marked by radius 0∼a, a∼2a and 2a∼3a. In the region 0∼a, the average reconstruction error of Ambisonics, the plane wave decomposition method, and the proposed method is below 10% across all cases. In other words, the area corresponding to error below 10% in the error map accounts for 100%. This indicates that all three methods can achieve accurate reconstruction within the region of 0∼a. However, the maximum error of the three methods is 9.4%, 1.5% and 1%, respectively, indicating that the proposed method has the higher reconstruction accuracy than the other two methods. In the region a∼2a, the average reconstruction error of Ambisonics exceeds 10% under all conditions. By comparison, the plane wave decomposition method and the proposed method achieve reconstruction error below 10% over 43% and 75% of the error map, respectively. In the region 2a∼3a, the three methods show larger reconstruction error, and the proportion of the area where the error is below 10% decreases to approximately 6.3% for the plane wave decomposition method and 18.8% for the proposed method. In addition, the average reconstruction error of the proposed method is obviously lower than the plane wave decomposition for all cases. In summary, the proposed method achieves higher reconstruction accuracy over a larger spatial region compared with both Ambisonics and the plane wave decomposition method. Specifically, the proposed method generally expands the effective reconstruction region to twice that of Ambisonics, and even up to three times under low-frequency and weakly reverberant conditions. Meanwhile, it improves the reconstruction accuracy (measured by the proportion of the area with error below 10%) by 32% and 12.5% within a region spanning two to three times the array radius, compared with the plane wave decomposition method.

Figure 7.

The average reconstruction error within regions marked by radius $0\sim a$ 、 $a\sim 2a$ and $2a\sim 3a$ .

The proposed method outperforms both Ambisonics and the plane wave decomposition method primarily because these conventional methods rely on spherical harmonics expansion, whose effective reconstruction region is inherently constrained by the microphone array.^35,36 Specifically, due to the finite number of microphones and discrete spatial sampling, the spherical harmonics expansion must be truncated, which limits the spatial region over which Ambisonics and plane wave decomposition can achieve accurate reconstruction. In contrast, the proposed method formulates the mixed wave model directly in the spatial domain rather than in the spherical harmonics domain. As a result, it avoids the truncation-related limitations imposed by the microphone array, enabling more accurate sound field reconstruction over an extended spatial region.

5. Finite element method-based virtual validation experiment

In this section, the performance of the proposed method is evaluated based on finite element analyses. ³⁷ A sound field model, as shown in Figure 8(a), is built in the COMSOL. The gray cuboid represents the computational area, with dimensions of $4 \mathrm{m}\times 1.5 \mathrm{m}\times 2.4 \mathrm{m}$ . The sphere inside the cuboid represents the rigid spherical microphone array (Brüel & Kjær Type 8608, 36-channel rigid array), and the symbol “★” denotes the sound source. The sound source position, boundary conditions, and reconstruction region are consistent with those in Section 4.2. Figure 8(b) shows the theoretical values of the real parts of the sound pressures at 1000 Hz, provided by COMSOL.

Figure 8.

Sound field model and the theoretical value of the real part of the sound pressure field at 1000 Hz provided by COMSOL.

Figure 9 presents the reconstruction results and reconstruction error of the Ambisonics method, the plane wave decomposition method, and the proposed method. It is clear that the proposed method continues to demonstrate the best reconstruction performance. It provides a larger region where reconstruction error is below 10%, effectively enabling accurate sound field reconstruction in reverberant environments. These findings are consistent with the conclusions in Section 4.2.

Figure 9.

Reconstruction results and error in the closed rectangular space based on finite element analyses.

6. Conclusions and perspectives

This paper proposes a method with rigid spherical microphone arrays for sound field reconstruction in reverberant environments. This method establishes a mixed wave model that represents incident and reflected sound fields using monopole sources and plane waves, and achieves sound field reconstruction via sparse constraint and ADMM. The proposed method generally doubles the effective reconstruction region of Ambisonics and improves reconstruction accuracy by up to 32% and 12.5% within a range of two to three times the array radius, compared with the plane wave decomposition method. Its superiority stems from operating in the spatial domain and avoiding limitations associated with spherical harmonics expansion.

As the frequency increases and reflections become stronger, the proposed method experiences some performance degradation. Future work will focus on addressing this limitation to enhance the method’s applicability. In addition, the real-world validation experiments are currently not possible due to the lack of suitable experimental conditions. We therefore plan to improve our experimental setup to provide validation under practical acoustic conditions in the future.