A convolutional plane wave model for sound field reconstruction

The true acoustic pressure $\mathbf{p}\in\mathbb{C}^{N}$ at frequency $f$ and $N$ positions $\mathbf{r}\in\mathbb{R}^{3}$ within a domain $\Omega$ is assumed to be modelled as a linear combination of basis functions

where $\mathbf{x}\in\mathbb{C}^{M}$ are the coefficients and $\mathbf{H}\in\mathbb{C}^{N\times M}$ contains the $M$ basis functions. For example, plane propagating waves $e^{-\mathrm{j}\mathbf{k}^{\sf T}\mathbf{r}}$ are often used to model reverberant or far-fields, where $\mathbf{k}$ is the wavenumber vector ($\|\mathbf{k}\|_{2}=2\pi/\lambda$, $\lambda$ is the wavelength). In the case of plane waves, the $n,m^{\text{th}}$ element of $\mathbf{H}$ is $e^{-\mathrm{j}\mathbf{k}_{m}^{\sf T}\mathbf{r}_{n}}$, with $\mathbf{r}_{n}$ denoting the $n^{\text{th}}$ position and $\mathbf{k}_{m}$ the wavenumber vector from the $m^{\text{th}}$ incidence direction.

For the equation to hold, $\mathbf{H}$ must span the observed sound pressure field $\mathbf{p}$ over the complete domain $\Omega$. An observation of the sound pressure field is

where $\mathbf{M}\in\{0,1\}^{N_{\mathrm{obs}}\times N}$ is a binary mask selecting the $N_{\mathrm{obs}}$ available observations from the sound field $\mathbf{p}$. $\mathbf{H}_{\mathrm{obs}}=\mathbf{MH}$ is the model at the observed positions and $\mathbf{n}$ is an error vector, that accounts for measurement noise and model error.

An optimal set of coefficients $\hat{\mathbf{x}}$ is found by inversion of Eq. (2). To arrive to a stable solution, regularization is necessary as $\mathbf{H}$ is typically ill-conditioned or rank-deficient.hansen1998a Typically, a structure in the coefficients is imposed, such as in

in which case a $\ell_{1}$-norm penalty is applied to promote a sparse structure in the coefficients and $\hat{\cdot}$ denotes an estimate. The sound field is then reconstructed as

where $\hat{\mathbf{p}}\in\mathbb{C}^{N}$ is the reconstructed sound pressure at $N$ reconstruction positions.

When reconstructing the sound field over large spatial domains (i.e. much larger than the acoustic wavelength), it is useful to partition the global domain $\Omega$ into smaller, overlapping subdomains $\Omega_{sub}$. Correspondingly, the sound field $\mathbf{p}$ can be described by a collection of $S$ subdomain sound fields $\mathbf{p}_{s}\in\mathbb{C}^{N_{s}}$

$\mathbf{R}_{s}\in\{0,1\}^{N_{s}\times N}$ is a binary extraction operator to select the $N_{s}$ positions contained in the $s^{\text{th}}$ subdomain $\Omega_{s}$. The partitioning is illustrated on the left side of Fig. 1. For simplicity, all subdomains within a sound field are considered to have the same extent, such that $N_{s}$ is constant. Specifically, this study considers subdomains of extent one wavelength $\lambda$ in each dimension of the aperture. Note that Eq. (5) yields a redundant representation of the sound field if the subdomains overlap ($N_{s}S\geq N$).

A sound field can then be reconstructed within each subdomain of the partitioned observations $\mathbf{P}_{\mathrm{obs}}$, for example by applying the procedure in Sec. II.1 and finding coefficients $\hat{\mathbf{x}}_{x}$ via Eq. (3) to estimate each $\hat{\mathbf{p}}_{s}$. This yields the collection of reconstructed subdomain sound fields [see center right in Fig. 1].

where $\mathbf{H}_{s}$ are the model functions at the desired positions within each local subdomain, for example $\mathbf{H}_{s}=\mathbf{R}_{1}\mathbf{H}$, and

the estimated local coefficients. The reconstructed field $\hat{\mathbf{p}}$ is reassembled from $\hat{\mathbf{P}}$ as the mean of overlapping subdomain representations [see right side of Fig. 1],

where the diagonal matrix $\mathbf{W}=\text{diag}(\sum_{s}\mathbf{R}_{s}^{T}\mathbf{1}_{N_{s}})^{-1}$ normalizes by the spatial overlap of the subdomains and $\mathbf{1}_{N_{s}}$ denotes a vector of $N_{s}$ ones.

Such partitioning approaches counteract model mismatch, which occurs if a global model is suboptimal. For example in the global sound field model of Eq. (1), the chosen model functions in $\mathbf{H}$ might not span observations of a sound field across a large spatial domain.

Estimating independent coefficients for each subdomain allows for arbitrarily different wave components in each subdomain. Such independent treatment of local representations disregards the similarity or even redundancy between sound fields in nearby or overlapping local subdomains.

The coefficients in the local subdomain sound field model can also be understood as a collection of coefficient maps $\mathbf{x}_{m}\in\mathbb{C}^{S}$, the rows in $\mathbf{X}=[\mathbf{x}^{\sf T}_{1}\cdots\mathbf{x}^{\sf T}_{M}]^{\sf T}$. The $m^{\text{th}}$ row in $\mathbf{X}$ contains a coefficient for the $m^{\text{th}}$ local model function across all $S$ subdomain locations. When plane waves are used, the coefficient map $\mathbf{x}_{m}$ contains the spatial distribution of coefficients for the $m^{\text{th}}$ plane wave across subdomain representations.

This study explores, how the spatial variations across $\mathbf{x}_{m}$, the rows of $\mathbf{X}$, can be taken into account for reconstruction, such that the global structure (or topology) of the sound field can be preserved. For the remainder of the paper, we consider the reconstruction positions on a regular grid (containing also the measured positions as a subset). Further, we assume subdomains of equal size and with full overlap, such that $S=N$, the number of subdomains equals the number of positions in $\mathbf{p}$ (circular boundary conditions).

The true sound field across an $L$-dimensional aperture can be rewritten as a sum of $M$ convolutions:

where $\mathbf{h}_{m}\in\mathbb{C}^{N_{s}}$ describes the $m^{\text{th}}$ local filter (the $m^{\text{th}}$ column of $\mathbf{H}_{s}$), for example a plane wave. $\circledast\overset{L}{\cdots}\circledast$ denotes a circular convolution along $L\in{1,2,3}$ spatial dimensions. For example for $L=1$, the $n^{\text{th}}$ element of $\mathbf{p}$ is,

where $n=1,\ldots N$ and $\%$ is the modulo operation.

Compared to the collection of local sound fields in Sec. II.2, this convolutional model reflects the spatial relations of coefficients and enables a joint analysis of all local representations in the global field. For example, each subdomain representation can take the coefficients in nearby subdomains into account. We propose to estimate the coefficients of Eq. (9) as

The $\ell_{1}$ penalty promotes sparse coefficients, notably applied on the global coefficient vector. A penalty on the spatial differences of the coefficients promotes smooth coefficient maps $\hat{\mathbf{x}}_{m}$, weighted by the regularization parameter $\mu$. Specifically, $\Delta_{l}\mathbf{x}_{m}$ are the first order finite differences of the $m^{\text{th}}$ coefficient map along the $l^{\text{th}}$ dimension. For example, when the considered aperture (and therefore $\mathbf{x}_{m}$) is one-dimensional, it is $\Delta\mathbf{x}_{m}=[x_{m1}-x_{m2},x_{m2}-x_{m3}\cdots x_{mS}-x_{m1}]^{\sf T}$. Smooth coefficient maps $\hat{\mathbf{x}}_{m}$ enforce similarity between nearby and overlapping representations, which seems particularly suitable for sound fields.

To reconstruct a sound field via the convolutional model, we solve Eq. (10) via the alternating direction method of multipliers (ADMM)boyd2010admm and rewrite Eq. (10) in matrix form

where $\mathbf{x}\in\mathbb{C}^{MN}$ are the stacked columns of $\mathbf{X}$ and $\mathbf{D_{*}}\in\mathbb{C}^{N\times MN}$ is convolutional dictionary matrix such that $\mathbf{D_{*}}\mathbf{x}=\sum_{m}\mathbf{h}_{m}\circledast\overset{L}{\cdots}\circledast\mathbf{x}_{m}$ (e.g.  for $L=1$, $\mathbf{D_{*}}$ is block-circulant with block $\mathbf{H}_{s}$). $\mathbf{G_{*}}_{l}\mathbf{x}_{m}$ calculates the first order finite differences of the $m^{\text{th}}$ coefficient map along the $l^{\text{th}}$ dimension. In the one-dimensional case, $\mathbf{G_{*}}$ is a circulant matrix with the first row $[1,-1,\mathbf{0}_{1\times S-2}]$, where $\mathbf{0}_{1\times S-2}$ is a row vector of $S-2$ zeros.

To solve Eq. (11), we split the variables and reformulate the joint problem wohlberg2016a ; wohlberg2017a ; wohlberg2018a as

The ADMM steps in the $k^{\text{th}}$ iteration are

where $\mathbf{u}=[\mathbf{u_{0}},\mathbf{u_{1}}]^{\sf T}$ is the dual variable. The upper index $k$ indicates the state before the $k^{\text{th}}$ iteration (omitted further on for readability). In spatial frequency domain, the convolutional matrices reduce to their transformed filters and Eq. (13) reduces to heide2015a

which can be efficiently solved via the Sherman-Morrison formula wohlberg-2014-efficient ; wohlberg2016b and where $\tilde{\cdot}$ indicates frequency domain quantities. Also, $\mathbf{\tilde{G}}^{\sf H}\mathbf{\tilde{G}}=\sum_{l}\mathbf{\tilde{G}}_{l}^{\sf H}\mathbf{\tilde{G}}_{l}$, where $\mathbf{\tilde{G}}_{l}$ is the frequency transformed finite difference matrix in the $l^{\text{th}}$ dimension.
To align dimensions of $\mathbf{\tilde{D}}$ to the sound field $\mathbf{\tilde{p}}$, zero-padding is typically applied to the local filters $\mathbf{h}_{m}$ (corresponding to $[\mathbf{H}_{s}^{\sf T},\mathbf{0}_{M\times N-N_{s}}]^{\sf T}$, i.e. the first $M$ columns of $\mathbf{D}_{*}$). Instead, we obtain $\mathbf{\tilde{D}}$ from plane waves functions evaluated over the complete sound field (i.e. the global plane wave expansion $\mathbf{H}$).
Equations (14) and (15) are separable, such that

where the $\mathbf{y_{0}}$ update Eq. (17) has an efficient closed-form solution in frequency domain

$\mathbf{y_{1}}$ is updated by soft-thresholding, separable along the elements of $\mathbf{y_{1}}$,

where $\mathcal{S}$ is the shrinkage operator

To assess a reconstructed pressure field $\hat{\mathbf{p}}$, it is compared to the true field $\mathbf{p}$ in terms of the normalized mean square error NMSE and the spatial similarity $C$. The NMSE is

The spatial similarity $C$ is assessed as

such that $C=0$ indicates no similiarity and $C=1$ means the fields are indistinguishable.

The first experiment reconstructs the sound field by a monopole at the end of a linear microphone array, see Fig. 2(a). The sound pressure is simulated radially across ten wavelengths ($\lambda$), spanning the near-field and the far-field of the monopole. The linear array consists of 31 microphones at $0.5\,\lambda$ to $10.5\,\lambda$ radial distance and with spacing of $\lambda/3$. Three reconstruction methods are compared:

1. i)

   global plane waves with least squares (no regularization in Eq. (3));

2. ii)

   local independent plane waves using compressive sensing, i.e. finding sparse representations via Eq. (3) for each local partition separately, then overlap and average;

3. iii)

   convolutional sparse plane waves with smooth coefficients via Eq. (10).

All three reconstruction methods use the same set of plane waves, with wavenumber $\|\mathbf{k}\|_{2}=\tfrac{2\pi}{\lambda}$ and 21 incidence angles equally spaced along a semicircle $[0\cdots\pi]$. This experiment interpolates from observations with resultion $\lambda/3$ to a grid spacing of $\lambda/24$. The aperture of length $10\lambda$ contains $N=241$ reconstruction positions. The local approaches operate use subdomains of size one wavelength. In this example, each subdomain contains $N_{s}=25$ reconstruction points and at most 4 observations. Both sides of the domain are padded with $N_{s}-1$ zeros, to reduce artifacts from circular wrapping. The zero-padded domain has size $N^{\prime}=N+2(N_{s}-1)=289$ and is partitioned in $S=N^{\prime}$ fully overlapping subdomains. After reconstruction, the sound field is cropped again to the original size $N$. For comparison, the local approaches determine $MS=21\times 289$ coefficients compared to $M=21$ in the global model.

The reconstructions $\hat{\mathbf{p}}$ are shown in Fig. 2(b) and the error $\hat{\mathbf{p}}-\mathbf{p}$ in Fig. 2(c). All methods yield good reconstructions. The error is lowest for the smooth convolutional model, where smooth coefficients (plane wave magnitude and direction) among neighbouring representations are enforced. Methods using locally variant plane wave coefficients are flexible enough to approximate all measurements (the error is zero at measurement positions). Still, they rely on sufficient measurements within a local partition to yield good predictions. The global plane wave model can not conform to all measurements across the array, because of the mismatch between the radial decay of the sound field and propagating plane waves.

The experiment is repeated for varying distances between the monopole source and the microphone array. The normalized spatial mean squared error (24) is shown in Fig. 3. When the array is close to the monopole, local representations improve reconstructions due to the field’s high curvature and strong decay with distance. The proposed smooth-convolutional approach gives the most accurate reconstructions up to 0.7$\lambda$. The error decreases with distance for all methods, due to the less pronounced magnitude decay in the field. The further the array is placed in far field, the more the true field approaches planar characteristics, such that also the global model fits to the observations well and yields the best predictions. It is to note that in an application scenario, the reconstruction quality depends on many factors, such as the aperture size, curvature of wavefronts within the aperture, number and distribution of measurements available and not at least the evaluation criteria.

The approach is tested for a two-dimensional aperture of size $(5\,\lambda)^{2}$ in the plane $z=0$. The sound field is generated by interference of a monopole at (0,0,$\lambda$/8) with a plane wave propagating in $\mathbf{k}/k=(k_{x},k_{y},k_{z})/k=(0.38,-0.76,0.52)$. Four plane wave reconstructions are tested:

1. i)

   global plane waves with ridge regression ($\ell_{2}$-norm regularization in Eq. (3), parameter $\beta$ via leave-one-out cross-validation);

2. ii)

   local independent, sparse plane waves (solving Eq. (3) via least-angle regression, regularization via leave-one-out cross-validation);

3. iii)

   convolutional sparse plane waves without smooth coefficients ($\mu=0$ in Eq. (10));

4. iv)

   convolutional sparse plane waves with smooth coefficients via Eq. (10).

All models use plane waves with propagation angles distributed in a fibonacci grid ($k_{z}\geq 0$ hemisphere). The global method i) uses $M$ = 1000 propagation angles, the local methods ii-iv) $M$ = 100 angles. The reconstruction grid is regular with spacing of $d=\lambda/10$ between positions ($N=51^{2}$). The local subdomains and filters in methods ii)-iv) have size $\lambda^{2}$ (i.e. $N_{s}=(\lfloor\lambda/d\rfloor+1)^{2}=11^{2}$ discrete positions). For methods iii) and iv), the domain is zeropadded (with $\sqrt{N_{s}}-1$ in each direction) to avoid artifacts from circular convolutions ($N^{\prime}=(\sqrt{N}+2(\sqrt{N_{s}}-1))^{2}=3481$). The regularization parameters are tuned to $\beta=1\times 10^{-5}$, $\mu=1\times 10^{-3}$ (0 for iii) ), $\rho=1\times 10^{-5}$, and the ADMM iterations are stopped after 500 iterations. Note that this demonstration showcase exhibits a high signal to noise ratio. In less favourable conditions, the regularization would likely need to be adjusted.

The results are shown in Fig. 4. For reconstructions from 100 microphones, the global and the proposed approach with smooth local coefficients capture the spatial phenomena better than the two other approaches ii) and iii), which only rely on local and global sparsity. The spatially invariant (global) and slowly varying (proposed) models prescribe the necessary spatial structure to reconstruct from few measurements. Both methods exploit measurements across a larger spatial range and capture the global structure of the sound field, which is critical in sparsely sampled scenarios (see Fig. 4(b) and (e) vs. (c-d). When more microphones are available, spatially variant (local) models benefit from their flexibility to model complex sound fields. The global approach exhibits artefacts around the monopole due to the strong decay and curvature of the wavefronts in this region. The proposed approach balances local flexibility with global structure and yield good reconstructions in both cases. Note that method iii) and the proposed method iv) apply the same shrinkage threshold $\beta/\rho$ to all local coefficients (see Eq. (22)). Dynamic regularization could further improve the results, as it was observed for the independent local approach ii) (which uses cross-validation to find an optimal $\beta$ for each subdomain).

The particle velocity and sound intensity xy-vector fields of the reconstructions from 300 microphones (second row Fig. 4) are shown in Fig. 5. Global representations can not conform to the drastic spatial variations close to the monopole, all particle velocity vectors point outwards from $(x,y)=(0,0)$. Local approaches recover the fine structure of particle velocity and intensity also around the monopole.

The same methods i)-iv) from Sec. III.2 (and parameters) are used to reconstruct the enclosed reverberant sound field in a classroom (DTU building 352, Lyngby, Denmark), shown in Fig. 6. The room dimensions are $(l_{x},l_{y},l_{z})=(6.63,9.45,2.97)$ m, the reverberation time approx. $T_{60}=0.5$ s, the Schröder frequency $f_{S}\approx 240$ Hz. The room is furnished and its walls are somewhat irregular, with wooden floor, scattering elements on the walls and absorbing ceiling. A loudspeaker (BM6, Dynaudio, Skanderborg, Denmark) placed in a room corner was used to excite the room with 10 s logarithmic sweeps from 20 Hz to 20 kHz. A total of 4761 frequency responses were measured using a robotic arm (UR5, Universal Robots, Odense, Denmark) with a 1/2 inch free field condenser microphone (Brüel&Kjäer, Nærum, Denmark). The positions are distributed over a $1.7\times 1.7$ m² planar aperture with $N=69^{2}$ positions on a regular grid with 2.5 cm spacing. We refer the reader to hahmann2021spatial, for more information on the room and measurements and to hahmann2021b, for the dataset.

The sound field in the classroom is reconstructed at 1000 Hz from 98 and 295 measurements, distributed with uniform probability (and a minimum distance of 7 cm) across the aperture. Reconstructions using i-iv) and measured reference (“true”) of the classroom sound field are shown in Fig. 7 for the two cases with $N_{\mathrm{obs}}$ = 98 and 295. To reconstruct from few measurements, it is necessary to capture the global structure of the sound field, as in the global and the proposed approach (Fig. 7(b) and (d)). Models based on local representations conform easier to many measurements due to their higher number of coefficients (Fig. 7(i-k)). Specifically in this test, subdomains of size $\lambda^{2}$ contain $N_{s}=24^{2}$ discrete reconstruction positions, such that local models use a total of ii) $MN=476100$ and with padding iii+iv) $MN^{\prime}=M(\sqrt{N}+2(\sqrt{N_{s}}-1))^{2}=1322500$ coefficients compared to 1000 in the global model. The proposed approach combines both local flexibility and joint global analysis. As a consequence, it yields the highest similarity and lowest reconstruction errors when compared to the measured field.

The benefit of smooth coefficients shows when comparing the two convolutional approaches. Both seek a sparse approximation of the measurements, but spatial continuity is required to reconstruct sound fields successfully via local representations. Also the local independent approach yields smooth reconstructions, namely by averaging of overlapping partitions. However, a joint analysis of nearby representations is needed to align nearby coefficients and hence, indirectly exploit nearby measurements. In this study, the sparsity constraint enables feasible reconstructions also when only few measurements are available within a local subdomain and the inverse problem is severely underdetermined. As such, it is not the goal to represent the sound field using the fewest number of coefficients, but local sparsity is a means of exploiting the available measurements.

The experiment in Fig. 7 is extended to other frequencies and reconstruct the sound field in the classroom from 500 Hz to 2 kHz using a fixed number of microphones. The NMSE results in Fig. 8 show that the proposed approach yields good reconstructions when average distance between measurements is lower than $\lambda/2$. The proposed approach yields significant improvement with errors close to (see Fig. 8(a)), or lower than global plane waves (for sufficient measurements, see Fig. 8(b,c)).