Adjoint-based sound reinforcement in the time domain
Lewin STEIN(1) , Florian STRAUBE(2) , Jörn SESTERHENN(1) , Stefan WEINZIERL(2) , Mathias LEMKE(1)(1) Computational Fluid Dynamics Group at Technische Universität Berlin, Germany(2) Audio Communication Group at Technische Universität Berlin, Germany

1 Jun 2019
Abstract
The determination of the electronic drive of line source arrays is an ill-posed inverse problem. The authors
present an adjoint-based approach in the time domain as an alternative to the commonly used frequency-based
methods. Therein, monopole sources in the E ULER equations are optimized to reinforce a prescribed target
sound field. In doing so, driving functions for the synthesis of the desired sound field, including a non-uniform
base flow, can be determined.

Keywords: loudspeaker arrays, sound reinforcement, acoustic time-reversal, monopole sources

1

INTRODUCTION

Sound reproduction, sound reinforcement, and source identification methods are typically based on frequency
domain approaches and restricted to the wave equation with a homogeneous background flow [1, 2, 3, 4]. In
contrast, we suggest an adjoint-based solution of the E ULER equations in the time domain for the following
reasons: The E ULER equations, as a generalized form of the wave equation, allow to consider inhomogeneousbackground flow and nonlinear sound effects. The adjoint-based approach (time-reversal) facilitates to simultaneously optimize all degrees of freedoms (control parameters), for example the speaker alignment, the signal
amplitude and the phase over the entire frequency range of arbitrary line arrays. The computational cost of aconventional parameter study (only forward in time) increases linearly with the degree of freedom, while thecomputational costs of the adjoint approach do not increase with the number of optimized parameters, instead,they remain constant.

2

METHODFor a better understanding of the results, a brief overview of the adjoint-based computational aero-acoustic(CAA) solver is given in the following. An introduction of the method with more details can be found in [5].
2.1 Basic idea of the adjoint approach
The adjoint equations can be derived in a continuous or discrete manner. Here, for the sake of simplicity, they
are introduced in a discrete version as in [6]. A matrix-vector notation is used. The vector space is the full
solution in space and time. For a complete discussion on the derivation, the reader is referred to [7, 8].
Adjoint equations arise by a so-called objective function J, which is defined by the product between a geometricweight g and the system state q.
J = gT q(1)
The system state q corresponds to the solution of the governing system
Aq = s(2)

with A as governing operator and s as source terms on the right-hand side. To optimize J by means of s, theequation has to be solved for all s. In order to reduce the computational effort, the adjoint equation can be

1162

used
with the adjoint variable q∗ . By

AT q∗ = g,

(3)

T

(4)

J = gT q = AT q∗

q = q∗ T Aq = q∗ T s

an expression is found, which enables the computation of J without the need to solve the system for every
discrete s at once. After calculating the adjoint equation, i.e. q∗ , the objective J can be determined by a simple
and computationally cheap scalar product (Eq. 4). Thus, the adjoint approach enables the efficient computationof J gradients with respect to s.
2.2 Incorporated CAA equations and optimization target
According to the intended application – the adjoint-based sound field generation in the time domain – theobjective function J is defined in space and time as
J=

1
2

ZZ

p0 − p0target

2

σ (xi ,t)dΩ.

(5)

The variable p0target (p = p0 + p0 ) denotes the desired sound field, resulting from quality characteristics for soundreinforcement tasks, such as a specific volume distribution or a smooth frequency response. The weight σ (xi ,t)defines where and when the objective is evaluated. Optimal sound reinforcement is realized if J reaches aminimum.
The minimum is to be achieved under the constraint that the E ULER equations   


ρui
ρ
0
0
∂t ρu j  + ∂xi ρui u j + pδi j  − ui ∂xi  0  = 0
(6)
ui pγ
p
p
s
γ−1
γ−1
abbreviated as E(q) q = s are satisfied, with q as state vector and s as a source term. The general goal is toobtain a solution which matches p0target (q) of Eq. 5 in an optimal sense by adapting the source term s, which
can be interpreted as multiple loudspeakers.
In order to use the adjoint approach for optimizing s, it is necessary to linearize the objective function and thegoverning system (6). This results in
δJ =

ZZ


q − qtarget σ δ qdΩ,
|
{z
}

(7)

=g

and

Elin (q0 )δ q = δ s(8)with a redefined weight g = (q − qtarget )σ . Combining the linearized system and the objective in a Lagrangianmanner by introducing the (adjoint) multiplier q∗ , leads to
δJ= gT δ q − q∗ T (Elin (q0 )δ q − δ s) ,
|
{z
}
∗T

= q δs+δq

T

(9)

=0

T
g − Elin(q0 )q∗ .

For the sake of simplicity, the integrals are not shown. The adjoint equations result inT
g − Elin(q0 )q∗ = 0 ,

(10)

1163

with q∗ as adjoint state vector. The change of the objective function reads
δ J = q∗ T δ s,(11)

similar to (4). Thus, the adjoint solution can be interpreted as a gradient of J with respect to the source term s∇s J = q∗ .

(12)

A more detailed derivation of the adjoint approach is discussed in [5, 8, 9].
2.3 The adjoint CAA solver
The driving functions for given speaker positions are determined and optimized iteratively:
1. Initially, guess the loudspeaker source terms s0 = 0.
2. Solve the governing E ULER equations (6) with si . Here, the exponent i = 0 is the first iteration step.
3. Compare with the desired target solution and break if a suitable convergence criterion is reached.
4. Solve the adjoint E ULER equations (10).
5. Compute the gradient ∇s J using the adjoint solution.
6. Update the current forcing si+1 = si + ∇s J.
7. Repeat, starting with 2.
The resulting optimal s is deconvolved with the considered primitive input signal defining p0target and the loudspeaker properties (e.g. electronic filters) resulting in an optimal driving for every considered loudspeaker. The
solutions shown in this work are the result of fifteen iteration steps.
Numerical specifications (parallelization, resolution, spatial derivative, time integration, filters, boundary conditions) of the adjoint CAA solver are listed in [5], i.e., originally in [7, 10].

3

TEST CASE

The three-dimensional test case studied here (Fig. 1) contains a line array with five monopoles speakers, i.e.
spherical loudspeakers, which are aligned in an x1 -line at x1 = (0.4, 0.6, 0.8, 1.0, 1.2) m, and constant x2 = 0.8 m
and x3 = 0.8 m.
The audience, i.e., the mask σ (Eq. 5), is located symmetrically around the line array at a minimal distance of0.175 m and a maximal distance of 0.4 m. The physical domain is evaluated at an equidistant rectangular gridwith a spacing of dx = 1.6/197 m and with a temporal sampling rate of 48 kHz.
As homogeneous background medium air is approximated as an ideal gas with a speed of sound c0 = 343 m/s,
a density of ρ0 = 1.21 kg/m3 , a pressure p0 = 101 325 Pa, an adiabatic index γ = 1.4 and a specific gas constantRs = 287.
As a test signal for all monopole speakers, a sine sweep with an exponential frequency increase and a constant
amplitude over time (see s̃CAA in Fig. 2) is used. Below 1.3 kHz and above 2.7 kHz the sweep has a fade in andfade out period. At the end of the sweep, there is 1 ms silence (the last 48 time steps). Throughout this paper
the frequency spectrum between 1.3 kHz and 2.7 kHz is evaluated. The initial transient is not discussed here. To
study a practically relevant spatially inhomogeneous SPL distribution the relative amplitudes and phases between
the signals differ among all speakers. With increasing x1 -position the signal amplitude decreases (cf. Fig. 1).

1164

p0 / Pa1
1.2
0.5

x2 / m

1
0.8

0

0.6

−0.5

0.4
0.4

0.6

0.8
x1 / m1

1.2−1Figure 1. Sound field radiated by a line array with five monopole sources. A snapshot of the pressure reference
field at a central x1 x2 -slice (x3 = 0.8 m) is shown, where all speakers are located. Isolines of the objective
weight function σ (Eq. 5) at values of 0.01, 0.5, and 0.99 are denoted by white, dashed lines.

4

RESULTSAs a reference solution, a complex-directivity point source (CDPS) model in the frequency domain is utilized[3, 4, 11]. The CDPS model provides, for a prescribed sine sweep s̃CDPS , the reference field p0target (Eq. 5)
for all spatial and temporal locations. With the reference field and the source locations as input, the CAA
framework iteratively recalculates the amplitude and phase of the driving sine sweep s. By comparing s withthe original signal s̃CDPS the adjoint-based approach is validated. We verify to which extent the adjoint approachis able to identify the optimal excitation signals for a given sound field.
In this work, the signals are optimized for predefined fixed speaker alignment. However, in general, the adjoint
approach is also capable to identify speaker locations [9].
4.1 Conversion and comparability of source terms with and without flow
To obtain comparable – physically equivalent – forms of source terms called s̃CAA and s̃CDPS , s of Eq. 6
needs to be converted to s̃CAA : While s̃ denotes a source term of the wave equation (second time derivative),s is a source term of the pressure equation (first time derivative). s̃CDPS of the CDPS model drives an ideal
monopole speaker, i.e., a Dirac delta function. In contrast s of the CAA solver consists of an entire Gaussian
distribution of monopoles, which taken all together approximate a single monopole speaker. A single Dirac delta
source would create instabilities (Runge phenomenon), while its smoothed alternative representation as Gaussian
distribution – applied here – is numerically more stable. We derive in the following a conversion rule betweens and s̃CAA , which will be applied later to compare s̃CAA and s̃CDPS in their physically equivalent form.
The linearized momentum and energy equation of the E ULER equations (Eq. 6) of the CAA solver (withoutflow) are!  !
ρ0 ∂t u0j
∂x j p00+ γ p0 ∇u0 =
.
(13)
∂t p0
s
γ−1
γ−1

1165

The time derivative and the combination of both equations gives the wave equation∂t2 p0 − c20 ∆p0 = (γ − 1)∂t s(y,t) ,(14)

where s(y,t) is the related force term. Since we assume an ideal gas, c20 = γ p0 /ρ0 holds.
The general solution of the inhomogeneous wave equation in free space is known as
p0 (x,t) =

γ −1
4π

Z
R3

∂τ s(y, τ)dy ,|x − y|(15)

with the delayed time τ = t − |x − y|/c0 . As discussed before, to ensure numerical stability, the CAA solver uses
an entire Gaussian shaped distribution of monopoless(y,t) =2 −2
1e−(ym −y0 ) µ δ (ym − y)sCAA (t)∑γ −1 m(16)

instead of a pure monopole source located at y0 . sCAA (t) is the time signal, which is driving all individual
monopoles of the Gaussian distribution. The monopole positions are numbered with the index m. The varianceof the Gauss is set as µ = 2∆x. Hence µ is eight times smaller than the shortest acoustic wavelength evaluatedhere (2.7 kHz).
Combining the latter equations provides
2 −2

p0 (x,t) = ∑ e−(ym −y0 ) µm∂τ sCAA (τm )
,
4π rm

(17)

being τm = t − rm /c0 the delayed time due to the distance rm = |x − ym | between the observer x and eachmonopole position ym .
In case of a homogeneous background flow in x1 -direction with a Mach number M, the solution follows
1 − M cos(θm )/am+0−(ym −y0 )2 µ −2 c0 M cos(θm )+sCAA (τm ) +p (x,t) = ∑ e∂ + sCAA (τm ) ,
(18)
2 a3
4π rm
4π rm (1 − M 2 ) am τmmbeing τm+ = t − rm (am − M cos θm )/(1 − M 2 )/c0 , a2m = 1 − M 2 sin2 (θm ), and cos(θm ) = (x1 − y1m )/rm , see [12]. In
case the mic is placed in the far-field downstream of the source with θm ∼= 0, this simplifies to2 −2

p0 (x,t) = ∑ e−(ym −y0 ) µm∂τ + sCAA (τm+ ).
4π rm (1 + M)

(19)

In contrast to Eq. 14, the CDPS model solves the wave equation∂t2 p0 − c20 ∆p0 = δ (y0 − y)s̃(t) ,(20)

with a different source type, denoted by the overlying tilde. Instead of a Gaussian monopole distribution, theCDPS model uses a single monopole at the center of the Gauss δ (y0 − y). The CDPS model sets s̃(t) = sCDPS
directly equal to the test signal, i.e., the exponential sine sweep. Please note here the difference between thesignals s̃(t) of Eq. 20 and s(t) of Eq. 14, especially the time derivative. The far-field solution of Eq. 20 is
p0 (x,t) =s̃(τ0+ ).
4π r0 (1 + M)

(21)

Since both sound fields p0 of Eq. 19 and Eq. 21 must be equal, s̃CAA and sCAA are related by+s̃CAA (τ0+ )2 −2 ∂τ + sCAA (τm )= ∑ e−(ym −y0 ) µ
.
r0
rm
m

1166

(22)

Amplitude / Pa s−2

8sCAA fs / Pa s−2 if M = 0
sCAA fs / Pa s−2 if M = 0.05 and θm = 0s̃CAA / Pa s−2 if M = 0

6
4
2
0

100

200

300

400

500

600

700

800

900

1,000

1,100

1,200

time step
Figure 2. Two different types of source terms sCAA and s̃CAA corresponding to the same exponential sine sweep
used as input signal throughout this work. While varying amplitude and phase this input signal is used for all
loudspeakers. For better visibility, only the upper half of the signal is shown. Furthermore, the peak envelope
for each signal is depicted, to highlight the different maximal amplitudes.
Thus, after converting sCAA to s̃CAA according to Eq. 22, s̃CAA and s̃CDPS are comparable.
For the present case without flow (M = 0) Fig. 2 shows the utilized exponential sine sweep, i.e., s̃CAA ∼sin(ω0 eα t t). α determines the exponent of the frequency increase over time. Besides an initial and final fadeout its amplitude is constant. The corresponding sCAA signal – perfectly related by Eq. 22 – is also shown in
Fig. 2.
The outer time derivative of a cos sweep gives a sine sweep (and vice versa), which explains the phase shift ofπ/2 between sCAA and s̃CAA . The inner time derivative of the exponential frequency dependency of the sweep
corresponds to the envelope curvature. In the case of homogeneous background flow with M c0 = 17.15 m/s, thetime delay between the different monopoles of the Gaussian distribution is either increased or reduced dependingon the observer position (θm ). However, this effect is restricted to high-frequency components (Fig. 2) and is
negligible (±0.1 dB) within the accuracy of the SPL errors discussed below. By increasing the ratio between
the Gauss variance and the shortest wavelength, this flow effect completely disappears.
Once Eq. 22 is evaluated for a given signal (here a sine sweep), the conversion factor is universal for all
speakers if – like in the present case – only the absolute phase and overall amplitude of the speaker signals
differ. Thus, also the spectral ratio F [s̃CAA ]( f )/F [sCAA ]( f ) needs to be determined only once. Therein, F isthe discrete Fourier transformation. Throughout the following spectral evaluation, this ratio will be applied toconvert sCAA ( f ) to s̃CAA ( f ).
4.2 Spectral comparison of the original excitation signal and the adjoint-based solution
Fig. 3 depicts the gains of s̃CAA provided by the CAA solver. Due to the used exponential sweep (cf. Fig. 2),
higher frequencies have smaller time shares reflected by a decaying gain with increasing frequency in Fig. 3.
The overall amplitude variation between the speakers causes the frequency independent offset.
Fig. 4 shows the phases of s̃CAA , featuring periodic 2π shifts as expected.
Fig. 5 and Fig. 6 show the gain and the phase differences between s̃CDPS (reference solution of the CDPS
model) and s̃CAA (adjoint CAA solution). The maximal gain deviation error at the frequency interval edges is±0.25 dB, and less than ±0.15 dB between 1.5 kHz and 2.4 kHz for all speakers (Fig. 5). Hence, also therelative pressure levels between all five speakers are correctly captured. The phase deviation between the CAA

1167

1

0

φopt / (2π)

Gopt / dB

0.8
−10
−20
−30

speaker 1
speaker 3
speaker 5

speaker 2
speaker 4

1.5

2
f / kHz

0.6
0.4
0.2

2.5

Figure 3. CAA solution of the input gains for all
speakers. The SPL is normalized to the loudest
speaker.

0

1.5

2
f / kHz

2.5

Figure 4. CAA solution of the input phase of theloudest speaker. All other speakers show qualitatively
similar behavior.

and the CDPS solution corresponds to a phase shift of less than half a percent (phase shift in rad/2π) within
the whole frequency interval (Fig. 6). Probably, the phase error is so small since a phase shift of p0 in Eq. 5
has a greater impact on the objective in contrast to a perturbation of the pressure magnitude.
To conclude, the adjoint-based CAA solver is able to provide driving functions for the line array, which arecomparable to the CDPS reference up to ±0.25 dB and a phase error of half a percent.

5

SUMMARY

An adjoint-based framework for the determination of optimal driving functions of monopole speakers for soundreinforcement applications is presented. The effects of the employed Gaussian monopole implementation is
discussed in detail. The required conversion function is derived analytically based on the analytic solution of
the inhomogeneous wave equation including background flow. In summary, it has been shown that the approachis able to determine optimal driving functions, i.e., gain and delay, for a representative line array arrangement.

ACKNOWLEDGEMENTS
We acknowledge the financial support by the Deutsche Forschungsgemeinschaft (DFG) within the projects LE
3888/2-1 and WE 4057/16-1.

REFERENCES[1] D. G. Meyer, “Computer simulation of loudspeaker directivity,” J. Audio Eng. Soc., vol. 32, pp. 294–315,
May 1984.
[2] G. W. J. van Beuningen and E. W. Start, “Optimizing directivity properties of DSP controlled loudspeakerarrays,” Proc. of the Inst. of Acoustics: Reproduced Sound, vol. 22, no. 6, pp. 17–37, 2000.
[3] P. Meyer and R. Schwenke, “Comparison of the directional point source model and BEM model for arrayed loudspeakers,” Proc. of the Institute of Acoustics, vol. 25, no. 4, 2003.
[4] S. Feistel, A. Thompson, and W. Ahnert, “Methods and limitations of line source simulation,” J. Audio
Eng. Soc., vol. 57, pp. 379–402, June 2009.

1168

[5] L. Stein, F. Straube, J. Sesterhenn, S. Weinzierl, and M. Lemke, “Adjoint-based optimization of soundreinforcement including non-uniform flow,” submitted to The Journal of the Acoustical Society of America,
2019.
[6] M. Giles and N. Pierce, “An introduction to the adjoint approach to design,” Flow, Turbulence and Combustion, vol. 65, pp. 393–415, 2000.
[7] M. Lemke, Adjoint based data assimilation in compressible flows with application to pressure determinationfrom PIV data. PhD thesis, Technische Universität Berlin, 2015.
[8] M. Lemke and J. Sesterhenn, “Adjoint-based pressure determination from PIV data in compressible flows– validation and assessment based on synthetic data,” European Journal of Mechanics - B Fluids, vol. 58,
pp. 29–38, 2016.
[9] M. Lemke, F. Straube, F. Schultz, J. Sesterhenn, and S. Weinzierl, “Adjoint-based time domain soundreinforcement,” in Audio Engineering Society Conference: 2017 AES International Conference on Sound Reinforcement – Open Air Venues, Aug 2017.
[10] L. Stein, Simulation and Modeling of a Helmholtz Resonator under Grazing Turbulent Flow. PhD thesis,
Technische Universität Berlin, 2019.
[11] F. Straube, F. Schultz, D. Bonillo, and S. Weinzerl, “An Analytical Approach for Optimizing the Curvingof Line Source Arrays,” Journal of the Audio Engineering Society, vol. 66, pp. 4–20, Feb. 2018.
[12] J. Delfs, “Basics of aeroacoustics,” in lecture notes, TU Braunschweig, 2018.
0.3
speaker 1
speaker 3
speaker 5

0.1
0
−0.1
−0.2
−0.3

speaker 1
speaker 3
speaker 5

4 · 10−3(φopt − φref ) / 2π

(Gopt − Gref ) / dB

0.2

6 · 10−3

speaker 2
speaker 4

speaker 2
speaker 4

2 · 10−3
0−2 · 10−3−4 · 10−3

1.5

2
f / kHz

2.5

Figure 5. Difference between the CAA and the CDPS
input gains of the speakers.

−6 · 10−3

1.5

2
f / kHz

2.5

Figure 6. Difference between the CAA and the CDPS
input phases of the speakers.

1169