This paper presents a numerical framework for designing diffuse fields in rooms of any shape and size, driven at arbitrary frequencies. That is, we aim at overcoming the Schroeder frequency limit for generating diffuse fields in an enclosed space. We formulate the problem as a Tikhonov regularized inverse problem and propose a low-rank approximation of the spatial correlation that results in significant computational gains. Our approximation is applicable to arbitrary sets of target points and allows us to produce an optimal design at a computational cost that grows only linearly with the (potentially large) number of target points. We demonstrate the feasibility of our approach through numerical examples where we approximate diffuse fields at frequencies well below the Schroeder limit.

## I. INTRODUCTION

^{1,2}depends on the chamber volume

*V*and absorption as per the relation

*c*is the phase speed, and

*T*

_{60}is the reverberation time required after source termination for the energy in the room to attenuate 60 dB. This expression provides an estimate of the frequency at which sufficient modal overlap first occurs. See Kuttruff

^{3}for a detailed derivation of (1).

Since $ f s 2$ is inversely proportional to the enclosure volume, the dimensions of the chamber determine the lowest frequency at which a diffuse field would naturally occur. Although there are means to improve or fully develop a diffuse field within a given chamber, such as splayed walls, rotating panels, and moving vanes,^{4} there are currently no known means to induce a spatially broad diffuse field in the room at frequencies below *f _{s}* where both the first and second order statistics of the field are satisfied.

Traditionally, random noise sources, such as large air horns, were used to generate the high amplitude sound pressure levels needed for testing within reverberation chambers. Recently, these sources have been augmented or replaced with concert-grade loudspeakers, allowing improved closed-loop control and the potential to achieve a wider range of test spectra.^{5,6} The addition of loudspeaker sources provides an opportunity for optimizing the acoustic field within a chamber, as will be shown in this paper. Synthesizing diffuse fields below the Schroeder frequency can lead to significant financial and time savings. For instance, physical dimensions of new chambers could be relaxed, and shipping of test articles to facilities with sufficiently low *f _{s}* could be avoided.

The design of diffuse acoustic fields has been explored, to some extent, in the open literature. Specific examples include the work of Maury and Bravo,^{7} who designed various types of acoustic random fields on the surface of a baffled thin plate, including diffuse fields, using a near-field 4 × 4 loudspeaker array. To that end, they solved a quadratic optimization problem in which they found optimal speaker drive signals that minimized the misfit between the predicted field and the spatial correlation of the target random field near the surface of the test body. In recent work, Alvarez Blanco *et al.*^{8} presented an approach to design controls for an array of loudspeakers to produce a spatially uniform sound pressure level during direct field acoustic tests. For this, they also used a pseudo-inverse strategy to obtain solutions of a quadratic optimization problem that directly provided the signal inputs. Although coherence and phase were not constrained, the resulting uniform pressure field could potentially be diffuse but was not guaranteed. For other relevant, recent work in this area, see Refs. 9–11.

As in the aforementioned work, we investigate the synthesis of diffuse fields at arbitrary frequencies (specifically below the Schroeder limit), although in contrast, the field is generated broadly within the room and without requiring that the acoustic sources be located in close proximity to the testing region. Hence, our main contribution is the development of a numerical approach for the design of acoustic diffuse fields in enclosed rooms. To that end, we put forward a partial differential equation (PDE)-constrained optimization formulation for the control problem. Although existing work has addressed the multiple input multiple output (MIMO) control problem,^{7,8} we provide further insight into Tikhonov regularization for these problems. Moreover, one of the numerical challenges in designing or producing samples of a diffuse field in simulations is the need to capture its spatial correlation. The latter may require fine discretizations (e.g., using finite elements) that lead to large dense matrices.

The computational challenge of factorizing the spatial correlation (in the context of diffuse fields) was recently addressed by Van Hoorickx and Reynders,^{12} who proposed to generate diffuse fields given by expansions on the eigenfunctions of the pressure correlation operator (i.e., using the Karhunen–Loève decomposition of the latter). When all target points lie on a line, the univariate correlation kernel has analytically known eigenfunctions (given by spheroidal wave functions) and, hence, may be accurately approximated in a separated-variable manner. For multidimensional sets of target points, the computational bottleneck of factorizing the dense correlation was avoided by replacing its kernel with the first term of its expansion given by the addition theorem [see Ref. 13, Eq. (10.60.2)]. Upon accepting the resulting truncation error for vectors not significantly exceeding the operating wavelength, this approximation recast the correlation as a product of univariate kernels, whose known eigenfunctions were then used.

Our approach avoids the above low-order series truncation by exploiting a more-accurate low-rank approximation of the correlation function, applicable to arbitrary sets of target points. This allows us to solve the least squares problem, producing an optimal design at a computational cost that grows only linearly with the (potentially large) number of target points (instead of quadratically if using directly the dense correlation matrix). In addition, this proposed treatment allows us to formulate simple and efficient regularized versions of the optimization problem, which cater to it being possibly ill-conditioned and remain computationally economical.

The rest of the paper is organized as follows. We provide a summary of the plane wave model for diffuse fields and the ensuing statistics. We then formulate the forward problem for modeling a random acoustic field in an enclosed room. Next, we provide the optimization formulation for the control problem and the description of our approach for the low-rank approximation of the spatial correlation. We also offer two Tikhonov regularization strategies for the control problem. Next, we demonstrate the feasibility of our proposed strategy using numerical examples of an enclosed room driven at frequencies below the Schroeder limit. Finally, we provide some conclusions and future directions.

## II. BACKGROUND

Here, we summarize existing results on the theoretical modeling of diffuse fields. We will adopt the plane wave model in which a random pressure field is conceived as the interaction of an infinite number of plane waves having randomized direction of propagation and phase.^{14} It is well known that this model leads to a Gaussian spatiotemporal pressure field, which is fully characterized by its mean and correlation function. In this work, for the sake of simplicity and without loss of generality, we will concentrate on pure-tone time-harmonic fields.

### A. Plane wave model

^{15}we model a pure-tone diffuse field as a random pressure field expressed as

*c*is the phase speed, $ D n$ are independent random vectors (uniformly distributed over the unit sphere) describing the direction of a plane wave,

*A*are independent random variables, and $ \Phi n$ are random phases. The vector

_{n}*θ*represents the collection of all the random variables in the model. Notice that, as a consequence of the central limit theorem, the random field $ P ( x , \omega ; \theta )$ is Gaussian.

^{15}Also, without loss of generality, we use in the sequel a constant amplitude $ A n : = p 0$.

### B. First and second order statistics of a diffuse field

^{16}

As previously noted, the random field $ P ( x , \omega ; \theta )$ is Gaussian. Hence, we can completely characterize a diffuse field through its mean and correlation. Moreover, the field is weakly isotropic (i.e., the mean is constant, and the correlation depends only on the distance between two points).^{17} A more extensive discussion on diffuse fields can be found in Refs. 15, 16, 18, and 19.

## III. DESIGN OF DIFFUSE FIELDS

As stated before, we seek to develop a numerical framework for designing approximately diffuse fields for rooms of any shape and size at arbitrary frequencies. To this end, we first introduce the acoustic equations and numerical approximations, followed by the design problem, and a formal regularized treatment of the ensuing ill-posed design problem.

As we will see later in this section, one of the key challenges in the design problem is the computational expense that arises from the discretization of the target correlation function. One of our contributions in this work is an efficient low-rank representation of the target correlation that renders the design problem tractable.

### A. Forward problem

*d*is the number of sources in the room and $ \Gamma = \Gamma r \u222a \Gamma N$. We assume in this work that all sources are located on the room walls; however, this is not a limitation of the method. The pressure field in the room satisfies

^{20}

*ω*is the angular frequency,

*ρ*

_{0}is the fluid density, $ \kappa = \omega / c$ is the wavenumber,

*s*is the amplitude of the normal component of acceleration over $ \Gamma N j$,

_{j}*Z*is the specific acoustic impedance, and $ n j$ is a unit vector normal to $ \Gamma j$ and directed out of the room. Using the finite element method

^{21}to obtain a discrete representation of (5), we arrive at

*N*is the number of finite element nodes, $ M \u2208 R N \xd7 N$ is the mathematical mass matrix, $ C \u2208 R N \xd7 N$ is the damping matrix, $ p \u2208 C N$ is the nodal pressure vector, and $ F \u2208 R N \xd7 d$ is a matrix that acts on the acceleration vector $ s \u2208 C d$, projecting it onto the finite element basis functions.

*m*target locations where we want the field to be diffuse. Now, let $ B \u2208 R m \xd7 N$ be the Boolean matrix such that $ p T = B p \u2208 C m$ collects the nodal pressures at the target locations. We then have

*R*is invertible for any frequency.

^{22}

*G*and the correlation of the sources, denoted as

*S*, at a given frequency is given as

*h*denotes complex conjugation and transposition. So, given

*S*, we can compute

*G*at given locations. Next, we develop a design methodology for diffuse fields based on the forward model presented in this section.

### B. Inverse or design problem

*G*(

*S*) solves (11) and $ \Vert G \Vert F$ is the Frobenius matrix norm, which is associated with an inner product: $ \Vert G \Vert F 2 = ( G , G ) F = : tr ( G h G )$. Then the optimal correlation,

*S*, can be obtained as

_{o}*a priori*be restricted to positive definite and Hermitian matrices

*S*, it turns out that the minimum-norm solution found by the unconstrained minimization (14) automatically satisfies those requirements (see Remark 3).

### C. Optimality condition and minimum-norm least squares solution

*S*is obtained as

_{o}*T*does not have full column rank

*d*), the above equation fails to provide a solution $ S o \u2208 C d \xd7 d$ that is invertible, let alone positive definite. The transfer matrix

*T*is therefore assumed henceforth to have rank

*d*.

*T*,

*d*nonzero singular values $ \sigma 1 \u2265 \sigma 2 \u2265 \cdots \sigma d \u2265 0$ of

*T*, while $ X \u2208 C m \xd7 d$ and $ Y \u2208 C d \xd7 d$ hold the

*d*left and right associated singular vectors (arranged columnwise), respectively. In particular, the matrices

*X*and

*Y*have the orthonormality properties $ X h X = I d$ and $ Y h Y = Y Y h = I d$, with

*I*the

_{d}*d*×

*d*identity matrix. On introducing the SVD (22) in (20) and using the latter properties of

*X*and

*Y*, the minimum-norm least squares solution

*S*is found as

_{o}*S*is in fact positive definite.

_{o}The expression (23) of *S _{o}* entails the evaluation of $ X h G \u0302 X$, whose $ O ( d m 2 + d 2 m / 2 )$ cost constitutes a potential computational bottleneck as $ G \u0302$ is a dense

*m*×

*m*matrix that may be large in realistic problems [e.g., $ m = O ( 10 4 )$ to $ O ( 10 6 )$]. To reduce the computational complexity in

*m*of the design solution method, we now introduce a low-rank approximation of $ G \u0302$.

### D. Low-rank approximation of $ G \u0302$

*P*terms, yields the best rank-

*P*approximation of $ G \u0302$, denoted as $ G \u0302 P$, in the sense of the Frobenius norm. Specifically, the relative truncation error is given as

The rate of decay of the eigenvalues of a correlation matrix depends on the correlation length. The latter is usually high, which allows a truncation order $ Q \u226a m$. However, setting up this approximation still entails obtaining a large enough number of eigenpairs of the *m* × *m* matrix $ G \u0302$ to achieve and verify a sufficiently low truncation error, and this remains often impractical.

*j*

_{0}is the spherical Bessel function of first kind and order zero, and $ r i j : = x i \u2212 x j$ is the position vector joining two generic target points. Now, for any $ z \u2208 R 3$, the function

*j*

_{0}admits the integral representation

*Q*-point quadrature rule with nodes $ \theta \u0302 q \u2208 S \u0302$ and positive weights

*w*( $ 1 \u2264 q \u2264 Q$), yielding

_{q}*Q*of the quadrature rule ensuring a desired (small enough) quadrature error depends on the oscillatory character of the integral (26) and, hence, on the magnitude of the argument $ | z |$ of

*j*

_{0}there. In this study, we have $ | z | \u2264 \kappa L$ with

*L*denoting the largest spatial separation between target points. The quadrature rule used in (27) is the product (in spherical coordinates) of a 2

*R*-point trapezoidal rule and an

*R*-point Gauss–Legendre rule (so that $ Q = 2 R 2$) and is known (see Ref. 23, Theorem 5.4) to integrate exactly all spherical harmonics $ Y l m$ of degree $ l \u2264 2 R \u2212 1$. Expressing $ j 0 ( | z | )$ in (26) using the Jacobi–Anger expansion [see Ref. 13, Eq. (10.60.7)], we can estimate $ \epsilon Q$ as $ | \epsilon Q | \u2264 C \u2211 l \u2265 2 R \u2211 m = \u2212 l l j l ( | z | ) \epsilon Q , l m$, with $ \epsilon Q , l m$ denoting the quadrature error for $ \u222b S \u0302 Y l m ( \theta \u0302 ) d S ( \theta \u0302 )$. By the large-order asymptotics of spherical Bessel functions [Ref. 13, Eq. (10.19.1) used in (10.47.3)], the coefficients $ j l ( | z | )$ are known to decay superexponentially fast with

*l*if $ l > e | z | / 2$, which suggests in the present context to set $ R \u2248 e \kappa L / 4$, i.e., $ Q \u2248 ( e \kappa L ) 2 / 8$. In particular,

*Q*does not depend on the number

*m*of target points if their maximum spatial separation

*L*is fixed. In the forthcoming examples, $ Q = O ( 10 2 )$, whereas $ m = O ( 10 4 )$, so that (29) accomplishes a low-rank approximation of $ G \u0302$, whose computation is moreover economical as the

*m*-vectors $ \varphi q$ are given explicitly.

*S*. Indeed, using (29) in (23) gives

_{o}*O*(

*m*) computational work and memory, down from $ O ( m 2 )$ if using the full matrix $ G \u0302$.

**Remark 1**. *Unlike in the expansion* $ G \u0302 = \u2211 q = 1 m \psi q \psi q T \lambda q$ *in terms of eigenpairs, the vectors* $ \varphi q$ *in* *(29)* *are not orthogonal*.

### E. Regularized least squares solution

*T*may be ill-conditioned (i.e., have a large condition number), with

^{hT}*T*still assumed to have full column rank. Let

*T*) the unique solution of the least squares problem

*S*given by (32), we have to solve at most

_{o}*Q*problems of the type (33). To cater for

*T*being possibly ill-conditioned, we add a regularization term in (33) to get the problem

^{h}T*S*are both given by finite sums of tensor products].

_{o}#### 1. Second regularization method

*S*is equal to the sum of the mean square acceleration of the sources, $ 1 / 2 \Vert S \Vert F = 1 / 2 tr ( S S H ) = 1 / 2 \u2211 i d | s i | 2$, and would be proportional to the total electric power supplied to the sources. The stationarity condition of (40) is

*T*, this becomes

*α*= 0 yields $ S o = Y H ( 0 ) Y h$ through (23) and that we have

**Remark 2**. *The first regularization yields* $ H i j ( \alpha ) = ( Z Z h ) i j / ( \sigma i + \alpha ) ( \sigma j + \alpha )$ *instead of* *(43)*, *with H_{ij} as in*

*(42)*.

*This shows that the two regularization approaches are not identical, although Remarks 3 and 4 show that they are similar in several ways*.

**Remark 3**.

*The minimum-norm solution*

*S*to problem_{o}*(14)*,

*as well as its regularized approximations*$ S o 1 ( \alpha )$

*and*$ S o 2 ( \alpha )$,

*are Hermitian and positive definite (and, hence, acceptable as correlation matrices), without those restrictions needing to be explicitly enforced (e.g., through constraints). In fact, any matrix $ S \u2208 C d \xd7 d$ can be additively decomposed into its Hermitian and skew-Hermitian parts, $ S = S 1 + S 2$ with $ S 1 h = S 1$ and $ S 2 h = \u2212 S 2$, and we have $ \Vert S \Vert F 2 = \Vert S 1 \Vert F 2 + \Vert S 2 \Vert F 2$. Moreover, it is easy to verify that $ G 1 : = A ( S 1 )$ and $ G 2 : = A ( S 2 )$ are, respectively, Hermitian and skew-Hermitian. For the objective functional*

*J*(*S*), this gives (since $ G \u0302$ is Hermitian)*T*has full (column) rank. A similar line of reasoning applies to the regularized versions of problem (14).

**Remark 4**. *From a computational complexity standpoint, both regularizations require $ O ( m d 2 ) + O ( Q d 2 ) + O ( d 3 )$ complex arithmetic operations, with $ m \u226b Q \u2265 d$ in the present context. The leading $ C \xd7 m d 2$ amount of arithmetic operations (where C is a method-dependent constant) arises from the decomposition of the transfer matrix T. Both regularization methods may (as explained) use the reduced SVD of T, in which case C = 6 (Ref. 24, Sec. 8.6). Alternatively, the first regularization may as easily be carried out using a thin QR factorization of $ [ T ; \alpha I ]$, resulting in C = 2 (Ref. 24, Sec. 5.2)*.

## IV. NUMERICAL RESULTS

### A. Problem description

In this section, we demonstrate how we can construct approximately diffuse fields in enclosed rooms at arbitrary frequencies (e.g., below the Schroeder frequency). That is, by exciting the room with signals drawn from a multivariate Gaussian vector with zero mean and an optimal source correlation, we can obtain a random field in the target region whose spatial correlation is close to (3) in the sense described by (14).

We consider a room with dimensions $ 5.75 \xd7 4.25 \xd7 3 m 3$. The specific acoustic impedance of the wall, ceiling, and floor is taken constant for all simulations and set to $ Z = 4.25 \xd7 10 5 kg / ( m 2 s$), while the speed of sound and mass density are $ c = 340 m / s$ and $ \rho 0 = 1 kg / m 3$, respectively. The side walls of the room contain uniformly spaced acoustic sources, each with an area of $ 0.25 \xd7 0.25 m 2$. The sources may be understood as each representative of an elastically rigid square planar loudspeaker diaphragm flush mounted in the room wall. We consider three different speaker configurations in this study: (1) 9 sources/wall (*d* = 36), (2) 16 sources/wall (*d* = 64), and (3) 25 sources/wall (*d* = 100). For each case, we compute the optimal correlation of the acceleration of the sources (i.e., loudspeaker diaphragms) over the frequency range $ [ 150 , 300 ]$ Hz in 10 Hz intervals and report detailed results for three frequencies: 150, 250, and 300 Hz. We chose a limited number of frequencies to demonstrate our approach since solutions at various frequencies are independent of each other. If a wider response band is needed, thanks to this independence, multiple runs can be conducted in a trivially parallel fashion, eliminating the need for communication and resulting in substantial computational efficiencies. A representative geometry for the room with 36 sources is shown in Fig. 1. The target region, also shown, has dimensions $ 2 \xd7 2 \xd7 2 m 3$ and is located in the center of the room.

For each of the studied cases, we first compute a low-rank approximation of the target correlation $ G \u0302$ as per (29), using $ Q = 2 R 2 = 200$. With reference to the discussion following Eq. (29), we have $ L \u2264 2 3$ m and (for the highest frequency) $ \kappa \u2248 5.54$ m^{−1}; hence, $ | z | \u2264 19.2$ in (26), which would result in *R* = 13 by our proposed adjustment rule if catering for the worst-case spatial target point separation. In practice, the used value *R* = 10 is consistent with the separation of most pairs of target points and results in a relative approximation error of less than 0.1% on $ G \u0302$ for all cases studied herein. For simplicity of implementation, the value of *Q* was not adjusted downward for the two lower frequencies used. We then solve the least squares problems (34) for $ 1 \u2264 q \u2264 Q$. Finally, the optimal source correlation *S _{o}* is obtained as per (36). The Tikhonov regularization parameter

*α*is determined using an L-curve approach.

^{25,26}To this end, we use the objective and regularization term in (40) for either of the regularization strategies described in Sec. III E. We point out that the computational cost of solving all the aforementioned optimization problems is negligible when compared to the computational cost of building the transfer matrix

*T*.

We use an in-house finite element code developed using the FEniCS library^{27} in conjunction with the parallel direct solver MUMPS for all the calculations presented herein. The models are meshed with four-node tetrahedral elements. All the results shown are generated with a mesh containing approximately 130 000 nodes and 750 000 elements, which is fine enough to achieve a low discretization error in all calculations. The number of points in the target region for all simulations was $ m \u2248 15 \u2009 000$.

^{28}

*α*is the random-incidence absorption coefficient

_{rand}^{20}

### B. Results

^{2}in all the examples. Hence, a useful metric to explore is how much the diagonal entries of $ G o = A S o$, the mean square pressure at the target locations, depart from unity. To this end, we define the error

^{15,16,18}and proved very useful for assessing the quality of our numerical solutions.

Table I contains a summary of the errors for all source configurations along with the Tikhonov parameter $ ( \alpha )$ for the three frequencies 150, 250, and 300 Hz studied in detail. We notice that for these frequencies, both errors, *ϵ _{G}* and

*ϵ*, decrease as the number of sources increases, as expected. Moreover, we investigated the behavior of these errors over the frequency range $ [ 150 , 300 ]$ Hz. The results are shown in Fig. 2. We see that as the number of sources increases, both errors decrease for all frequencies. It will be shown below that all errors are related to satisfactory approximations of a diffuse field when we explore in more depth at three frequencies of interest.

_{ms}Sources . | 150 Hz . | 250 Hz . | 300 Hz . | ||||||
---|---|---|---|---|---|---|---|---|---|

ϵ
. _{ms} | ϵ
. _{G} | α
. | ϵ
. _{ms} | ϵ
. _{G} | α
. | ϵ
. _{ms} | ϵ
. _{G} | α
. | |

36 | 0.17 | 0.33 | $ 2 \xd7 10 \u2212 3$ | 0.30 | 0.49 | $ 5 \xd7 10 \u2212 3$ | 0.43 | 0.64 | $ 2 \xd7 10 \u2212 3$ |

64 | 0.09 | 0.17 | $ 1 \xd7 10 \u2212 3$ | 0.10 | 0.26 | $ 5 \xd7 10 \u2212 3$ | 0.18 | 0.37 | $ 1 \xd7 10 \u2212 3$ |

100 | 0.07 | 0.15 | $ 1 \xd7 10 \u2212 4$ | 0.08 | 0.23 | $ 1 \xd7 10 \u2212 4$ | 0.08 | 0.24 | $ 1 \xd7 10 \u2212 4$ |

Sources . | 150 Hz . | 250 Hz . | 300 Hz . | ||||||
---|---|---|---|---|---|---|---|---|---|

ϵ
. _{ms} | ϵ
. _{G} | α
. | ϵ
. _{ms} | ϵ
. _{G} | α
. | ϵ
. _{ms} | ϵ
. _{G} | α
. | |

36 | 0.17 | 0.33 | $ 2 \xd7 10 \u2212 3$ | 0.30 | 0.49 | $ 5 \xd7 10 \u2212 3$ | 0.43 | 0.64 | $ 2 \xd7 10 \u2212 3$ |

64 | 0.09 | 0.17 | $ 1 \xd7 10 \u2212 3$ | 0.10 | 0.26 | $ 5 \xd7 10 \u2212 3$ | 0.18 | 0.37 | $ 1 \xd7 10 \u2212 3$ |

100 | 0.07 | 0.15 | $ 1 \xd7 10 \u2212 4$ | 0.08 | 0.23 | $ 1 \xd7 10 \u2212 4$ | 0.08 | 0.24 | $ 1 \xd7 10 \u2212 4$ |

Interestingly, the top two plots in Fig. 2 show that the errors increase with frequency in a non-monotonic fashion. This trend could be explained, at least qualitatively, through the sources of error involved in the computation of *G _{o}*. There are three main sources of error in our numerical approximation of a diffuse field in a room below the Schroeder limit, namely, the number

*d*and location of the sources, the properties of the room (both shape and wall impedance), and the numerical discretization error in the finite element approximation of the governing equations. From our mesh convergence studies, we concluded we could ignore the discretization errors and focus on the sources and the room.

*d*due to the rank of

*T*. Using the best approximation error $ G \u0302 d$ in Eq. (24), we can obtain a lower bound $ E ( d )$ for the relative residual of our inverse problem,

*d*. Comparing this plot to the plots of

*ϵ*and

_{G}*ϵ*in the top of the figure, we see that the lower bound is well below these errors. This behavior indicates the discrepancy in error is likely due to how well the spatial distribution of the sources with an optimized correlation can excite room modes to produce fields that are nearly diffuse.

_{ms}Another interesting feature we observe in Fig. 2 is that there are frequencies where the errors have local minima. For example, at 240 Hz, the 64-source configuration produces errors almost as low as the 100-source case. This trend could be explained as an excitation of relatively more favorable room modes for the 64-source configuration near this frequency, which in turn would lead to a very interesting conclusion: There exist optimal locations for the sources in terms of diffuse field creation. This suggests the addition of a computational strategy involving source location optimization within our approach could increase the efficiency and quality of diffuse field creation. Such addition lies outside the scope of the current work and will be investigated in a future paper. Moreover, the foregoing explanations require careful investigation, which we also plan to do as part of our future work.

To better illustrate the quality of our solutions, we plot the spatial correlation field with respect to the center of the target region for different combinations of sources and frequencies as shown in Fig. 3. We observe that the optimized correlation closely resembles that of a diffuse field. Furthermore, we compare the target and estimated fields along a diagonal across the target region in Fig. 4. We again observe a close match of the computed and target correlations for different frequency and source combinations. It is important to point out that we are showing only the target region at the center of the room in these images. As we will show later, the acoustic field departs from being purely diffuse at locations near the walls, as expected.

Recall that the mean square pressure (i.e., diagonal entries of the spatial correlation matrix) is spatially constant in a diffuse field (as captured by the metric *ϵ _{ms}* shown in Table I). We now provide more global representations of this behavior. Figure 5(a) shows the mean square pressure along a line through the center of the room in the

*x*direction for all cases. We can see that the field is close to unity and constant in a region slightly larger than the target domain and departs from pure diffusion close to the walls. Furthermore, for comparison purposes, we show in Fig. 5(b) the mean square pressure obtained using random realizations from uncorrelated sources assuming a standard normal Gaussian distribution. Notice that indeed uncorrelated sources cannot produce the desired constant field at this given frequency, reinforcing the success of the proposed optimization approach.

Last, Fig. 6 shows the mean square pressure field for the three studied cases. Again, we notice in all cases the field is constant over a region larger than the target one but departs from diffuse behavior away from the target region, as expected.

## V. SUMMARY AND CONCLUSIONS

We presented how we can use an optimization approach to design sources that yield diffuse fields in enclosed rooms at arbitrary frequencies (even below the Schroeder limit). To this end, we characterized a diffuse field solely by its mean and correlation as these are stochastic homogeneous and isotropic Gaussian processes. Then we postulated an optimization problem in which we sought the correlation structure of input sources that minimized the distance between the output correlation and that of a diffuse field over a target region. We addressed the large computational cost that arises from the discretization of the target spatial correlation using a low-rank expansion based on the integral representation of spherical Bessel functions. Moreover, we formally showed how to regularize the ensuing ill-posed inverse problem.

Our results demonstrated that it is possible to obtain approximate diffuse fields in enclosed rooms even at frequencies below the Schroeder limit by driving correlated sources in an optimal way. Also, we found that there is a limitation in the quality of the approximation that depends on the number of sources and the frequencies of interest. As frequency increases, a larger number of sources seems to be needed to maintain a given level of error in the diffuse field approximation. However, we also observed that lower errors could be obtained by placing sources in an optimal way.

A direction for future work is to study the influence of source location on the approximation error. Furthermore, it is possible to devise optimization algorithms to find such locations. Also, incorporating uncertainty in boundary conditions, material properties, etc. in the optimization formulation would be highly desirable. Finally, experimental validation of the method will be performed in future work.

## ACKNOWLEDGMENTS

This work was supported by the Laboratory Directed Research and Development program at Sandia National Laboratories, a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy's National Nuclear Security Administration under Contract No. DE-NA-0003525.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

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

## REFERENCES

*The NIST Handbook of Mathematical Functions*

*The Diffuse Sound Field: Statistical Considerations Concerning the Reverberant Field in the Steady State*

*Stochastic Systems: Uncertainty Quantification and Propagation*

*Theoretical Acoustics*

*Finite Element Analysis of Acoustic Scattering*

*Spherical Harmonics and Approximations on the Unit Sphere: An Introduction*

*Matrix Computations*

*Discrete Inverse Problems: Insight and Algorithms*