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.

Many aerospace structures, satellites, and internal components experience acoustic loads that are diffuse in nature. That is, these loads are composed of a large number of waves having random amplitude, propagation direction, and phase. Ground-based qualification testing of these structures is typically performed in acoustic reverberation chambers within which a diffuse field arises naturally. The minimum frequency beyond which an acoustic field in an enclosure is (naturally) diffuse, termed the Schroeder frequency,1,2 depends on the chamber volume V and absorption as per the relation
f s = c 3   T 60 4 ln ( 10 ) V .
(1)
Here, c is the phase speed, and T60 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 Kuttruff3 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 fs 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 fs 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.

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.

Following Jacobsen,15 we model a pure-tone diffuse field as a random pressure field expressed as
P ( x , ω ; θ ) = lim N 1 N n = 1 N A n e i ( κ D n · x + Φ n ) ,
(2)
where κ = ω / c is the wavenumber, c is the phase speed, D n are independent random vectors (uniformly distributed over the unit sphere) describing the direction of a plane wave, An are independent random variables, and Φ n are random phases. The vector θ 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 , ω ; θ ) is Gaussian.15 Also, without loss of generality, we use in the sequel a constant amplitude A n : = p 0.
Here, we summarize well-known results on the statistics of the stochastic field shown in (2). It is straightforward to show that
E [ P ( x , ω ; θ ) ] = 0 x , ω ,
where E [ · ] denotes expectation.
Let r q s = x q x s be the distance between any two points x q and x s. Then the spatial correlation can be shown to be given by16 
G ( r q s , ω ) = E [ P ( x q , ω ; θ ) P ¯ ( x s , ω ; θ ) ] = p 0 2 2 sin ( κ ( ω ) r q s ) κ ( ω ) r q s ,
(3)
where the overbar denotes complex conjugation.

As previously noted, the random field P ( x , ω ; θ ) 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.

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.

We model a room as a bounded domain Ω R 3 with boundary (walls) Γ. The region of the walls having a prescribed acoustic impedance is denoted as Γ r, while that occupied by the sources is denoted as Γ N = j Γ N j , j = 1 , , d, where d is the number of sources in the room and Γ = Γ r Γ 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 satisfies20 
2 p + κ 2 p = 0 in   Ω ,
(4)
p · n + i Z 1 ω ρ 0 p = 0 on   Γ r ,
(5)
p · n j ρ 0 s j = 0 on   Γ N j for   j = 1 , , d ,
(6)
where ω is the angular frequency, ρ0 is the fluid density, κ = ω / c is the wavenumber, sj is the amplitude of the normal component of acceleration over Γ N j, Z is the specific acoustic impedance, and n j is a unit vector normal to Γ j and directed out of the room. Using the finite element method21 to obtain a discrete representation of (5), we arrive at
R p = F s,
(7)
where R : = K κ 2 M + i ω C, K R N × N is the acoustic stiffness, N is the number of finite element nodes, M R N × N is the mathematical mass matrix, C R N × N is the damping matrix, p C N is the nodal pressure vector, and F R N × d is a matrix that acts on the acceleration vector s C d, projecting it onto the finite element basis functions.
Let T : = { x 1 , , x m } be a finite set of m target locations where we want the field to be diffuse. Now, let B R m × N be the Boolean matrix such that p T = B p C m collects the nodal pressures at the target locations. We then have
p T = B p
(8)
= B R 1 F s
(9)
= T s ,
(10)
where we have introduced the transfer matrix T : = B R 1 F C m × d. Also, we point out that if Z 0, R is invertible for any frequency.22 
For random sources, the relationship between the spatial correlation of the target pressures denoted as G and the correlation of the sources, denoted as S, at a given frequency is given as
G = E [ p T p T h ] = T E [ s s h ] T h = T S T h ,
(11)
where the superscript 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.
Our design problem can be described as follows: Given a spatial correlation of the pressure at a set T of target nodes, determine the correlation of the acceleration of the input sources. To this end, let x i , x j T be two target locations. Then, from (3), the component of the target spatial correlation matrix at a given frequency, G ̂ i j, is given as
G ̂ i j = p o 2 2 sin ( κ r i j ) κ r i j .
(12)
We first introduce an un-regularized inverse problem for the sake of simplicity. This formulation leads to a decomposition of the inverse problem that strongly reduces computational cost, as will be shown. After this, a regularized version of the problem follows naturally. Define an objective function as
J ( S ) : = 1 2 G ( S ) G ̂ F 2 ,
(13)
where G(S) solves (11) and G F is the Frobenius matrix norm, which is associated with an inner product: G F 2 = ( G , G ) F = : tr ( G h G ). Then the optimal correlation, So, can be obtained as
S o = arg min S C d × d J ( S ) .
(14)
While the search space for the above minimization should 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).
To simplify our derivations, define a linear operator A : C d × d C m × m as
A S : = T S T h .
(15)
Substituting this expression into the objective (13), we get
J ( S ) : = 1 2 A S G ̂ F 2 .
(16)
The first-order optimality condition for problem (14) is that the gradient of the objective be zero at the minimizer. Using the inner product associated with the Frobenius norm, the directional derivative of the objective at S C d × d in an arbitrary direction H C d × d is obtained as
J ( S ) , H = d d ϵ J ( S + ϵ H ) | ϵ = 0
(17)
= Re ( A S G ̂ , A H ) F
(18)
= Re ( A A S A G ̂ , H ) F ,
(19)
where A denotes the adjoint of A, defined by ( W , A V ) F = ( A W , V ) F for any W C m × m and V C d × d. Then, from the above expression, the first-order optimality condition J ( S o ) = 0 verified by So is obtained as
A A S o A G ̂ = 0 ,
(20)
which is in essence the normal equation for the least squares problem (14). Then, using (15) in (20) and simplifying, we get
T h T S o T h T = T h G ̂ T .
(21)
If T h T is not invertible (i.e., if T does not have full column rank d), the above equation fails to provide a solution S o C d × d that is invertible, let alone positive definite. The transfer matrix T is therefore assumed henceforth to have rank d.
To solve Eq. (20), and also to later address regularized versions of problem (14), it is convenient and computationally reasonable (under the present operating conditions) to introduce and use the reduced singular value decomposition (SVD) of T,
T = X Σ Y h ,
(22)
where Σ R d × d is a diagonal matrix holding the d nonzero singular values σ 1 σ 2 σ d 0 of T, while X C m × d and Y C d × 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 Id the 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 So is found as
S o = Y Σ 1 ( X h G ̂ X ) Σ 1 Y h
(23)
and is clearly Hermitian and positive definite. Moreover, for any s C d, we have s h S o s = ( X Σ 1 Y h s ) h   G ̂   ( X Σ 1 Y h s ) > 0 by virtue of the positive definiteness of G ̂, showing that So is in fact positive definite.

The expression (23) of So entails the evaluation of X h G ̂ X, whose O ( d m 2 + d 2 m / 2 ) cost constitutes a potential computational bottleneck as G ̂ 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 ̂.

The target spatial correlation G ̂ being symmetric and positive definite, we have G ̂ = k = 1 m ψ k ψ k h λ k, where ψ k , λ k are the eigenpairs of G ̂ numbered so that λ 1 λ 2 λ m > 0. By the Eckart–Young theorem, this expansion, truncated to its first P terms, yields the best rank-P approximation of G ̂, denoted as G ̂ P, in the sense of the Frobenius norm. Specifically, the relative truncation error is given as
E 2 ( P ) : = G ̂ G ̂ P F 2 G ̂ F 2 = k = P + 1 m λ k 2 j = 1 m λ j 2 .
(24)

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 m. However, setting up this approximation still entails obtaining a large enough number of eigenpairs of the m × m matrix G ̂ to achieve and verify a sufficiently low truncation error, and this remains often impractical.

We therefore propose an alternative strategy for deriving low-rank approximations of G ̂. It is based on observing that the generic entry G ̂ i j of G ̂ [see (12)] is, in fact, equivalently given as
G ̂ i j = p o 2 2 j 0 ( κ | r i j | ) ,
(25)
where j0 is the spherical Bessel function of first kind and order zero, and r i j : = x i x j is the position vector joining two generic target points. Now, for any z R 3, the function j0 admits the integral representation
j 0 ( | z | ) = 1 4 π S ̂ e i z · θ ̂   d S ( θ ̂ ) ,
(26)
where S ̂ is the unit sphere (spanned by unit vectors θ ̂). For instance, expressing the above integral using spherical angular coordinates reduces it to the one-dimensional integral representation formula (10.54.1) given in Ref. 13. Let the above integral be approximated by a Q-point quadrature rule with nodes θ ̂ q S ̂ and positive weights wq ( 1 q Q), yielding
j 0 ( | z | ) = ( 1 4 π q = 1 Q w q e i z · θ ̂ q   ) + ε Q ,
(27)
ε Q being the quadrature error. Setting z = κ ( x i x j ), we, thus, approximate G ̂ i j as
G ̂ i j p o 2 8 π q = 1 Q w q e i x i · θ ̂ q   e i x j · θ ̂ q ,
(28)
a result that, in turn, yields, upon application to all pairs of target points, the following (approximate) decomposition of the target correlation matrix G ̂:
G ̂ Φ Φ h = q = 1 Q ϕ q ϕ q h , Φ = [ ϕ 1 , , ϕ q ] C m × Q , Φ j q = ( ϕ q ) j = w q e i x j · θ ̂ q .
(29)
The size 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 j0 there. In this study, we have | z | κ 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 2R-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 2 R 1. Expressing j 0 ( | z | ) in (26) using the Jacobi–Anger expansion [see Ref. 13, Eq. (10.60.7)], we can estimate ε Q as | ε Q | C l 2 R m = l l j l ( | z | ) ε Q , l m, with ε Q , l m denoting the quadrature error for S ̂ Y l m ( θ ̂ )   d S ( θ ̂ ). 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 e κ L / 4, i.e., Q ( e κ 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 ̂, whose computation is moreover economical as the m-vectors ϕ q are given explicitly.
The low-rank approximation (29) greatly reduces the computational load in evaluating So. Indeed, using (29) in (23) gives
S o = ( Y Σ 1 Z ) ( Y Σ 1 Z ) h , Z : = X h Φ C d × Q ,
(30)
whose evaluation needs only O(m) computational work and memory, down from O ( m 2 ) if using the full matrix G ̂.

Remark 1. Unlike in the expansion G ̂ = q = 1 m ψ q ψ q T λ q in terms of eigenpairs, the vectors ϕ q in (29) are not orthogonal.

We now address the case where ThT may be ill-conditioned (i.e., have a large condition number), with T still assumed to have full column rank. Let
s q = Y Σ 1 X h ϕ q .
(31)
Then the minimum-norm solution to (14) is given by
S o = q = 1 Q s q s q h .
(32)
We notice that s q given by (31) is (by our assumption on T) the unique solution of the least squares problem
s q = arg min u C d 1 2 T u ϕ k 2 2 ,
(33)
where · 2 is the Euclidean norm. Therefore, to obtain So given by (32), we have to solve at most Q problems of the type (33). To cater for ThT being possibly ill-conditioned, we add a regularization term in (33) to get the problem
s q ( α ) = arg min u C d 1 2 ( T u ϕ q 2 2 + α u 2 2   )
(34)
(where α > 0 is a regularization parameter), whose unique minimizer is given in closed form as
s q ( α ) = Y ( Σ + α   I ) 1 X h ϕ q .
(35)
The resulting input correlation matrix is then given by
S o 1 ( α ) : = q = 1 Q s q ( α ) s q h ( α ) = ( Y ( Σ + α   I ) 1 Z )   ( Y ( Σ + α   I ) 1 Z ) h .
(36)
We can show that S o 1 ( α ) converges to the minimum-norm solution S o 1 of our original problem (14) as α 0. Indeed, from (31) and (35), we have
s q ( α ) s q = Y [ ( Σ + α   I ) 1 Σ 1 ] X h ϕ q
(37)
= Y ( Σ + α   I ) 1 [ I ( Σ + α   I ) Σ 1 ] X h ϕ q
(38)
= α Y ( Σ + α   I ) 1 Σ 1 X h ϕ q .
(39)
Therefore, s q ( α ) s q 0 and S o 1 ( α ) S o F 0 as α 0 [since S o 1 ( α ) and So are both given by finite sums of tensor products].

1. Second regularization method

Alternatively, we can consider the regularized version
S o 2 ( α ) : = arg min S C d × d J α ( S ) , J α ( S ) : = 1 2 A S G ̂ F 2 + α 2 S F 2
(40)
of the original least squares problem (14). In the regularization term, one-half the Frobenius norm of S is equal to the sum of the mean square acceleration of the sources, 1 / 2 S F = 1 / 2 tr ( S S H ) = 1 / 2 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 h T S o 2 ( α ) T h T + α S o 2 ( α ) = T G ̂ T h .
(41)
Invoking again the reduced SVD (22) of T, this becomes
Σ 2 H ( α ) Σ 2 + α H ( α ) = Σ H ̂ Σ , with   H : = Y h S o 2 ( a ) Y .
(42)
Since Σ is diagonal, the above equation decouples into component-wise scalar equations, whereby
H i j ( α ) = σ i σ j σ i 2 σ j 2 + α ( Z Z h ) i j , 1 i , j d ,
(43)
and S o 2 ( α ) = Y H ( α ) Y h is readily found once H ( α ) is evaluated using the above formula. Moreover, it is easy to verify that (43) with α = 0 yields S o = Y H ( 0 ) Y h through (23) and that we have
H i j ( α ) H i j ( 0 ) = α σ i σ j ( σ i 2 σ j 2 + α ) ( Z Z h ) i j 1 i , j d .
(44)
Consequently, this second regularization approach also verifies S o 2 ( α ) S o 0 as α 0.

Remark 2. The first regularization yields H i j ( α ) = ( Z Z h ) i j / ( σ i + α ) ( σ j + α ) instead of (43), with Hij 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 So to problem (14), as well as its regularized approximations S o 1 ( α ) and S o 2 ( α ), 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 C d × 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 = S 2, and we have S F 2 = S 1 F 2 + S 2 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 ̂ is Hermitian)
2 J ( S ) = G 1 G ̂ F 2 + G 2 F 2 ,
(45)
so that optimality implies G 2 = 0; hence, S 2 = 0, since by assumption, 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 Q d in the present context. The leading C × 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 ; α I ], resulting in C = 2 (Ref. 24, Sec. 5.2).

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 × 4.25 × 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 × 10 5   kg / ( m 2   s), while the speed of sound and mass density are c = 340   m / s and ρ 0 = 1   kg / m 3, respectively. The side walls of the room contain uniformly spaced acoustic sources, each with an area of 0.25 × 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 × 2 × 2   m 3 and is located in the center of the room.

FIG. 1.

(Color online) Room geometry, shown with 36 wall sources.

FIG. 1.

(Color online) Room geometry, shown with 36 wall sources.

Close modal

For each of the studied cases, we first compute a low-rank approximation of the target correlation G ̂ as per (29), using Q = 2 R 2 = 200. With reference to the discussion following Eq. (29), we have L 2 3 m and (for the highest frequency) κ 5.54 m−1; hence, | z | 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 ̂ 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 q Q. Finally, the optimal source correlation So 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 library27 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 15 000.

The Schroeder frequency for this room is calculated from (1). The reverberation time is estimated using the Norris–Eyring relation,28 
T 60 = 24 V ln ( 10 ) A   c ln ( 1 α rand ) ,
(46)
where A = 2 ( L x L y + L x L z + L y L z ) is the room surface area, and αrand is the random-incidence absorption coefficient20 
α rand = 1 0 π / 2 | Z cos ( θ ) ρ 0 c Z cos ( θ ) + ρ 0 c | 2 sin ( 2 θ )   d θ .
(47)
Using the above expressions in (1), we obtain f s = 1001 Hz. Hence, we point out that the frequencies used in the examples (i.e., [ 150 , 300 ] Hz) are well below the Schroeder limit for this room.
We define the relative error in spatial correlation (or residual) as
ϵ G : = A S o G ̂ F G ̂ F .
(48)
Without loss of generality, we use p o 2 / 2 = 1 Pa2 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
ϵ m s : = i | G o i i 1 | 2 m ,
(49)
indicative of how much the mean square pressure (again, normalized to unity in our case) departs from being spatially constant. This is a widely used metric to judge the level of sound diffusion in laboratory experiments15,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 ( α ) for the three frequencies 150, 250, and 300 Hz studied in detail. We notice that for these frequencies, both errors, ϵG and ϵms, 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.

TABLE I.

Results summary.

Sources 150 Hz 250 Hz 300 Hz
ϵms ϵG α ϵms ϵG α ϵms ϵG α
36  0.17  0.33  2 × 10 3  0.30  0.49  5 × 10 3  0.43  0.64  2 × 10 3 
64  0.09  0.17  1 × 10 3  0.10  0.26  5 × 10 3  0.18  0.37  1 × 10 3 
100  0.07  0.15  1 × 10 4  0.08  0.23  1 × 10 4  0.08  0.24  1 × 10 4 
Sources 150 Hz 250 Hz 300 Hz
ϵms ϵG α ϵms ϵG α ϵms ϵG α
36  0.17  0.33  2 × 10 3  0.30  0.49  5 × 10 3  0.43  0.64  2 × 10 3 
64  0.09  0.17  1 × 10 3  0.10  0.26  5 × 10 3  0.18  0.37  1 × 10 3 
100  0.07  0.15  1 × 10 4  0.08  0.23  1 × 10 4  0.08  0.24  1 × 10 4 
FIG. 2.

(Color online) Error behavior with frequency.

FIG. 2.

(Color online) Error behavior with frequency.

Close modal

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 Go. 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.

To this end, recall that the rank of the approximate spatial correlation G o = T S o T h is at most d due to the rank of T. Using the best approximation error G ̂ d in Eq. (24), we can obtain a lower bound E ( d ) for the relative residual of our inverse problem,
G ̂ G o ( d ) F G ̂ F E ( d ) k = d + 1 m λ k 2 j = 1 m λ j 2 .
(50)
We see that the lower bound depends on the decay of the eigenvalue spectrum, which in turn depends on frequency byway of the spatial correlation (3). In the bottom of Fig. 2, the lower bound or truncation error is plotted as a function of frequency for different d. Comparing this plot to the plots of ϵG and ϵms 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.

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.

FIG. 3.

(Color online) Spatial correlation with respect to the center of the target region.

FIG. 3.

(Color online) Spatial correlation with respect to the center of the target region.

Close modal
FIG. 4.

(Color online) Spatial correlation along a diagonal across the target region.

FIG. 4.

(Color online) Spatial correlation along a diagonal across the target region.

Close modal

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.

FIG. 5.

(Color online) Mean square pressure in the entire room along x axis. Results for the uncorrelated sources are normalized by dividing by the maximum mean square pressure of the 36-source optimal solution along the diagonal.

FIG. 5.

(Color online) Mean square pressure in the entire room along x axis. Results for the uncorrelated sources are normalized by dividing by the maximum mean square pressure of the 36-source optimal solution along the diagonal.

Close modal

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.

FIG. 6.

(Color online) Mean square pressure in the entire room.

FIG. 6.

(Color online) Mean square pressure in the entire room.

Close modal

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.

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.

The authors have no conflicts to disclose.

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

1.
M.
Schroeder
, “
Die statistischen Parameter der Frequenzkurven von grossßen Räumen” (“The statistical parameters of the frequency curves of large rooms”)
,
Acustica
4
,
594
600
(
1954
).
2.
M. R.
Schroeder
, “
Statistical parameters of the frequency response curves of large rooms
,”
J. Audio Eng. Soc.
35
(
5
),
299
306
(
1987
).
3.
H.
Kuttruff
,
Room Acoustics
, 5th ed. (
Spon
,
London
,
2009
).
4.
T. J.
Schultz
, “
Diffusion in reverberation rooms
,”
J. Sound Vibr.
16
(
1
),
17
28
(
1971
).
5.
P. A.
Larkin
and
D. O.
Smallwood
, “
Control of an acoustical speaker system in a reverberant chamber
,” Sandia Technical Memo SAND2003-3008C (
Sandia National Laboratories
,
Albuquerque, NM
,
2003
).
6.
F. W.
Grosveld
and
S. A.
Rizzi
, “
Controlled reverberant acoustic excitation capabilities at NASA Langley Research Center
,” in
Proceedings of the 43rd AIAA Aerospace Sciences Meeting
, Reno, NV (January 10–13,
2005
) (
AIAA
,
Reston, VA
).
7.
C.
Maury
and
T.
Bravo
, “
The experimental synthesis of random pressure fields: Practical feasibility
,”
J. Acoust. Soc. Am.
120
,
2712
2723
(
2006
).
8.
M.
Alvarez Blanco
,
P.
Van Vlierberghe
,
M.
Rossetti
,
K.
Janssens
,
B.
Peeters
, and
W.
Desmet
, “
Pre-test analysis to reproduce random pressure fields with multi-channel acoustic control
,”
Mech. Syst. Signal Proc.
163
(
May 2021
),
108103
(
2022
).
9.
S. J.
Elliott
and
J.
Cheer
, “
Modeling local active sound control with remote sensors in spatially random pressure fields
,”
J. Acoust. Soc. Am.
137
(
4
),
1936
1946
(
2015
).
10.
S.
Zhao
,
Q.
Zhu
,
E.
Cheng
, and
I. S.
Burnett
, “
A room impulse response database for multizone sound field reproduction (L)
,”
J. Acoust. Soc. Am.
152
(
4
),
2505
2512
(
2022
).
11.
A. G.
de Miguel
,
M.
Alvarez Blanco
,
E.
Matas
,
H.
Bériot
,
J.
Cuenca
,
O.
Atak
,
K.
Janssens
, and
B.
Peeters
, “
Virtual pre-test analysis for optimization of multi-channel control strategies in direct field acoustic testing
,”
Mech. Syst. Signal Proc.
184
,
109652
(
2023
).
12.
C.
Van Hoorickx
and
E. P. B.
Reynders
, “
Numerical realization of diffuse sound pressure fields using prolate spheroidal wave functions
,”
J. Acoust. Soc. Am.
151
(
3
),
1710
1721
(
2022
).
13.
F. W. J.
Olver
,
D. W.
Lozier
,
R. F.
Boisvert
, and
C. W.
Clark
,
The NIST Handbook of Mathematical Functions
(
Cambridge University Press
,
New York
,
2010
).
14.
R. V.
Waterhouse
, “
Statistical properties of reverberant sound fields
,”
J. Acoust. Soc. Am.
43
(
6
),
1436
1444
(
1968
).
15.
F.
Jacobsen
,
The Diffuse Sound Field: Statistical Considerations Concerning the Reverberant Field in the Steady State
(
Technical University of Denmark
,
Kongens Lyngby, Denmark
,
1979
).
16.
B.
Rafaely
, “
Spatial-temporal correlation of a diffuse sound field
,”
J. Acoust. Soc. Am.
107
(
6
),
3254
3258
(
2000
).
17.
M.
Grigoriu
,
Stochastic Systems: Uncertainty Quantification and Propagation
(
Springer
,
London
,
2012
).
18.
F.
Jacobsen
and
T.
Roisin
, “
The coherence of reverberant sound fields
,”
J. Acoust. Soc. Am.
108
(
1
),
204
210
(
2000
).
19.
H.
Nélisse
and
J.
Nicolas
, “
Characterization of a diffuse field in a reverberant room
,”
J. Acoust. Soc. Am.
101
(
6
),
3517
3524
(
1997
).
20.
P. M.
Morse
and
K. U.
Ingard
,
Theoretical Acoustics
(
Princeton University
,
Princeton, NJ
,
1986
).
21.
F.
Ihlenburg
,
Finite Element Analysis of Acoustic Scattering
(
Springer
,
New York
,
1998
).
22.
L.
Demkowicz
, “
Asymptotic convergence in finite and boundary element methods: Part 1: Theoretical results
,”
Comput. Math. Appl.
27
(
12
),
69
84
(
1994
).
23.
K.
Atkinson
and
W.
Han
,
Spherical Harmonics and Approximations on the Unit Sphere: An Introduction
(
Springer
,
New York
,
2012
).
24.
G. H.
Golub
and
C. F.
Van Loan
,
Matrix Computations
, 4th ed. (
Johns Hopkins University Press
,
Baltimore
,
2013
).
25.
P. C.
Hansen
, “
Analysis of discrete ill-posed problems by means of the L-curve
,”
SIAM Rev.
34
(
4
),
561
580
(
1992
).
26.
P. C.
Hansen
,
Discrete Inverse Problems: Insight and Algorithms
(
Society for Industrial and Applied Mathematics
,
Philadelphia
,
2010
).
27.
M.
Alnaes
,
J.
Blechta
,
J.
Hake
,
A.
Johansson
,
B.
Kehlet
,
A.
Logg
,
C.
Richardson
,
J.
Ring
,
M. E.
Rognes
, and
G. N.
Wells
, “
The FEniCS project version 1.5
,”
Arch. Numer. Softw.
3
(
100
),
9
23
(
2015
).
28.
C. F.
Eyring
, “
Reverberation time in ‘dead’ rooms
,”
J. Acoust. Soc. Am.
1
(
2A
),
217
241
(
1930
).