Sound field synthesis (SFS) includes all methods that try to generate a defined
sound field in an extended area that is surrounded by loudspeakers. This page
focuses on those methods that provide analytical solutions to the underlying
mathematical problem, namely WFS, NFC-HOA, and the SDM.

The Toolboxes provide you with the implementation of the underlying mathematics.
You can make numerical simulations of the resulting sound fields and can even
create binaural simulations of the same sound fields. This enables you to listen
to large loudspeaker arrays, even if you don’t have one in your laboratory or at
home. In addition, you can easily plug-in your own algorithms in order to test
or compare them.

The SFS Toolbox project is structured in the following three sub-projects.

Most of the figures in this page are directly created by the SFS Toolbox for
Python. All of them display the corresponding code for creating them directly
before the actual figure. In order to recreate them, you have to execute the
following code first:

Figure shows the coordinate system that is used in the following
chapters. A vector \(\x\) can be described by its position
\((x,y,z)\) in space or by its length, azimuth angle
\(\phi \in [0,2\pi[\), and elevation
\(\theta \in \left[-\frac{\pi}{2},\frac{\pi}{2}\right]\).
The azimuth is measured counterclockwise and elevation is positive for
positive \(z\)-values.

Fig. 1 Coordinate system used in this document. The vector \(x\) can also be
described by its length, its azimuth angle \(\phi\), and its elevation
\(\theta\).

Fig. 2 Illustration of the geometry used to discuss the physical fundamentals of
sound field synthesis and the single-layer potential.

The problem of sound field synthesis can be formulated after as follows. Assume
a volume \(V \subset \mathbb{R}^n\) which is free of any sources and sinks,
surrounded by a distribution of monopole sources on its surface \(\partial
V\). The pressure \(P(\x,\w)\) at a point \(\x\in V\) is then given
by the single-layer potential (compare p. 39 in [ColtonKress1998])

where \(G(\x-\x_0,\w)\) denotes the sound propagation of the source at
location \(\x_0 \in \partial V\), and \(D(\x_0,\w)\) its weight,
usually referred to as driving function. The sources on the surface are called
secondary sources in sound field synthesis, analogue to the case of acoustical
scattering problems. The single-layer potential can be derived from the
Kirchhoff-Helmholtz integral [Williams1999]. The challenge in sound field
synthesis is to solve the integral with respect to \(D(\x_0,\w)\) for a
desired sound field \(P = S\) in \(V\). It has unique solutions which
[ZotterSpors2013] explicitly showed for the spherical case and [Fazi2010]
(Chap.4.3) for the planar case.

In the following the single-layer potential for different dimensions is
discussed. An approach to formulate the desired sound field \(S\) is
described and finally it is shown how to derive the driving function
\(D\).

The integral equation (3) states a Fredholm equation of first
kind with a Green’s function as kernel. This type of equation can be solved in a
straightforward manner for geometries that have a complete set of orthogonal
basis functions. Then the involved functions are expanded into the basis
functions \(\psi_n\) as [MorseFeshbach1981], p. (940)

where \(\tilde{G}_n, \tilde{D}_n, \tilde{S}_n\) denote the series expansion
coefficients, \(n \in \mathbb{Z}\), and \(\langle\psi_n,
\psi_{n'}\rangle = 0\,\) for \(n \ne n'\).
If the underlying space is not compact the equations will involve an integration
instead of a summation

This means that the Fredholm equation (3) states a convolution.
For geometries where the required orthogonal basis functions exist, (10)
follows directly via the convolution theorem [ArfkenWeber2005], eq. (1013).
Due to the division of the desired sound field by the spectrum of the Green’s
function this kind of approach has been named SDM[AhrensSpors2010]. For
circular and spherical geometries the term NFC-HOA is more common due to the
corresponding basis functions. “Near-field compensated” highlights the usage of
point sources as secondary sources in contrast to Ambisonics and HOA that
assume plane waves as secondary sources.

The challenge is to find a set of basis functions for a given geometry.
In the following paragraphs three simple geometries and their widely
known sets of basis functions will be discussed.

where \(P_n^{|m|}\) are the associated Legendre functions. Note that
this function may also be defined in a slightly different way, omitting
the \((-1)^m\) factor, see for example [Williams1999], eq. (6.20).

The complex conjugate of \(Y_n^m\) is given by negating the degree
\(m\) as

For a spherical secondary source distribution with a radius of \(R_0\) the
sound field can be calculated by a convolution along the surface. The driving
function is then given by a simple division as [Ahrens2012], eq. (3.21) [2]

where \(\breve{S}_n^m\) denote the spherical expansion coefficients of the
source model, \(\theta_\text{s}\), \(\phi_\text{s}\), and
\(r_\text{s}\) its directional dependency, and \(\breve{G}_n^0\) the
spherical expansion coefficients of a secondary monopole source located at
the north pole of the sphere \(\x_0 = (\frac{\pi}{2},0,R_0)\). For a point
source this is given as [SchultzSpors2014], eq. (25)

For a circular secondary source distribution with a radius of \(R_0\) the
driving function can be calculated by a convolution along the surface of the
circle as explicitly shown by [AhrensSpors2009a] and is then given as

where \(\breve{S}_m\) denotes the circular expansion coefficients for the
source model, \(\phi_\text{s}\), and \(r_\text{s}\) its directional
dependency, and \(\breve{G}_m\) the circular expansion coefficients for a
secondary monopole source. For a line source located at \(\x_0 = (0,R_0)\)
this is given as

where \(\breve{S}\) denotes the planar expansion coefficients for the source
model, \(y_\text{s}\) its positional dependency, and \(\breve{G}\) the
planar expansion coefficients of a secondary point source with
[SchultzSpors2014], eq. (49)

For the planar and the following linear geometries the Fredholm equation is
solved for a non compact space \(V\), which leads to an infinite and
non-denumerable number of basis functions as opposed to the denumerable case for
compact spaces [SchultzSpors2014].

For an infinitely long secondary source distribution located on the
\(x\)-axis the driving function for \({\mathbb{R}}^2\) can be calculated
by a convolution along this axis as [Ahrens2012], eq. (3.73)

where \(\breve{S}\) denotes the linear expansion coefficients for the source
model, \(y_\text{s}\), \(z_\text{s}\) its positional dependency, and
\(\breve{G}\) the linear expansion coefficients of a secondary line source
with

The single-layer potential (3) satisfies the homogeneous
Helmholtz equation both in the interior and exterior regions \(V\) and
\(V^* {\mathrel{\!\mathop:}=}{\mathbb{R}}^n \setminus (V \cup \partial V)\,\).
If \(D(\x_0,\w)\) is continuous, the pressure \(P(\x,\w)\) is
continuous when approaching the surface \(\partial V\) from the inside and
outside. Due to the presence of the secondary sources at the surface
\(\partial V\), the gradient of \(P(\x,\w)\) is discontinuous when
approaching the surface. The strength of the secondary sources is then given by
the differences of the gradients approaching \(\partial V\) from both sides
as [FaziNelson2013]

where \(\partial_\n{\mathrel{\mathop:}=}\scalarprod{\nabla}{\n}\) is
the directional gradient in direction \(\n\) – see Fig. 2.
Due to the symmetry of the problem the solution for an infinite planar boundary
\(\partial V\) is given as

where the pressure in the outside region is the mirrored interior pressure given
by the source model \(S(\x,\w)\) for \(\x\in V\). The integral
equation resulting from introducing (28) into (3) for a
planar boundary \(\partial V\) is known as Rayleigh’s first integral
equation. This solution is identical to the explicit solution for planar
geometries (21) in \({\mathbb{R}}^3\) and for linear
geometries (25) in \({\mathbb{R}}^2\).

A solution of (27) for arbitrary boundaries can be found by applying
the Kirchhoff or physical optics approximation[ColtonKress1983], p. 53–54.
In acoustics this is also known as determining the visible elements for the
high frequency boundary element method [Herrin2003]. Here, it is assumed that
a bent surface can be approximated by a set of small planar surfaces for
which (28) holds locally. In general, this will be the case if the wave
length is much smaller than the size of a planar surface patch and the position
of the listener is far away from the secondary sources. [3] Additionally, only
one part of the surface is active: the area that is illuminated from the
incident field of the source model.

The outlined approximation can be formulated by introducing a window function
\(w(\x_0)\) for the selection of the active secondary sources
into (28) as

In the SFS Toolbox we assume convex secondary source distributions, which
allows to formulate the window function by a scalar product with the normal
vector of the secondary source distribution. In general, also non-convex
secondary source distributions can be used with WFS – compare the appendix in
[LaxFeshbach1947] and [4].

One of the advantages of the applied approximation is that due to its local
character the solution of the driving function (28) does not depend on
the geometry of the secondary sources. This dependency applies to the direct
solutions presented in Solution for Special Geometries: NFC-HOA and SDM.

The single-layer potential (3) is valid for all \(V \subset
{\mathbb{R}}^n\). Consequentially, for practical applications a two-dimensional
(2D) as well as a three-dimensional (3D) synthesis is possible. Two-dimensional
is not referring to a synthesis in a plane only, but describes a setup that is
independent of one dimension. For example, an infinite cylinder is independent
of the dimension along its axis. The same is true for secondary source
distributions in 2D synthesis. They exhibit line source characteristics and are
aligned in parallel to the independent dimension. Typical arrangements of such
secondary sources are a circular or a linear setup.

The characteristics of the secondary sources limit the set of possible sources
which can be synthesized. For example, when using a 2D secondary source setup it
is not possible to synthesize the amplitude decay of a point source.

For a 3D synthesis the involved secondary sources depend on all dimensions and
exhibit point source characteristics. In this scenario classical secondary
sources setups would be a sphere or a plane.

Fig. 3 Sound pressure in decibel for secondary source distributions with different
dimensionality all driven by the same signals. The sound pressure is color
coded, lighter color corresponds to lower pressure. In the 3D case a planar
distribution of point sources is applied, in the 2.5D case a linear
distribution of point sources, and in the 2D case a linear distribution of
line sources.

In practice, the most common setups of secondary sources are 2D setups,
employing cabinet loudspeakers. A cabinet loudspeaker does not show the
characteristics of a line source, but of a point source. This dimensionality
mismatch prevents perfect synthesis within the desired plane. The combination of
a 2D secondary source setup with secondary sources that exhibit 3D
characteristics has led to naming such configurations 2.5D synthesis[Start1997]. Such scenarios are associated with a wrong amplitude decay due to
the inherent mismatch of secondary sources as is highlighted in
Fig. 3. In general, the amplitude is only
correct at a given reference point \(\xref\).

For a circular secondary source distribution with point source characteristic
the 2.5D driving function can be derived by introducing expansion coefficients
for the spherical case into the driving function (17). The equation
is than solved for \(\theta = 0{^\circ}\) and \(r_\text{ref} = 0\). This
results in a 2.5D driving function given as [Ahrens2012], eq. (3.49)

For a linear secondary source distribution with point source characteristics the
2.5D driving function is derived by introducing the linear expansion
coefficients for a monopole source (43) into
the driving function (25) and solving the equation for \(y =
y_\text{ref}\) and \(z = 0\). This results in a 2.5D driving function given
as [Ahrens2012], eq. (3.77)

A driving function for the 2.5D situation in the context of WFS and arbitrary
2D geometries of the secondary source distribution can be achieved by applying
the far-field approximation \(\Hankel{2}{0}{\zeta} \approx
\sqrt{\frac{2\i}{\pi\zeta}} \e{-\i\zeta}\) for
\(\zeta \gg 1\) to the 2D Green’s function [Williams1999], eq. (4.23).
Using this the following relationship between the 2D and 3D Green’s functions
can be established.

Knowing the pressure field of the desired source \(S(\x,\w)\) is
required in order to derive the driving signal for the secondary source
distribution. It can either be measured, i.e. recorded, or modeled. While the
former is known as data-based rendering, the latter is known as model-based
rendering. For data-based rendering, the problem of how to capture a complete
sound field still has to be solved. [Avni2013] et al. discuss some influences
of the recording limitations on the perception of the reproduced sound field.
This document will consider only model-based rendering.

Frequently applied models in model-based rendering are plane waves, point
sources, or sources with a prescribed complex directivity. In the following the
models used within the SFS Toolbox are presented.

xs=0,0,0# position of sourceomega=2*np.pi*800# frequencygrid=sfs.util.xyz_grid([-1.75,1.75],[-1.75,1.75],0,spacing=0.02)p=sfs.mono.source.point(omega,xs,[],grid)normalization=4*np.pisfs.plot.soundfield(normalization*p,grid);

xs=0,0,0# position of sourcens=sfs.util.direction_vector(0)# direction of sourceomega=2*np.pi*800# frequencygrid=sfs.util.xyz_grid([-1.75,1.75],[-1.75,1.75],0,spacing=0.02)p=sfs.mono.source.point_dipole(omega,xs,ns,grid)sfs.plot.soundfield(p,grid);

The source model for a three dimensional dipole source is given by the
directional derivative of the three dimensional Green’s function with respect to
\({\n_\text{s}}\) defining the orientation of the dipole source.

xs=0,0,0# position of sourceomega=2*np.pi*800# frequencygrid=sfs.util.xyz_grid([-1.75,1.75],[-1.75,1.75],0,spacing=0.02)p=sfs.mono.source.line(omega,xs,None,grid)normalization=np.sqrt(8*np.pi*omega/sfs.defs.c)*np.exp(1j*np.pi/4)sfs.plot.soundfield(normalization*p,grid);

In the following, driving functions for NFC-HOA and SDM are derived
for spherical, circular, and linear secondary source distributions. Among the
possible combinations of methods and secondary sources not all are meaningful.
Hence, only the relevant ones will be presented. The same holds for the
introduced source models of plane waves, point sources, line sources and focused
sources. [AhrensSpors2010] in addition have considered SDM
driving functions for planar secondary source distributions.

For NFC-HOA, temporal-domain implementations for the
2.5D cases are available for a plane wave and a point source as source models. The
derivation of the implementation is not explicitly shown here, but is described
in [Spors2011].

For a spherical secondary source distribution with radius \(R_0\) the
spherical expansion coefficients of a plane
wave (37) and of the Green’s function for a
point source (14) are inserted into (13) and yield
[SchultzSpors2014], eq. (A3)

For a circular secondary source distribution with radius \(R_0\) the
circular expansion coefficients of a plane
wave (38) and of the Green’s function for a
line source (18) are inserted into (17) and yield
[AhrensSpors2009a], eq. (16)

For a circular secondary source distribution with radius \(R_0\) and point
source as Green’s function the 2.5D driving function is given by inserting the
spherical expansion coefficients for a plane
wave (37) and a point
source (42) into (30) as

For an infinite linear secondary source distribution located on the
\(x\)-axis the 2.5D driving function is given by inserting the linear
expansion coefficients for a point source as Green’s
function (43) and a plane
wave (39) into (31) and exploiting
the fact that \((\wc )^2 - k_{x_\text{s}}\) is constant.
Assuming \(0 \le |k_{x_\text{s}}| \le |\wc |\) this results in
[AhrensSpors2010], eq. (17)

The advantage of this result is that it can be implemented by a simple weighting
and delaying of the signal, plus one convolution with \(h(t)\). The same
holds for the driving functions of WFS as presented in the next section.

xs=0,2.5,0# position of sourceomega=2*np.pi*1000# frequencyR0=1.5# radius of secondary sourcesx0,n0,a0=sfs.array.circular(200,R0)grid=sfs.util.xyz_grid([-1.75,1.75],[-1.75,1.75],0,spacing=0.02)d=sfs.mono.drivingfunction.nfchoa_25d_point(omega,x0,R0,xs)p=sfs.mono.synthesized.generic(omega,x0,n0,d*a0,grid,source=sfs.mono.source.point)normalization=20sfs.plot.soundfield(normalization*p,grid);sfs.plot.secondarysource_2d(x0,n0,grid)

For a spherical secondary source distribution with radius \(R_0\) the
spherical coefficients of a point
source (42) and of the Green’s
function (14) are inserted into (13) and yield
[Ahrens2012], eq. (5.7) [2]

For a circular secondary source distribution with radius \(R_0\) and point
source as secondary sources the 2.5D driving function is given by inserting the
spherical coefficients (42)
and (14) into (30). This results in [Ahrens2012],
eq. (5.8)

For an infinite linear secondary source distribution located on the
\(x\)-axis and point sources as secondary sources the 2.5D driving function
for a point source is given by inserting the corresponding linear expansion
coefficients (43) and (26)
into (31). Assuming \(0 \le |k_x| < |\wc |\) this
results in [Ahrens2012], eq. (4.53)

For a spherical secondary source distribution with radius \(R_0\) the spherical
coefficients of a line source (48) and of
the Green’s function (14) are inserted into (13) and
yield [Hahn2015], eq. (20)

For a circular secondary source distribution with radius \(R_0\) and point
sources as secondary sources the 2.5D driving function is given by inserting the
spherical coefficients (48)
and (14) into (30) as [Hahn2015], eq. (23)

For an infinite linear secondary source distribution located on the
\(x\)-axis and line sources as secondary sources the driving function is
given by inserting the linear coefficients (50)
and (26) into (25) as

Focused sources mimic point or line sources that are located inside the audience
area. For the single-layer potential the assumption is that the audience area is
free from sources and sinks. However, a focused source is neither of them. It
represents a sound field that converges towards a focal point and diverges
afterwards. This can be achieved by reversing the driving function of a point or
line source in time which is known as time reversal focusing [Yon2003].

Nonetheless, the single-layer potential should not be solved for focused sources
without any approximation. In the near field of a source, evanescent waves
appear for spatial frequencies \(k_x > |\wc |\)[Williams1999],
eq. (24). They decay exponentially with the distance from the source. An exact
solution for a focused source is supposed to include these evanescent waves
around the focal point. That is only possible by applying very large amplitudes
to the secondary sources, compare Fig. 2a in [SporsAhrens2010]. Since the
evanescent waves decay rapidly and are hence not influencing the perception,
they can easily be omitted. For corresponding driving functions for focused
sources without the evanescent part of the sound field see [SporsAhrens2010]
for SDM and [AhrensSpors2009b] for NFC-HOA.

In the SFS Toolbox only focused sources in WFS are considered at the moment.

In the following, the driving functions for WFS in the frequency and temporal
domain for selected source models are presented. The temporal domain functions
consist of a filtering of the source signal and a weighting and delaying of the
individual secondary source signals. This property allows for a very efficient
implementation of WFS driving functions in the temporal domain. It is one of the
main advantages of WFS in comparison to most of the NFC-HOA, SDM solutions
discussed above.

It has the advantage that its temporal domain version could again be implemented
as a simple weighting- and delaying-mechanism.

To reach at 2.5D for a point source, we will start in 3D and apply stationary
phase approximations instead of directly using (34) – see discussion
after [Schultz2016], (2.146). Under the assumption of \(\frac{\omega}{c}
(|\x_0-\xs| + |\x-\x_0|) \gg 1\) it then follows [Schultz2016], eq. (2.137),
[Start1997], eq. (3.10, 3.11)

whereby \(\xref\) is a reference point at which the synthesis is correct.
A second stationary phase approximation can be applied to reach at
[Schultz2016], eq. (2.131, 2.141), [Start1997], eq. (3.16, 3.17)

which is the traditional formulation of a point source in WFS as given by eq.
(2.27) in [Verheijen1997][6]. Now \(d_\text{ref}\) is the distance of a
line parallel to the secondary source distribution and \(d_\text{s}\) the
shortest possible distance from the point source to the linear secondary source
distribution.

The default WFS driving functions for a point source in the SFS Toolbox are
(75) and (76). Transferring both to the temporal
domain via an inverse Fourier transform (2) it follows

xs=0,2.5,0# position of sourceomega=2*np.pi*1000# frequencyx0,n0,a0=sfs.array.circular(200,1.5)grid=sfs.util.xyz_grid([-1.75,1.75],[-1.75,1.75],0,spacing=0.02)d=sfs.mono.drivingfunction.wfs_2d_line(omega,x0,n0,xs)a=sfs.mono.drivingfunction.source_selection_line(n0,x0,xs)twin=sfs.tapering.tukey(a,.3)p=sfs.mono.synthesized.generic(omega,x0,n0,d*twin*a0,grid,source=sfs.mono.source.point)normalization=7sfs.plot.soundfield(normalization*p,grid);sfs.plot.secondarysource_2d(x0,n0,grid)

where \(|\n_\text{s}| = 1\). For a 2D or 2.5D secondary source setup and
a line source orientation perpendicular to the plane where the
secondary sources are located this automatically simplifies to \(\vec{v} =
\x_0 - \xs\).

By inserting the source model for a line source (47) into (28)
and (34) and calculating the derivate of the Hankel function after eq.
(9.1.20) in [AbramowitzStegun1972] it follows

As mentioned before, focused sources exhibit a field that converges in a focal
point inside the audience area. After passing the focal point, the field becomes
a diverging one as can be seen in Fig. 13. In order
to choose the active secondary sources, especially for circular or spherical
geometries, the focused source also needs a direction \(\n_\text{s}\).

The driving function for a focused source is given by the time-reversed
versions of the driving function for a point source (75) and
(76) as

where \(d_\text{ref}\) is the distance of a line parallel to the secondary
source distribution and \(d_\text{s}\) the shortest possible distance from
the focused source to the linear secondary source distribution.

Transferred to the temporal domain via an inverse Fourier transform (2) it
follows

The reproduction accuracy of WFS is limited due to practical aspects. For the
audible frequency range the desired sound field can not be synthesized
aliasing-free over an extended listening area, which is surrounded by a discrete
ensemble of individually driven loudspeakers. However, it is suitable for
certain applications to increase reproduction accuracy inside a smaller (local)
listening region while stronger artifacts outside are permitted. This approach
is termed LSFS in general.

The implemented Local Wave Field Synthesis method utilizes focused sources as a
distribution of virtual loudspeakers which are placed more densely around the
local listening area. These virtual loudspeakers can be driven by conventional
SFS techniques, like e.g. WFS or NFC-HOA. The results are similar to
band-limited NFC-HOA, with the difference that the form and position of the
enhanced area can freely be chosen within the listening area.

The set of focused sources is treated as a virtual loudspeaker distribution and
their positions \({\x_\text{fs}}\) are subsumed under
\(\mathcal{X}_{\mathrm{fs}}\). Therefore, each focused source is driven
individually by \(D_\text{l}({\x_\text{fs}}, \w)\), which in principle
can be any driving function for real loudspeakers mentioned in previous
sections. At the moment however, only WFS and NFC-HOA driving functions are
supported. The resulting driving function for a loudspeaker located at
\(\x_0\) reads

which is superposition of the driving function
\(D_{\mathrm{fs}}(\x_0,{\x_\text{fs}},\w)\) reproducing a single focused
source at \({\x_\text{fs}}\) weighted by \(D_\text{l}({\x_\text{fs}},
\w)\). Former is derived by replacing \(\xs\) with
\({\x_\text{fs}}\) in the WFS driving functions and for focused sources.
This yields

The solution mentioned by [LaxFeshbach1947] assumes that the listener is
far away from the radiator and that the radiator is a physical source not a
notional one as the secondary sources. In this case the selection criterion
has to be chosen more carefully, incorporating the exact position of the
listener and the virtual source. See also the related discussion.

Note that Williams defines the Fourier transform with transposed signs as
\(F(\w) = \int f(t) \e{\i\w t}\). This leads also to changed signs in
his definitions of the Green’s functions and field expansions.