Directional sound source modeling using the adjoint Euler equations in a finite-difference time-domain approach

: An adjoint-based approach for synthesizing complex sound sources by discrete, grid-based monopoles in ﬁnite-difference time-domain simulations is presented. Previously, Stein, Straube, Sesterhenn, Weinzierl, and Lemke [( 2019 ). J. Acoust. Soc. Am. 146 (3), 1774–1785] demonstrated that the approach allows one to consider unsteady and non-uniform ambient conditions such as wind ﬂow and thermal gradient in contrast to standard methods of numerical sound ﬁeld simulation. In this work, it is proven that not only ideal monopoles but also realistic sound sources with complex directivity characteristics can be synthesized. In detail, an oscillating circular piston and a real two-way near-ﬁeld monitor are modeled. The required number of monopoles in terms of the sound pressure level deviation between the directivity of the original and the synthesized source is analyzed. Since the computational effort is independent of the number of monopoles used for the synthesis, also more complex sources can be reproduced by increasing the number of monopoles utilized. In contrast to classical least-square problem solvers, this does not increase the computational effort, which makes the method attractive for predicting the effect of sound reinforcement systems with highly directional sources under difﬁcult acoustic boundary conditions.


I. INTRODUCTION
Accurate modeling of the directivity of complex sound sources is a prerequisite for many applications in technical acoustics, sound reinforcement, sound field synthesis, and auralization. If a mechanical model of the source is available, the problem can be solved by finite element method (FEM) simulations of the entire system (Karjalainen and Savioja, 2001). More complex composite sources can be designed by a combination of given generic elements and filters, once an underlying speaker model is known (Feistel et al., 2009;van Beuningen and de Vries, 1997).
Grid-based approaches, on the other hand, do not require any a priori knowledge of the source but try to fit the coefficients of a synthetic source such that the resulting sound field optimally reproduces a given reference. This can be considered as an inverse optimization approach, where, typically, some kind of least-squares problem needs to be solved.
Traditionally, these problems are formulated in the frequency domain (Mehra et al., 2014;Ochmann, 1992). As a least-square fitting algorithm, a pseudo-inverse can be calculated using a singular value decomposition (SVD) (Wang and Wu, 1997). In doing so-inverting an ill-conditioned system-special care such as a Tikhonov regularisation has to be taken (Kim and Nelson, 2004).
Alternatively, grid-based least-squares problems of source synthesis can be based on a finite-difference timedomain (FDTD) method and solved in the time domain, for example, by determining a minimum eigenvector with an Arnoldi type algorithm (Takeuchi et al., 2019) or by obtaining a least-square fit by calculating a Moore-Penrose pseudo-inverse (Bilbao et al., 2019;Escolano et al., 2007). Other comparable finite-differences sound source modeling techniques are the pseudospectral time-domain method (Georgiou and Hornikx, 2016).
In contrast to previous approaches of grid-based source synthesis, the present method is based on the Euler equation as a more general form of the wave equation. This generalized form allows one to model the behavior of acoustic sources in situations with complex background flow and thermal stratification, relevant for applications such as ventilation systems with thermal gradients, sound reinforcement in open spaces such as train stations or sports stadiums, or noise pollution in urban environments. In a previous study (Stein et al., 2019) we have demonstrated how to model and optimize sound sources, for example, for open-air concerts with a non-uniform velocity and temperature background. While in that study only idealized spherical sources with uniform directivity were used, the present study shows that the adjoint FDTD a) Electronic mail: stein@cfd.tu-berlin.de, ORCID: 0000-0002-4298-2001.
Whereas other least-square fitting algorithms such as the SVD, the Arnoldi iteration, or the Moore-Penroseinverse scale quadratically or worse with the number of synthetic source coefficients to be optimized, the cost of the proposed adjoint-based optimization method depends on the domain size only but not on the number of coefficients, i.e., the number of tuned grid-based monopoles inside the domain. An adjoint-based source modeling therefore not only has a wider range of applications concerning the acoustic conditions, but it can also be more efficient as the complexity of the source directivity increases, i.e., when more monopoles need to be used.
The paper is structured as follows. Section II introduces the adjoint-based concept of source modeling. In Sec. III, we set up the underlying FDTD simulation utilized in the analysis. Section IV yields the source synthesis of a circular piston, while in Sec. V, a realistic speaker is modeled based on its measured directivity. Section VI draws conclusions and discusses their impact.

II. ADJOINT-BASED SOURCE MODELING
In the following, we present an adjoint-based optimization approach as a novel method to model complex sound sources. The aim is to synthesize acoustic monopole sources, which reproduce the sound field of the original source in an optimal sense.

A. Governing equations
The governing equations for the analyzes are the Euler equations (Landau and Lifshitz, 1987) in the time domain expanded by a source term s p ðx j ; tÞ on the right side of the pressure (energy) equation to model sound sources. Thus, on each grid point of the computational domain, a monopole source is modeled. In general, mass and momentum sources could also be used. They are not considered here but are the target of future investigations. The variable . denotes the density, u i the velocity in x i -direction, p the pressure, and c the heat capacity ratio. For the indices i and j the summation convention applies over 1 to 3. For details on the formulation of the energy equation in terms of pressure, see Lemke et al. (2014). The Euler equations were chosen as the most general governing equations for acoustic problems, including nonuniform base flows and, in perspective, reflecting boundaries. However, in the following only free space propagation examples in an environment at rest are presented. To consider more complex boundary conditions, such as wall impedance, the direct and the adjoint solver need to be extended by corresponding boundary conditions which is beyond the scope of this work.
The source term sðx j ; tÞ is used to modify the solution of the governing Euler equations such that a given reference sound field is reproduced in space and time. By this, an optimal source distribution is synthesized.
The optimum is defined in terms of the objective function in which pðx j ; tÞ denotes the pressure part of the solution of the governing equations and p ref ðx j ; tÞ a given or desired (sound) field, while p 0 and p 0 ref denote acoustic fluctuations in the sense of p 0 ¼ p À p 1 . The integration is processed for the whole computational domain in space and time with dX ¼ dx 1 dx 2 dx 3 dt. The weight rðx j Þ controls where the objective is evaluated. An optimal source synthesis is found if J reaches a minimum.
To minimize the objective function the highdimensional gradient is required. The degrees of freedom of s p correspond to the number of grid points multiplied by the number of considered time steps. Different methods for computation of the gradient are available. A computationally efficient technique to determine the gradient is an adjoint-based approach well known from different fields, e.g., meteorology (Courtier et al., 1994) and fluid dynamics (Giles and Pierce, 2000).

B. Adjoint approach
In the following, the adjoint approach is introduced by a simple example in a discrete manner. The section is based on Giles and Pierce (2000) and Lemke (2015). A matrixvector notation is used.
The adjoint equation arises by the scalar-valued objective function J, such as Eq. (2). In the matrix-vector notation, it is given by the scalar product between a user-defined weight vector g and a system state vector q (4) The system state vector q is the solution of the governing system A q ¼ s; with A as governing operator and s as right-hand side forcing. To optimize J by means of s in terms of finite differences the governing equation has to be solved multiple times corresponding to the number of entries in s.
To reduce the computational effort, the adjoint equation can be introduced with the adjoint variable q Ã . Thus, a formulation is found, which allows the computation of the objective J without solving the governing system for different s. Instead, the objective can be calculated by a scalar product using the solution of the adjoint equation (6).

C. Adjoint Euler equations
In the following the derivation and use of the adjoint Euler equation are shown. Therefore, the set of governing equations (1) is abbreviated by Therein, the variable qðx j ; tÞ comprises all primitive variables of Eq. (1), i.e., q ¼ ½.; u 1 ; u 2 ; u 3 ; p, and s the source term s p by means of s ¼ ½0; 0; 0; 0; s p .
To derive the corresponding adjoint, the equation and the objective function has to be linearized, The latter is expressed in a more general form using q. However, only the pressure part is evaluated in the following. Combining the linearized governing equations and objective in a Lagrangian manner leads to with subsequent partial integration from Eq. (11) to Eq. (12) to factor out dq. For the sake of simplicity, the integrals a T b :¼ Ð abdX are not shown. Instead, the scalar products indicated by the transpose reminds of the integration. The partial integration gives rise to the adjoint boundary and initial conditions which are discussed in Lemke (2015) in detail. Here, the initial adjoint solution is given by q Ã init ðx j Þ ¼ 0, and non-reflecting adjoint boundary conditions are applied.
The adjoint equation results from demanding in Eq. (12) to be zero. By this dJ simplifies to similar to Eq. (7). The required solution of the adjoint equation can be interpreted as a high-dimensional gradient of J with respect to the source terms s Expanding Eq. (13) the full adjoint Euler equations read

D. Iterative framework
The adjoint-based gradient [Eq. (15)] is used with an iterative procedure, see Fig. 1. First, the governing equations 1 are solved, with s ½0 ¼ 0 as an initial guess. The superscript in rectangular brackets stands for the iteration step. Another initialization for s is also possible. Subsequently, the adjoint equations (13) are solved. Based on the latter solution, the required gradient r s J is determined and used to update the source distribution s ½n s nþ1 with a denoting an appropriate step size, n the iteration number, and hðx j Þ as spatial weight defining the permitted source volume. The procedure is repeated until a user-defined convergence criterion is reached, here 15 loops. Measurement errors included in the reference q ref but contradicting the physics described by the Euler equations are not reproduced by s. By construction of the adjoint-based approach-insertion of Eq. (8) in Eq. (11)-these artificial contributions are filtered out naturally. This build-in filtering by the governing equations is another advantage of the approach.
The presented procedure optimizes to local extrema only. The identification of a global optimum is not ensured. The computational costs scale linearly with the spatial and temporal resolution of the computational domain only, as the adjoint is solved for the whole computational domain and time. In particular, the costs are independent of the number of sources and their arrangement, which is beneficial if numerous grid-based monopoles are employed for the source synthesis. No information about the underlying loudspeaker is needed, but only the reference sound field p ref . The computational problem is fully parallelizable. The approach works equally well if only experimentally measured signals at specific locations are available (Gray et al., 2017; Lemke, 2015; Lemke and Stein, 2020).

III. FDTD SETUP AND AMBIENT CONDITIONS
Our adjoint-based source modeling method (Sec. II) relies on the FDTD solution of Eqs. (1) and (13). For both the same discretization is used. We discretize spatial derivatives with a compact scheme of sixth order (Lele, 1992). For time-wise integration, we employ the explicit Runge-Kuttascheme of fourth order. We consider an unbounded domain with non-reflective boundaries, realized by characteristic boundary conditions (Poinsot and Lele, 1992;Stein, 2018) and a quadratic sponge layer (Mani, 2012). Throughout this work, a cubical domain is resolved by equidistantly distributed points with seven gridpoints per 3 kHz wave. It includes the entire space from the loudspeaker at the center at ½0; 0; 0 m to the most distant microphone position of interest. We are assuming an ideal gas with a heat capacity ratio of c ¼ 1:4. The ambient density corresponds to a speed of sound of 343 m/s and a pressure p 1 of 101 325 Pa.
The masking function of the monopole source region h [Eq. (17)] and the mask of the objective area r [in Eqs. (2) and (10)] are composed of error instead of Heaviside functions, to avoid steep gradients in the excitation of the governing and the corresponding adjoint equations; covers the objective area. The r region starts at a radius r J and is limited outwards by a cube. x i;s and x i;e mark the effective start and end of the cube. l ¼ 20 mm is the selected transition length of the mask edges; masks the monopole source region with a radius r src . Figure 2 presents the entire three-dimensional isosurfaces of the r objective and h source masks according to Eqs. (18) and (19), respectively.

IV. RESULTS PART 1: CIRCULAR PISTON SYNTHESIS
As the first result, we present the reproduction of an ideal circular piston. This serves as the first proof, that the adjoint-based source modeling concept can capture general directivity characteristics.

A. Definition of the circular piston reference
A circular piston in a rigid baffle is an idealized model of a combined woofer and midrange speaker of a two-way system. Its resulting far-field p 0 ref can be calculated analytically with a standard complex directivity point source model (Straube, 2019). Later it serves as reference field p ref ¼ p 1 þ p 0 ref of our adjoint-based optimization in Eq.
(2). The loudspeaker model consists of the following parts. First, the Green's function describes the wave propagation in free three-dimensional space. r is the distance between the origin of the point source and the observer position. In the present work, the origin of the speaker is always at ½0; 0; 0 m. k ¼ 2p f =c is the angular wavenumber for the frequency f, being c the ambient speed of sound and i the imaginary unit. Second, the directivity of the point source according to Skudrzyk (2013) [Chap. 26, Eq. (42) The piston surface is orientated in the x 2 À x 3 plane, being x 1 the main radiation axis. J 1 is the Bessel function of first kind and first order. b is the angle between the main radiation axis (x 1 ) and the line connecting the observer position and the origin of the speaker (see Fig. 3). On the main radiation axis (b ¼ 0), the directivity is lim x!0 2 J 1 ðxÞ=x ¼ 1. The radius of the speaker membrane, i.e., the piston radius, is set as r mem ¼ 82 mm. Convolution of the total transfer function with the test signal (audio file) a(t) provides the resulting sound field at a given observer location and time Our first test signal is a sine-sweep, with a constant amplitude and a frequency linearly increasing from 1 to 3 kHz. Below 1 kHz and above 3 kHz, the sweep amplitude is smoothly faded to zero. Excluding the sponge layer zone, the effective domain-considered in the following-is one cubic meter. The total underlying resolution is 197 Â 197 Â99 gridpoints in space and 1000 steps in time with a sampling frequency of 48 kHz. As radius r src of the source region Eq. (19) we select the same size as the radius r mem ¼ 82 mm of the circular piston model. Equation (18) defines the objective area with r J ¼ 175 mm; x i;s ¼ À0:4 m, and x i;e ¼ 0:4 m.
After 15 iteration loops (refer to Fig. 1), the difference between the reproduced sound field and the reference field in the objective area reached convergence. Per loop, the local 16-core workstation (AMD EPYC 7351) took 67 min on average. Below, we spectrally analyze the time-series measured by the virtual microphones introduced by Fig. 3: Fig. 4 depicts the narrowband gain at the three microphone positions resulting from the adjoint approach. We normalize the gain at the loudest frequency of the first microphone (at b ¼ 0 ), On the main, i.e., loudest, radiation axis (b ¼ 0 ), the amplitude of the linear sweep is almost perfectly constant within its entire frequency interval from 1 to 3 kHz. As expected, with increasing b, the gain decreases among the three curves of the circular piston model. With increasing frequency, more side-lobes also appear at smaller b values. The gain drop reflects this at 3.6 kHz for the b ¼ 45 microphone.
To evaluate the reproduced sound field p 0 opt further, we compare it to the reference field p 0 ref . Figure 5 presents the sound pressure level (SPL) differences at the microphone positions. Throughout this work narrowband SPL differences are The maximum SPL deviation is less than 61 dB near 1 kHz and 60.5 dB between 1.2 and 3 kHz for all three microphones. These deviations are inside the expected range, comparable to Stein et al. (2019). A possible explanation of the slightly larger deviation at the low-frequency edge of 1 kHz is the corresponding wavelength of 0.343 m. This big wavelength has the same order of magnitude as the computational domain and such might suffer from cut-off effects and a less functional sponge layer boundary treatment. Figure 6 depicts the phase differences of the sound pressure signal at the microphones. The most significant difference corresponds to a phase shift of less than 0.8% (phase shift in rad = 2p) within the considered frequency interval.
So, besides the magnitude also the phase between all microphone signals is correctly captured by our optimization algorithm.

C. Directivity characteristics of the circular piston
This section discusses mono-frequent test signals as the input of the circular piston model to study frequency components individually. The setup is the same as in Sec. IV B. The only difference is that the test signal is resolved by 800 instead of 1000 time steps, reducing the average computation time per loop from 67 to 55 min.
First, we demonstrate that the adjoint-based source synthesis reproduces similar directives like the reference. Due to the rotational symmetry around the x 1 -axis and a mirror symmetry at the x 2 À x 3 plane, it is sufficient to present here x 2 À x 1 slices for positive x 1 only. Figure 7 displays the relative SPL directivity for three different test cases calculated with a pure harmonic 2 and 3 kHz signal, respectively. At 2 kHz we find a maximum difference of À0.3 dB rel (b ¼ 0 ); at 3 kHz the greatest deviation is þ0.3 dB rel (b ¼ 90 ).
Quantitatively these deviations of directivity fully comply with practical sound reinforcement or wavefield synthesis needs. The error is comparable or even smaller to other synthesis approaches (Hamilton and Bilbao, 2017).

D. Spectral source patterns of the circular piston
An interesting question is: How did the adjoint-based monopole synthesis manage to reproduce p 0 ref ? To answer this question, we study the underlying source distribution s p of monopoles, which was created by our optimization algorithm. Also, we want to investigate how the size of the source region (h mask) and the frequency affects the quality of the reproduction.
The following discussion is based on 6 Â 9 different cases, in which the extension of the spherical source region and the frequency of the test signal are varied. Six different radii r src ¼ ½5; 6; …; 10 Dx 1 ¼ ½41; 49; 57; 65; 73; 82 mm of the source region (h of Fig. 2) and nine different frequencies f ¼ ½1; 1:25; 1:5; …; 3 kHz are considered. Figure 8 shows a matrix of these 6 Â 9 cases. For each case, the source distribution hs p i of an x 1 -x 2 region with 204 Â 204 mm (25 Â 25 gridpoints) is depicted. The radius r src ¼ 10 Dx 1 of the least compact source region is nearly as big as the smallest considered wavelength (3 kHz): r src =ðc=f Þ ¼ 0:71.
Like p 0 , s p Dt is also a small scale fluctuation around zero in Pa [see Eq. (1) for the unit]. A simple temporal average of absolute s p values discards phase information of the monopole distribution. Therefore, we define hs p i ¼ hsign s p ð40 mm; 0 mm; 0 mm; tÞ Â Ã s p ð:; :; :; tÞi t : Normalization with the s p -sign at ½40; 0; 0 mm ensures that at this gridpoint the argument inside the brackets of Eq. (25) is positive for all times. The algorithm successfully transferred the symmetry of the sound field to the source distribution of Fig. 8, as expected. At low frequencies near 1 kHz, the source distribution is Gaussian (circular) with a uniform phase. Around 2 kHz, the radiation in the x 2 direction (b ¼ 90 ) is reduced, by placing some monopoles with anti-phase along the x 2 axis. Near 3 kHz, the presence of the first side lobe is reflected by the separation of two visually distinguishable vertically aligned ellipsoids (black, in anti-phase to the main lobe at b ¼ 0 ). The adjoint-based optimization always utilizes the entire available region provided by the hðr src Þ sphere to reach the optimization target. Thus, the spatial distribution of active monopoles increases cubically with r src .
To assess the overall reproduction quality with a single scalar number, the use of the objective function J [Eq.
(2)] is evident. We use J ½1 to normalize the objective function, Strictly speaking,J ½n is only a measure of how strong the algorithm reduces the objective function in n loops. A real convergence criterion should be proportional to the local minimum, which ideally can be reached. To this end, we definẽ J conv: ¼J n conv: ½ being n conv: the first n wherẽ J n conv: À1 ½ ÀJ n conv: ½ < 0:007: (27) Figure 9 shows theJ conv: values for the same 6 Â 9 radius-frequency matrix as before. Larger r src lengths guarantee the best results, especially important in case of high frequencies. Near 1 kHz small source regions of r src ¼ 5 Dx 1 are sufficient. The monopol-like directivity below 1 kHz is simply captured. The greatest reduction to 1% is reached at 3 kHz and r src ¼ 10 Dx 1 of Fig. 9. The least reduction to 11% is visible at 3 kHz and r src ¼ 5 Dx 1 . With a smaller number of available monopole sources, the more detailed side-lobe pattern at high frequencies is harder to capture.
In summary, we found that a spherical source region with a radius of ten gridpoints is adequate to reproduce the directivity of the circular piston reference (Sec. IV A) up to 3 kHz, using our adjoint approach. Near or below 1 kHz, already a radius of five gridpoints is adequate to model the circular piston. The largest appearing SPL error is 1.5 dB (see Fig. 7), while above 1.2 kHz, the error is always below 0.5 dB (see Fig. 5).

V. RESULTS PART 2: SYNTHESIS OF REAL-LIFE SOURCES
In the following, we demonstrate that our source synthesis method can model the directivity of a Genelec 8020c 2-way near-field monitor as an example for real-life sources. In contrast to Sec. IV, its magnitude and phase directivity is less smooth and holds no clear symmetry.

A. Origin of the experimental reference
The impulse response of the Genelec 8020c monitor was measured in the semi-anechoic chamber at RWTH Aachen on a 2 Â 2 equi-angular sampling grid using an exponential sweep of order 14 and a G.R.A.S. 40AF half-inch free-field microphone placed at a distance of 2 m from the speaker. In post-processing, the common propagation delay was removed and a subsonic high-pass filter was applied (fourth order Butterworth À3 dB with 30 Hz). The impulse responsesmeasured with 44 100 Hz-were truncated to 9.5 ms to discard reflections by applying a 2 ms Hann window fade in/out. A discrete spherical Fourier transform according to Rafaely (2015) [Eq. (3.34)] with a spherical harmonics order of 20 was applied to the magnitude and unwrapped phase spectra. The data are publicly available in Asp€ ock et al. (2019).
To obtain the impulse response Hðx 1 ; x 2 ; x 3 ; f Þ at all necessary locations we did a range extrapolation (Brinkmann and Weinzierl, 2017) of the spherical harmonic coefficients soft-limited by 40 dB following Bernsch€ utz et al. (2011). Thus, for a given test signal a(t), the sound field of the Genelec speaker is Here, G 3D and D of Eq. (22) became the impulse response H.

B. Genelec monitor synthesis with linear sweep test signal
As a test signal, we use the same 1-3 kHz sine-sweep as in Sec. IV B. To avoid errors caused by a huge range extrapolation, we increased the utilized computational domain to 2.6 m in all room directions. Subsequently, we can place the objective mask Eq. (18)-the area where the difference between the reference and the reproduced field is minimized-at a minimum radial distance of r J ¼ 0:9 m with respect to the source center at ½0; 0; 0 m. x i;s ; x i;e are set 61.2 m, respectively. To enclose the Genelec 8020c dimensions of ½143; 242; 151 mm we confine the region utilized for the monopole synthesis with r src ¼ 0:17 m [Eq. (19)]. In total the spatial domain is resolved with 197 Â 197 Â 197 points and the temporal domain with 974 steps at 44 100 Hz. With this resolution, the 16-core workstation (AMD EPYC 7351) provided one loop every 190 min.
To compare our results with the reference, we pick three microphone positions inside the objective mask r [Eq. (18)]. The position r, b in the x 1 -x 2 plane with x 3 ¼ 0 (r, b are defined as in Fig. 3) are given in the legend of the following graphs. Figures 10 and 11 show the SPL and the phase differences between our model and the measurement (Sec. V A).
The maximum SPL deviation for the microphones is 1.6 dB. At low frequencies, e.g., at 1 kHz, the deviation decreases to 0.9 dB. The phase deviation displays shifts of less than eight percent in the whole frequency interval [phase shift in radian=ð2pÞ]. Below 2.5 kHz the phase shift is just 1%. At low frequencies, the reproduction is more accurate since the sound field is less complex. Please note, 3 kHz is the bass/treble crossover frequency of the monitor.

C. Directivity characteristics of the Genelec monitor
Next, we repeat the setup of Sec. V B but replace the sweep with mono-frequent test signals. This allows us to reduce the number of timesteps to 674, corresponding to 113 min of computation time per loop. We analyze the spectral components of the directivity over all b angles, instead of selecting three microphone positions only. The following polar plots are all located in the x 1 -x 2 plane at a radial distance of r ¼ 1 m. We normalized both the SPL and the phase by the maximum value. Figure 12 depicts the SPL directivity. At 2 and 3 kHz, the maximal difference between the experimental reference and our reconstruction is 1.8 and 10.4 dB, respectively. These maximal differences, however, appear only on the rear side of the monitor where absolute SPL levels are 15 dB below the maximum SPL value. The adjoint approach gives less prominence to regions, where relatively small dB values contribute only minimally to the overall objective function Eq. (2). As mentioned before, 3 kHz is especially critical due to the interference of the bass/treble drivers at the crossover frequency. In the main radiation direction of b 2 ½À45 ; 45 the maximal SPL deviations are 1.6 and 2.4 dB for 2 and 3 kHz, respectively. Figure 13 shows the normalized phase in radian over 2p. Similarly to the SPL directivity, the maximum phase deviations arise at the rear side, too, with 8% at 2 kHz up to 114% at 3 kHz. In the most relevant listening region (b 2 ½À45 ; 45 ) the maximal deviation is less than 8% for both frequencies.
To illustrate the accuracy of the reconstruction of the directivity as a weighted average over all radiation directions, we have calculated the directivity index for six monofrequent signals in octaves from 250 Hz to 4 kHz, So, in addition to the polar radiation patterns in the x 1 -x 2 plane (Sec. V C) also the directivity index-as a threedimensional integral measure-is accurately modeled by our optimization algorithm.
Overall the occurring deviations are comparable to other grid-based source synthesis approaches without Euler equations (Bilbao et al., 2019;Feistel et al., 2009;Karjalainen and Savioja, 2001;van Beuningen and de Vries, 1997). To further improve the reproduction quality of the presented approach a finer numerical resolution, yielding more grid-based monopoles within the same volume, could be used.

VI. CONCLUSION
A grid-based monopole synthesis method based on the Euler equation as a more general form of the wave equation is introduced and validated on the examples of an analytic piston model (Sec. IV) and a measured two-way monitor (Sec. V). We demonstrated that our source synthesis method can model the directivity of practically relevant sound sources in amplitude and phase. The reproduction quality is best for the radiation angles that contribute most to the overall SPL.
As discussed in Sec. II D, the synthesized source is corrected for measurement errors contradicting the Euler equations. Therefore, some deviations to the measured reference are not numerical errors but a "physically filtered" version of the measurements.
The computational effort of the utilized adjoint-based optimization algorithm depends solely on the FDTD solution of the equation system. They do not increase with the number of monopole sources within the same finite differences domain. This is particularly advantageous when very complex sources need to be reproduced with high accuracy using numerous synthetic sources. The method autonomously seeks an optimal representation in terms of a local minimum of the objective function without the need for prior knowledge about the synthesized source itself.