Amplitude Matching for Multizone Sound Field Control

A multizone sound field control method, called amplitude matching, is proposed. The objective of amplitude matching is to synthesize a desired amplitude (or magnitude) distribution over a target region with multiple loudspeakers, whereas the phase distribution is arbitrary. Most of the current multizone sound field control methods are intended to synthesize a specific sound field including phase or to control acoustic potential energy inside the target region. In amplitude matching, a specific desired amplitude distribution can be set, ignoring sound propagation directions. Although the optimization problem of amplitude matching does not have a closed-form solution, our proposed algorithm based on the alternating direction method of multipliers (ADMM) allows us to accurately and efficiently synthesize the desired amplitude distribution. We also introduce the differential-norm penalty for a time-domain filter design with a small filter length. The experimental results indicated that the proposed method outperforms current multizone sound field control methods in terms of accuracy of the synthesized amplitude distribution.

Abstract-A multizone sound field control method, called amplitude matching, is proposed. The objective of amplitude matching is to synthesize a desired amplitude (or magnitude) distribution over a target region with multiple loudspeakers, whereas the phase distribution is arbitrary. Most of the current multizone sound field control methods are intended to synthesize a specific sound field including phase or to control acoustic potential energy inside the target region. In amplitude matching, a specific desired amplitude distribution can be set, ignoring sound propagation directions. Although the optimization problem of amplitude matching does not have a closed-form solution, our proposed algorithm based on the alternating direction method of multipliers (ADMM) allows us to accurately and efficiently synthesize the desired amplitude distribution. We also introduce the differential-norm penalty for a time-domain filter design with a small filter length. The experimental results indicated that the proposed method outperforms current multizone sound field control methods in terms of accuracy of the synthesized amplitude distribution.
Index Terms-Multizone sound field control, personal audio, pressure matching, amplitude matching.

I. INTRODUCTION
S YNTHESIZING a desired sound field inside a target region using multiple loudspeakers (secondary sources) has various applications, such as virtual/augmented reality, spatial noise cancellation, and personal sound-zone generation. Sound field synthesis/control methods can be classified into two categories. One includes methods based on boundary integral equations analytically derived from the Helmholtz equation, such as wave field synthesis and higher-order ambisonics [1], [2], [3], [4], [5], [6]. The other includes methods based on numerical optimization to minimize error between synthesized and desired sound fields inside the target region, such as pressure matching (PM) and (weighted) mode matching [2], [7], [8], [9], [10]. Since the numerical-optimization-based methods enable us to The authors are with the Graduate School of Information Science and Technology, The University of Tokyo, Tokyo 113-8656, Japan (e-mail: takumi-abe@g.ecc.u-tokyo.ac.jp; koyama.shoichi@ieee.org; natsuki.ueno@ieee.org; hiroshi_saruwatari@ipc.i.u-tokyo.ac.jp).
Digital Object Identifier 10.1109/TASLP.2022.3231715 generate complex sound fields with a flexible array geometry of loudspeakers, they have a broad range of practical applications.
Multizone sound field control [11], [12], [13], [14], [15], which is one particular problem in sound field control, aims at generating different sound fields inside multiple target regions for personal audio applications [16], [17]. For example, planewave fields with different propagation angles are synthesized inside two regions. However, this problem is sometimes impractical because its physical feasibility becomes significantly low depending on the desired propagation angles of plane-waves.
Meanwhile, in some applications, it is desired to generate sound fields of certain acoustic power levels rather than synthesizing specific sound fields such as plane-wave fields of particular propagation angles. For example, the acoustic power should be high in one region, but suppressed in another region, focusing on whether the sound is audible or not. In previous studies, this type of problem has been addressed as an acoustic contrast control (ACC) problem [18], [19], where the ratio of acoustic potential energy in one region to that in another region is maximized to generate acoustic bright and dark zones. However, the power distribution inside the target region cannot be controlled by ACC because only the total energy is taken into account. Furthermore, ACC suffers from undesired distortion of filters due to discontinuities between frequencies when applied in the frequency domain. When ACC is applied in the time domain, a flat frequency response cannot be guaranteed.
We consider the problem of synthesizing specific amplitude (or magnitude) distributions inside the target regions by using secondary sources, which is referred to as amplitude matching. For example, a flat amplitude distribution over a target region is set as a desired distribution but its phase distribution is left arbitrary. Such an optimization problem does not have closedform solution; therefore, iterative algorithms are necessary [20]. However, because of the indifferentiability of the cost function of amplitude matching, general gradient methods, such as gradient descent and (quasi-)Newton's methods, cannot be applied in a strict sense. Even if these algorithms work in practice, putting the theoretical validity aside, they basically require a high computational load for computing the inverse of a Hessian matrix (or its proxy) or searching step size parameter at each iteration.
In our previous study [21], we proposed amplitude matching based on the majorization-minimization (MM) algorithm [22], [23], where the optimization problem of amplitude matching is efficiently solved in the frequency domain with a guarantee This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ of a monotonic nonincrease in the cost function. However, this algorithm tends to be stagnated at local minima or saddle points. Furthermore, solving the amplitude matching problem at each frequency can lead to discontinuities in phase between frequency bins as in ACC. As a result, in a broadband case, the time-domain filter to obtain driving signals of the secondary sources can be unnecessarily long.
We propose amplitude matching based on the alternating direction method of multipliers (ADMM) [24] to achieve a more accurate synthesis of the amplitude distribution with an efficient procedure. We also introduce the differential-norm penalty to prevent discontinuities between frequency bins in the broadband case. The proposed amplitude matching method is evaluated by numerical experiments and compared with current multizone sound field control methods and amplitude matching based on MM algorithm. In addition, results of experiments using real data obtained in a practical environment are also presented.
The rest of this paper is organized as follows. In Section II, prior works on multizone sound field control methods are briefly summarized. The problem statement on amplitude matching is described with current algorithms for solving it in Section III. In Section IV, the proposed ADMM-based algorithm is developed. Differential-norm penalty for a time-domain filter design is introduced in Section V. Experimental results in the narrowband and broadband cases are presented in Section VI. Finally, Section VII concludes this paper.

A. Notation
Italic letters denote scalars, lower case boldface italic letters denote vectors, and upper case boldface italic letters denote matrices. The sets of real and complex numbers are denoted by R and C, respectively. Subscripts of scalars, vectors, and matrices indicate their indexes. To illustrate, x i,j is the (i, j)th entry of the matrix X. The imaginary unit and Napier's constant are denoted by j (j 2 := −1) and e, respectively. The complex conjugate, transpose, conjugate transpose, and inverse are denoted by superscripts (·) * , (·) T , (·) H , and (·) −1 , respectively. The absolute value of a scalar x is denoted by |x|. The argument of complex-valued scalar x is denoted by arg(x). The signs | · | and arg(·), and exponential e are also applied to vectors for each element. For example, |x| e j arg(y) means the vector consisting of |x i |e j arg(y i ) with the Hadamard product . The Euclidean norm of a vector x is denoted by x .
The angular frequency, sound velocity, and wavenumber are denoted by ω, c, and k = ω/c, respectively. The harmonic time dependence e −jωt with the time t is assumed according to conventions.

II. PRIOR WORKS ON SOUND FIELD CONTROL
Sound field control is aimed at synthesizing a desired sound field inside a target region. Driving signals of L secondary sources are controlled so that the synthesized sound field corresponds to the desired one inside the target region Ω. The secondary sources are arbitrarily placed in the region outside Ω, as shown in Fig. 1. The target region Ω can consist of multiple non-overlapped regions. The synthesized sound pressure u syn (r, ω) of angular frequency ω at position r ∈ Ω is represented by a linear combination of transfer functions of the secondary sources as where d l (ω) and g l (r, ω) are the driving signal and the transfer function (Green's function) at r for the lth secondary source (l ∈ {1, . . . , L}), respectively. Hereafter, the angular frequency ω is omitted for notational simplicity. The objective of general sound field control is to obtain the driving signals d = [d 1 , . . . , d L ] T ∈ C L so that the following cost function F (d) is minimized: where u des (r) is the desired sound field. When Ω is separated into multiple regions, and different u des (r) values are set inside each region, this minimization problem is specifically called multizone sound field control. Since minimizing (2) is difficult to solve owing to the regional integral, several approximation methods have been proposed.

A. Pressure Matching
Pressure matching is a widely used sound field control method since its practical implementation is simple and flexible. The cost function F (d) is approximated by discretizing the target region Ω into control points [7], [10], [11]. The number of control points and the position of the mth control point are denoted by M and r m ∈ Ω (m ∈ {1, . . . , M}), respectively. We define the vector u des ∈ C M consisting of u des (r m ) and the matrix G ∈ C M ×L consisting of g l (r m ). Thus, the optimization problem for the pressure matching is formulated as where γ is a constant parameter. The first term is the square error between the synthesized and desired sound pressures at the control points, and the second term is the regularization term to prevent an excessive amplitude of d. To make the problem balanced or overdetermined, the number of control points M is generally set to be larger than or equal to L. This optimization problem is solved in a closed form as the regularized leastsquares solution:d The pressure-matching method can be directly applied to the multizone sound field control. Instead of discretizing the target region Ω, the sound field is represented by spherical wave function expansion in modematching methods [2]. The driving signal d is obtained so that the expansion coefficients of u syn (r) correspond to those of u des (r). The weighting factors for each expansion coefficient are introduced in the weighted mode-matching [9] for accurate approximation of the original cost function (2) without empirical truncation of the expansion order. The mode-matching methods are also extended to multizone sound field control [12], [13].

B. Acoustic Contrast Control
Generating regions of high-and low-acoustic power distributions is one particular example of the multizone sound field control. In ACC [18], [19], the target region Ω is separated into The regions Ω B and Ω D are the so-called (acoustic) bright and dark zones, respectively. The goal is to synthesize high acoustic potential energy inside Ω B and low acoustic potential energy inside Ω D .
The optimization problem of ACC is defined as the maximization of the ratio of acoustic potential energy inside Ω B to that inside Ω D as where G B and G D are the transfer function matrices at the control points inside Ω B and Ω D , respectively. The optimal d can be obtained as the eigenvector that corresponds to the largest eigenvalue of the matrix where μ is the regularization parameter. This solution has indefiniteness for the scaling and phase rotation of d. Since only the total acoustic energy inside the target region is considered in (5), it is not possible to synthesize specific power distributions. Therefore, a large variation of acoustic power inside the bright zone can be generated. Furthermore, solving (5) for each frequency leads to discontinuities in phase between frequency bins.

III. AMPLITUDE MATCHING
Pressure matching is applicable to synthesizing various desired sound fields inside the target region. However, in some applications, it is necessary to synthesize a sound field of the desired amplitude inside the target region, whereas the phase of the desired sound field is arbitrary. For example, the sound field of high-and low-acoustic-power distributions should be synthesized similarly to that in ACC, which will be useful in car audio and sound navigation systems. We formulate the optimization problem to achieve sound field control only with the amplitude constraint, i.e., amplitude matching, as where λ is the constant parameter. The first term is the square error between the synthesized and desired amplitude distributions at the control points, and the second term is the regularization term of d. In contrast to the cost function of ACC (5), (7) is designed to measure the error of spatial amplitude distributions. Thus, various amplitude distributions can be set as the desired sound field in the amplitude matching.
Since the optimization problem of the amplitude matching does not have closed-form solution owing to the nonconvexity and indifferentiability of the cost function J(d), iterative algorithms are necessary to solve (7). Gradient methods, such as gradient descent and (quasi-)Newton's method, are widely used nonlinear optimization algorithms. Since J(d) is not differentiable at an arbitrary point, the gradient methods are not applicable to solving (7) in a strict sense, but they can work in practice. The gradient of J(d) can be obtained as where Gd is assumed not to contain zero-valued components. By using (8), we can construct the gradient methods for solving (7), such as gradient descent and quasi-Newton's methods. In addition to the lack of theoretical validity, the computational complexity of these algorithms can be high, as discussed later.

A. MM Algorithm
An amplitude matching method based on the MM algorithm has been proposed by Koyama et al. [21]. In the MM algorithm, a surrogate function for a nonconvex objective function, which can be simply minimized, is constructed. The monotonic nonincrease in the objective function can be guaranteed by alternately updating the variable of the surrogate function and the variable to be optimized. The resulting update rules can be computationally efficient by appropriate construction of the surrogate function.
In [21], the surrogate function of J(d) is obtained as where v (i) is the auxiliary variable defined with the iteration index i as Here, d (i) is d at the ith iteration. This surrogate function satisfies where the equality holds for d = d (i) . See Appendix A for the proof. Thus, the driving signal d is obtained by alternately updating v (i) and d (i) as Algorithm 1: MM Algorithm for Amplitude Matching.
The monotonic nonincrease in the objective function can be confirmed as The iterative update rules are computed until a stopping condition is satisfied, for example, by setting a threshold for the variation of J(d (i) ) or d (i) . The MM algorithm for the amplitude matching is summarized in Algorithm 1. Note that (G H G + λI) −1 G H in the fourth line can be calculated before the iteration. Thus, the computational cost for each iteration is O(LM ). This algorithm corresponds to the Griffin-Lim algorithm [25] and Gerchberg-Saxton algorithm [26] for the phase retrieval problem, whose cost function is similar to that in (7).

IV. ALTERNATING DIRECTION METHOD OF MULTIPLIERS FOR AMPLITUDE MATCHING
The MM algorithm presented in Section III-A is computationally efficient and simple for implementation. However, in some cases, the cost function does not sufficiently decrease, as later shown in the experiments. We propose the amplitude matching algorithm based on ADMM. For convex optimization problems, the convergence to the optimal solution can be guaranteed by ADMM [24]. It has been experimentally validated that ADMM effectively works for some nonconvex optimization problems [24], [27], [28]. In particular, Liang et al. [27] have shown that ADMM performs well in phase retrieval problems.
Since ADMM cannot be directly applied to solve (7), we reformulate the optimization problem by using auxiliary variables of amplitude and phase defined as Gd = a e jθ as where a ≥ 0 means that each element of a is equal to or larger than 0. Then, the augmented Lagrangian function for (15) can be defined as where w ∈ C M is the Lagrange multiplier, [·] represents the real part of complex value, and ρ > 0 is the penalty parameter. In ADMM, each variable is alternately updated, starting with initial values as The augmented Lagrangian function for each update is minimized for one or two variables while fixing the other variables.
To obtain the update rules of a and θ, we modify the subproblem (17) as where h (i) = Gd (i) + w (i) /ρ. In the second line, the terms involving a e jθ are combined into one term. The terms not related to the optimization with respect to a and θ are omitted in the last line. The phase θ is included only in the second term, and choosing θ (i+1) according to (21) is optimal, regardless of a.
Then, a is updated by minimizing the quadratic of a as The update rule of d is obtained by solving where the same modification as that in (20) is used. This quadratic of d can also be simply minimized as Algorithm 2: ADMM Algorithm for Amplitude Matching.
We summarize the proposed ADMM-based amplitude matching algorithm in Algorithm 2. The computational cost of the proposed algorithm for each iteration is O(LM ), which is identical to that of the MM-based algorithm. Among the gradient methods, the computational cost of the gradient descent method is relatively small, which requires O(LM ) for computing the gradient (8) and the cost for the line search of the step-size parameter. However, its convergence speed is generally low. The quasi-Newton's methods, such as the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method, are practically useful because of their fast convergence and low computational cost by approximating the inverse of the Hessian matrix. The computational cost for each iteration is O(L 2 + LM ) with the line search of the step-size parameter. Therefore, the proposed method is computationally efficient as the MM-based algorithm. Further investigation is given in the experiments.

V. DIFFERENTIAL-NORM PENALTY FOR TIME-DOMAIN FILTER DESIGN
In the previous section, the optimal driving signals for amplitude matching are obtained in the frequency domain. In practical implementations, the frequency-domain driving signals of target frequency bins are transformed into the time domain by inverse Fourier transform to construct a single-input-multiple-output (SIMO) finite impulse response (FIR) filter, where the input is the source signal for playback and the output is the driving signals of the secondary sources. The procedure to obtain broadband driving signals is equivalent to solving the following optimization problem: where the index of the frequency bin k (∈ {1, . . . , K}) is introduced in the subscript for each frequency-dependent variable.
Here, |u des | is assumed to be independent of frequency to obtain a flat amplitude response filter for each control point. Since the loss term and the regularization term of 2 -norm for each frequency are simply added, (25) can be solved at every frequency as in Section IV. However, since the driving signals are independently determined for each frequency, the resulting solution can have discontinuities between frequencies. These discontinuities can lead to an unnecessarily large FIR filter length. To achieve amplitude matching with a small FIR filter length, we propose the use of a penalty term for the discontinuities of d k , which is called the differential-norm penalty. We define the differential norm as Thus, the optimization problem of the amplitude matching with the differential-norm penalty is described as The optimization problem (28) can be solved by ADMM as in the single-frequency case. First, (28) is reformulated as The augmented Lagrangian function is defined with the La- The update rule of each variable is obtained as with The update rules ofθ,ā, andw are efficiently computed. On the other hand, the update ofd in (34) requires large computational complexity owing to M −1 . Since the computation of M −1Ḡ H remains unchanged for each iteration, which can be performed before the iteration, which requires O((M + L)L 2 K 3 ), then the computational cost of the update ofd becomes O(LM K 2 ). However, these computational complexities can be reduced by using the property of the block tridiagonal matrix of M . We apply the block birecurrence method [29]

VI. EXPERIMENTS
We conducted experiments to evaluate our proposed method. First, experimental results in the narrowband case are shown to compare the performance of the amplitude matching algorithms. Next, the broadband case is investigated to validate the effectiveness of the differential norm. Finally, we demonstrate experimental results using real data.

A. Narrowband Case
We conducted numerical experiments using 2D free-field simulations in the frequency domain, focusing on the comparison of the amplitude matching (AM) algorithms described in Section III with the conventional pressure matching (PM) and acoustic contrast control (ACC). As shown in Fig. 2, a circular loudspeaker array was placed, whose radius was 1.5 m and the center was at the origin. The array consisted of 48 loudspeakers that were equiangularly aligned. Each loudspeaker was assumed to be the 2D point source, whose transfer function is defined as For the desired amplitude distributions, we set 1.0 in Ω U and 0.0 in Ω L ; therefore, Ω U and Ω L are supposed to be acoustic bright and dark zones, respectively. Hence, ACC was applied to synthesizing bright and dark zones in each region. The amplitude of the driving signal for ACC was determined so that the average amplitude inside Ω U becomes 1.0.
To solve the optimization problem of AM, we compared the proposed ADMM-based algorithm (AM-ADMM), MM-based algorithm (AM-MM), and gradient method using BFGS (AM-BFGS). PM and ACC were also evaluated. Since the experiments were conducted in the frequency domain, we used the frequencydomain ADMM algorithm (Algorithm 2) in this section. The initial value of d for the three AM algorithms was the solution of PM. The regularization parameter λ in (7) for the AM algorithms and γ in (5) for ACC were set as σ max (G) × 10 −3 , where σ max (G) means the largest singular value of G. The parameter ρ in Algorithm 2 was set to 1.0.
As an evaluation measure, we define the mean square error (MSE) between synthesized and desired amplitude distributions as where |u syn (ω)| is the amplitude distribution of the synthesized sound field. We also evaluate with acoustic contrast (AC) defined as AC(ω) = 10 log 10 u syn where u syn Ω U and u syn Ω L are the synthesized pressure in Ω U and Ω L , respectively. The experiments were conducted using a 3.00 GHz 8-Core Intel Core i7-9700F processor and implemented with Python.
An appropriate setting of u des in (3), i.e., the desired phase distribution, for PM is not obvious, especially when the desired Fig. 3. MSE with respect to plane wave direction φ pw of the desired field in Ω U at 1400 Hz for narrowband case. Each PM solution was used as the initial value of AM-ADMM. The case when u des is set as the desired amplitude distribution with zero phase is also shown, which is denoted as PM (zero phase) and AM-ADMM (zero phase).
amplitude distribution and arrangement of the target regions are complicated, e.g., the case of non-uniform |u des |. However, in this setting, it can be inferred that the desired sound field in Ω U for PM should be a plane wave propagating in the direction of 0 • or 180 • . The setting of u des for PM also has an impact on the AM methods because the solution of PM is used for the initial value in this experiment. Fig. 3 shows the MSE of PM with respect to the plane wave direction φ pw of the desired sound field and AM-ADMM with the initial value of each PM solution. The case that u des is set as the desired amplitude distribution with zero phase, i.e., 1.0 in Ω U and 0.0 in Ω L , was also investigated, which is denoted as PM (zero phase) and AM-ADMM (zero phase). Although synthesizing a uniform pressure distribution of zero phase is impractical, this setting will be useful for initial value of the AM methods when an appropriate setting of u des is difficult. Note that the results of PM (zero phase) and AM-ADMM (zero phase) are not dependent on φ pw . As expected, the lowest MSE of PM was achieved at φ pw = 0 • and 180 • , whose value was −43.3 dB. The MSE of PM (zero phase) was −4.2 dB, which was larger than the highest MSE of PM (−12.0 dB at φ pw = 90 • and 270 • ). In AM-ADMM, the lowest MSE was −40.2 dB at φ pw = 0 • and 180 • , which was slightly larger than that of PM at the same φ pw . However, the MSE of AM-ADMM was maintained at around −36 dB for φ pw ∈ [45 • , 135 • ] and [225 • , 315 • ], which was close to the MSE of AM-ADMM (zero phase), i.e., −36.4 dB. AM-ADMM attained low MSE even when the initial value was the PM solution with the high MSE value; therefore, PM (zero phase) is a reasonable compromise for the initial value of the AM methods. Table I shows the MSE and AC of each method at 1400 Hz. PM (φ pw = 0 • ) and AM-ADMM (φ pw = 0 • ) indicate PM with the desired field of the plane wave with the optimal angle φ pw = 0 • and AM-ADMM with the initial value of the solution of PM (φ pw = 0 • ), respectively. In the other AM methods, PM (zero phase) was used for the initial value. The synthesized pressure distribution of each method is shown in Fig. 4. In ACC, the highest AC was achieved with a low acoustic potential energy in Ω L , but the amplitude distribution in Ω U was not uniform. Thus, even though the MSE of ACC was very low in Ω L , that of ACC was high in Ω U . The plane-wave-like pressure distribution was synthesized by AM-ADMM. In AM-MM and AM-BFGS,  the synthesized pressure distributions were still close to that in PM (zero phase), i.e., the initial value. The total MSE was the lowest in AM-ADMM among the methods not using the optimal plane wave angle φ pw = 0 • .
Figs. 5 and 6 show the MSE and AC of each method with respect to the frequency from 100 to 4000 Hz at intervals of 50 Hz. Since the optimization problem of AM was solved for each frequency, the MSEs of the AM algorithms largely varied at each frequency. However, they were significantly lower than those of PM (zero phase) and ACC. In contrast, the highest AC was achieved by ACC because this measure corresponds to the objective function of ACC. Among the three AM algorithms, the MSEs of AM-MM were relatively higher than those of AM-ADMM and AM-BFGS, especially at the frequencies below 1400 Hz. This means that AM-MM did not reach the optimal   solution owing to its initial value dependence at these frequencies. The performance of AM-ADMM was comparable to that of AM-BFGS, but high MSEs appeared at several frequencies. Nevertheless, their general performance was better than that of the other methods.
Next, we investigated the case when the desired amplitude distribution is not uniform. The desired amplitude |u des | in Ω U is set to be a truncated Gaussian function with its mean at the center of Ω U and variance 0.04. The desired field u des for PM is difficult to set, and ACC is not applicable to synthesize a specific phase distribution; therefore, PM (zero phase) and AM methods were compared. In PM (zero phase), the desired pressure distribution was set to be |u des | with zero phase. Fig. 7 shows the MSE of each frequency. The desired amplitude distribution was accurately synthesized by the AM methods although their general MSEs were relatively high, compared to those in Fig. 5.
Finally, we demonstrate the case that the desired amplitude distributions are more complex. We set three circular target regions of 0.3 m radius with their centers at We denote these upper left, Fig. 8. MSE with respect to φ pw at 1400 Hz for narrowband case when three target regions were set, where φ pw denotes the plane wave angle of the desired field in Ω UL . Each PM solution was used as the initial value of AM-ADMM. The case when u des is set as the desired amplitude distribution with zero phase is also shown, which is denoted as PM (zero phase) and AM-ADMM (zero phase). upper right, and lower center target regions as Ω UL , Ω UR , and Ω LC , respectively. The desired amplitude distributions were uniform: 1.0 in Ω UL , 0.5 in Ω UR , and 0.0 in Ω LC . ACC was applied by separately obtaining the driving signals for synthesizing a single bright zone in Ω UL or Ω UR ; then, these signals were normalized and summed. Since the desired amplitude distributions are uniform, plane wave fields will be a reasonable setting for the desired field of PM; however, it is difficult to infer physically feasible plane wave angles because of complex arrangement of the target regions. Fig. 8 shows the MSE with respect to the plane wave angle φ pw in Ω UL at 1400 Hz when the desired field in Ω UR was fixed to be a plane wave field of φ pw = 210 • . The direction φ pw = 210 • in Ω UR was taken from the propagation direction of the sound field synthesized by AM-ADMM. The MSE of AM-ADMM was lower than that of PM for all φ pw . AM-ADMM (zero phase) attained almost the same or even lower MSE of AM-ADMM. The MSE of each method at 1400 Hz is shown in Table II. The MSE with respect to frequency is shown in Fig. 9.
Since AC cannot be defined for three regions, we only evaluated with MSE. Again, the lowest MSE was achieved by AM-ADMM and AM-BFGS at most of frequencies. The synthesized pressure distributions at 1400 Hz are shown in Fig. 10. These results indicate that our proposed algorithm performs well even when the desired amplitude distributions are complex. Fig. 11 shows the cost function value of the AM algorithms in terms of computation time. The decreasing speed of the cost function was high in AM-MM and AM-ADMM, but it stagnated at around 6.7 in AM-MM. The cost function decreased to 3.2 in Fig. 9. MSE with respect to frequency for narrowband case when three target regions were set.   AM-ADMM and 2.7 in AM-BFGS. The computation time until the cost function reached 3.5 was 56.8 ms in AM-BFGS and 7.9 ms in AM-ADMM. A similar tendency can be seen at other frequencies. Therefore, the proposed ADMM algorithm is effective both in control accuracy and computational time.

B. Broadband Case
Next, we show experimental results in the broadband case, focusing on the effectiveness validation of the differential-norm penalty for AM. Again, numerical simulation is conducted under the 2D free-field assumption. The placements of loudspeakers, target regions, and control points, and the settings of the desired amplitude were the same as those in Section VI-A (see Fig. 2).
We compare AM with the differential-norm penalty (28) and 2 -norm regularization (25), which are denoted as AM w/ Diff.-norm and AM w/ 2 -norm, respectively. The parameter λ in (28) was set to 10.0, and the other parameters were the same as those in AM-ADMM in the narrowband case. After generating the frequency-domain driving signals using the ADMM-based algorithm, time-domain FIR filters were obtained by their inverse fast Fourier transform (FFT). The sampling frequency was 16 kHz, and the number of frequency bins, i.e., FFT length, was 32768. For evaluation measures in the broadband case, we define average MSE and AC for all the frequency bins as MSE = 10 log 10 where F is the number of frequency bins, f is its index, and u syn Ω B and u syn Ω D are synthesized pressure in bright and dark zones, respectively. The synthesized pressures in the time domain at the evaluation points were transformed into the frequency domain to compute these evaluation measures. Fig. 12 shows the FIR filters of AM w/ Diff.-norm and 2 -norm before truncating the filter length. Apparently, the decrease in the filter amplitude from the peak of the AM w/ Diff.-norm is much sharper than that of 2 -norm. When the optimal filter length is defined so that the decrease from the peak amplitude is 30 dB, that of AM w/ 2 -norm was 12888. On the other hand, the optimal filter length of AM w/ Diff.-norm was 162. The magnitude and phase responses of the time-domain filters of AM w/ 2 -norm and AM w/ Diff.-norm are shown in Figs. 13 and 14, respectively. The discontinuities in the frequency response were suppressed by the differential norm penalty. Table III shows the MSE and AC of each method when the filter length is truncated to 220. This filter length is defined so  that the decrease in filter amplitude from the peak is 40 dB in AM w/ Diff.-norm. MSE in the frequency domain is also plotted in Fig. 15. MSE of AM w/ Diff.-norm was particularly lower than that of 2 -norm inside Ω U . The frequency-domain MSE of AM w/ 2 -norm largely fluctuated for each frequency bin at high frequencies, whereas that of Diff.-norm was relatively smooth.

C. Experiments Using Real Data
To investigate the performance of the proposed method in a practical environment, we conducted experiments using the   impulse response dataset presented in [30]. Configurations of loudspeakers and control points are shown in Fig. 16. We set two 2D square target regions of acoustic bright and dark zones,    First, the filters generated by AM w/ 2 -norm and Diff.-norm are compared. The parameters for ADMM-based algorithms were the same as those in Section VI-B. Fig. 17 shows the filters before truncation obtained by AM w/ 2 -norm and Diff.-norm. Since the impulse responses of the loudspeakers include reverberation, the difference in the necessary filter length for these two methods was much larger than that in Section VI-B.  Fig. 20. In Ω D , the magnitude distribution of the proposed method was low even at the positions other than the control points, compared with those of the other methods.

VII. CONCLUSION
We proposed the amplitude matching method for multizone sound field control. The amplitude matching aims to synthesize a desired amplitude distribution over a target region, whereas the phase distribution is arbitrary, using multiple secondary sources. Since the cost function of amplitude matching is neither linear, convex, nor differentiable, it is important to develop efficient algorithms for solving the amplitude matching problem. An algorithm based on ADMM is formulated both in the frequency and time domains. To avoid an unnecessarily large time-domain filter length in the broadband case, the differential-norm penalty for adjacent frequency bins is also introduced. In the experiments, for k = k b − 1, k b − 2, . . . , 1, and Takumi Abe received the B.E. degree in mathematical engineering and information physics in 2021 from the University of Tokyo, Tokyo, Japan, where he is currently working toward the master's degree. His research interests include haptics and ultrasonic wave. Natsuki Ueno (Member, IEEE) received the B.E. degree in engineering from Kyoto University, Kyoto, Japan, in 2016, and the M.S. and Ph.D. degrees in information science and technology from the University of Tokyo, Tokyo, Japan, in 2018 and 2021, respectively. He is currently a Project Assistant Professor with Tokyo Metropolitan University, Tokyo, Japan. His research interests include spatial audio and acoustic signal processing.
Hiroshi Saruwatari (Member, IEEE) received the B.E., M.E., and Ph.D. degrees from Nagoya University, Nagoya, Japan, in 1991, 1993, and 2000, respectively. In 1993, he joined SECOM IS Laboratory, Tokyo, Japan, and in 2000, Nara Institute of Science and Technology, Ikoma, Japan. Since 2014, he has been a Professor with The University of Tokyo, Tokyo, Japan. His research interests include statistical audio signal processing, blind source separation, and speech enhancement. He has put his research into the world's first commercially available independentcomponent-analysis-