Wave-domain acoustic energy difference maximization method for multi-zone sound field reproduction

Abstract

Multi-zone sound field reproduction aims to create personal sound zones within distinct spatial regions of a shared environment using a loudspeaker array. This paper proposes a wave-domain acoustic energy difference maximization (WDAEDM) method for reproducing multi-zone sound field. The proposed method first represents sound fields and transfer functions with spatial harmonic expansion, then constructs an acoustic energy difference maximization (AEDM) model in the wave domain using sound field coefficients and their relationship with loudspeaker weights, and finally solves this model through eigenvalue decomposition. Simulation results of free-field and automotive cabin acoustic environments have demonstrated that the proposed WDAEDM method presents a more uniform energy distribution within zones compared with the existing spatial-domain methods, and also eliminates the parameter selection issues encountered in existing wave-domain methods.

Introduction

Multi-zone sound field reproduction utilizes a loudspeaker array to deliver well-isolated sound zones in different regions of a shared space,^1–5 which has broad application prospects in smart car cabins, smart mobile devices, and advanced audio-visual systems. The essence of this technology is to design loudspeaker weights. With the effects of loudspeaker weights and transfer functions, bright zones and dark zones are formed in the space. In bright zones, the high-energy signal is reproduced, while in dark zones, the signal is attenuated as much as possible.

The existing reproduction methods are classified into the energy-based methods and the matching methods. The representatives are acoustic contrast control (ACC)^6–8 and pressure matching (PM),^9–11 respectively. ACC obtains the loudspeaker weights by maximizing the ratio of the energy between bright zones and dark zones, that is, maximizing acoustic contrast (AC). PM achieves weight estimation by minimizing the difference between the reproduced and target signal. ACC method shows higher AC but involves matrix inversion, which leads to high sensitivity to noise when the inverted matrix is severely ill-conditioned.^12–16 To mitigate the effects caused by the inversion of ill-conditioned matrices, regularization methods have been adopted.^17,18 It is found that ACC is highly sensitive to the regularization parameter, and the optimal regularization parameter varies significantly across different frequencies, system configurations, and noise conditions. Therefore, it is a challenging task to select an appropriate regularization parameter. Shin et al. proposed the acoustic energy difference maximization (AEDM) method, which avoids matrix inversion and thereby circumvents regularization. ¹⁹ Unlike ACC, AEDM obtains weights by maximizing the energy difference between bright and dark zones and also provides sufficiently high AC.

Both ACC and AEDM are based on discrete spatial point control. However, the discreteness and the limited number of control points lead to performance degradation away from these points, resulting in uneven energy distribution within bright and dark zones. To overcome this limitation, Han et al. proposed wave-domain acoustic contrast control (WDACC),^20,21 which represents sound fields through spatial harmonic expansion and uses sound field coefficients to establish the ACC model in the wave domain. Recently, Wen et al. ²² also exploited harmonic expansion to formulate a wave-domain multi-objective optimization problem to achieve a tradeoff between AC, signal distortion and array effort. Compared with spatial-domain methods, these wave-domain methods present a more uniform energy distribution within zones. However, both methods require regularization parameters and/or Lagrange multipliers, whose values strongly depend on frequency and substantially affect performance, thereby making their determination time-consuming and labor-intensive.

Inspired by the wave-domain methods and AEDM, this paper proposes wave-domain acoustic energy difference maximization (WDAEDM). The proposed method first performs spatial harmonic expansion on the sound field and transfer functions with cylindrical harmonic functions. It then uses sound field coefficients and their relationship with loudspeaker weights to construct a model that maximizes the acoustic energy difference between bright zones and dark zones, and finally applies eigenvalue decomposition to obtain the solutions. By incorporating harmonic expansion and the AEDM concept into the model formulation, the proposed method inherits the advantages of both wave-domain approaches and AEDM, namely, achieving a more uniform energy distribution within zones while avoiding the need to select regularization parameters and Lagrange multipliers. In this paper, the environment is assumed to be stable. The paper is organized as follows: In Theory, the proposed WDAEDM method is elaborated. In Simulation, the performance is evaluated. The conclusions are summerized in the last section.

Theory

Model construction

Figure 1 shows a multi-zone sound field reproduction system. A circular loudspeaker array with a radius of $r_{\mathrm{L}}$ is used to create a bright zone and a dark zone within the region. The number of loudspeakers is L, and the radii of the bright zone and the dark zone are $r_{\mathrm{B}}$ and $r_{\mathrm{D}}$ , respectively. $r_{\mathrm{B}}$ and $r_{\mathrm{D}}$ are smaller than $r_{\mathrm{L}}$ . $M_{\mathrm{B}}$ and $M_{\mathrm{D}}$ control points are placed in the bright zone and the dark zone. Each loudspeaker is regarded as a point source and placed in a free field, and thus the transfer function between the $m\hspace{0.17em}\text{th}\left(m=1,2,\cdots ,M_{\mathrm{B}\left(\mathrm{D}\right)}\right)$ control point and the $l\hspace{0.17em}\text{th}\left(l=1,2,\cdots ,L\right)$ loudspeaker is

$g_{l,m}^{\mathrm{B}\left(\mathrm{D}\right)}\left(\omega \right)=\frac{{\mathrm{e}}^{-\mathrm{i}k\left|{\mathbf{r}}_m-{\mathbf{r}}_l\right|}}{4\mathrm{\pi }\left|{\mathbf{r}}_m-{\mathbf{r}}_l\right|},$

(1)

where ${\mathbf{r}}_m$ and ${\mathbf{r}}_l$ are positions of the $m\hspace{0.17em}\text{th}$ control point and the $l\hspace{0.17em}\text{th}$ loudspeaker, $\mathrm{i}=\sqrt{-1}$ is the imaginary unit, $k=\omega /c$ is the wave number, $\omega$ is the circular frequency, and $c$ is the speed of sound.

Figure 1.

Multi-zone sound field reproduction system.

The weight of the $l\hspace{0.17em}\text{th}$ loudspeaker is denoted as $q_l$ , and the pressure at each control point within the bright zone and the dark zone is expressed as

$p_m^{\mathrm{B}\left(\mathrm{D}\right)}={\displaystyle \sum }_l^L{g}_{l,m}^{\mathrm{B}\left(\mathrm{D}\right)}{q}_l.$

(2)

For brevity, $\omega$ is omitted in this and subsequent equations. Let ${\mathbf{p}}^{\mathrm{B}}={\left[p_1^{\mathrm{B}},p_2^{\mathrm{B}},\cdots ,p_{M_{\mathrm{B}}}^{\mathrm{B}}\right]}^{\mathrm{T}}$ $\in {\mathbb{C}}^{M_{\mathrm{B}}\times 1}$ and ${\mathbf{p}}^{\mathrm{D}}={\left[p_1^{\mathrm{D}},p_2^{\mathrm{D}},\cdots ,p_{M_{\mathrm{D}}}^{\mathrm{D}}\right]}^{\mathrm{T}}\in {\mathbb{C}}^{M_{\mathrm{D}}\times 1}$ denote the pressure vectors of the bright zone and the dark zone, respectively. The superscript “^T” represents the transport operator. According to equation (2), ${\mathbf{p}}^{\mathrm{B}}$ and ${\mathbf{p}}^{\mathrm{D}}$ can be written as

$\{\begin{array}{l}{\mathbf{p}}^{\mathrm{B}}={\mathbf{G}}^{\mathrm{B}}\mathbf{q}\\ {\mathbf{p}}^{\mathrm{D}}={\mathbf{G}}^{\mathrm{D}}\mathbf{q}\end{array},$

(3)

where $\mathbf{q}={\left[q_1,q_2,\cdots ,q_L\right]}^{\mathrm{T}}\in {\mathbb{C}}^{L\times 1}$ is the weight vector, and ${\mathbf{G}}^{\mathrm{B}}\in {\mathbb{C}}^{M_{\mathrm{B}}\times L}$ and ${\mathbf{G}}^{\mathrm{D}}\in {\mathbb{C}}^{M_{\mathrm{D}}\times L}$ are the transfer function matrices which take equation (1) as elements.

Acoustic energy difference maximization

AEDM method estimates the loudspeaker weights by maximizing the energy difference between the bright zone and the dark zone. The energy within the bright zone and the dark zone is

$\{\begin{array}{l}{E}^{\mathrm{B}}=\frac{{\left({\mathbf{p}}^{\mathrm{B}}\right)}^{\mathrm{H}}{\mathbf{p}}^{\mathrm{B}}}{M_{\mathrm{B}}}={\mathbf{q}}^{\mathrm{H}}{\mathbf{R}}^{\mathrm{B}}\mathbf{q}\\ E^{\mathrm{D}}=\frac{{\left({\mathbf{p}}^{\mathrm{D}}\right)}^{\mathrm{H}}{\mathbf{p}}^{\mathrm{D}}}{M_{\mathrm{D}}}={\mathbf{q}}^{\mathrm{H}}{\mathbf{R}}^{\mathrm{D}}\mathbf{q}\end{array},$

(4)

where ${\mathbf{R}}^{\mathrm{B}}={\left({\mathbf{G}}^{\mathrm{B}}\right)}^{\mathrm{H}}{\mathbf{G}}^{\mathrm{B}}/M_{\mathrm{B}}\in {\mathbb{C}}^{L\times L}$ and ${\mathbf{R}}^{\mathrm{D}}={\left({\mathbf{G}}^{\mathrm{D}}\right)}^{\mathrm{H}}{\mathbf{G}}^{\mathrm{D}}/M_{\mathrm{D}}\in {\mathbb{C}}^{L\times L}$ are the spatial auto-correlation matrices. The superscript “^H” represents the Hermitian operator. To maximize the energy difference between the bright zone and the dark zone, the cost function of AEDM is formulated as

$\underset{\mathbf{q}}{\max }\hspace{0.17em}J\left(\mathbf{q}\right)=\frac{E^{\mathrm{B}}-\xi E^{\mathrm{D}}}{{\mathbf{q}}^{\mathrm{H}}\mathbf{q}}=\frac{{\mathbf{q}}^{\mathrm{H}}\left({\mathbf{R}}^{\mathrm{B}}-\xi {\mathbf{R}}^{\mathrm{D}}\right)\mathbf{q}}{{\mathbf{q}}^{\mathrm{H}}\mathbf{q}},$

(5)

where $E^{\mathrm{L}}={\mathbf{q}}^{\mathrm{H}}\mathbf{q}$ is the energy of the array input, which is used to normalize the energy difference. $\xi \left(\xi >0\right)$ is a tuning factor. With a small $\xi$ , the cost function tends to maximize the energy within the bright zone. On the contrary, when $\xi$ is relatively large, the cost function primarily aims to suppress the energy in the dark zone.

To obtain $\widehat{\mathbf{q}}$ that maximizes equation (5), we first multiply ${\mathbf{q}}^{\mathrm{H}}\mathbf{q}$ on the both sides of equation (5), and the cost function becomes

$J\left(\mathbf{q}\right){\mathbf{q}}^{\mathrm{H}}\mathbf{q}={\mathbf{q}}^{\mathrm{H}}\left({\mathbf{R}}^{\mathrm{B}}-\xi {\mathbf{R}}^{\mathrm{D}}\right)\mathbf{q}.$

(6)

Then, by taking the partial derivatives of equation (6) with regard to $\mathbf{q}$ , we obtain

$\frac{\partial J\left(\mathbf{q}\right)}{\partial \mathbf{q}}{\mathbf{q}}^{\mathrm{H}}\mathbf{q}+J\left(\mathbf{q}\right)\frac{\partial \left({\mathbf{q}}^{\mathrm{H}}\mathbf{q}\right)}{\partial \mathbf{q}}=\frac{\partial \left({\mathbf{q}}^{\mathrm{H}}\left({\mathbf{R}}^{\mathrm{B}}-\xi {\mathbf{R}}^{\mathrm{D}}\right)\mathbf{q}\right)}{\partial \mathbf{q}}.$

(7)

When $J\left(\mathbf{q}\right)$ takes its extremum, $\partial J\left(\mathbf{q}\right)/\partial \mathbf{q}=0$ holds. Thus, equation (7) becomes

$J\left(\mathbf{q}\right)\frac{\partial \left({\mathbf{q}}^{\mathrm{H}}\mathbf{q}\right)}{\partial \mathbf{q}}=\frac{\partial \left({\mathbf{q}}^{\mathrm{H}}\left({\mathbf{R}}^{\mathrm{B}}-\xi {\mathbf{R}}^{\mathrm{D}}\right)\mathbf{q}\right)}{\partial \mathbf{q}}.$

(8)

Since $\partial \left({\mathbf{q}}^{\mathrm{H}}\mathbf{q}\right)/\partial \mathbf{q}=2\mathbf{q}$ , and $\partial \left({\mathbf{q}}^{\mathrm{H}}\left({\mathbf{R}}^{\mathrm{B}}-\xi {\mathbf{R}}^{\mathrm{D}}\right)\mathbf{q}\right)/\partial \mathbf{q}=2\left({\mathbf{R}}^{\mathrm{B}}-\xi {\mathbf{R}}^{\mathrm{D}}\right)\mathbf{q}$ ,

$J\left(\mathbf{q}\right)\mathbf{q}=\left({\mathbf{R}}^{\mathrm{B}}-\xi {\mathbf{R}}^{\mathrm{D}}\right)\mathbf{q}.$

(9)

Obviously, the extremum of $J\left(\mathbf{q}\right)$ is the largest eigenvalue of the matrix ${\mathbf{R}}^{\mathrm{B}}-\xi {\mathbf{R}}^{\mathrm{D}}$ , and the optimal loudspeaker weights ${\widehat{\mathbf{q}}}_{\text{AEDM}}$ is the eigenvector corresponding to the largest eigenvalue.

Wave-domain acoustic energy difference maximization

AEDM achieves multi-zone sound reproduction based on discrete spatial control points. However, due to the discreteness and the limited number of these points, the performance of AEDM degrades in the regions away from the control points, resulting in uneven energy distribution within the zones. To overcome this limitation, spatial harmonic expansion is introduced into AEDM, leading to WDAEDM.

As shown in Figure 1, there is no source within the zones. Therefore, the sound pressure at any point within these zones can be expanded with respect to the zone center $O^{B\left(D\right)}$ using cylindrical harmonic functions,

$p^{\mathrm{B}\left(\mathrm{D}\right)}\left(\stackrel{⌢}{\mathbf{r}}\right)={\displaystyle \sum }_{n=-\infty }^{\infty }{a}_n^{\mathrm{B}\left(\mathrm{D}\right)}{J}_n\left(k\stackrel{⌢}{r}\right){\mathrm{e}}^{\mathrm{i}n\stackrel{⌢}{\varphi }},$

(10)

where $\stackrel{⌢}{\mathbf{r}}=\left(\stackrel{⌢}{r},\varphi \right)$ denotes the position of this point measured from $O^{\mathrm{B}\left(\mathrm{D}\right)}$ , $J_n(·)$ is the first Bessel function with the order of n, and $a_n^{\mathrm{B}\left(\mathrm{D}\right)}$ are the sound field coefficients. Similarly, transfer functions represent the responses of unit-amplitude sources at given points, and can be also represented by

$g_l^{\mathrm{B}\left(\mathrm{D}\right)}\left(\stackrel{⌢}{\mathbf{r}}\right)={\displaystyle \sum }_{n=-\infty }^{\infty }{b}_{n,l}^{\mathrm{B}\left(\mathrm{D}\right)}{J}_n\left(k\stackrel{⌢}{r}\right)e^{\mathrm{i}n\stackrel{⌢}{\varphi }},$

(11)

where $g_l^{\mathrm{B}\left(\mathrm{D}\right)}\left(\stackrel{⌢}{\mathbf{r}}\right)$ is the transfer function between the $l\hspace{0.17em}\text{th}$ loudspeaker and the position $\stackrel{⌢}{\mathbf{r}}$ . $b_n^{\mathrm{B}\left(\mathrm{D}\right)}$ are the sound field coefficients of transfer functions. $b_n^{\mathrm{B}\left(\mathrm{D}\right)}$ can be obtained by sampling the sound pressures at evenly spaced control points on a circle, ²²

$b_n^{\mathrm{B}\left(\mathrm{D}\right)}=\frac{1}{J_n\left(k{r}_{\mathrm{B}\left(\mathrm{D}\right)}\right)}{\displaystyle \sum }_{m=1}^{M_{\mathrm{B}\left(\mathrm{D}\right)}}\frac{g_l^{\mathrm{B}\left(\mathrm{D}\right)}\left({\stackrel{⌢}{\mathbf{r}}}_m\right)}{M_{\mathrm{B}\left(\mathrm{D}\right)}}{\mathrm{e}}^{-\mathrm{i}2\mathrm{\pi }nm/M_{\mathrm{B}\left(\mathrm{D}\right)}}.$

(12)

where $g_l^{\mathrm{B}\left(\mathrm{D}\right)}\left({\stackrel{⌢}{\mathbf{r}}}_m\right)$ denotes the transfer function between the $l\hspace{0.17em}\text{th}$ loudspeaker and the $m\hspace{0.17em}\text{th}$ control point. Combining equation (2), equation (10), and equation (11), we could build the relationship between $a_n^{\mathrm{B}\left(\mathrm{D}\right)}$ and $b_{n,l}^{\mathrm{B}\left(\mathrm{D}\right)}$

$a_n^{\mathrm{B}\left(\mathrm{D}\right)}={\displaystyle \sum }_{l=1}^L{b}_{n,l}^{\mathrm{B}\left(\mathrm{D}\right)}{q}_l.$

(13)

In the wave domain, the energy in the bright zone and the dark zone is expressed as

$\{\begin{array}{l}{E}^{\mathrm{B}}=\frac{1}{S_{\mathrm{B}}}{\iint }_{S_{\mathrm{B}}}{\left|p^{\mathrm{B}}\left(\stackrel{⌢}{\mathbf{r}}\right)\right|}^2\stackrel{⌢}{\mathbf{r}}\mathrm{d}\stackrel{⌢}{\mathbf{r}}\mathrm{d}\stackrel{⌢}{\varphi }=\frac{2}{r_{\mathrm{B}}^2}{\displaystyle \sum }_{n=-\infty }^{\infty }{\left(a_n^{\mathrm{B}}\right)}^2{\int }_0^{r_{\mathrm{B}}}{J}_n^2\left(k\stackrel{⌢}{r}\right)\stackrel{⌢}{r}\mathrm{d}\stackrel{⌢}{r}={\displaystyle \sum }_{n=-\infty }^{\infty }{\left(a_n^{\mathrm{B}}\right)}^2{w}_n\left(k{r}_{\mathrm{B}}\right)\\ E^{\mathrm{D}}=\frac{1}{S_{\mathrm{D}}}{\iint }_{S_{\mathrm{D}}}{\left|p^{\mathrm{D}}\left(\stackrel{⌢}{\mathbf{r}}\right)\right|}^2\stackrel{⌢}{\mathbf{r}}\mathrm{d}\stackrel{⌢}{\mathbf{r}}\mathrm{d}\stackrel{⌢}{\varphi }=\frac{2}{r_{\mathrm{D}}^2}{\displaystyle \sum }_{n=-\infty }^{\infty }{\left(a_n^{\mathrm{D}}\right)}^2{\int }_0^{r_{\mathrm{D}}}{J}_n^2\left(k\stackrel{⌢}{r}\right)\stackrel{⌢}{r}\mathrm{d}\stackrel{⌢}{r}={\displaystyle \sum }_{n=-\infty }^{\infty }{\left(a_n^{\mathrm{D}}\right)}^2{w}_n\left(k{r}_{\mathrm{D}}\right)\end{array},$

(14)

where $S_{\mathrm{B}}$ and $S_{\mathrm{D}}$ are the areas of the bright zone and the dark zone, and $w_n\left(k{r}_{\mathrm{B}\left(\mathrm{D}\right)}\right)=2/r_{\mathrm{B}\left(\mathrm{D}\right)}^2{\int }_0^{r_{\mathrm{B}\left(\mathrm{D}\right)}}{J}_n^2\left(k\stackrel{⌢}{r}\right)\stackrel{⌢}{r}\mathrm{d}\stackrel{⌢}{r}=J_n^2\left(k{r}_{\mathrm{B}\left(\mathrm{D}\right)}\right)-J_{n-1}\left(r_{\mathrm{B}\left(\mathrm{D}\right)}\right)J_{n+1}\left(r_{\mathrm{B}\left(\mathrm{D}\right)}\right)$ is the weight. ²¹ Equation (4) indicates that, in the spatial domain, the energy within a sound zone is estimated using the sound pressure at discrete control points. Equation (14) shows that, in the wave domain, the energy is obtained by integrating over the entire continuous region. This is the key to the WDAEDM method to overcome the problem of uneven energy distribution.

Previous studies^16,22–24 have shown that, when n exceeds $N\approx \lceil k{r}_{\mathrm{B}\left(\mathrm{D}\right)}\rceil$ , where $\lceil ·\rceil$ is the ceiling function, $w_n\left(k{r}_{\mathrm{B}\left(\mathrm{D}\right)}\right)$ drops to zeros. That is, only the basis functions $J_n\left(k\stackrel{⌢}{r}\right){\mathrm{e}}^{\mathrm{i}n\stackrel{⌢}{\varphi }}$ with indices in the set [-N, N] contribute significant energy to the sound zones. Therefore, we truncate the order to $N\approx \lceil k{r}_{\mathrm{B}\left(\mathrm{D}\right)}\rceil$ . Define ${\mathbf{a}}^{\mathrm{B}\left(\mathrm{D}\right)}={\left[a_{-N}^{\mathrm{B}\left(\mathrm{D}\right)},\cdots ,a_N^{\mathrm{B}\left(\mathrm{D}\right)}\right]}^{\mathrm{T}}\in {\mathbb{C}}^{\left(2N+1\right)\times 1}$ , ${\mathbf{b}}_l^{\mathrm{B}\left(\mathrm{D}\right)}={\left[b_{-N,l}^{\mathrm{B}\left(\mathrm{D}\right)},\cdots ,b_{N,l}^{\mathrm{B}\left(\mathrm{D}\right)}\right]}^{\mathrm{T}}\in {\mathbb{C}}^{\left(2N+1\right)\times 1}$ , ${\mathbf{B}}^{\mathrm{B}\left(\mathrm{D}\right)}=\left[{\mathbf{b}}_1^{\mathrm{B}\left(\mathrm{D}\right)},{\mathbf{b}}_2^{\mathrm{B}\left(\mathrm{D}\right)},\cdots ,{\mathbf{b}}_L^{\mathrm{B}\left(\mathrm{D}\right)}\right]\in {\mathbb{C}}^{\left(2N+1\right)\times L}$ , ${\mathbf{W}}^{\mathrm{B}\left(\mathrm{D}\right)}=$ $\text{Diag}\left(\left[w_{-N}\left(k{r}_{\mathrm{B}\left(\mathrm{D}\right)}\right),\cdots ,w_N\left(k{r}_{\mathrm{B}\left(\mathrm{D}\right)}\right)\right]\right)\in {\mathbb{C}}^{\left(2N+1\right)\times \left(2N+1\right)}$ , where $\text{Diag}(·)$ constructs a diagonal matrix using the vector in parenesis as its diagonal. We rewrite equation (13) and equation (14) into a matrix form,

${\mathbf{a}}^{\mathrm{B}\left(\mathrm{D}\right)}={\mathbf{B}}^{\mathrm{B}\left(\mathrm{D}\right)}\mathbf{q}.$

(15)

$\{\begin{array}{l}{E}^{\mathrm{B}}={\left({\mathbf{a}}^{\mathrm{B}}\right)}^{\mathrm{H}}{\mathbf{W}}^{\mathrm{B}}{\mathbf{a}}^{\mathrm{B}}\\ E^{\mathrm{D}}={\left({\mathbf{a}}^{\mathrm{D}}\right)}^{\mathrm{H}}{\mathbf{W}}^{\mathrm{D}}{\mathbf{a}}^{\mathrm{D}}\end{array}.$

(16)

By using the energy expression in equation (16) and the relationship in equation (15), we construct the AEDM model,

$\underset{\mathbf{q}}{\max }\hspace{0.5em}J\left(\mathbf{q}\right)=\frac{E_{\mathrm{B}}-\xi E_{\mathrm{D}}}{{\mathbf{q}}^{\mathrm{H}}\mathbf{q}}=\frac{{\left({\mathbf{a}}^{\mathrm{B}}\right)}^{\mathrm{H}}{\mathbf{W}}^{\mathrm{B}}{\mathbf{a}}^{\mathrm{B}}-\xi {\left({\mathbf{a}}^{\mathrm{D}}\right)}^{\mathrm{H}}{\mathbf{W}}^{\mathrm{D}}{\mathbf{a}}^{\mathrm{D}}}{{\mathbf{q}}^{\mathrm{H}}\mathbf{q}}=\frac{{\mathbf{q}}^{\mathrm{H}}\left({\overline{\mathbf{R}}}^{\mathrm{B}}-\xi {\overline{\mathbf{R}}}^{\mathrm{D}}\right)\mathbf{q}}{{\mathbf{q}}^{\mathrm{H}}\mathbf{q}},$

(17)

where ${\overline{\mathbf{R}}}^{\mathrm{B}}={\left({\mathbf{B}}^{\mathrm{B}}\right)}^{\mathrm{H}}{\mathbf{W}}^{\mathrm{B}}{\mathbf{B}}^{\mathrm{B}}$ and ${\overline{\mathbf{R}}}^{\mathrm{D}}={\left({\mathbf{B}}^{\mathrm{D}}\right)}^{\mathrm{H}}{\mathbf{W}}^{\mathrm{D}}{\mathbf{B}}^{\mathrm{D}}$ are the auto-correlation matrices in the wave domain. The steps used to solve equation (5) are used here to solve equation (17), and we get

$J\left(\mathbf{q}\right)\mathbf{q}=\left({\overline{\mathbf{R}}}^{\mathrm{B}}-\xi {\overline{\mathbf{R}}}^{\mathrm{D}}\right)\mathbf{q}$

(18)

Similarly, the extremum of $J\left(\mathbf{q}\right)$ is the largest eigenvalue of the matrix ${\overline{\mathbf{R}}}^{\mathrm{B}}-\xi {\overline{\mathbf{R}}}^{\mathrm{D}}$ , and the optimal loudspeaker weights ${\widehat{\mathbf{q}}}_{\text{WDAEDM}}$ is the eigenvector corresponding to the largest eigenvalue.

Simulations

Free-field simulation

This subsection examines the performance of the proposed method based on the free-field simulation. This simulation is performed in MATLAB 2022. Figure 2 shows the layout of a multi-zone sound reproduction system. A circular array with 30 loudspeakers and a radius of 2 m is placed at the origin. Inside the circular array, two circular zones with a radius of 0.5 m are defined as the bright zone and the dark zone. The centers of the two zones are located at (0.6 m, 0.5 m) and (−0.5 m, −0.7 m), respectively. 22 control points are evenly distributed along the circumferences of the bright zone and the dark zone. Multiple evaluation points are uniformly distributed within two zones, with a spacing of 0.01 m between the adjacent evaluation points.

Figure 2.

The layout of the multi-zone sound reproduction system used in the simulation. (▼ loudspeakers,● control points, █ evaluation points).

Evaluation metrices

We use the pressure at the evaluation points to calculate acoustic contrast (AC) and audible gain (AG). AC is defined as the ratio of the energy between the bright zone and the dark zone,

$\text{AC}=10\hspace{0.17em}{\log }_{10}\frac{E^{\mathrm{B}}}{E^{\mathrm{D}}}\approx 10\hspace{0.17em}{\log }_{10}\left(\frac{{\stackrel{\sim }{M}}_{\mathrm{D}}}{{\stackrel{\sim }{M}}_{\mathrm{B}}}\frac{{\mathbf{q}}^{\mathrm{H}}{\left({\stackrel{\sim }{\mathbf{G}}}^{\mathrm{B}}\right)}^{\mathrm{H}}{\stackrel{\sim }{\mathbf{G}}}^{\mathrm{B}}\mathbf{q}}{{\mathbf{q}}^{\mathrm{H}}{\left({\stackrel{\sim }{\mathbf{G}}}^{\mathrm{D}}\right)}^{\mathrm{H}}{\stackrel{\sim }{\mathbf{G}}}^{\mathrm{D}}\mathbf{q}}\right),$

(19)

where ${\stackrel{\sim }{M}}_{\mathrm{B}}$ and ${\stackrel{\sim }{M}}_{\mathrm{D}}$ are the number of evaluation points in the bright zone and the dark zone, and ${\stackrel{\sim }{\mathbf{G}}}^{\mathrm{B}}$ and ${\stackrel{\sim }{\mathbf{G}}}^{\mathrm{D}}$ represent the matrices whose elements are the transfer functions between the loudspeakers and the evaluation points. AG is defined as the ratio of the energy between the bright zone and the loudspeaker array input,

$\text{AG}=10\hspace{0.17em}{\log }_{10}\frac{E^{\mathrm{B}}}{E^{\mathrm{L}}}\approx 10\hspace{0.17em}{\log }_{10}\left(\frac{1}{{\stackrel{\sim }{M}}_{\mathrm{B}}}\frac{{\mathbf{q}}^{\mathrm{H}}{\left({\stackrel{\sim }{\mathbf{G}}}^{\mathrm{B}}\right)}^{\mathrm{H}}{\stackrel{\sim }{\mathbf{G}}}^{\mathrm{B}}\mathbf{q}}{{\mathbf{q}}^{\mathrm{H}}\mathbf{q}}\right)$

(20)

A high AC indicates strong isolation between the bright zone and the dark zone, while a high AG reflects efficient sound energy radiation by the loudspeaker array. The larger the two metrices, the better the performance.

The influence of the tuning factor

According to equation (14), the tuning factor $\xi$ simultaneously affects AC and AG. This subsection analyzes the influence of $\xi$ on the performance of the proposed WDAEDM method. $\xi$ is logarithmically sampled within the range [10⁻², 10⁵], with nine data points selected per order of magnitude, resulting in a total of 63 sampling points. The frequency ranges from 20 Hz to 1200 Hz. Figure 3 shows the AC and AG at different tuning factors across various frequencies.

Figure 3.

AC and AG under different tuning factors and frequencies.

As shown in Figure 3(a), when the frequency is fixed, AC increases with $\xi$ . For each frequency, AC reaches its peak and stops increasing once $\xi$ exceeds 1000. The black dashed curve in this figure indicates the variation of $\xi$ with frequency when AC is fixed at 20 dB. This curve shows that, within the frequency range of 20–1100 Hz, $\xi$ fluctuates between 0 and 60 to maintain an AC of 20 dB. When $\xi$ is set to 60, the WDAEDM method achieves an AC above 20 dB across the entire 20–1100 Hz frequency range. Figure 3(b) shows that, below 200 Hz, AG decreases as $\xi$ increases. However, above 200 Hz, $\xi$ has little impact on AG.

Performance comparison

This subsection compares the performance of the AEDM, ¹⁹ WDACC,^20,21 WDAEDM, and Lagrange method (LM). ²² Multi-zone sound field reproduction is performed at 200 Hz and 1200 Hz. For WDAEDM and AEDM, $\xi$ is set to 60. For WDACC, the regularization parameter corresponding to the optimal AC is selected using the traversal method (hereafter referred to as the optimal regularization parameter). For LM, its Lagrange multipliers are selected to satisfy the constraints on the energy in the dark zone and array effort. Figure 4 shows the energy distribution maps of the reproduced sound fields. At 200 Hz (Figure 4(a)–(d)), all these methods exhibit good sound zone separation. AEDM, WDACC, and WDAEDM achieve an AC above 28 dB and exhibit comparable performance in suppressing energy in the dark zone, whereas LM yields an AC slightly below 28 dB. The AG of the four methods is −18.61 dB, −17.15 dB, −16.77 dB, and −19.49 dB, respectively, indicating that the WDACC and WDAEDM methods achieve higher radiation efficiency in the bright zone compared to AEDM and LM. In addition, the maps reveal that WDACC, WDAEDM, and LM do not exhibit the dropouts present in the AEDM method within the bright zone, resulting in a more uniform energy distribution. As frequency increases to 1200 Hz, LM fails, and the performance of the other three methods decreases to varying degrees. The AC of AEDM drops to 12.12 dB, while that of WDACC and WDAEDM decreases to approximately 19.5 dB, indicating that the latter two methods perform similarly in terms of AC control and outperform AEDM. The underlying reason is that AEDM fails to sufficiently suppress energy in the dark zone, leading to uneven energy distribution. At this frequency, the AG of the three methods is −24.34 dB, −23.55 dB, and −23.70 dB, respectively. The WDACC and WDAEDM methods demonstrate higher radiation efficiency in the bright zone compared to AEDM. In summary, WDACC and WDAEDM outperform AEDM and LM.

Figure 4.

The energy distribution maps at 200 Hz and 1200 Hz.

In this part, we further compare the performance of WDACC and WDAEDM. The frequency ranges from 20 Hz to 1200 Hz. For WDAEDM, $\xi$ is set to 1000 at each frequency to achieve optimal AC. For WDACC, two strategies are adopted to determine the regularization parameter: (1) the traversal method and (2) the empirical maximal singular value-based method (SV). ²⁵

Figure 5 shows AC, AG, and sound pressure levels (SPLs) in the dark zone of each method. It is interesting to find that: (1) WDAEDM and WDACC (with the optimal regularization parameter) show comparable AC, AG, and SPL in the dark zone. However, it is worth noting that the optimal regularization parameter for WDACC varies significantly across frequencies, as shown in Figure 6. To achieve an optimal AC at each frequency, WDACC requires substantial manual effort and time to determine the optimal parameter frequency-by-frequency. In contrast, WDAEDM avoids this issue, as a fixed tuning factor can yield performance equivalent to WDACC across most frequencies, which is a key advantage of WDAEDM. (2) The WDACC method using the SV-based regularization fails to effectively suppress energy in the dark zone at several frequencies, resulting in SPLs exceeding 50 dB, which further leads to a marked decline in AC at those frequencies. The reason is that the empirically chosen regularization parameter does not reliably mitigate the impact of ill-conditioned matrix inversion on the performance. Therefore, WDAEDM has better performance than WDACC.

Figure 5.

AC, AG, and SPL in dark zone versus frequency.

Figure 6.

The optimal regularization parameter of WDACC versus frequency.

To summarize, compared with the other three methods, the proposed WDAEDM achieves higher AC and AG, resolves the issue of uneven energy distribution, and avoid spending extensive effort and time on regularization. Overall, WDAEDM delivers the best performance among all approaches discussed above.

Finite element-based simulation

This subsection conducts finite element-based simulations to evaluate the performance of the proposed method in a car cabin acoustic environment. The acoustic model of the car cabin, as shown in Figure 7, is built on the COMSOL Multiphysics platform. The car windows, dashboards, and doors are modeled with constant absorption coefficients ( $\alpha$ = 0.005, 0.01, and 0.01, respectively). ²⁶ The complex impedance of leather seats is set based on experimental data, ²⁷ and the ceiling and floor are defined using porous layers under impedance boundary conditions. The depths of ceiling and carpet are 0.7 cm and 0.5 cm. The flow resistance rates are 20,000 $\text{Pa}·\mathrm{s}/{\mathrm{m}}^2$ and 10,000 $\text{Pa}·\mathrm{s}/{\mathrm{m}}^2$ . After completing the mesh independence verification, a grid size of approximately 6 cm was finalized for the finite element model. The head positions of the driver and front passenger are defined as the bright zone and the dark zone, respectively. The linear loudspeaker array is used. The control points (microphones) are positioned at ear level, with the center located 0.55 m away from the loudspeaker array along the normal direction. 30 evaluation points (microphones) are densely distributed within the bright zone and the dark zone to assess performance.

Figure 7.

The finite element acoustic model of the car cabin.

AEDM, WDACC, WDAEDM, and LM are used for multi-zone reproduction. AC, AG, and SPLs in the dark zone (SPL_d) at 500 Hz for each method are listed in Table 1. The corresponding energy distribution maps of the reproduced sound fields are shown in Figure 8. For AEDM, the performance metrics calculated using the pressure at the evaluation points are generally inferior to those calculated from the control points. This is because AEDM is based on discrete point control, providing good performance only near those control points, with degradation observed farther away. For the WDACC and WDAEDM, they show comparable performance, and their performance metrics based on evaluation points are better than those based on control points, due to its ability to control the sound field over continuous regions. In addition, WDACC and WDAEDM outperform LM in term of AC and SPLd. However, the proposed WDAEDM does not need regularization parameters.

Table 1.

The comparison of performance metrics.

	AC/(dB)		AG/(dB)		SPLd/(dB)
Method	Control point	Evaluation point	Control point	Evaluation point	Control point	Evaluation point
AEDM	19.75	17.08	−34.19	−34.39	38.54	40.09
WDACC	21.85	22.21	−31.02	−31.01	39.71	39.39
WDAEDM	19.90	21.98	−35.05	−34.37	38.07	36.40
LM	17.13	17.10	−39.94	−34.93	40.21	40.12

(Note. “Control points” refer to results calculated based on the sound pressure at control points, while “Evaluation points” refer to results calculated based on the sound pressure at evaluation points).

Figure 8.

The energy distribution maps of the reproduced sound field in the car cabin.

Conclusion

This paper proposes the WDAEDM method for multi-zone sound field reproduction based on spatial harmonic expansion and acoustic energy difference maximization. Its performance is examined through the free-field and the car cabin acoustic environment simulations. Results show that compared with spatial-domain AEDM, WDAEDM yields a more uniform sound energy distribution. In comparison with other wave-domain methods, WDAEDM outperforms LM, and attains performance comparable to WDACC, while avoiding selecting a regularization parameter.

It is worth mentioning that the computational complexity of the proposed method increases with frequency because a higher order is required. Consequently, more loudspeakers are needed, resulting in increased cost. Therefore, the proposed method is recommended for use primarily at low frequencies. In addition, this paper does not account for environmental variations, and thus WDAEDM is currently suitable only for relatively stable acoustic environments. Future work will focus on developing robust WDAEDM against dynamic environments, but this lies beyond the scope of this paper.