2.5D Local wave field synthesis of a virtual plane wave
using a time domain representation of spherical harmonics expansion

Nara HAHN; Fiete WINTER; Sascha SPORS
Institute of Communications Engineering
University of Rostock, Rostock, Germany
nara.hahn@uni-rostock.de

ABSTRACT
Wave field synthesis is a well established sound reproduction technique, where the wavefront of a desired soundfield is physically recreated within an extended target area by using a large number of individually driven loudspeakers. For a non-continuous loudspeaker distribution, the synthesized sound field typically exhibits spatial
aliasing artifacts which occur above an upper frequency limit. If the synthesis in a relatively small area is concerned, called local wave field synthesis, the physical accuracy of the sound field can be improved significantly.
This was achieved in a recently proposed method, where the target sound field is expanded into a finite numberof cylindrical harmonic components with respect to a properly chosen expansion center. The present paper deals
with the time domain realization of this method for a 2.5-dimensional setup. The driving function for a virtual
plane wave is derived based on the time domain version of the spherical harmonics expansion. Due to the low
computational cost, the proposed approach is well suited for real-time implementation. The spatial, temporal,
and spectral properties of the synthesized sound fields are evaluated by numerical simulations.

Keywords: local wave field synthesis, spherical harmonics expansion, time domain

1

INTRODUCTION

Wave Field Synthesis (WFS), sometimes called wavefront synthesis, is a sound field synthesis method which
is based on the high frequency approximation of the Kirchhoff-Helmholtz integral equation [1, 2]. While thehigh spatial fidelity of WFS is consistently supported by listening experiments [3–6], the timbral degradationis still a challenge [7, 8]. This is mainly attributed to the discrete distribution of secondary sources, whichcauses spectral distortions, called spatial aliasing artifacts. The artifacts are pronounced at high frequencies,
where, roughly speaking, the half wavelength is smaller than the distance between the loudspeakers. Unless acontinuous (or at least quasi-continuous in a practical sense) loudspeaker array is developed, spatial aliasing isinevitable and should be accepted as a physical limitation.
Although the spatial aliasing artifacts cannot be removed completely, it is possible to spatially redistribute themin a controlled way. If the requirement for an accurate synthesis is relaxed to a smaller area, the performancecan be improved at the expanse of stronger artifacts occurring elsewhere. Such a technique is called local WFS
and a couple of methods were proposed in the past years [9–11]. The present study focuses on the method
using spatial band limitation [11], where the target sound field is expanded in terms of spherical/cylindricalharmonics with respect to the target position, and the size of the aliasing-free area is controlled by varying the
maximum order of the expansion. The driving function of this method was initially derived in the temporalfrequency domain, which demands intense computations due to the evaluation of the Bessel functions at all
frequency bins. Two time domain realizations are introduced in [12, 13], one of which exploits the time domain
representation of the cylindrical harmonics expansion [13]. While it exhibits significantly lower computations,some tweaks are needed to cope with temporal aliasing and excessive amplitude.

1132

In this paper, an improved version of the time domain driving function is proposed for local WFS. The proposed driving function is derived in Sec. 2. The so called 2.5-dimensional (2.5D) setup is considered, meaningthat a 2D target sound field is synthesized by means of secondary point sources distributed in the horizontal
plane. Similar to [13], the target sound field is initially described in terms of cylindrical harmonics with alimited spatial bandwidth, i.e. finite number of cylindrical harmonics. The spatial windowing associated withthe directional gradient and secondary source selection is performed in the cylindrical harmonics domain, andthe resulting coefficient is converted into the spherical harmonics domain, followed by amplitude correction andequalization. In Sec. 3, the proposed driving function is used in numerical simulations in order to evaluate theaccuracy of the synthesized sound fields. The paper closes with a brief conclusion in Sec. 4.

1.1 Nomenclature
The following notational conventions are used throughout this paper. Frequency domain sound fields are denotedby uppercase S(x, ω) and time domain sound fields by lowercase s(x, t), where ω = 2πf denotes the angularfrequency and t the time. Position vectors x = (r, θ, φ) are denoted by lowercase boldface. Unless statedotherwise, the spherical coordinates are used, where θ denotes the colatitude and φ the azimuth. A scalar
product of two vectors is denoted by h·, ·i. The imaginary unit is denoted by i and the speed of sound by c.
The cylindrical harmonics expansion of a 2D interior sound field reads [14, Eq. (4.49)]
S(x, ω) =∞X

S̊m (ω) i−m Jm ( ωc r sin θ) eimφ ,

(1)

m=−∞

where Jm (·) denotes the cylindrical Bessel function of the first kind, and S̊m (ω) the expansion coefficient. The
spherical harmonics expansion of an interior sound field reads [14, Eq. (6.21)]
S(x, ω) =n∞ X
XS̆nm (ω) i−n jn ( ωc r) Ynm (θ, φ)(2)

n=0 m=−nwhere jn (·) denotes the spherical Bessel function of the first kind, and S̆nm (ω) the expansion coefficient. The
spherical harmonic of degree n and order m is defined asq(n−m)! mimφYnm (θ, φ) = 2n+1,
(3)
4π (n+m)! Pn (cos θ)ewhere Pnm (·) denotes the associated Legendre polynomial.

2

DRIVING FUNCTIONS

2.1 Conventional WFS
The conventional 2.5D driving function for a virtual plane wave propagating in the direction npw = (1, π2 , φpw )reads [15, Eq. (2.177)]qp
ω
DWFS (x0 , ω) = 8πkxref − x0 k i ωc max{hnpw , n0 (x0 )i, 0} e−i c hnpw ,xi
(4)
|
{z
}
a(x0 )

in the frequency domain, and [15, Eq. (2.178)]p
hn ,xi 
dWFS (x0 , t) = 8πkxref − x0 k h(t) ∗t a(x0 ) δ t − pwc(5)in the time domain, where δ(t) denotes the Dirac delta function and ∗t the convolution in the time domain. The
position of the secondary sources is denoted by x0 = (r0 , π2 , φ0 ), and the inward surface normal vector at x0 by

1133

p
n0 (x0 ) = (1, π2 , φn0 ). The amplitude errors occurring in 2.5D configurations can be corrected by kxref − x0 k
for a reference position xref . Compensating the dimensionality mismatch
between the 2D sound field and the 3Dpsecondary sources (point sources) also requires the equalization i ωc which exhibits a slope of +3 dB/octavepand a constant phase of π4 . In the time domain, this is performed by the filter h(t) = F −1 { i ωc }. The spatial
window a(x0 ) is a rectified cosine function which activates the secondary sources only if hnpw , n0 (x0 )i > 0.
It also applies a cosine weight hnpw , n0 (x0 )i > 0 which comes from the direction gradient of the sound field.
The interested readers are referred to [2, 15] for more detailed treatments.

2.2 Local WFS
For local WFS of a plane wave, a spatial band limitation is applied to the sound field in the cylindrical harmonics domain. In the horizontal plane (θ = π2 ), this reads [16, Sec. 2.6]
M
X

S(x, ω) =−imφpw −mimφω,|e {z } i Jm ( c r)e(6)m=−M S̊ (φ )mpw

where M denotes the spatial bandwidth of the expansion. Notice that the coefficient S̊m (φpw ) = e−imφpw also
describes the 2D plane wave expansion of the sound field [16, Sec. 2.6],Z 2π XMω(7)
S(x, ω) =S̊m (φpw )eimφ̄ e−i c r cos(φ−φ̄) dφ̄,
0

m=−M

|

{z

S̄(φ̄)}

where the sound field is represented as a superposition of plane waves. The direction of the individual planewaves is denoted by n̄ = (1, θ̄, φ̄) and the corresponding amplitude by S̄(φ̄). A spatial band limitation constitutes a rectangular windowing of S̊m (ω) in the m domain, which corresponds to a circular convolution witha periodic sinc function in the φ̄ domain. For a single plane wave, i.e. S̄(φ̄) = δ(φ̄ − φpw ), a spatial band
limitation smears the spatial structure of the sound field.
For a given secondary source position x0 , the spatial windowa(x0 ) = max{hn̄, n0 (x0 )i, 0}(8)

has to be applied to every plane wave in (7). The window function can be approximated as a truncated Fourier
series expansion (equivalent to cylindrical harmonics expansion),a(x0 ) ≈MaXåm (φn0 ) eimφ̄ ,

(9)

m=−Ma

where Ma denotes the modal bandwidth and åm (φn0 ) the Fourier coefficients,
(−1)m/2 −imφn0,em evenZ ππ(1 − m2 )1
a(x0 ) e−imφ̄ dφ̄ = 1 −imφn
åm =
0,

2π −π
e
m = ±1



4
0,
else.

(10)

The spatial windowing is performed in the cylindrical harmonic domain which constitutes a convolution (denotedby ∗m ) of the coefficients,D̊m (φpw , φn0 ) = S̊m (φpw ) ∗m åm (φn0 ) :=µmax
Xµ=µmin1134

S̊µ (φpw ) åm−µ (φn0 ),m = −(M + Ma ), . . . , (M + Ma ), (11)

where µmin = max{−M, −Ma + m} and µmax = min{M, Ma + m}. The coefficient D̊m (φpw , φn0 ) exhibits aspatial bandwidth of Md = M + Ma . Notice that D̊m (φpw , φn0 ) describes the spatially windowed version of
the sound field which corresponds to the last two terms in the right-hand side of (4) or (5). The amplitude
correction and equalization are not applied yet.
At this point, D̊m (φpw , φn0 ) can be used to represent the driving function as a 2D plane wave expansion orequivalently as a cylindrical harmonics expansion. However, as discussed in [13], the analytical time domain
modeling of the radial function Jm ( ωc r) causes temporal aliasing which requires some heuristic adjustments. In
the present work, this is avoided by replacing the cylindrical harmonics expansion with a spherical harmonics
expansion. For 2D sound fields, the coefficient can be converted as [17, Eq. (14)]∗D̆nm (φpw , φn0 ) = 4πYnm( π2 , 0) D̊m (φpw , φn0 ).

(12)Note that, an additional term im−n is present in [17, Eq. (14)], as the definition of the coefficients is slightly
different. The resulting spherical harmonics expansion readsD̃(x0 , ω) =

Md
X

∞
X

∗
4π Ynm( π2 , 0) D̊m (φpw , φn0 ) i−n jn ( ωc r0 ) Ynm ( π2 , φ0 )(13)∗
4π Ynm( π2 , 0) D̊m (φpw , φn0 ) i−n jn ( ωc r0 ) Ynm ( π2 , φ0 ).

(14)m=−Md n=|m|

≈

Md
X

N
X

m=−Md n=|m|While (13) is equivalent to the cylindrical harmonics expansion, (14) is an approximation of the former where
the degree n is truncated to N . The radial part is now described by spherical Bessel functions jn ( ωc r0 ) which
exhibit a strong roll-off at high frequencies, therefore, the time discretization causes significantly less temporalaliasing.
The only frequency dependent term in (14) is jn ( ωc r0 ) whose inverse Fourier transform reads [18, Eq. (10.59.1)]
F −1 {i−n jn ( ωc r0 )} = 2rc0 Pn ( rct0 ),(15)

with Pn (·) denoting the Legendre polynomial. The time domain representation of the spherical harmonics expansion thus reads [19]NMdX∗
˜ 0 , t) = c
d(x
4πYnm
( π2 , 0)F −1 {D̊m (φpw , φn0 )} ∗t Pn ( rct0 ) Ynm ( π2 , φ0 )2r0 nm

(16)

=

NM
c Xd
∗
4πYnm
( π2 , 0)D̊m (φpw , φn0 )Pn ( rct0 ) Ynm ( π2 , φ0 )2r0 nm

(17)

=

NM
c Xd
(2n + 1)Pn ( rct0 )D̊m (φpw , φn0 )Knm eimφ0 ,2r0 nm(18)

where a shorthand notation is used for the double summation,

PN Md
nm

:=

PMd

m=−Md

PN

n=|m| .

The frequency

independence of D̆nm (φpw , φn0 ) is exploited in the second equality, i.e.
F −1 {D̊m (φpw , φn0 )} = D̊m (φpw , φn0 ) · δ(t).
The coefficient Knm in the third equality is defined as [20, p. 783]
( (n−m−1)!!(n+m−1)!!
2,(n−m)!
(n+m)!!(n−m)!!
mm
Kn := (n+m)! Pn (0) =0,n + m even(19)n + m odd,1135

with (·)!! denoting the double factorial. It can be computed efficiently by exploiting K00 = 1 and the recurrence
relations,m−1
Knm = n+m−1
n+m Kn−1 ,m+1
Knm = n−m−1
n−m Kn−1 .

(20)

Notice from (19) that the spherical harmonics expansion of a 2D sound field requires only the coefficients witheven n+m. Finally, the time domain driving function is obtained by multiplying the amplitude correction factor
and applying the equalization filter (cf. (5)),
dLWFS (x0 , t) =p
8πkxref − x0 k h(t) ∗tNM
c Xd(2n + 1)Pn ( rct0 ) D̊m (φpw , φn0 ) Knm eimφ0 .
2r0 nm(21)

Once the target position xref is determined, the coordinate system can be translated so that it coincides with theorigin, i.e. x0 = x − xref . The driving function is then simplified as
s
N
Md
X
0
2π
0
(22)
dLWFS (x0 , t) = ch(t)∗(2n + 1)Pn ( rct0 ) D̊m (φpw , φn0 ) Knm eimφ0 ,
t
0
r00
nm
where x00 = x0 − xref = (r00 , π2 , φ00 ) denotes the secondary source position in the new coordinate system.

3

EVALUATION

The local synthesis of a virtual plane wave is simulated by using the proposed driving function (22). A totalof 60 ideal secondary point sources are uniformly distributed on a circle with radius of 1.5 m (spacing of
15.7 cm). The virtual plane wave propagates in the −y direction, i.e. npw = (1, π2 , − π2 ). Four different target
positions (x, y, z) are chosen, where the spatial aliasing artifacts should be suppressed:‘Center’ : (0, 0, 0),‘Right’ : (0.75, 0, 0),‘Front’ : (0, 0.75, 0),‘Back’ : (0, −0.75, 0).

In each case, the target position is used as the expansion center of the driving function. While the spatial
bandwidth of the sound field is varied M = 15, 23, a fixed bandwidth Ma = 20 is used for the spatial windowa(x0 ). The maximum order of the spherical harmonics expansion is thus Md = M + Ma = 35, 43, and the
maximum degree n is set equal to the maximum order N = Md . The discrete-time driving function is obtainedby sampling (22) at integer multiples of the sampling interval. A max-r E window is applied to the driving
function [21, Eq. (10)]. The simulations are performed at a sampling rate of fs = 88.2 kHz, but the frequency
responses are evaluated only for audio band f ≤ 22 kHz. The pre-equalization filter is realized as a linear phase
FIR filter (513 taps), therefore the constant phase π4 is omitted. The speed of sound is set to c = 343 m/s.
The driving functions are shown in Fig 1, where the equalization filter is not applied to better visualize thetemporal structure. The first thing to notice is the high amplitude pulses along the curve t = rc0 cos(φ0 − φpw )for φ0 ∈ [0, π), which is almost independent to the target positions. This shows that, in the low frequencies, thelocal WFS driving function coincides with the conventional one (5). The fine structure of the driving functions,representing the high frequency components, depends on the target position with varying spatial and temporalcharacteristics. The non-causality of the driving functions is due to the definition of the plane wave which
passes the origin at time t = 0. In practice, a pre-delay has to be added.
The synthesized sound fields on the xy-plane are visualized in Fig. 2. At each target position (indicated by‘+’), a strong impulse is received which corresponds to the wavefront of the virtual plane wave. In the bottomrow, the driving functions are band limited in the temporal frequency domain, as indicated at the top-left cornerof each plot. For a bandwidth of 500 Hz (leftmost), the wavefront is almost perfectly synthesized within theentire area surrounded by the loudspeakers. A phase distortion can be seen which is due to the linear phasedesign of the equalization filter h(t). As the temporal bandwidth increases, the synthesis is correct only within

1136

9

dB
Center

Right

Front

0

Back

t / ms

6
3

−30

0
−60

−3
−6
−9

0.0

0.5

1.0
φ0 / π

1.5

2.0

0.0

0.5

1.0
φ0 / π

1.5

2.0

0.0

0.5

1.0
φ0 / π

1.5

2.0

0.0

0.5

1.0
φ0 / π

1.5

−90

2.0

Figure 1. Proposed local WFS driving functions (φpw = − π2 , M = 23, max-r E window, without equalization).
dB
Center

Right

Front

Back

0

y/m

1

0

−15

−1

< 500 Hz

< 1 kHz< 2 kHz−30< 4 kHz

y/m

1

−45

0

−1
−1

0
x/m

1

−1

0
x/m

1

−1

0
x/m

1

−1

0
x/m

1

−60

Figure 2. Synthesized sound fields (φpw = − π2 , M = 23, max-r E window). Top: full-band sound fields
for different target positions. Bottom: temporally band limited sound fields for target position ‘Right’. The
secondary point sources are indicated by the loudspeaker symbols and the target positions by ‘+’.
a circular area, the radius of which scales down. For the rest of the listening area, the wavefront is smeared andsound waves propagating in wrong directions are superimposed. The former is due to spatial band limitation
whereas the latter is caused by spatial aliasing.
The received time signals at the target positions are shown in Fig. 3a (left). The phase distortion mentioned
above is clearly seen. In Fig. 3a (right), the impulse responses at the same receiver positions are evaluated,but the driving function for the target position ‘Right’ is used. The temporal smearing can be seen at off-target
positions.
In Fig. 3b, the frequency responses at the target positions are depicted by coloured lines. Above 300 Hz,
the magnitude responses are almost flat with moderate fluctuations (< 3 dB). The roll-off in the lower end
is attributed to the high-frequency approximation employed in WFS. Since the low frequency deviations arealmost constant, the overall system response can be equalized, e.g. by a shelving filter. The spectral variations
around the target positions are also evaluated, which might occur due to listener’s movements or head rotations.

1137

Off target

Amplitude

On target
Center/Center

Right/Center

Right/Right

Right/Right

Front/Front

Right/Front

Back/Back−2

Magnitude / dB

0

−1

0
t / ms1

Right/Back
−2

2

(a)

M = 23−1

0
t / ms1

2

M = 15Center

−30

Right

−60

Front

−90

Back
102

103
f / Hz

104

102

(b)

103
f / Hz104

Figure 3. (a) Impulse responses (M = 23) for different target/receiver positions. (b) Frequency responses(M = 23, 15) at on- and off-target positions. The thin gray curves show the frequency responses on a circle
centered at each target position (radius of 8.5 cm, 20 uniform angles).
Twenty equiangular points are chosen on a circle with radius 8.5 cm centred at each target position. The radiuscorresponds to the half wavelength of 2 kHz. The frequency responses at those positions are plotted with thingray curves in Fig. 3b. High frequency (> 5 kHz) attenuation is observed in most cases, which is a typicaleffect of a spatial band limitation. However, in Fig. 3b (left, M = 23), the frequency responses around ‘Front’
exhibit high frequency boost which is caused by spatial aliasing. This agrees with the observation in [22], thatthe spatial aliasing is more pronounced at positions those are closer to the virtual sound source. The influence of
spatial aliasing can be reduced to some extent by decreasing the spatial bandwidth. As shown in 3b (right), thespectral distortion at ‘Front’ is significantly reduced, but the high frequency attenuation at other target positions
becomes more severe. The spatial bandwidth has to be chosen in such a way that the spectral distortion due tothe spatial band limitation and the spatial aliasing artifacts are in an acceptable level within the listening area.
Alternatively, the spatial bandwidth can be adaptively changed according to the tracked listener position. It canbenefit from a geometric prediction model of spatial aliasing artifacts, such as [22].

4

CONCLUSION

In this paper, a time domain driving function for local WFS was introduced which is represented as a sphericalharmonics expansion. The significant reduction of spatial aliasing around the target positions is demonstrated bynumerical simulations The size of the aliasing-free area can be controlled by varying the spatial bandwidth. No
noticeable artifact due to temporal aliasing is found which was problematic in the earlier approach [13]. The
computational efficiency is still maintained, due to the recurrence relation of the coefficients. This makes the
proposed method suited for real time and dynamic scenarios, where the listener movements are tracked and thedriving functions are updated accordingly.

1138

REFERENCES[1] A. J. Berkhout, D. de Vries, and P. Vogel, “Acoustic Control by Wave Field Synthesis,” J. Acoust. Soc. Am. (JASA), vol. 93,
no. 5, pp. 2764–2778, 1993.
[2] S. Spors, R. Rabenstein, and J. Ahrens, “The Theory of Wave Field Synthesis Revisited,” in Proc. 124th Audio Eng. Soc.
(AES) Conv., 2008.
[3] H. Wierstorf, Perceptual Assessment of Sound Field Synthesis. Doctoral thesis, Technical University of Berlin, 2014.
[4] P. Vogel, Application of Wave Field Synthesis in Room Acoustics. PhD thesis, Delft University of Technology, 1995.
[5] E. Verheijen, Sound Reproduction by Wave Field Synthesis. PhD thesis, Delft University of Technology, 1997.
[6] E. W. Start, Direct Sound Enhancement by Wave Field Synthesis. PhD thesis, Delft University of Technology, 1997.
[7] H. Wierstorf, C. Hohnerlein, S. Spors, and A. Raake, “Coloration in Wave Field Synthesis,” in Proc. 55th Audio Eng. Soc.
(AES) Int. Conf. Spatial Audio, Audio Engineering Society, 2014.
[8] H. Wittek, Perceptual Differences between Wavefield Synthesis and Stereophony. PhD thesis, University of Surrey Surrey,
UK, 2007.
[9] S. Spors and J. Ahrens, “Local Sound Field Synthesis by Virtual Secondary Sources,” in Proc. 40th Audio Eng. Soc. (AES)Int. Conf. Spatial Audio, (Tokyo, Japan), Oct. 2010.
[10] S. Spors, K. Helwani, and J. Ahrens, “Local Sound Field Synthesis by Virtual Acoustic Scattering and Time-reversal,” inProc. 131th Audio Eng. Soc. (AES) Conv., (New York, USA), Oct. 2011.
[11] N. Hahn, F. Winter, and S. Spors, “Local Wave Field Synthesis by Spatial Band-limitation in the Circular/Spherical Harmonics Domain,” in Proc. 140th Audio Eng. Soc. (AES) Conv., (Paris, France), June 2016.
[12] F. Winter, N. Hahn, and S. Spors, “Time-Domain Realisation of Model-Based Rendering for 2.5D Local Wave Field Synthesis
Using Spatial Bandwidth-Limitation,” in Proc. Eur. Signal Process. Conf. (EUSIPCO), (Kos, Greece), 2017.
[13] N. Hahn, F. Winter, and S. Spors, “Synthesis of a Spatially Band-Limited Plane Wave in the Time-Domain Using Wave Field
Synthesis,” in Proc. Eur. Signal Process. Conf. (EUSIPCO), (Kos, Greece), 2017.
[14] E. G. Williams, Fourier Acoustics: Sound Radiation and Nearfield Acoustical Holography. Academic press, 1999.
[15] F. Schultz, Sound Field Synthesis for Line Source Array Applications in Large-Scale Sound Reinforcement. Doctoral thesis,
University of Rostock, 2016.
[16] A. Kuntz, Wave Field Analysis using Virtual Circular Microphone Arrays. Verlag Dr. Hut, 2008.
[17] N. Hahn and S. Spors, “Sound Field Synthesis of Virtual Cylindrical Waves using Circular and Spherical Loudspeaker
Arrays,” in Proc. 138th Audio Eng. Soc. (AES) Conv., Audio Engineering Society, 2015.
[18] F. W. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, NIST Handbook of Mathematical Functions Handback. Cambridge
University Press, 2010.
[19] N. Hahn and S. Spors, “Time-Domain Representations of a Plane Wave with Spatial Band-Limitation in the Spherical Harmonics Domain,” in Proc. 45th German Annu. Conf. Acoust. (DAGA), (Rostock, Germany), Mar. 2019.
[20] G. B. Arfken and H. J. Weber, Mathematical Methods for Physicists. Academic press, 2005.
[21] F. Zotter and M. Frank, “All-Round Ambisonic Panning and Decoding,” J. Audio Eng. Soc, vol. 60, no. 10, pp. 807–820,
2012.
[22] F. Winter, F. Schultz, G. Firtha, and S. Spors, “A Geometric Model for Prediction of Spatial Aliasing in 2.5D Sound Field
Synthesis,” IEEE/ACM Trans. Audio, Speech, and Language Process., vol. 27, pp. 1031–1046, June 2019.

1139