Technical University of Berlin
Faculty I Humanities and Educational Sciences
Institute for Speech and Communication
Audio Communication Group
Master’s thesis in the program
“Audio Communication and Technology” (M. Sc.)
Adjoint-Based Monopole Synthesis
of Sound Sources with Complex Directivities
Arne Birk Hölter
July 2020
Supervision:
Prof. Dr. Stefan Weinzierl
Dr. - Ing. Mathias Lemke
Dr. rer. nat. Florian Straube
Abstract
Numerical time domain methods are continuously improving and gain in importance in virtual
acoustics. The accurate modeling of directive sound sources is a prerequisite for acoustical
simulations. In contrast to frequency domain methods, the directivity pattern of sound
sources cannot be analytically implemented in the time domain directly. The sound sources
are rather approximated by the superposition of monopoles around the directive sound source.
For that purpose, an adjoint-based monopole synthesis method is discretized in a finite
differences time domain (FDTD) scheme. Therein, the full non-linear Euler equations are
solved by means of computational aeroacoustics (CAA).
To efficiently compute reference sound fields, an analytical complex directivity point source
(CDPS) model is embedded into the existing architecture of the CAA solver, which enables
a decomposition of the computing domain to parallelize the computation. In order to avoid
unfavorable interferences between the monopoles in FDTD, its spatial expansion is analyzed.
Finally the adjoint-based monopole synthesis method is considered for a dipole, a quadrupole
and a (complex) circular piston model. The results are evaluated by graphical representations
and technical measures, e.g., 3D directivity pattern figures. The analysis is not only limited
to the reproduction of the reference directivity patterns, as it is as well the intention of this
thesis to examine the procedure of the adjoint-based monopole synthesis method.
I
Zusammenfassung
Numerische Methoden im Zeitbereich sind ein aktiver Forschungsbereich in der virtuellen
Akustik. Die Modellierung von richtungsabhängigen Schallquellen ist eine Vorraussetzung
für akustische Simulationen. Im Gegensatz zu den Methoden im Frequenzbereich können
Richtcharakteristiken von Schallquellen nicht direkt analytisch im Zeitbereich implementiert
werden. Vielmehr werden die Schallquellen durch eine Superposition von Monopolen um
die richtungsabhängige Schallquelle approximiert. Zu diesem Zweck wird eine adjungierten-
basierte Monopolsynthese mithilfe eines finite Differenzenschema im Zeitbereich (FDTD)
diskretisiert. Darin werden die kompletten nicht-linearen Euler Gleichungen mittels numeri-
scher Strömungsakustik (CAA) gelöst.
Zur effizienten Berechnung von Referenzschallfeldern wird ein analytisches Complex Directi-
vity Point Source (CDPS) Modell in die bestehende Architektur des CAA Lösers eingebet-
tet, die eine Zerlegung des Rechengebiets zur Parallelisierung der Rechnungen ermöglicht.
Um unerwünschte Interferenzen zwischen den Monopolen in FDTD zu vermeiden, wird die
räumliche Ausdehnung der Monopole untersucht. Abschließend wird die adjungierten-basierte
Monopolsynthese für einen Dipol, einen Quadrupol und einen (komplexen) Rundkolbenstrah-
ler (Circular Piston) angewendet. Die Ergebnisse werden mittels graphischer Darstellung und
technischer Beurteilungen, wie beispielsweise 3D Abbildungen von Richtcharakteristiken, aus-
gewertet. Die Analyse ist nicht nur auf die Nachbildung der Richtcharakteristiken beschränkt,
sondern konzentriert sich auch auf die Vorgehensweise der adjungierten-basierten Monopol-
synthese.
II
Contents
1 Introduction 1
1.1 StateoftheArt ................................ 2
1.2 Objective.................................... 3
2 Methods 5
2.1 Complex Directivity Point Source Model . . . . . . . . . . . . . . . . . . . 5
2.1.1 DirectivityPattern........................... 6
2.1.2 Free Space Green’s Function . . . . . . . . . . . . . . . . . . . . . . 8
2.2 Adjoint-Based Monopole Synthesis . . . . . . . . . . . . . . . . . . . . . . 8
2.2.1 AdjointEquations ........................... 9
2.2.2 Adjoint Sound Field Generation . . . . . . . . . . . . . . . . . . . . 9
2.2.3 Monopole Sound Source . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2.4 SourceUpdating............................ 13
2.2.5 Objective Area Mask Function . . . . . . . . . . . . . . . . . . . . . 14
2.2.6 IterativeProcess............................ 15
3 Simulation Settings and Setup 17
4 Evaluation Methods 20
5 Testcases 24
5.1 Case (I): Implementation of the CDPS Solver . . . . . . . . . . . . . . . . . 24
5.2 Case (II): Nearest Distance of Neighbouring Monopoles . . . . . . . . . . . 27
5.3 Case(III):Dipole................................ 31
5.4 Case(IV):Quadrupole............................. 35
5.5 Case (V): (Complex) Circular Piston . . . . . . . . . . . . . . . . . . . . . 40
6 Conclusion 48
A Objective Functions and Source Positions 54
III
List of Figures
1 Normalized Gaussian distribution on a 1D grid at xs≈0.5m of two different
grid discretization schemes. . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2 1D objective area mask function σ(x, t)with x0= 0.5m, rin = 0.2m, rout =
0.3m and l= 0.05 m. ............................. 14
3 Iterative process of the adjoint-based monopole synthesis method. . . . . . . 15
4 2D slice of the considered grid in the x1-x2-plane at x3= 0.5m........ 17
5 Virtual microphone positions (expanded view of Fig. 4). . . . . . . . . . . . 21
6 Computation time over the parallelism cores of two grid discretization schemes. 26
7 Relative objective funtion J / Jmax over distance between two monopole sources
ingridpoints(nodes).............................. 29
8 Gaussian distribution of two monopoles with 11 grid nodes distance to each
other....................................... 29
9 Polar directivity patterns of two spatially distant monopoles with variable xdist
at2kHz..................................... 30
10 Forcing signals of the dipole. . . . . . . . . . . . . . . . . . . . . . . . . . 32
11 Source positions xs,1-xs,3ofthedipole. ................... 32
12 SPL and Phase deviation of the reference and optimized sound field at the
virtual microphone positions of the dipole. . . . . . . . . . . . . . . . . . . 33
13 Amplitude balloon plots of the dipole. . . . . . . . . . . . . . . . . . . . . . 33
14 Polar directivity patterns of the dipole. . . . . . . . . . . . . . . . . . . . . 34
15 Forcing signals of the quadrupole. . . . . . . . . . . . . . . . . . . . . . . . 37
16 Determined source positions of the quadrupole. . . . . . . . . . . . . . . . . 37
17 SPL and phase deviation of the quadrupole. . . . . . . . . . . . . . . . . . 38
18 Amplitude balloon plots of the quadrupole. . . . . . . . . . . . . . . . . . . 38
19 Polar directivity patterns of the quadrupole. . . . . . . . . . . . . . . . . . . 39
20 Forcing signals of the circular piston. . . . . . . . . . . . . . . . . . . . . . 42
21 Amplitude frequency response deviation of sp,2and sp,4to sp,1. ....... 43
22 SPL and phase deviation at the virtual microphone positions of the circular
piston. ..................................... 43
23 Polar directivity patterns of the circular piston. . . . . . . . . . . . . . . . . 44
24 Amplitude balloon plots of the circular piston. . . . . . . . . . . . . . . . . 45
25 Evolution process of the optimized circular piston at 2.5 kHz. . . . . . . . . 45
26 SPL and phase deviation at the virtual microphone positions of the complex
circular piston (K= 0.1)............................ 46
27 Amplitude balloon plots of the complex circular piston (K= 0.1). . . . . . . 46
28 SPL and phase deviation at the virtual microphone positions of the complex
circular piston (K= 1.0)............................ 47
29 Relative objective function Jn/J0over the iteration loops n.......... 54
IV
List of Tables
1 Settings of the CDPS and CAA simulations. . . . . . . . . . . . . . . . . . 19
2 Computation time of the Matlab and Fortran implementation. . . . . . . . . 26
3 Considered monopole distances. . . . . . . . . . . . . . . . . . . . . . . . . 28
4 Monopole source positions of the cases (III) – (Va) determined by the adjoint-
based monopole synthesis method. . . . . . . . . . . . . . . . . . . . . . . 55
5 Monopole source positions of the cases (Vb) –(Vc) determined by the adjoint-
based monopole synthesis method. . . . . . . . . . . . . . . . . . . . . . . 55
V
List of Symbols
Symbol Definition
Latin
Agoverning operator
cspeed of sound
Ddriving function of a sound source
EEuler equations
ffrequency
fccut-off frequency
fssampling frequency
fin input / test signal
ggeometric weight
G0,3D 3D Green’s function, free field acoustic transfer function
Hdirectivity pattern
j imaginary unit, j2=−1
Jobjective function
Kphase scaling parameter of the complex circular piston
ltransition length of the objective area mask function
lssource expansion diameter
Lpsound pressure level
Mmol molar mass
ppressure field in the time domain
p0ambient pressure
p0sound pressure
p∗adjoint pressure
PATF sound pressure transfer function in the frequency domain
Psound field in the frequency domain (complex frequency spectrum)
qsystem state
VI
List of Symbols (Continuation)
Symbol Definition
q∗adjoint system state
r0distance from the source to the receiver position, |x−x0|
rin effective start of the objective area mask function
rout effective end of the objective area mask function
rssource region radius
Rgas constant
Rsspecific gas constant
ssource term
sfsource forcing signal
sgGaussian distribution for the spatial source activation
spsource term regarding the pressure
ttime
Ttotal time
T0reference temperature
utot total velocity, consists of |u|+c
xdiscretized computation area (grid), consists of the three spatial direc-
tions [x1, x2, x3]T
x0reference sound source location
xgnumber of the discretized grid node
xs,moptimized sound source locations
VII
List of Symbols (Continuation)
Symbol Definition
Greek
αactive radiating factor (ARF)
αsline search step width
βangle from the source to the receiver, contains the azimuth ϕand the
elevation ϑ
γisentropic exponent
λwavelength
Λyheight of the loudspeaker
ϕazimuth angle
Ψmask function of the source region
%0reference density of the air
σevaluation mask function
σsd standard deviation
ϑelevation angle
Θradius of the circular piston
VIII
1 Introduction
1 Introduction
To improve the acoustical behaviour of environments such as rooms, concert venues etc.,
acoustical optimizations are unavoidable. In most cases experimental optimizations are ex-
pensive and very time-consuming. The refined way is to apply acoustical simulations.
Acoustical simulations have been broadly developed and studied in the recent years, where
mainly two categories of simulation methods have been formed: geometrical- and wave-based
methods. Geometrical-based methods assume the sound propagation as a ray, while wave-
based methods numerically approximate the solution of the wave equation (Takeuchi et al.,
2019). Both allow a prediction of sound fields in a specific area emitted by any number of
sound sources and receiver positions. To approximate real environments, such as rooms or
open air venues, the actual simulation strategies still come up against limiting factors, e.g.,
the implementation of boundary conditions or non-uniform flow (Stein et al., 2019). Among
these limiting factors, directive sound sources have to be highlighted as this thesis will be
focused on them. Currently, in frequency domain approaches the frequency-dependent di-
rectivity is already considered, e.g., in the complex directivity point source (CDPS) model
(Meyer, 1984; Feistel et al., 2009) or at least with simple source directivities in geometrical
room acoustic simulations (Poirier-Quinot et al., 2017). For the most common time domain
approaches, that mostly make use of finite difference time domain (FDTD) schemes, no
method is available to date that is able to model complex directivities in an adequate way
(see Sec. 1.1), i.e., directivity pattern impulse responses with changing amplitude and phase
in terms of direction and frequency (Takeuchi et al., 2019).
More general, directivity means that different sound pressure amplitudes and phases are ob-
tained in different directions on equidistant evaluation positions around the source, i.e., on
a spherical surface with the source at the origin. If a sound field cannot be considered as
diffuse, the obtained sound field is strongly influenced by directive sound sources. Therefore
an implementation method of directive sound sources in FDTD is in need.
This thesis makes use of an adjoint-based approach which is able to solve the full non-linear
Euler equations and the corresponding adjoint in the time domain by means of computa-
1
1 Introduction
tional aeroacoustic (CAA) techniques (see Sec. 2). Stein et al. (2019) demonstrated that the
method is able to optimize driving functions of sound sources. Moreover, the ability of the
method to find optimal monopole source locations and the corresponding monopole weights
to reconstruct directive sound sources will be analyzed.
1.1 State of the Art
The reproduction of sound sources with 3D audio systems in the frequency domain is well
described in the literature, e.g., Wave Field Synthesis (WFS) or Near Field Compensated
Higher Order Ambisonics (NFC - HOA) (Berkhout et al., 1992; de Vries et al., 1994; Slavik
and Weinzierl, 2008; Ahrens, 2012). Usually, the sound sources are inserted as monopoles,
i.e., point sources, because they are easily implemented and they may be described analyti-
cally well.
The complex directivity point source (CDPS) model is an analytical method to compute
sound field predictions in the free field and the frequency domain (Meyer, 1984; van Beunin-
gen and Start, 2000; Meyer and Schwenke, 2003; Feistel et al., 2009). Further, the model is
able to deal with modeled and measured loudspeaker directivity data (Straube, 2019). The
sound field is separately computed for every considered frequency in the CDPS model.
As the available computational power is increasing significantly in the recent time, time do-
main approaches become increasingly important. In this thesis a CAA method in a FDTD
scheme is considered. It is based on the wave equation in the time domain instead of the
Helmholtz equation. If the wave equation is expanded to the full non-linear Euler equations,
an easier treatment of, i.a., non-uniform flow, boundary conditions and heat stratifications
is possible (Stein et al., 2019). Instead of the complex directivity in the frequency domain
approaches, the directivity is up to date modeled as a decomposition of the source into spa-
tially located monopoles.
Source modeling in FDTD has already approached with many different methods. In general
two different modeling strategies have been mainly considered: superposition of point sound
sources and expansion into spherical harmonics or multipoles (Takeuchi et al., 2019).
2
1 Introduction
First, the superposition of point sound sources approximates the directive source with sec-
ondary sources in the frequency domain. Subsequently, a set of appropriate coefficients must
be determined for the predefined point sources (Redlich, 2017; Opdam et al., 2016), which
may be solved, e.g., through the least-squares-method (Escolano et al., 2007). However,
these methods must convert the results into the time domain, which may degrade the final
accuracy.
Spherical-harmonics- or multipole-based methods consider the directivity in the time domain
directly (Takeuchi et al., 2019). In the spherical-harmonics-based method a spatial Gaussian
pulse is multiplied with the spherical harmonics, which approximates the desired directivity
(Sakamoto and Takahashi, 2013). The method was already fitted to measured loudspeaker
data in the work of Bilbao et al. (2019). The multipole-based method uses spatial derivatives
of Dirac delta functions to express a desired directivity by their linear combination (Bilbao
and Hamilton, 2018). The mentioned time domain methods do not consider frequency-wise
directivity because they consider the directional pattern for all frequencies simultaneously.
Thus, no method exists which is able to control (time-directional) frequency-wise directivity
and can be implemented directly in the time domain (Takeuchi et al., 2019).
1.2 Objective
Recent research has shown that the numerical time domain methods are continuously improv-
ing and gain in importance in acoustics as well as in the research area of virtual acoustics. As
outlined above, the implementation of directive sound sources is still an active research area.
A novel grid-based monopole synthesis approach was presented by Stein et al. (2020) using
an adjoint-based CAA method which will be applied in this thesis too. It treats all grid nodes
in a specific source region as monopole sources, which can be activated if necessary. But
the method does not provide information of the activation process of the multiple monopoles
directly. In contrast to the work of Stein et al. (2020), this thesis aims at synthesizing the
directive sound source by adding gradually single monopoles in a predefined source region.
Also, the objective of this procedure is to understand the optimization process of the adjoint-
3
1 Introduction
based monopole synthesis approach. Thus, directive sound sources shall be synthesized by a
small number of monopole sources providing a result in the optimal sense.
The following sections are divided as follows: Sec. 2 highlights the implemented complex
directivity point source (CDPS) model as well as the adjoint-based CAA method. Ensuing,
Sec. 3 describes the simulation settings and the general setup of both methods. After docu-
menting the computational tools, Sec. 4 gives an overview of the evaluation and visualization
methods. The evaluated results and the discussion of the regarded testcases are given in
Sec. 5. Finally, Sec. 6 compares the findings with the study of Stein et al. (2020) and draws
a conclusion of the intended adjoint-based monopole synthesis method.
4
2 Methods
2 Methods
As described in Sec. 1, this thesis will focus on the implementation of directive sound sources
in a FDTD scheme. To achieve this, an adjoint-based time domain method will be considered.
An analytical sound field simulation which is based on the complex directivity point source
(CDPS) model will serve as a reference for the validation during the CAA computation and
in the post-processing.
2.1 Complex Directivity Point Source Model
The analytic CDPS model is widely used for sound field prediction in virtual acoustics, sound
reproduction and sound reinforcement. It models sound sources as point sources, i.e., the
loudspeakers’ acoustical centers, with a complex directivity in the frequency domain. Subse-
quently, the propagation of all sources is computed and the produced partial sound pressure
fields are added.
Formulated in equations, the fundamental equation of the CDPS model (Meyer, 1984,
Eq. (5)); (van Beuningen and Start, 2000, Eq. (3-5)); (Meyer and Schwenke, 2003, Sec. 1.1);
(Feistel et al., 2009, Eq. (11)) reads
PATF(x, ω) =
M
X
m=1
(H(β(x,x0,m), ω)◦G0,3D(x,x0,m, ω)) D(x0,m, ω),(1)
where H(β(x,x0,m), ω)are the directivity patterns of the loudspeakers with β(x,x0,m)as
the angle from the source position x0,mto the receiver positions x, i.e., to all discretized grid
nodes in the three spatial directions x1,x2and x3, at the angular frequencies ω. Note that β
contains both spherical angles, the azimuth ϕand the elevation ϑ.G0,3D(x,x0,m, ω)denotes
the three dimensional free space Green’s function, i.e., the ideal point source (Williams, 1999,
Eq. (8.41), p. 265) and ◦is the Hadamard product, i.e., element-wise matrix multiplication.
D(x0,m, ω)are the loudspeakers’ driving functions and the output PATF(x, ω)is the sound
pressure transfer function at the receiver position xat the angular frequency ω.
5
2 Methods
As this thesis aims to reconstruct directivity patterns, the driving functions of the loudspeakers
D(x0,m, ω)will be considered as uniformly driven sources
D(x0,m, ω)=1 ∀mand ∀ω. (2)
Since the adjoint-based CAA solver is a time domain solver, the output PATF(x, ω)of the
CDPS model has to be transformed into the time domain using an inverse discrete (fast)
Fourier transform
p0(x, t) = F−1
t(PATF(x, ω)) (3)
and is excited by a sound signal fin(t)at every receiver position x
p0
out(x, t) = p0(x, t)∗tfin(t),(4)
whereby the asterisk ∗tdenotes the convolution with respect to time.
2.1.1 Directivity Pattern
First, variations of analytic directivity patterns H(β(x,x0,m), ω)in the frequency domain
are considered in this section. They characterize the generated reference sound field in the
free field. Subsequently, the reference sound field is reproduced by a direct implementation
in the time domain using the adjoint-based monopole synthesis method (see Sec. 2.2).
The variations range from simple directivity patterns as monopoles, dipoles and quadrupoles
to more complex directivity patterns like the baffled circular piston model. The dipole and
quadrupole will only be described by cosine functions. If an implementation of multipoles
of a higher order is required—the dipole has the order 1 and the quadrupole the order 2—a
spherical harmonics representation is recommended (Ahrens, 2012, Sec.2).
The simple and not frequency dependent directivity patterns may be described as follows:
the monopole radiates to all directions with the same amplitude. Since only the directivity
6
2 Methods
is of interest, it reads
Hmono (β(x,x0), ω)=1.(5)
The dipole may be described by a cosine (Sarradj, 2016, Eq. (10.30)) and thus the directivity
reads
Hdi (β(x,x0), ω) = cos (β(x,x0)).(6)
The formula of quadrupoles is given by Sarradj (2016, Eq.(10.41)) and depends on two
directions. Hence, a longitudinal and a lateral quadrupole can be described
Hquad,long (β(x,x0), ω) = cos (β(x,x0))2,(7)
Hquad,lat (β(x,x0), ω) = cos (β(x,x0)) cos β(x,x0) + π
2.(8)
Note that Eq. (5) - (8) are reduced to the directivity pattern with unique amplitude.
A frequency dependent directivity pattern is the baffled circular piston model with radius Θ
and a constant surface velocity. Its formula is given by (Skudrzyk, 1971, Eq. (26.42))
Hcirc(β(x,x0), ω) = 2J1ω
cΘ sin(β(x,x0))
ω
cΘ sin(β(x,x0)) ,(9)
denoting the cylindrical Bessel function of first kind of first order as J1(Olver et al., 2010,
Eq. (10.2.2)). Depending on the loudspeaker height Λyand the active radiating factor (ARF)
αof the loudspeaker, the radius Θof the circular piston may be calculated by (Schultz et al.,
2015)
Θ = αΛy
2.(10)
Note that the 3D circular piston has a circular symmetry around the main radiation axis, while
the dipole and quadrupole directivity are plane mirrored on the x1-x2-plane. The circular
piston model is widely used to model woofer and midrange speaker.
7
2 Methods
2.1.2 Free Space Green’s Function
The three dimensional Green’s function of point sources reads (Howe, 2002, Eq. (3.2.6))
G0,3D(x,x0, w) = e−jω
c|x−x0|
4π|x−x0|(11)
and describes the sound propagation of an ideal point source from the source position x0
to the receiver position xat the angular frequency ωand is therefore called the acoustical
transfer function (ATF). Further, air absorption is neglected in Eq.(11) and the speed of
sound is given by (Möser, 2012, Eq. (2.18))
c=sγ R T0
Mmol
=sγ p0
%0
.(12)
For diatomic gases the isentropic exponent amounts to γ≈1.4, the density of air will be
assumed with %≈1.2kg/m3and the atmospheric pressure is 101325 Pa. Hence, the assumed
speed of sound amounts to c≈343 m/s.
2.2 Adjoint-Based Monopole Synthesis
This section is based on Stein et al. (2019, Sec. 2). In Lemke (2015) the adjoint equations
are derived and discussed in more detail, e.g., the adjoint Euler equations with initial and
boundary conditions as well as the adjoint compressible Navier-Stokes equations.
In contrast to the commonly used frequency domain approaches for sound field generation
(see Sec. 2.1), the adjoint-based method uses a representation of wave propagation in the
time domain. Here, the full non-linear Euler equations are solved forward in time at the direct
computation and the corresponding adjoint Euler equations are solved backward in time at
the adjoint computation. Hereinafter, only the adjoint equations are described.
8
2 Methods
2.2.1 Adjoint Equations
The adjoint equations are introduced in a discrete version as in Giles and Pierce (2000).
Moreover, the vector space is the full solution in space and time.
The adjoint equations arise by an objective function J, which is defined by the product
between a geometric weight gand the system state q
J=gTq,g,q∈Rn,(13)
where qcorresponds to the solution of the governing system, which reads
Aq=s, A ∈Rn×n,s∈Rn.(14)
Therein, Adenotes the governing operator and sthe sources. The computational effort of
computing Jcan be reduced by the use of the adjoint equation
ATq∗=g(15)
with q∗as the adjoint variable. Combining Eq. (13), Eq.(14) and Eq. (15) gives
J=q∗Ts.(16)
After solving the adjoint equation, the objective Jcan be determined by a scalar product.
Therefore gradients of the objective Jcan be computed efficiently.
2.2.2 Adjoint Sound Field Generation
The adjoint-based equations of Sec. 2.2.1 referred to sound field generation are given in this
section. The total pressure pis given by
p=p0+p0,(17)
9
2 Methods
where the ambient pressure is p0= 101325 Pa and p0denotes the sound pressure. The
adjoint-based equations arise by a so-called objective function
J=1
2ZZ (p0(x, t)−p0
ref(x, t))2σ(x, t)dΩ,(18)
where p0(x, t)is now the sound pressure in the CAA model, p0
ref(x, t)is the reference target
sound pressure provided by the CDPS model described in Sec.2.1 and σ(x, t)is an additional
weight that defines the evaluation time and location of the objective function Jdefined in
space and time with dΩ = dxidt. Optimal sound field generation is met, when Jreaches
a minimum regarding the sound pressure with subject to the Euler equations including the
source forcing terms s(x, t):
min
pJ(19)
subject to E(q(x, t)) = s(x, t).
The constraint is the governing system abbreviating the Euler equations
∂t
%
%uj
p
γ−1
+∂x
%ui
%uiuj+pδij
uipγ
γ−1
−ui∂x
0
0
p
=
0
0
sp
(20)
with the velocity u(x, t), the density %(x, t)and γas the heat capacity ratio. Optimizing
the sound sources s(x, t)requires a linearization of Eq. (18) and the Euler equations in
Eq. (19):
δJ =ZZ (p0(x, t)−p0
ref(x, t))σ(x, t)
|{z }
=gp(x,t)
δp(x, t)dΩ(21)
10
2 Methods
and
Elin(q0(x, t))δq(x, t) = δs(x, t)(22)
with g(x, t) = [0,0,∆p(x, t)]Tσ(x, t). The combination of Eq.(21) and Eq. (22) in La-
grangian manner with an (adjoint) multiplier q∗(x, t)=[%∗(x, t), u∗
j(x, t), p∗(x, t)]Tleads
to
δJ =ZZ gT(x, t)δq(x, t)dΩ−ZZ q∗T(x, t) (Elin(q0(x, t))δq(x, t)−δs(x, t))
| {z }
=0
dΩ
=ZZ q∗T(x, t)δs(x, t)dΩ + ZZ δqT(x, t)g(x, t)−ET
lin(q0(x, t))q∗(x, t)dΩ.(23)
By demanding
g(x, t)−ET
lin(q0(x, t))q∗(x, t)=0,(24)
the second term of the right hand side of Eq. (23) vanishes and the change of the objective
function reads
δJ =ZZ q∗T(x, t)δsdΩ.(25)
The change of the objective function in Eq. (25) may be interpreted as a gradient of Jwith
respect to the sources s(x, t):
∆sJ=q∗(x, t).(26)
2.2.3 Monopole Sound Source
In contrast to the monopole implementation in frequency domain methods (e.g., Eq.(5)),
the implementation in the present FDTD scheme is approximated by Gaussian distributions.
11
2 Methods
The equation of the Gaussian distribution reads
sg(x) = 2π−3
2(σsd∆x)−3·e−0.5
(σsd∆x)2(x−xs)2,(27)
where the standard deviation σsd = 2 is used for all cases in this thesis. ∆xis the increment
of the discretized grid. Only cartesian grids are considered and ∆xis chosen equally for all
three spatial directions. The activation of the spatial expanded monopoles is attained by the
multiplication with a forcing signal sf(t)and thus defined in space and time:
sp(x, t) = sg(x)·sf(t).(28)
Fig. 1 depicts a normalized 1D Gaussian distribution using Eq. (27) for two different grid
discretization schemes. Small ∆xcause less absolute spatial expansion of the distribution.
One more way to reduce the spatial expansion might be the reduction of σsd. In that case
the Gaussian distribution will have less supporting points (nodes) to generate the waves.
Following Tam and Webb (1993) more than four nodes per wavelength are required to
ensure a satisfactory CAA transmission behaviour. The disregard would directly reduce
the CAA transmission behaviour and could cause parasitic, physically undesired waves. A
corresponding analysis is omitted here.
0.40 0.45 0.50 0.55 0.60
x1in m
0.0
0.2
0.4
0.6
0.8
1.0
Amplitude
(a) ∆x≈7.87 ·10−3m, i.e., 128 grid points /m
0.40 0.45 0.50 0.55 0.60
x1in m
0.0
0.2
0.4
0.6
0.8
1.0
Amplitude
(b) ∆x≈3.92 ·10−3m, i.e., 256 grid points /m
Fig. 1: Normalized Gaussian distribution on a 1D grid at xs≈0.5m of two different grid
discretization schemes.
12
2 Methods
2.2.4 Source Updating
To minimize the objective function Jiteratively, the sound sources sp(x, t), which generate
the sound field p0(x, t), are adapted after the adjoint computation. The intended procedure
to update the source distribution is separated into two steps.
At certain iteration loops na monopole source is added to the existing source distribution.
The source location of the added sound source xsis determined at the maximum time
averaged gradient of J:
xs=x"arg max
x 1
T
T
X
t
(p∗(x, t))2!Ψ(x)!#.(29)
The variable Ψ(x)is a mask function that spatially restricts the source region, e.g., the
evaluation region of Eq. (29). p∗(x, t)is the adjoint solution arised by the gradient of the
objective function Jwith respect to the source forcing sf(t). It provides the optimal position
to minimize the objective (Lemke, 2015, Sec. 3.6.1). Thus, a minimization of the objective J
may be reached in optimal sense by placing an additional monopole where p∗(x, t)2spatially
maximizes. The source region is defined as a sphere around the directive point source x0of
the CDPS model. Thus, Ψ(x)reads
Ψ(x) =
1if |x−x0| ≤ rs
0else,
(30)
where rsdenotes the source region radius.
Subsequently, the force term of the source distribution is updated at each iteration loop
based on the gradient of the objective J, given in Eq. (26), by
s[n+1]
f(t) = s[n]
f(t) + αs∆sJ, (31)
where αsdenotes a step width and nthe loop iteration number. The pressure term of the
spatial expanded source s[n+1]
pis then constructed by Eq. (28). An appropriate choice of αsis
13
2 Methods
essential for the convergence behaviour of the steepest descent approach. Further information
on the deployed line search technique can be found in Lemke (2015, Sec. 3.6.3).
2.2.5 Objective Area Mask Function
The objective area mask function σ(x, t)defines where and when the objective function
(Eq. 18 and 21) is evaluated in space and time. The objective Jis evaluated at every time
step but spatially only within a predefined area. Therefore, σ(x, t)only varies on the spatial
axes in this work:
σ(x, t) = erf h√π
l(r0−rin)i−erf h√π
l(r0−rout)i
2,r0=|x−x0|.(32)
Therein, rin denotes the effective start and rout the effective end of the evaluation area, where
σ(x, t)>0.5.lis a transition length of the mask edges and r0contains the distance between
all grid nodes and the acoustical center of the sound source defined in the CDPS model.
The objective area mask function σ(x, t)is illustrated in Fig. 2 for x0= 0.5m, rin = 0.2m,
rout = 0.3m and l= 0.05 m on a 1D grid.
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
xin m
0.0
0.2
0.4
0.6
0.8
1.0
σ(x, t)
Fig. 2: 1D objective area mask function σ(x, t)with x0= 0.5m, rin = 0.2m, rout = 0.3m
and l= 0.05 m.
14
2 Methods
2.2.6 Iterative Process
A general overview of the iterative processing chain of the adjoint-based monopole synthesis
method is depicted in Fig. 3.
add source at
xs,m=
x(arg max(|p∗|))
difference to ref
∆p=p−pref
p, ∆p
p
adjoint solution
ET
lin(q0)q∗=
[0,0,∆p]T
(direct)
Euler solution
E(q) = s[n]
update source
forcing
s[n+1]
p−s[n]
p∝p∗
- line search
p∗
s[n+1]
initial guess of
source s[0] = 0
directive sound
source
analytic solution
(CDPS model)
reference sound
field pref =p0+p0
ref
optimal repro-
duced source
sopt for m
monopole sources
if m=Mmonopole sources
and/or convergence min(∆p)
Fig. 3: Iterative process of the adjoint-based monopole synthesis method. Computationally
intensive steps including CAA methods are marked in gray.
The process starts with an initial guess of the sound sources s[0]
p= 0. The governing equations
are solved forward in time. Subsequently, the results are compared with the reference state
p0
ref(x, t)provided by the CDPS model. The obtained difference ∆p(x, t)drives the adjoint
equations, which are solved backward in time. Further, the adjoint solution provides the
gradient of ∆s[n]Jto add a sound source and update the actual forcing of the sound sources
s[n+1]. The process is repeated until the maximum amount of sound sources and convergence
in terms of the objective function Jis reached. If Jincreases during the iterative process,
15
2 Methods
the loop iteration number neval
neval =arg min
n(J(n)) (33)
is used and subsequently refined by a lower step width αsuntil convergence is achieved. Note
that in Eq. (33) Jis defined over the loop iteration numbers nas each loop produces one
value for Jthrough Eq. (18). In the evaluation of the testcases in Sec.5 only n=neval is
mentioned, but the refinement is always applied.
16
3 Simulation Settings and Setup
3 Simulation Settings and Setup
This section contains the settings of the CDPS solver and the CAA solver. An overview of
the general settings is given in Tab. 1. Note that the case relevant CDPS settings are not
listed in Tab. 1, but are described in the Secs.5.2 -5.5 of the testcases. Using the parameters
of the grid, objective area, source region and the sponge area given in Tab. 1, a 2D slice
of the x1-x2-plane (x3= 0.5m) is depicted in Fig. 4. Claiming at least six grid points per
wavelength, the cut-off frequency of ∆xis fc,∆x=c/(6 ·∆x)≈7341 Hz.
0.0 0.2 0.4 0.6 0.8 1.0
x1in m
0.0
0.2
0.4
0.6
0.8
1.0
x2in m
0.0
0.2
0.4
0.6
0.8
1.0
σ(x, t)
Fig. 4: 2D slice of the considered grid in the x1-x2-plane at x3= 0.5m. The darkgrey painted
area surrounding the boundaries represents the sponge area, the lightgrey circle area
is the source region and the colorbar running from blue to yellow indicates the mask
function values of the objective area σ(x, t).
CAA Solver
The CAA solver is discretized in a finite differences scheme in the time domain (FDTD). It
means that the computation area is represented by a grid with a finite number of nodes. The
17
3 Simulation Settings and Setup
underlying differential equation is replaced by an equivalent approximation (finite difference)
at each node (Sesterhenn, 2018, p.8). The necessary discretizations and settings to keep
the error as low as possible are shortly described: the spatial discretization is applied with a
sixth order accurate compact symmetric derivative stencil (Lele, 1992). The time discretiza-
tion is employed by an explicit fourth-order Runge-Kutta scheme for time-wise integration
(Sesterhenn, 2018, Sec. 5.4). All boundaries are handled as non-reflective using so-called
”characteristic boundary conditions” (Poinsot and Lele, 1992; Stein, 2019). To further sup-
press reflections, a quadratic sponge layer is additionally added to all boundaries, acting on
a side margin with a width of 0.2 m (Mani, 2012). The stability of the system is ensured
with an implicit filter of sixth order applied at each time step (Gaitonde and Visbal, 2000).
The CFL-number, which states the progression of a cell per time step (Courant et al., 1928),
amounts to
CFL =utot∆t
∆x=(
=0
z}|{
|u|+c)∆t
∆x=qγp0
%0
1
fs
∆x≈0.908 (34)
for the deployed constant node distance ∆x≈7.87 ·10−3m. To ensure the stability of the
method, the condition CFL <1should hold.
18
3 Simulation Settings and Setup
Tab. 1: Settings of the CDPS and CAA simulations.
variable / description solver value
xboth 0m< xi<1m for i=1,2,3
xgboth uniform, resolution 128 ×128 ×128, except case (I)
∆xboth ≈7.87 ·10−3m, equidistant
fsboth 48 kHz
%0both ≈1.21 kg / m3
p0both 101325 Pa
γboth 1.4
Rsboth 287
D(x, ω)CDPS 1
x0CDPS [0.5,0.5,0.5]Tm, except case (II)
spatial discretization CAA sixth order accurate compact symmetric derivative stencil
CFL - number CAA ≈0.908
time discretization CAA explicit fourth-order Runge-Kutta scheme
boundary condition CAA characteristic boundary condition: non-reflective
sponge CAA quadratic sponge with 0.2 m width
stability CAA implicit sixth order filter at each time step
objective area CAA rin = 0.175 m, rout = 0.275 m and l= 0.025 m
source area CAA rs= 0.15 m
19
4 Evaluation Methods
4 Evaluation Methods
The evaluation of the resulting sound fields p0
ref(x, t)and p0
opt(x, t)in Sec. 5.2 - 5.5 will be
realized in the frequency domain. Thus, a discrete (fast) Fourier transform is applied
Pref(x, ω) = Ft(p0
ref(x, t)) and Popt(x, ω) = Ft(p0
opt(x, t)).(35)
Subsequently, a normalization of each sound field is employed at f= 2 kHz and
x= [0.8,0.5,0.5]Tm, which corresponds to ϕ= 0 rad, ϑ= 0 rad and r= 0.3m in spherical
coordinates with the origin at x0= [0.5,0.5,0.5]Tm. This enables a comparison of Pref(x, ω)
and Popt(x, ω)by various evaluation measures, which are described below.
Amplitude and Phase Spectrum
The amplitude and phase spectrum will be analyzed at three specific virtual microphone
positions on the x1-x2-plane (x3= 0.5m). Fig. 4 is expanded by the microphone positions
in Fig. 5 by red dots at xmic,1 = [0.80,0.50,0.50]Tm, xmic,2 = [0.76,0.65,0.50]Tm and
xmic 3 = [0.21,0.58,0.50]Tm. The microphone positions expressed in polar coordinates on
the x1-x2-plane (x3= 0.5m) amount for microphone 1: rmic 1 = 0.3m, ϕmic 1 = 0 rad,
for microphone 2: rmic 2 = 0.3m, ϕmic 2 =π/6rad and for microphone 3: rmic 3 = 0.3m,
ϕmic 3 = 11π/12 rad.
The amplitude spectrum |P(x, ω)|is the absolute value
|P(x, ω)|=q<{P(x, ω)}2+={P(x, ω)}2(36)
and the phase spectrum is the angle ]P(x, ω) = φ(x, ω)
]P(x, ω) = φ(x, ω) = arctan ={P(x, ω)}
<{P(x, ω)}!(37)
20
4 Evaluation Methods
0.0 0.2 0.4 0.6 0.8 1.0
x1in m
0.0
0.2
0.4
0.6
0.8
1.0
x2in m
xmic 1
xmic 2
xmic 3
0.0
0.2
0.4
0.6
0.8
1.0
σ(x, t)
Fig. 5: Expanded view of Fig. 4: The virtual microphone positions are marked with red dots.
For a further explanation see Fig. 4.
of the complex frequency spectrum P(x, ω). The deviation of the reference and optimized
sound fields will then be visualized. The amplitude figures are depicted as sound pressure
level (SPL) deviation in dBrel
(Lp,opt −Lp,ref) = (20 ·log10 (|Popt(x, ω)|)−20 ·log10 (|Pref(x, ω)|)) dBrel
= 20 ·log10 |Popt(x, ω)|
|Pref(x, ω)|!dBrel (38)
and finally smoothed by one neighbouring data point. The phase will also be shown as the
deviation
(φopt(x, ω)−φref(x, ω)) rad.(39)
The phase deviation in percentage can be calculated by
φopt(x, ω)−φref(x, ω)
2π·100 %.(40)
21
4 Evaluation Methods
Typically, the tolerance range of the phase deviation is smaller than |10 %|or
|π/5|rad ≈ |0.628|rad.
2D Amplitude Directivity Pattern
The 2D directivity patterns are suitable to compare the reference and optimized sound source
at a certain frequency f. Here, the x1-x2-plane (x3= 0.5m) is considered. The directivity
pattern is always determined at a circle radius of r= 0.3m around the reference sound
source x0= [0.5,0.5,0.5]Tm. Thus, the circle radius is only applied at the azimuth angle
ϕ. The 2D directivity patterns are logarithmically scaled
Lp(x, ω) = 20 ·log10 (|P(x, ω)|)dBrel (41)
and subsequently linearly smoothed by one neighbouring data point.
3D Amplitude Directivity Pattern
The 3D amplitude directivity figures, also known as balloon plots, may be used for a quali-
tative analysis. Thus, the circle radius r= 0.3m is applied for both spherical coordinates,
the azimuth ϕand the elevation ϑ. As the balloon plots shall suite for a qualitative analysis,
no axis are employed in the figures and the absolute sound pressure with a subsequent linear
smoothing of two neighbouring data points is shown.
FDTD Source Data
The determined source information of the CAA method are the forcing of the sources sf,m(t)
and the source positions within the source region xs,m. Note that the forcing signals sf,m(t)
will only be evaluated at the monopole center positions xs,m, where the normalized Gaussian
22
4 Evaluation Methods
distribution equals to one and thus
sp,m(x, t) = sf,m(t)·sg,m(x)
|{z }
=1
=sf,m(t)for x=xs,m.(42)
Considering both information may give a comprehension of the optimization process of the
adjoint-based monopole synthesis method.
23
5 Testcases
5 Testcases
5.1 Case (I): Implementation of the CDPS Solver
To speed up the workflow and to reduce the susceptibility to errors in the conversion between
the CDPS solver in Matlab and the CAA solver written in the Fortran programming language,
a CDPS solver is embedded into the architecture of the CAA solver. Hence, the CDPS solver
could be easily parallelized by taking advantage of the existing Open MPI architecture (Gabriel
et al., 2004) of the CAA solver. This enables an efficient computation of the reference sound
fields.
Architecture
The CDPS module consists of a main routine called „sound field prediction” (sfp) that
reads and allocates the parameters given in the parameter.dat-file. Subsequently, the
main routine calls subroutines which are structured as follows:
1. Initialization of the frequency vector.
2. Computation of the driving functions. Note that the sources are uniformly driven in
this thesis with D(x0,m, ω) = 1, but included in the code as a segment for other
possible applications.
3. for-loop over all sources x0,m:
(a) Computation of the angle β(x,x0,m)between the considered sound source and
the receiver positions.
(b) Computation of the directivity pattern given in Sec.2.1.1.
(c) Computation of the sound field transfer function P(x, ω)and excitation of the
sound field by the input signal fin(t)as described in Sec. 2.1.
4. Addition of all partial sound fields to p0
out(x, t).
24
5 Testcases
5. Saving of the computed data in the .h5-format.
Comparison of the Implemetation in Matlab and Fortran
In this section the optimization of the computing time of the CDPS model in Fortran in
contrast to the existing implementation in Matlab will be investigated. Moreover, the effi-
ciency of the parallelization in Fortran is analyzed. The setup is as follows: One sound source
with a circular piston directivity pattern is employed in a 3D grid of 32 ×32 ×32 (A) and
64 ×64 ×64 (B) grid nodes. The sound field is stimulated by a sine function of f= 2 kHz
with 1000 time steps. The sampling frequency fs= 48 kHz is considered, which amounts to
a physical time span of approximately 20.8ms. The simulation is then executed once with
Matlab on one core and with Fortran parallelized on [1, 2, 4, 8, 16, 32, 64] cores.
The computing time of the considered cases is shown in Tab. 2. The CDPS solver in Fortran
is approximately twice as fast as the Matlab solver (both executed on one core). Moreover,
the execution time required by the Fortran program scales approximately linear with the used
CPU cores. The computing time over the CPU parallelism cores is shown in a double loga-
rithmic frame in Fig. 6. The deviations from the optimal linear scaling line are small. Besides,
the visually larger deviation at 64 parallelism cores is tolerable due to the short computation
time.
The Fortran CDPS solver makes use of an efficient implementation of the convolution with
respect to time in the the time domain (see Eq. (4)), while the most efficient implementation
in Matlab is to apply twice a Fourier transform and a multiplication:
p0(x, t)∗tfin(t) = F−1
t{Ft{p0(x, t)}·Ft{fin(t)}} .(43)
Note that even though the transfer function PATF(x, ω)is computed in the frequency domain,
it has to be transformed into the time domain due to the adapation to the test signal.
As outlined above, the computation time of the CDPS solver is eventually improved by the
factor 2 for one CPU core. The parallelization linearly scales up to the considered 64 CPU
cores.
25
5 Testcases
Tab. 2: Computation time of the Matlab and Fortran implementation. Grid A consists of
32 ×32 ×32 and Grid B of 64 ×64 ×64 grid points (nodes).
Programming Language Parallelism Cores Time (Grid A) Time (Grid B)
Matlab 1 44.82 s 438.05 s
Fortran 1 23.35 s 181.90 s
Fortran 2 11.26 s 94.88 s
Fortran 4 6.39 s 52.04 s
Fortran 8 3.21 s 25.49 s
Fortran 16 1.59 s 12.69 s
Fortran 32 0.80 s 6.39 s
Fortran 64 0.46 s 4.11 s
100101
Parallelism Cores
100
101
Computation Time in s
(a) 32 ×32 ×32 grid points
100101
Parallelism Cores
101
102
Computation Time in s
(b) 64 ×64 ×64 grid points
Optimal Linear Scaling CDPS Solver (Fortran)
Fig. 6: Computation time over the parallelism cores of two grid discretization schemes.
26
5 Testcases
5.2 Case (II): Nearest Distance of Neighbouring Monopoles
Sec. 2.2.3 describes that the considered monopole sound sources in the FDTD scheme are
spatially expanded. If more than one monopole has to be implemented, as it is necessary
in monopole synthesis approaches, it may lead to problematic interferences between sound
sources if they are positioned too close to each other. This case shall ensure that the source
positions xsdo not interfere with each other in the following monopole synthesis cases in
Sec. 5.3 - 5.5. It means that only the reproduction of the reference sound field is regarded,
but not the directivity pattern itself.
Specific Settings and Setup
Two monopole sound sources are placed in a specific distance to each other. Step by step,
the distance of the sources to each other is reduced until the centers of the monopole sound
sources are at neighbouring grid nodes, i.e., one grid node distance to each other. The
considered distances are listed in Tab. 3. The sources are stimulated with a band limited
white noise (0.8 - 4.0 kHz) with 1000 time steps. Additionally, the test signals are in phase
opposition to each other. Only one optimization loop is applied at each simulation of this
case.
Results
The relative objective function J/Jmax is depicted in Fig. 7 over the distance difference in
grid points (nodes). It is here defined by a normalization over the cases with Jmax as the
maximum value of all cases. A minimum is present at 11 grid nodes distance. Polar figures
of grid node distances between 7 and 25 nodes at 2kHz are shown in Fig. 9. It can be seen
that until 10 grid nodes in Fig. 9(a-d) the deviation between Lp,ref and Lp,opt is below 1 dB,
but increases with reduced grid node distance rapidly (see Fig. 9(e-f)). Considering the 1D
Gaussian distributions of the 11 grid nodes distance case in Fig. 8, an overlapping of the
27
5 Testcases
Tab. 3: Considered monopole distances.
∆grid nodes xdist in 10−3m∆grid nodes xdist in 10−3m
25 ≈196.85 10 ≈78.74
21 ≈165.35 9≈70.87
17 ≈133.85 7≈55.12
15 ≈118.11 5≈39.37
13 ≈102.36 3≈23.62
12 ≈94.49 1≈7.87
11 ≈86.61
Gaussian distributions is only observed at the edge of the Gaussian distribution where less
than 10 % of the maximum amplitude is present.
Evaluation and Discussion
The increasing deviation between pref(x, t)and popt(x, t)for less than 10 grid nodes differ-
ence is comprehensible due to the increasing overlapping of the Gaussian distributions of the
monopole sources. As the two input signals are in phase opposition, the cancellation of the
Gaussian distributions to each other enlarges with decreasing xdist.
For more than 11 grid nodes distance, the reference directivity pattern shows more constric-
tions, but also a larger source expansion, which probably causes more near field shares in
the directivity pattern. As a conclusion of the case, one may find that the distance between
the centers of the monopole sources should be at least 11 grid nodes to ensure unfavorable
influence to each other. Consequently, Eq. (30) is modified to
Ψ(x) =
1if |x−x0| ≤ rsand |x−xs,ex| ≥ 11∆x
0else,
(44)
28
5 Testcases
where xs,ex denotes the already existing monopole source center positions. Hence, Eq. (44)
is implemented at the following cases in Sec. 5.3 - 5.5.
0 5 10 15 20 25 30
∆xin grid points
0.92
0.94
0.96
0.98
1.00
J / Jmax
Fig. 7: Relative objective funtion J / Jmax over distance between two monopole sources in
grid points (nodes).
0.400 0.425 0.450 0.475 0.500 0.525 0.550 0.575 0.600
x1in m
0.0
0.2
0.4
0.6
0.8
1.0
Amplitude
xs= [0.46,0.50,0.50]Tmxs= [0.54,0.50,0.50]Tm
Fig. 8: Gaussian distribution of two monopoles with 11 grid nodes distance to each other.
29
5 Testcases
0°
45°
90°
135°
180°
225°
270°
315°
−10
−5
−2
0
dBrel
(a) 25 grid nodes
0°
45°
90°
135°
180°
225°
270°
315°
−10
−5
−2
0
dBrel
(b) 13 grid nodes
0°
45°
90°
135°
180°
225°
270°
315°
−10
−5
−2
0
dBrel
(c) 11 grid nodes
0°
45°
90°
135°
180°
225°
270°
315°
−10
−5
−2
0
dBrel
(d) 10 grid nodes
0°
45°
90°
135°
180°
225°
270°
315°
−10
−5
−2
0
dBrel
(e) 9 grid nodes
0°
45°
90°
135°
180°
225°
270°
315°
−10
−5
−2
0
dBrel
(f) 7 grid nodes
Lp,ref Lp,opt
Fig. 9: Polar directivity patterns of two spatially distant monopoles with variable xdist at
2 kHz. The polar figures are taken on the x1-x2-plane (x3= 0.5m) at the radius of
0.3 m from the center of the source at [0.5,0.5,0.5]Tm.
30
5 Testcases
5.3 Case (III): Dipole
This case deals with the adjoint-based monopole synthesis method described in Sec. 2. The
loop circle of the method (see Fig. 3) is repeated for n= 30 iteration loops while a monopole
sound source is added every third loop.
The sound field is stimulated by a band limited white noise between 0.8 kHz and 4kHz with
2000 time steps, which amounts to a physical time span of approximately 41.6ms using
the given parameters of Tab. 1. A similar case with slightly different settings was already
presented in Lemke et al. (2020, Sec. „Klassische Monopolsynthese (A)“).
Evaluation and Discussion
The results are shown for the iteration loop n= 20 as it shows the lowest relative objective
function value Jn/J0in Fig. 29 in App.A. Starting with the discussion of the forcing signals
sp,m, it can be seen in Fig. 10 that the forcing signals sp,1-sp,3are more than one order of
magnitude larger than sp,4-sp,7. An additional look at the relative objective function Jn/J0
in Fig. 29 shows that the graph of the dipole is rapidly decreasing while adding the first
three monopole sources. The following four monopole sources subsequently refine the dipole
source. Further, the source forcing of sp,2and sp,3is approximately in phase opposition to
sp,1. Including the source positions in Fig. 11, the phase opposition enables a reinforcement
in x1-direction and a reduction in x2- and x3-direction. The source positions of xs,4-xs,7are
located on the x2-x3-plane (x1=0.5 m) and refine the reduction in both directions as their
forcing signal is in phase opposition to the signals of xs,2and xs,3.
Regarding the resulting sound fields, the deviation at the three virtual microphone positions
in Fig. 12 is absolutely less than 2dBrel between 1.1 kHz and 2.9 kHz. In that range, the
phase deviation is likewise low with a maximum of absolutely 0.2 rad. The 2D directivity
patterns in Fig. 14 as well as the 3D directivity patterns in Fig. 13 confirm the behaviour
that the dipole is well modeled in the range mentioned above. The discussion of the cut-off
frequencies of the simulation settings used in the cases (III) to (V) can be found in Sec. 5.4.
Comparing the visualization of the 2D directivity pattern at 1.5 kHz with Lemke et al. (2020,
31
5 Testcases
Fig. 2), almost equal results are obtained. This supports the assumption from above that
the monopole sources xs,4-xs,7only refine the dipole model, while the xs,1-xs,3mainly
reproduce the dipole characteristic.
0246810
tin ms
−1.0
−0.5
0.0
0.5
1.0
Amplitude
sp,1
sp,2
sp,3
(a) force signals from sp,1to sp,3
0246810
tin ms
−1.0
−0.5
0.0
0.5
1.0
Amplitude
sp,4
sp,5
sp,6
sp,7
(b) force signals from sp,4to sp,7
Fig. 10: Extract between t= 0 ms and t= 10 ms of the force signals sp,m of the determined
monopole sources xs,mof the dipole.
0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75
x1in m
0.35
0.40
0.45
0.50
0.55
0.60
0.65
x2in m
xs,1
xs,2xs,3
xs,m
x0
Fig. 11: Determined source positions of xs,1-xs,3on the x1-x2-plane (x3= 0.5m) are
marked by an octagon. The green star marks the position of the reference sound
source x0. The red painted circles mark the spatial expansion of the monopole
sources and the blue circle is the source region with an radius of rs= 0.15 m. The
coordinate grid corresponds to the spatial discretization.
32
5 Testcases
1000 1500 2000 2500 3000
fin Hz
−4
−3
−2
−1
0
1
2
(Lp,opt −Lp,ref) in dBrel
(a) SPL deviation
1000 1500 2000 2500 3000
fin Hz
−0.2
0.0
0.2
(φopt −φref) in rad
(b) Phase deviation
mic 1 (r= 0.3 m, φ = 0 rad) mic 2 (r= 0.3 m, φ =π/6rad) mic 3 (r= 0.3 m, φ =11π/12 rad)
Fig. 12: SPL and Phase deviation of the reference and optimized sound field at the vir-
tual microphone positions. Microphone 1 is at [0.8,0.5,0.5]Tm, microphone 2 at
[0.76,0.65,0.50]Tm and microphone 3 at [0.21,0.58,0.50]Tm.
(a) 1.5 kHz (b) 2.0 kHz (c) 2.5 kHz (d) 3.0 kHz
(e) 1.5 kHz (f) 2.0 kHz (g) 2.5 kHz (h) 3.0 kHz
Fig. 13: Amplitude balloon plots of the reference (a-d) and optimized (e-h) dipole source.
The plots are taken at the circle radius r= 0.3m from the center of the source at
[0.5,0.5,0.5]Tm.
33
5 Testcases
0°
45°
90°
135°
180°
225°
270°
315°
−20
−10
−5
0
dBrel
(a) 1.5 kHz
0°
45°
90°
135°
180°
225°
270°
315°
−20
−10
−5
0
dBrel
(b) 2.0 kHz
0°
45°
90°
135°
180°
225°
270°
315°
−20
−10
−5
0
dBrel
(c) 2.5 kHz
0°
45°
90°
135°
180°
225°
270°
315°
−20
−10
−5
0
dBrel
(d) 3.0 kHz
Lp,ref Lp,opt
Fig. 14: Polar directivity patterns of the reference and optimized dipole source. The fig-
ures are taken at the circle radius r= 0.3m from the center of the source at
[0.5,0.5,0.5]Tm on the x1-x2-plane (x3= 0.5m).
34
5 Testcases
5.4 Case (IV): Quadrupole
The modeling of a lateral quadrupole sound source is pursued in this case, which is given
in Eq. (8). It is applied as a frequency and phase constant model, but more complex in
comparation to the dipole in Sec. 5.3, caused by the multiple constrictions of the quadrupole.
The input signal is equal to Sec. 5.3 (white noise, 0.8- 4.0kHz, 2000 time steps).
Evaluation and Discussion
The loop iteration number n= 27 is evaluated because it provides the lowest relative
objective value Jn/J0(see Fig. 29). The results of the evaluation methods described in
Sec. 4 are depicted in Fig. 15 - 19.
Analyzing the source positions and their forcing in Fig. 15 and Fig.16, a recurring algorithm
of placing and forcing the source terms can be determined. The first monopole source is
placed by the monopole synthesis method at the same position as the reference directive
source and has the strongest forcing signal. The four following added monopole sources
surround the first monopole source by
ϕ=i·π
2+π
4with i= 0,1,2,3(45)
and r≈0.089 m which corresponds to approximately 11 grid nodes. Their forcing signals
are nearly identical (see Fig. 15 (a)) and are approximately in phase opposition to sp,1(similar
to Sec. 5.3). The same behaviour may be observed for the next four monopole sources too,
where ϕ=i·(π/2) with i= 0,1,2,3and a wider radius of 16 grid nodes is present (except
xs,6with 18 grid nodes).
Thus, it is concluded that the adjoint optimization process of a quadrupole can be reduced
by only locating the first source of each source circle around xs,1or x0. The positions of the
following three monopole sources can be calculated by Eq. (45) and the nearest possible radius
to the first monopole source at the acoustic center. Hence, the amount of possible source
circles depends on the grid discretization and the given source region radius rs. Analyzing
35
5 Testcases
Fig. 29 in the App. A, the relatively steepest decrease is at n= 1,13,25, that corresponds to
the loop iteration number nwhen a circle is ”completely filled” by monopole sources. Thus,
the number of deployed monopole sources should be
m= 1 + 4i, i ≥0∈Z.(46)
Considering the resulting sound fields, it can be noticed in Fig. 17(a) that the SPL deviation
of microphone 1 between 1.7 kHz and 2.7 kHz is below 1dBrel, but increases outside that
range. A possible explanation of the large deviation at the low-frequency edge of 1 kHz is
the corresponding wavelength of 0.343 m. This wavelength has the same order of magnitude
as the computational domain and such might suffer from cut-off effects and a less functional
sponge boundary (Stein et al., 2020). An explanation of the deviation above 2.5 kHz might
be the disregard of the far-field condition suggested by Möser (2012, p. 109)
r
ls
>ls
λ,(47)
where lsis the expansion of the sound source. Assuming ls= 0.2m, the cut-off frequency
is fc≈2570 Hz. The large deviation of microphone 2 accounts to the strong constriction
of the quadrupole source at 30 deg with approximately −12 dB and is only in a small range
between 2.3 kHz and 2.7 kHz below 1dBrel. At 2.0 kHz the deviation is 2.0dBrel, but the
qualitative reproduction in Fig. 19 (b) is satisfactory at 30deg. Thus, larger deviations at the
constrictions of the directivity pattern are tolerable, e.g., if the reference states −∞ dB, a
reproduction with −15 dB is adequate (see Fig. 19(b)). The phase deviation in Fig. 17(b) is
relatively low with a maximum of ±0.2rad (exluding the collapse at 1 kHz) or expressed in
percentage (0.2/2π)·100 % ≈3.2 %.
The 2D directivity patterns in Fig. 19 confirm the observations of Fig. 17, that the reproduc-
tion of the directivy pattern at 2.0 kHz and 2.5 kHz is satisfactory, while larger deviation may
be observed at 1.5 kHz and 3.0 kHz. This behaviour is further qualitatively confirmed by the
3D directivity patterns in Fig. 18.
36
5 Testcases
0246810
tin ms
−1.0
−0.5
0.0
0.5
1.0
Amplitude
sp,1
sp,2
sp,3
sp,4
sp,5
(a) force signals from sp,1to sp,5
0246810
tin ms
−1.0
−0.5
0.0
0.5
1.0
Amplitude
sp,6
sp,7
sp,8
sp,9
(b) force signals from sp,6to sp,9
Fig. 15: Extract between t= 0 ms and t= 10 ms of the force signals sp,m of the determined
monopole sources xs,m.
0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80
x1in m
0.35
0.40
0.45
0.50
0.55
0.60
0.65
0.70
x2in m
xs,1
xs,2
xs,3
xs,4
xs,5
xs,6
xs,7
xs,8
xs,9
xs,m
x0
Fig. 16: Determined source positions of loop iteration number n= 27 in the x1-x2-plane
(x3= 0.5m) are marked by an octagon. The green star markes the position of the
reference sound source x0. The red painted circles mark the spatial expansion of
the monopoles and the blue circle is the source region with an radius of rs= 0.15 m.
The coordinate grid corresponds to the spatial discretization.
37
5 Testcases
1000 1500 2000 2500 3000
fin Hz
−6
−5
−4
−3
−2
−1
0
1
2
3
4
5
6
(Lp,opt −Lp,ref) in dBrel
(a) SPL deviation
1000 1500 2000 2500 3000
fin Hz
−0.6
−0.4
−0.2
0.0
0.2
(φopt −φref) in rad
(b) Phase deviation
mic 1 (r= 0.3 m, φ = 0 rad) mic 2 (r= 0.3 m, φ =π/6rad) mic 3 (r= 0.3 m, φ =11π/12 rad)
Fig. 17: SPL and phase deviation of the reference and optimized sound field. Microphone 1
is at [0.8,0.5,0.5]Tm, microphone 2 at [0.76,0.65,0.50]Tm and microphone 3 at
[0.21,0.58,0.50]Tm.
(a) 1.5 kHz (b) 2.0 kHz (c) 2.5 kHz (d) 3.0 kHz
(e) 1.5 kHz (f) 2.0 kHz (g) 2.5 kHz (h) 3.0 kHz
Fig. 18: Amplitude balloon plots of the reference (a-d) and optimized (e-h) quadrupole. The
plots are taken at the circle radius r= 0.3m from the center of the source at
[0.5,0.5,0.5]Tm.
38
5 Testcases
0°
45°
90°
135°
180°
225°
270°
315°
−20
−10
−5
0
dBrel
(a) 1.5 kHz
0°
45°
90°
135°
180°
225°
270°
315°
−20
−10
−5
0
dBrel
(b) 2.0 kHz
0°
45°
90°
135°
180°
225°
270°
315°
−20
−10
−5
0
dBrel
(c) 2.5 kHz
0°
45°
90°
135°
180°
225°
270°
315°
−20
−10
−5
0
dBrel
(d) 3.0 kHz
Lp,ref Lp,opt
Fig. 19: Polar directivity patterns of the reference and optimized quadrupole. The figures are
taken at the circle radius r= 0.3m from the center of the source at [0.5,0.5,0.5]Tm
on the x1-x2-plane (x3= 0.5m).
39
5 Testcases
5.5 Case (V): (Complex) Circular Piston
The circular piston model is widely used to model an idealized 2-way woofer and midrange
speaker. Its formula is given in Eq. (9), where the amplitude of the model is frequency
dependent, but the phase is constant. Usually, realistic, measured loudspeaker directivity
patterns are measured as impulse responses in the time domain and subsequently Fourier
transformed into a complex transfer function in the frequency domain. Therefore, Eq. (9)
is subsequently expanded by a frequency and a location dependent complex exponential
function:
Hccp(β(x,x0), ω) = 2J1ω
cΘ sin(β(x,x0))
ω
cΘ sin(β(x,x0)) ·ejω
cK·ejβ(x,x0)K,(48)
where Kis a scaling parameter. The chosen loudspeaker parameters are the loudspeaker
height Λy= 0.2m and the active radiating factor α= 0.82. According to Eq. (10), the
radius of the circular piston amounts to Θ=0.082 m. The complex exponential functions
result in a phase rotation in terms of the variable frequency and the radiation angle, since a
complex exponential function can be expressed by Euler’s formula
ejβ= cos(β) + jsin(β).(49)
The real circular piston (Case Va) given in Eq. (9) may be described by Eq.(48) with K= 0.0
as well. Subsequently, a complex circular piston will be analyzed for K= 0.1(Case Vb) and
for a faster phase rotation by K= 1.0(Case Vc).
Again, similar to Sec. 5.3, the test signal is a band-limited white noise (0.8 -4.0 kHz, 2000
time steps).
Evaluation and Discussion
Fig. 29 shows that the minimum of the relative objective function is present at the loop iter-
ation number n= 23 for the real and the complex circular piston and is therefore evaluated.
40
5 Testcases
First, considering the forcing signals of the real circular piston in Fig.20 together with the
source positions xs,min Tab. 4, a similar behaviour to the previous sections is observed: sp,1
is the „main” monopole source, while sp,2and sp,3have a reduced and anti-phased signal to
sp,1. Further, xs,1is surrounded by xs,2and xs,3.xs,4-xs,8are positioned on the x2-x3-plane
(x1= 0.5m) and refine the source modeling. As the circular piston is a frequency dependent
model, a closer look at the amplitude frequency spectrum of the forcing signals sp,1,sp,2and
sp,4might be of interest and is depicted as the deviation of sp,2and sp,4to sp,1in Fig. 21. It
demonstrates that almost only sp,1is activated at low frequencies around 1 kHz. Here, the
directivity pattern of the circular piston has a monopole-like character. With increasing fre-
quency the other monopole sources are continuously activated and show only approximately
3 dB difference to sp,1at 3 kHz where the first side lobes begin to spread. To clarify the evo-
lution process, Fig. 25 depicts the 3D directivity patterns for n= [3,6,9,12,15,18,21,23]
at f= 2.5kHz. As described, the first three monopoles at n= 9 (Fig. 25 (c)) already shape
the dipole character of the model at f= 2.5kHz. Following the directivity pattern is refined
until n= 23.
Analyzing the reproduction of the reference directivity pattern at the virtual microphone po-
sitions, Fig. 22 shows that they are correctly captured about the whole analyzed frequency
range with a maximum deviation of ±2dB. The polar plots in Fig. 23 confirm this behaviour,
but show deviations at the side lobe at π/2rad and 3π/2rad at f= 3 kHz. The adjoint-
based optimization gives less prominence to regions with small SPL values and can therefore
not capture the small side lobes correctly, if a low amount of monopoles is used. In this
context, it can be assumed that more monpoles are necessary. It is considered in the work of
Stein et al. (2020, Fig. 8) that includes—as mentioned in Sec. 1.2—a grid- and adjoint-based
monopole synthesis method. It means, that all grid nodes in a specific source regions are as-
sumed as monopole source, i.e., Dirac functions, with the condition that the Euler equations
hold. Comparing the results with the findings of Stein et al. (2020, Fig. 8), the observed
deviations of the polar plots are throughout larger. Due to the variable location and the
higher number of monopoles, the grid-based method is superior to the considered method in
this thesis, that uses spatially expanded monopoles. Besides a higher spatial discretization
41
5 Testcases
was applied in the work of Stein et al. (2020).
Analyzing the complex circular piston model given in Eq. (48) for K= 0.1, the amplitude
and phase spectra are shown in Fig. 26. The deviation is almost equal to the results of the
completely real circular piston case (K= 0.0). The 3D directivity patterns of the complex
circular piston in Fig. 27 show almost the same shapes as Fig. 24. Different to K= 0.1,
the amplitude and phase spectra of the complex circular piston with K= 1.0in Fig. 28
shows deviations up to approximately 4.3dBrel and 1.4rad. Thus, it can be concluded that
the method is able to reproduce complex directivity pattern for slow phase variations in the
same manner as completely real directivity pattern. However, fast phase variations obviously
require more monopoles similar to the findings with respect to the amplitude variations. In
summary, the more and stronger the variations in terms of phase and amplitude are, the
more monopoles are necessary for satisfactory reproduction of the directive sound source.
0246810
tin ms
−1.0
−0.5
0.0
0.5
1.0
Amplitude
sp,1
sp,2
sp,3
sp,4
(a) force signals from sp,1to sp,4
0246810
tin ms
−1.0
−0.5
0.0
0.5
1.0
Amplitude
sp,5
sp,6
sp,7
sp,8
(b) force signal from sp,5to sp,8
Fig. 20: Extract between t= 0 ms and t= 10 ms of the force signals sp,m of the determined
monopole sources xs,mof the real circular piston.
42
5 Testcases
1000 1500 2000 2500 3000
fin Hz
−50
−40
−30
−20
−10
0
Lpin dBrel
sp,2−sp,1
sp,4−sp,1
Fig. 21: Amplitude frequency response deviation of sp,2and sp,4to sp,1.
1000 1500 2000 2500 3000
fin Hz
−3
−2
−1
0
1
2
3
(Lp,opt −Lp,ref) in dBrel
(a) SPL deviation
1000 1500 2000 2500 3000
fin Hz
0.0
0.2
0.4
(φopt −φref) in rad
(b) Phase deviation
mic 1 (r= 0.3 m, φ = 0 rad) mic 2 (r= 0.3 m, φ =π/6rad) mic 3 (r= 0.3 m, φ =11π/12 rad)
Fig. 22: SPL and phase deviation of the reference and optimized sound field at the virtual
microphone positions of the circular piston (K= 0.0).
43
5 Testcases
0°
45°
90°
135°
180°
225°
270°
315°
−20
−10
−5
0
dBrel
(a) 1.5 kHz
0°
45°
90°
135°
180°
225°
270°
315°
−20
−10
−5
0
dBrel
(b) 2.0 kHz
0°
45°
90°
135°
180°
225°
270°
315°
−20
−10
−5
0
dBrel
(c) 2.5 kHz
0°
45°
90°
135°
180°
225°
270°
315°
−20
−10
−5
0
dBrel
(d) 3.0 kHz
Lp,ref Lp,opt
Fig. 23: Polar directivity patterns of the reference and optimized circular piston. The fig-
ures are taken at the circle radius r= 0.3m from the center of the source at
[0.5,0.5,0.5]Tm on the x1-x2-plane (x3= 0.5m).
44
5 Testcases
(a) 1.5 kHz (b) 2.0 kHz (c) 2.5 kHz (d) 3.0 kHz
(e) 1.5 kHz (f) 2.0 kHz (g) 2.5 kHz (h) 3.0 kHz
Fig. 24: Amplitude balloon plots of the reference (a-d) and optimized (e-h) circular piston
sources. The plots are taken at the circle radius r= 0.3m from the center of the
source at [0.5,0.5,0.5]Tm.
(a) n= 3 (b) n= 6 (c) n= 9 (d) n= 12
(e) n= 15 (f) n= 18 (g) n= 21 (h) n= 23
Fig. 25: Evolution process of the optimized circular piston at 2.5 kHz. The plots are taken
at the circle radius r= 0.3m from the center of the source at [0.5,0.5,0.5]Tm.
45
5 Testcases
1000 1500 2000 2500 3000
fin Hz
−3
−2
−1
0
1
2
3
(Lp,opt −Lp,ref) in dBrel
(a) SPL deviation
1000 1500 2000 2500 3000
fin Hz
0.0
0.2
0.4
(φopt −φref) in rad
(b) Phase deviation
mic 1 (r= 0.3 m, φ = 0 rad) mic 2 (r= 0.3 m, φ =π/6rad) mic 3 (r= 0.3 m, φ =11π/12 rad)
Fig. 26: SPL and phase deviation of the reference and optimized sound field at the virtual
microphone positions of the complex circular piston with K= 0.1.
(a) 1.5 kHz (b) 2.0 kHz (c) 2.5 kHz (d) 3.0 kHz
(e) 1.5 kHz (f) 2.0 kHz (g) 2.5 kHz (h) 3.0 kHz
Fig. 27: Amplitude balloon plots of the reference (a-d) and optimized (e-h) complex circular
piston with K= 0.1. The plots are taken at the circle radius r= 0.3m from the
center of the source at [0.5,0.5,0.5]Tm.
46
5 Testcases
1000 1500 2000 2500 3000
fin Hz
−4
−3
−2
−1
0
1
2
3
4
5
(Lp,opt −Lp,ref) in dBrel
(a) SPL deviation
1000 1500 2000 2500 3000
fin Hz
−1.0
−0.5
0.0
0.5
1.0
1.5
(φopt −φref) in rad
(b) Phase deviation
mic 1 (r= 0.3 m, φ = 0 rad) mic 2 (r= 0.3 m, φ =π/6rad) mic 3 (r= 0.3 m, φ =11π/12 rad)
Fig. 28: SPL and phase deviation of the reference and optimized sound field at the virtual
microphone positions of the complex circular piston with K= 1.0.
47
6 Conclusion
6 Conclusion
The present thesis introduced a novel method of monopole synthesis using an adjoint-based
CAA solver in a finite differences time domain (FDTD) discretization scheme. In addition,
a complex directivity point source (CDPS) solver in the frequency domain was implemented
into the existing environment of the computational aeroacoustics (CAA) solver to efficiently
compute reference sound fields as demonstrated in Sec. 5.1. Sec. 5.2 investigated the spatial
expansion of the monopole sources. Following, the method was evaluated for a dipole, a
quadrupole and a (complex) circular piston model.
In general, the method provided satisfactory results in the frequency range between 1.5 kHz
and 2.5 kHz for the employed simulation setup. Also, it was demonstrated that the method
is able to reproduce complex directivity patterns in the same manner as directivity patterns
consisting only of real values. However, a higher complexity of the directive source in terms of
amplitude and phase variations results in an inferior reproduction when the same amount of
monopoles is used. It can be assumed that the correctly reproduced frequency range can be
enlarged by (i) a larger computing domain, i.e., a larger spherical receiver radius to avoid near
field shares on the evaluation sphere and (ii) a higher grid resolution, i.e., more grid nodes
per meter in the computing domain to reduce the source expansion or increase the possible
amount of monopole sources. A disadvantage of the method is that positioned monopole
sources remain on their positions for the whole computation and a spatial adjustment is not
possible.
The ability of the adjoint-based monopole synthesis method to reproduce sound sources with
a high complexity, e.g., real world loudspeakers, was already demonstrated by Stein et al.
(2020). Therein, the deviations are lower to the reference directivity pattern due to the
variable location and the higher number of monopoles, or rather, Dirac pulses. In fact, the
focus of this thesis was also to show how the adjoint-based monopole synthesis method
proceeds to reproduce the directivity pattern. It principally locates monopoles around a
main monopole in the center and transfers the forcing signal of the main monopole to
the surrounding monopoles. Subsequently, the amplitude and the phase is adjusted. The
48
6 Conclusion
presented method is also able to efficiently reproduce simple sources with a low number
of monopoles, such as the dipole or the quadrupole. Moreover it was recognized that the
objective function (see Fig. 29) begins to increase at a certain iteration loop due to the
constant step width αsat the line search process. To find a variable and optimized αs, the
Gaussian distribution has to be adjoint as well. This was not investigated since the grid-based
method by Stein et al. (2020) provides more promising results.
In future, the method could be tested in a larger computing domain with a higher grid
resolution (smaller ∆x) to reproduce a reference directivity pattern with greater precision
and even synthesize very complex sources. As the available computing power is increasing
rapidly in the recent time, more accurate calculations will be possible. A grid refinement of
the source region—while the free space around the source region remain unrefined—could
improve the directivity pattern reproduction as well.
49
References
References
Ahrens, J. (2012): Analytic Methods of Sound Field Synthesis. T-Labs Series in Telecom-
munication Services. Berlin, Heidelberg: Springer-Verlag.
Berkhout, A. J.; de Vries, D.; and Vogel, P. (1992): “Wave Front Synthesis: A New Direction
in Electroacoustics.” In: Proc. of the 93rd Audio Eng. Soc. Conv.
Bilbao, S.; Ahrens, J.; and Hamilton, B. (2019): “Incorporating Source Directivity in Wave-
Based Virtual Acoustics: Time-Domain Models and Fitting to Measured Data.” In: The
Journal of the Acoustical Society of America,146(4), pp. 2692–2703.
Bilbao, S. and Hamilton, B. (2018): “Directional Sources in Wave-Based Acoustic Simula-
tion.” In: IEEE/ACM Transactions on Audio, Speech, and Language Processing,27(2),
pp. 415–428.
Courant, R.; Friedrichs, K.; and Lewy, H. (1928): “Über die partiellen Differenzengleichungen
der mathematischen Physik.” In: Mathematische Annalen,100, pp. 32–74.
de Vries, D.; Start, E. W.; and Valstar, V. G. (1994): “The Wave Field Synthesis Concept
Applied to Sound Reinforcement: Restrictions and Solutions.” In: Proc. of the 96th Audio
Eng. Soc. Conv.
Escolano, J.; Lopez, J.; and Pueo, B. (2007): “Directive Sources in Acoustic Discrete-Time
Domain Simulations Based on Directivity Diagrams.” In: JASA Express Letters,121, pp.
256–262.
Feistel, S.; Thompson, A.; and Ahnert, W. (2009): “Methods and Limitations of Line Source
Simulation.” In: J. Audio Eng. Soc.,57(6), pp. 379–402.
Gabriel, E.; et al. (2004): “Open MPI: Goals, Concept, and Design of a Next Generation
MPI Implementation.” In: Proceedings, 11th European PVM/MPI Users’ Group Meeting.
Budapest, Hungary, pp. 97–104.
50
References
Gaitonde, D. V. and Visbal, M. R. (2000): “Pade-Type Higher-Order Boundary Filter for the
Navier-Stokes Equations.” In: AIAA Journal,38, pp. 2103 – 2112.
Giles, M. B. and Pierce, N. A. (2000): “An Introduction to the Adjoint Approach to Design.”
In: Flow, Turbulence and Combustion,65(3), pp. 393–415.
Howe, M. S. (2002): Theory of Vortex Sound. Cambridge Texts in Applied Mathematics.
Cambridge University Press.
Lele, S. K. (1992): “Compact Finite Difference Schemes with Spectral-Like Resolution.” In:
Journal of Computational Physics,103(1), pp. 16 – 42.
Lemke, M. (2015): Adjoint Based Data Assimilation in Compressible Flows with Application
to Pressure Determination from PIV Data. Ph.D. thesis, Technische Universität Berlin.
Lemke, M.; Stein, L.; Hölter, A.; Straube, F.; and Weinzierl, S. (2020): “Synthese komplexer
Richtcharakteristiken für eine Schallfeldoptimierung im Zeitbereich.” In: Fortschritte der
Akustik: Tagungsband d. 46. DAGA. pp. 1192–1195.
Mani, A. (2012): “Analysis and Optimization of Numerical Sponge Layers as a Nonreflective
Boundary Treatment.” In: Journal of Computational Physics,231(2), pp. 704–716.
Meyer, D. G. (1984): “Computer Simulation of Loudspeaker Directivity.” In: J. Audio Eng.
Soc.,32(5), pp. 294–315.
Meyer, P. and Schwenke, R. (2003): “Comparison of the Directional Point Source Model
and BEM Model for Arrayed Loudspeakers.” In: Proc. of the Inst. of Acoustics,25(4).
Möser, M. (2012): Technische Akustik. 9. Berlin, Heidelberg, New York: Springer.
Olver, F. W. J.; Lozier, D. W.; Boisvert, R. F.; and Clark, C. W. (2010): NIST Handbook
of Mathematical Functions. 1. Cambridge University Press.
Opdam, R.; Latzke, K.; Behler, G.; and Vorländer, M. (2016): “Directivity as Monopole
Decomposition for Active Noise Reduction.” In: Fortschritte der Akustik: Tagungsband d.
42. DAGA. pp. 461–463.
51
References
Poinsot, T. J. and Lele, S. K. (1992): “Boundary Conditions for Direct Simulations of
Compressible Viscous Flows.” In: Journal of Computational Physics,101(1), pp. 104 –
129.
Poirier-Quinot, D.; Katz, B.; and Noisternig, M. (2017): “EVERTims: Open Source Frame-
work for Real-Time Auralization in Architectural Acoustics and Virtual Reality.” In: Proc.
of the 20th Int. Conf. on Digital Audio Effects. pp. 323–328.
Redlich, J. (2017): Zur Rekonstruktion gerichteter Schallquellen mittels Monopolsynthese.
Master’s thesis, Technische Universität Berlin.
Sakamoto, S. and Takahashi, R. (2013): “Directional Sound Source Modeling by Using
Spherical Harmonic Functions for Finite-Difference Time-Domain Analysis.” In: Proceed-
ings of Meetings on Acoustics, vol. 19. p. 015129.
Sarradj, E. (2016): Kurzskript zur Vorlesung „Luftschall-Grundlagen”. Technische Universität
Berlin. Unpublished Lecture Notes.
Schultz, F.; Straube, F.; and Spors, S. (2015): “Discussion of the Wavefront Sculpture
Technology Criteria for Straight Line Arrays.” In: Proc. of the 138th Audio Eng. Soc.
Conv., Warsaw, #9323.
Sesterhenn, J. (2018): Skript zur Vorlesung „Numerische Strömungsakustik”. Technische
Universität Berlin. Unpublished Lecture Notes.
Skudrzyk, E. (1971): The Foundations of Acoustics. New York: Springer.
Slavik, K. M. and Weinzierl, S. (2008): “Wiedergabeverfahren, 3D Audio.” In: Stefan
Weinzierl (Ed.) Handbuch der Audiotechnik, chap. 11. Berlin, Heidelberg, New York:
Springer, pp. 657–681.
Stein, L. (2019): Simulation and Modeling of a Helmholtz Resonator under Grazing Turbulent
Flow. Ph.D. thesis, Technische Universität Berlin.
52
References
Stein, L.; Straube, F.; Sesterhenn, J.; Weinzierl, S.; and Lemke, M. (2019): “Adjoint-Based
Optimization of Sound Reinforcement Including Non-Uniform Flow.” In: The Journal of
the Acoustical Society of America,146(3), pp. 1774–1785.
Stein, L.; Straube, F.; Weinzierl, S.; and Lemke, M. (2020): “Adjoint-Based FDTD Direc-
tional Sound Source Modeling Using Euler Equations.” In: The Journal of the Acoustical
Society of America. Submitted.
Straube, F. (2019): Optimized Geometric and Electronic Wavefront Shaping with Line Source
Arrays for Large-Scale Sound Reinforcement. Ph.D. thesis, Technische Universität Berlin.
Takeuchi, D.; Yatabe, K.; and Oikawa, Y. (2019): “Source Directivity Approximation for
Finite-Difference Time-Domain Simulation by Estimating Initial Value.” In: The Journal
of the Acoustical Society of America,145(4), pp. 2638–2649.
Tam, C. K. W. and Webb, J. C. (1993): “Dispersion-Relation-Preserving Finite Difference
Schemes for Computational Acoustics.” In: Journal of Computational Physics,107, pp.
262–281.
van Beuningen, G. W. J. and Start, E. W. (2000): “Optimizing Directivity Properties of DSP
Controlled Loudspeaker Arrays.” In: Proc. of the Inst. of Acoustics: Reproduced Sound,
22(6), pp. 17–37.
Williams, E. G. (1999): Fourier Acoustics – Sound Radiation and Nearfield Acoustical Holog-
raphy. Academic Press.
53
A Objective Functions and Source Positions
A Objective Functions and Source Positions
Fig. 29 shows the relative objective function Jn/J0over the iteration loop numbers nfor the
cases (III) - (V). The source positions xs,mup to the minimum of Jn/J0are given in Tab. 4
and Tab. 5.
1 4 7 10 13 16 19 22 25 28
loop iteration number n
10−1
100
Jn/J0
Case (III): Dipole
Case (IV): Quadrupole
Case (Va): Real Circular Piston (K= 0.0)
Case (Vb): Complex Circular Piston (K= 0.1)
Case (Vc): Complex Circular Piston (K= 1.0)
Fig. 29: Relative objective function Jn/J0over the iteration loops n. A new monopole source
is added at every tick of the loop iteration number n.
54
A Objective Functions and Source Positions
Tab. 4: Monopole source positions of the cases (III) –(Va) determined by the adjoint-based
monopole synthesis method.
Case
Source (III): Dipole (IV): Quadrupole (Va): Circ. Piston (K= 0.0)
xs,1in m [0.50,0.50,0.50]T[0.50,0.50,0.50]T[0.46,0.50,0.50]T
xs,2in m [0.42,0.50,0.50]T[0.57,0.57,0.50]T[0.42,0.50,0.50]T
xs,3in m [0.59,0.50,0.51]T[0.44,0.44,0.50]T[0.59,0.50,0.51]T
xs,4in m [0.50,0.50,0.42]T[0.44,0.57,0.50]T[0.50,0.50,0.42]T
xs,5in m [0.50,0.58,0.46]T[0.57,0.44,0.50]T[0.50,0.58,0.46]T
xs,6in m [0.50,0.42,0.46]T[0.50,0.65,0.50]T[0.50,0.42,0.46]T
xs,7in m [0.50,0.58,0.55]T[0.50,0.38,0.50]T[0.50,0.57,0.56]T
xs,8in m –[0.63,0.50,0.50]T[0.50,0.49,0.59]T
xs,9in m –[0.38,0.50,0.50]T–
Tab. 5: Monopole source positions of the cases (Vb) –(Vc) determined by the adjoint-based
monopole synthesis method.
Case
Source (Vb): Cmp. Circ. Piston (K= 0.1) (Vc): Cmp. Circ. Piston (K= 1.0)
xs,1in m [0.50,0.50,0.50]T[0.47,0.50,0.50]T
xs,2in m [0.59,0.50,0.51]T[0.38,0.50,0.50]T
xs,3in m [0.42,0.50,0.50]T[0.58,0.50,0.50]T
xs,4in m [0.50,0.50,0.42]T[0.48,0.57,0.59]T
xs,5in m [0.50,0.58,0.46]T[0.48,0.58,0.54]T
xs,6in m [0.50,0.42,0.46]T[0.49,0.57,0.46]T
xs,7in m [0.50,0.57,0.56]T[0.48,0.43,0.54]T
xs,8in m [0.50,0.49,0.59]T[0.48,0.50,0.42]T
55
