© The Author(s), Published by EDP Sciences, 2024

1 Introduction

Nowadays, modern large-scale sound reinforcement systems often use line-source loudspeaker arrays to provide high-quality sound for the largest parts of the audience area [1]. The optimization of these line sources is a widely discussed topic in the literature [2–11].

For surround sound reinforcement, it has been shown that reproduction at off-center listening positions can be improved by line-source loudspeakers [12–18]. In particular, optimal surround sound reinforcement considers a dual-target design of the direct-sound: 0 dB per distance doubling (dod) to preserve the mixing balance and −3 dB/dod to preserve the envelopment at off-center listening positions. Considering a live concert scenario, single-target line-array designs between 0 dB/dod and −1.5 dB/dod for both, direct sound objects and enveloping parts, achieve similar good results as the proposed dual-target design, also at off-center listening positions [18]. However, the line-array sound pressure design is limited to a line in the audience that lies in the axial direction of the loudspeaker arrangement. By contrast, planar arrays allow the sound pressure to be controlled across the depth and the width of the audience area. Therefore planar arrays should be investigated to find out whether further improvement can be achieved for surround sound reinforcement at off-center listening positions.

Planar arrays are known from antenna theory and telecommunication applications, in which a beam is electronically formed by applying amplitude and/or phase weights [19–23]. Antennas with co-secant squared patterns are used to track moving targets at constant height, such as airplanes, where the received power should stay independent of the position of the object [24–26]. This design approach is similar to that of loudspeaker arrays that aim for 0 dB/dod.

In acoustics, loudspeaker arrays are often used to reconstruct acoustic wave fields with desired wavefront properties [27, 28]. To control the direct sound, Wave Front Synthesis [29, 30] has been introduced that applies amplitude and phase load (driving function) to linear and planar loudspeaker arrays based on the Kirchhoff-Helmholtz and Rayleigh integrals. While the Rayleigh integral applies as an exact formulation to the planar point-source arrays of three-dimensional Wavefield Synthesis (3D WFS), the stationary-phase approximation of the Kirchhoff-Helmholz integral is required to obtain the driving function in the case of 2.5D WFS and its linear loudspeaker arrangements [31–36]. Due to practical restrictions, most applications have employed linear loudspeaker arrays as opposed to planar loudspeaker arrays, so far. The WFS theory has been shown to be applicable, more generally, to curved linear arrays in sound reinforcement to design their optimal curvature [37].

Recently, Holoplot¹ has released planar loudspeaker arrays, also referred to as matrix arrays, as commercial products for professional sound reinforcement applications, e.g. targeting 3D audio at live events [38, 39]. In particular, this application relies on FIR filters, gains, and delays to synthesize a virtual source emitting an optimal wavefront.

Although the available array literature in sound reinforcement outlines and summarizes the basic principles of planar arrays, it appears to lack a detailed definition of how to design the delay progression over an ideal continuous planar source in order to produce an adjustable sound pressure level decay across the audience area.

Based on the ideas that have been developed for continuously curved and phased line sources [11], recently, this paper presents a design equation for the continuous delay length profile of a phased planar source. The goal is to obtain an equalized direct-sound pressure level decay of −6 ⋅ β dB/dod. β is a design parameter whose value ranges typically from 0 and 0.5 yielding sources with distance decays between 0 dB/dod and −3 dB/dod. This essentially yields a (near-field) delay-and-sum beamformer design. Similar to a design based on linear arrays [11], broadband targets can be met when the planar array is tall and wide enough, and its drivers are spaced sufficiently close to each other.

The paper is structured as follows: Section 2.2 utilizes the stationary-phase approximation to get an estimate of the sound pressure level from the two-dimensional integral of the delayed Green’s function over a plane. Constrained to target the roll-off r^−β, its result defines a second-order non-linear differential equation describing a continuous phasing across the plane in Section 2.4 that is expressed as a delay length and solved in Algorithms 1–3. Section 3.1 shows simulations of this continuously phased planar source as a proof of concept, and Section 3.3 presents simulations of a discrete source geometry with discrete phase load to discuss the effects of discretization and spatial aliasing. Section 4 presents the applicability to practical applications by comparing the results of the simulations with measurements using on a small prototype array.

2 Continuous phasing

Our objective is to evaluate the sound pressure p of a planar source fed by a progressive time delay τ given as length w = c τ. The integral of a delayed Green’s function $G\left(r\right)=\frac{e^{-\mathrm{i} k(r+w)}}{4 \pi r}$ over the horizontal source coordinate h and the vertical source coordinate v describes the resulting sound pressure,

$$p=\iint \mathrm{}\frac{e^{-i k \left(r+w\right)}}{4 \pi r} \mathrm{d}h \mathrm{d}v.$$(1)

The imaginary unit yields i² = −1, the wave number is $k=\frac{2 \pi f}{c}$, the speed of sound is c = 343 m/s, and f denotes the frequency.

2.1 Source geometry

The geometry of the planar source x_s is described by

$${\mathit{x}}_{\mathrm{s}}={\mathit{x}}_0+h{\mathit{e}}_{\mathrm{h}}+v{\mathit{e}}_{\mathrm{v}},$$(2)

where e_h = [0 1 0]^⊺ is a unit-vector in horizontal direction and the source is symmetrical with regard to the x-z plane, $h\in [-\frac{H}{2},\frac{H}{2}]$. v defines the source length V in vertical direction, i.e. v ∈ [0, −V] and η in e_v = [sin η 0 cos η]^⊺ is the vertical inclination of the source, cf. Figures 1 and 2. x₀ = [0 0 z₀]^⊺ is a point on the source that defines the height of the source, i.e. the vertical distance z₀ from the highest point of the source to the audience area that is assumed to lie on x-y plane, i.e. z_r = 0. The distance to the receiver at x_r = [x_r y_r 0]^⊺ is

$$r=\Vert {\mathit{x}}_{\mathrm{s}}-{\mathit{x}}_{\mathrm{r}}\Vert =\sqrt{({\mathit{x}}_{\mathrm{s}}-{\mathit{x}}_{\mathrm{r}}{)}^{⊺}({\mathit{x}}_{\mathrm{s}}-{\mathit{x}}_{\mathrm{r}})}.$$(3)

Figure 1 Geometry of a planar source (gray rectangle) with horizontal extent H and vertical extent V showing the unit vectors eh in horizontal direction and ev in vertical direction.

Figure 2 Side view of a planar source (solid gray) inclined by η showing the stationary-phase direction u, cf. Section 2.3 and the beamforming angles ${\vartheta }_0=-\arcsin \frac{\partial w}{\partial v}{|}_{v =\mathrm{ }0, h = 0}$ on top and ${\vartheta }_V=-\arcsin \frac{\partial w}{\partial v}{|}_{v = -V, h = 0}$ on the bottom of the source for a desired depth coverage between xr0 and xrV, cf. Section 2.3.

Figure 3 Discretization of geometry of a phased planar source.

Figure 4 Side view of simulated sound pressure level maps at y = 0 for phased planar arrays in 0 dB/dod configuration comparing the continuous source to a source with discretized geometry for f = 1 kHz and f = 6 kHz.

Figure 5 Directivity index of phased planar sources with b = 0 (solid lines) and b = 0.2531 (dashed lines) for different decay parameters β.

Figure 6 Delay lengths of a continuous phased planar source mounted straight on the wall (left, η = 0°) compared to the delay lengths of a continuous phase planar source inclined by η = 6.5° with β = 0 and b = 0.3324.

Figure 7 Experimental setup using line array enclosures of [42] mounted as 7 × 8 element planar array.

2.2 Stationary-phase approximation

At high frequencies, the integrand of equation (1) oscillates rapidly. Therefore, the stationary-phase approximation is appropriate. Typically, designs should have a single stationary-phase point on the source plane which is characterized by a minimum of delay $\tau =\frac{w}{c}$ plus time of flight to the receiver. Equivalently, the first-order derivatives of the corresponding distances r + w must vanish, $\frac{\partial r}{\partial h}+\frac{\partial w}{\partial h}=0$ and $\frac{\partial r}{\partial v}+\frac{\partial w}{\partial v}=0$. Evaluation at this stationary-phase source point defines the approximate sound pressure

$$p\approx \frac{e^{-\mathrm{i} k r-\mathrm{i} \frac{\pi }{4} \mathrm{sig}\left\{\mathit{H}\right(r+w\left)\right\}}}{2 k}\frac{1}{\sqrt{r^2\left|\mathit{H}\right(r+w\left)\right|}},$$(4)

where sig{H(·)} denotes the signature of the Hessian matrix and |H(·)| its determinant,

$$\left|\mathit{H}\left(r+w\right)\right|=\left|\begin{array}{ll}\frac{{\partial }^{2}r}{\partial h^2}+\frac{{\partial }^{2}w}{\partial h^2}& \frac{{\partial }^{2}r}{\partial h\partial v}+\frac{{\partial }^{2}w}{\partial h\partial v}\\ \frac{{\partial }^{2}r}{\partial h\partial v}+\frac{{\partial }^{2}w}{\partial h\partial v}& \frac{{\partial }^{2}r}{\partial v^2}+\frac{{\partial }^{2}w}{\partial v^2}\end{array}\right|=\left(\frac{{\partial }^{2}r}{\partial h^2}+\frac{{\partial }^{2}w}{\partial h^2}\right)\left(\frac{{\partial }^{2}r}{\partial v^2}+\frac{{\partial }^{2}w}{\partial v^2}\right)-$$(5)

$$-{\left(\frac{{\partial }^{2}r}{\partial h\partial v}+\frac{{\partial }^{2}w}{\partial h\partial v}\right)}^2.$$(6)

The derivatives of the distance are

$$\begin{array}{cc}\frac{\partial r}{\partial h}=\frac{{\mathit{e}}_{\mathrm{h}}^{⊺}({\mathit{x}}_{\mathrm{s}}-{\mathit{x}}_{\mathrm{r}})}{r}& \frac{\partial r}{\partial v}=\frac{{\mathit{e}}_{\mathrm{v}}^{⊺} ({\mathit{x}}_{\mathrm{s}}-{\mathit{x}}_{\mathrm{r}})}{r},\end{array}$$(7)

$$\begin{array}{cc}\frac{{\partial }^{2}r}{\partial h^2}=\frac{1-{\left(\frac{\partial r}{\partial h}\right)}^2}{r}& \frac{{\partial }^{2}r}{\partial v^2}=\frac{1-{\left(\frac{\partial r}{\partial v}\right)}^2}{r},\end{array}$$(8)

$$\frac{{\partial }^{2}r}{\partial h\partial v}=-\frac{{\mathit{e}}_{\mathrm{h}}^{⊺}({\mathit{x}}_{\mathrm{s}}-{\mathit{x}}_{\mathrm{r}})\frac{\partial r}{\partial v}}{r^2}=-\frac{\frac{\partial r}{\partial h}\frac{\partial r}{\partial v}}{r}.$$(9)

2.3 Stationary-phase direction and distance

The direction $\mathit{u}=-\frac{{\mathit{x}}_{\mathrm{s}}-x_{\mathrm{r}}}{r}={\left[\begin{array}{ccc}{u}_x& u_y& u_z\end{array}\right]}^{⊺}$ points from a position (h, v) on the source to a receiver at z_r = 0. With the local coordinates (Eq. (2)) and derivatives (Eq. (7)), we can define its entries using $\frac{\partial r}{\partial h}=-{\mathit{e}}_{\mathrm{h}}^{⊺}\mathit{u}=-u_y$ and $\frac{\partial r}{\partial v}=-{\mathit{e}}_{\mathrm{v}}^{⊺}\mathit{u}=-\sin \eta u_x-\cos \eta u_z$. Furthermore, u_z, u_x follow from the quadratic equation $u_x^2+u_y^2+u_z^2=1$ as solved in Appendix A,

$$\mathit{u}=\left[\begin{array}{l}\sqrt{1-u_y^2-u_z^2}\\ -\frac{\partial r}{\partial h}\\ -\frac{\partial r}{\partial v} \cos \eta -\sin \eta \sqrt{1-(\frac{\partial r}{\partial h}{)}^2-(\frac{\partial r}{\partial v}{)}^2}\end{array}\right].$$(10)

The stationary-phase conditions $\frac{\partial r}{\partial h}+\frac{\partial w}{\partial h}=0$ and $\frac{\partial r}{\partial v}+\frac{\partial w}{\partial v}=0$ yield a stationary-phase direction from (h, v)

$$\mathit{u}=\left[\begin{array}{l}\sqrt{1-u_y^2-u_z^2}\\ \frac{\partial w}{\partial h}\\ \frac{\partial w}{\partial v}\cos \eta -\sin \eta \sqrt{1-(\frac{\partial w}{\partial h}{)}^2-(\frac{\partial w}{\partial v}{)}^2}\end{array}\right]$$(11)

targeting a receiver at z_r = 0 or [0 0 1][x_s + r u]=0. Therefore, the stationary-phase distance is

$$r=\frac{v\cos \eta +z_0}{-\frac{\partial w}{\partial v}\cos \eta +\sin \eta \sqrt{1-(\frac{\partial w}{\partial h}{)}^2-(\frac{\partial w}{\partial v}{)}^2}}.$$(12)

2.4 Target coverage

Our goal is to find the delay lengths w = c τ for which the sound pressure

$$p|\propto r^{-\beta }$$(13)

decays by −6 ⋅ β dB per doubling to the distance on the audience area that lies at z_r = 0. In typical designs, values between 0 and 0.5 are used for the design parameter β. We introduce a gain g and desire for the equalized magnitude square,

$$4{ k}^2 |p{|}^2=\frac{1}{r^2\left|\mathit{H}\right(r+w\left)\right|}={\left(\frac{g}{r^{\beta }}\right)}^2$$(14)

with $r=\frac{v\cos \eta +z_0}{-\frac{\partial w}{\partial v}\cos \eta +\sin \eta \sqrt{1-(\frac{\partial w}{\partial h}{)}^2-(\frac{\partial w}{\partial v}{)}^2}}$.

With stationary-phase conditions inserted in equation (6), this yields a design equation for the delay length w,

$$\begin{array}{c}{g}^{-2}{r}^{2\beta }=\left[1-{\left(\frac{\partial w}{\partial h}\right)}^2+r\frac{{\partial }^{2}w}{\partial h^2}\right]\\ \cdot \left[1-{\left(\frac{\partial w}{\partial v}\right)}^2+r\frac{{\partial }^{2}w}{\partial v^2}\right]\\ -{\left[-\frac{\partial w}{\partial h}\frac{\partial w}{\partial v}+r\frac{{\partial }^{2}w}{\partial h\partial v}\right]}^2.\end{array}$$(15)

For a given stationary-phase direction and distance specified by $\frac{\partial w}{\partial h},\frac{\partial w}{\partial v},r$, this differential equation is a single condition for three curvature parameters $\frac{{\partial }^{2}w}{\partial h^2},\frac{{\partial }^{2}w}{\partial v^2},\frac{{\partial }^{2}w}{\partial h\partial v}$. There are many solutions to this underdetermined design equation.

However, it is reasonable to introduce some manual separation and control. The first substitution considers the control of the horizontal coverage. In typical practical applications, w and $\frac{\partial w}{\partial h}$ should be 0 for h = 0. At v = 0, w must increase to the sides of the source |h| > 0 for a convex wave front covering a horizontal target range. Replacing $r\frac{{\partial }^{2}w}{\partial h^2}$ with ${\left(\frac{\partial w}{\partial h}\right)}^2+r a$ removes $\frac{{\partial }^{2}w}{\partial h^2}$ from the first bracket of equation (15) and defines a separate equation

$$\frac{{\partial }^{2}w}{\partial h^2}=a+\frac{1}{r}{\left(\frac{\partial w}{\partial h}\right)}^2,$$(16)

which yields a parabola $w\approx \frac{a^{2}h}{2}$ for $(\frac{\partial w}{\partial h}{)}^2\approx 0$. A positive a is used to control the horizontal coverage.

We consider a second substitution for the mixed derivative $\frac{{\partial }^{2}w}{\partial h \partial v}$. It controls the amount by which the vertical inclination changes up or down or with the lateral coordinate h. Replacing $r \frac{{\partial }^{2}w}{\partial h \partial v}$ with $\frac{\partial w}{\partial h} \frac{\partial w}{\partial v}+r b$ removes $\frac{{\partial }^{2}w}{\partial h \partial v}$ from the third bracket of equation (15) and defines the equation

$$\frac{{\partial }^{2}w}{\partial h\partial v}=b+\frac{1}{r}\frac{\partial w}{\partial h}\frac{\partial w}{\partial v}.$$(17)

Hereby, the off-diagonals of the Hessian are controlled by the parameter b, which influences the twist of the delay profile in the h-v plane. If b = 0, the mixed derivative $\frac{{\partial }^{2}r}{\partial h \partial v}$ of equation (9) and $\frac{{\partial }^{2}w}{\partial h \partial v}$ are equal, resulting in observation points for v = 0 that build a conical section on the listening area that is close to a circle. For b > 0 this arc curves faster towards x = 0. For b < 0 it curves less towards x = 0.

In the design equation (15), these substitutions turn the first bracket into 1 – r a and the last one into r²b². In addition to equations (16) and (17) a differential equation for the vertical curvature remains,

$$\frac{{\partial }^{2}w}{\partial v^2}=\frac{1}{r}\left[\frac{g^{-2}{r}^{2\beta }+r^2{b}^2}{1+\mathrm{ra}}+{\left(\frac{\partial w}{\partial v}\right)}^2-1\right]$$(18)

with $r=\frac{v\cos \eta +z_0}{-\frac{\partial w}{\partial v}\cos \eta +\sin \eta \sqrt{1-{\left(\frac{\partial w}{\partial h}\right)}^2-(\frac{\partial w}{\partial v}{)}^2}}$.

2.5 Optimal tilting and depth range

Typical applications do not apply horizontal beamforming at h = 0 which implies that the first order derivative of w with regard to h must vanish, $\frac{\partial w}{\partial h}{|}_{h =\mathrm{ }0}=0$. The parameters of the horizontal progression are initially defined at this point, so we also specify the initial values for the vertical progression of the delay length here.

By choosing the farthest observation point x_r0, i.e. the last seat row, and the height of the source measured from its highest point downwards to the listening area at z_r0, the tilt angle on top of the source is defined by

$$\nu =\arctan \left(\frac{z_{\mathrm{r}0}}{x_{\mathrm{r}0}}\right).$$(19)

This tilt angle is composed linearly by the inclination of the source η and a beamforming angle ϑ₀ on top of the source, ν = η + ϑ₀. The beamforming angle ϑ₀ is related to the tangent of the delay length on top of the source via the sine that defines the initial condition for the first derivative of w with regard to v, cf. Figure 2,

$$\frac{\partial w}{\partial v}{|}_{h =\mathrm{ }0, \mathrm{v }=\mathrm{ }0}=-\sin (\nu -\eta ).$$(20)

2.6 Start gain parameter

Typical practical applications require a convex progression of the delay length. The zero point of curvature at the highest point of the source determines the maximum value of the gain parameter for a convex delay length profile,

$$g_{\max }^2=\frac{r^{2\beta }}{\left[1-{\left(\frac{\partial w}{\partial v}\right)}^2\right](1+r a)-r^2{b}^2}.$$(21)

g defines how straight the starting curvature $\frac{{\partial }^{2}w}{\partial v^2}$ should be. A large parameter yields small curvature and requires a longer curve length to hit the closest observation point x_rV. To find this parameter automatically for a given source length V, a fixed-point iteration has already been proposed in [40] to find the gain parameter g for a line source, fitting the curvature between a given nearest x_rS and farthest observation point x_r0. In our case, we solve the differential equation (18) for h = 0 and calculate the x-coordinate of the closest observed point in the audience ${\widehat{x}}_{\mathrm{rS}}$ using the stationary-phase direction x_r = x_s + r u evaluated at h = 0 and v = −V,

$${\widehat{x}}_{\mathrm{rS}}=x_{\mathrm{s}}+r\sqrt{1-u_z^2}{|}_{h = 0,v = -V}.$$(22)

x_s is the x coordinate of the source and u_z denotes the z coordinate of the stationary-phase direction u. The algorithm starts with $g^2=0.99\mathrm{ }{g}_{\max }^2$ and compares the resulting depth coverage $\widehat{X}=x_{\mathrm{r}0}-{\widehat{x}}_{\mathrm{rS}}$ to the desired depth coverage X = x_r0 − x_rS, yielding an exponent parameter ${\alpha }_1=\frac{X}{\widehat{X}}$ that is used to update the parameter g,

$$g^2\leftarrow g_{\max }^2{\left(\frac{g^2}{g_{\max }^2}\right)}^{{\alpha }_1}.$$(23)

Whenever the exponent parameter α₁ is smaller than unity, it increases the gain parameter g resulting in a smaller vertical curvature which happens for an audience range $\widehat{X}$ that is too long. α₁ becomes large, if the resulting closest observation point is farther away. Therefore g² must be decreased. Furthermore, we introduce κ₁ to accelerate convergence in ${\alpha }_1={\left(\frac{X}{\widehat{X}}\right)}^{{\kappa }_1}$. We found that κ₁ = 3 is sufficiently fast to meet X within ±1% in a little number of iterations.

2.7 Horizontal coverage and concaveness parameter

The same fixed-point iteration is applied to find optimal values for the parameters a and b for a desired horizontal coverage, defined by an outermost point with x-coordinate X and y-coordinate Y. Its corresponding stationary-phase point on the source is at $h=\frac{H}{2}$ and v = 0. Typical designs require $\frac{\partial w}{\partial h}=0$ at h = 0 defining a positive non-zero value a > 0 in $\frac{{\partial }^{2}w}{\partial h^2}=\frac{1}{r}\frac{\partial w}{\partial h}+a$. As the goal is to reach an appropriate horizontal curvature, a_min is set to 0.3. The algorithm starts with a = 1.2 · a_min = 0.36 and compares the resulting width $\tilde{Y}$,

$$\tilde{Y}=y_{\mathrm{s}}+r\frac{\partial w}{\partial h}{|}_{h = \frac{H}{2},v = 0},$$(24)

to the desired width Y, where y_s is the y coordinate of the source. The comparison of the desired width to the observed width yields the parameter ${\alpha }_2={\left(\frac{Y}{\tilde{Y}}\right)}^{{\kappa }_2}$ that is used to update the parameter a,

$$a\leftarrow a_{\min }{\left(\frac{a}{a_{\min }}\right)}^{{\alpha }_2}.$$(25)

The same procedure is applied to update the parameter b. The observation points build a conical section for b = 0. This section gets warped when choosing other values for b or modifying the outermost observation point for $h=\frac{H}{2}$ and v = 0. For the fixed-point algorithm, we chose a minimum value of b_min = 0.01 and an initial value of b = 1.2 · b_min = 0.012. The corresponding x coordinate of the outermost observation point is

$$\tilde{X}=r\sqrt{1-u_y^2+u_z^2}{|}_{h = \frac{H}{2},v = 0},$$(26)

and the comparison to the desired x coordinate yields ${\alpha }_3={\left(\frac{\tilde{X}}{X}\right)}^{{\kappa }_3}$ to update the parameter b,

$$b\leftarrow b_{\min }{\left(\frac{b}{b_{\min }}\right)}^{{\alpha }_3}.$$(27)

Furthermore, we use the constants κ₂ = 0.2 and κ₃ = 0.2 to reach convergence in an acceptable number of iterations meeting X and Y within ±1%.

3 Simulation study

Similarly to the theory of line sources [11], we apply step-wise updates to solve the differential equation (18) numerically as outlined in Algorithm 3. As we assume symmetry, the procedure is applied for $h\in [0,\frac{H}{2}]$ and the solution is then mirrored. The step sizes Δv and Δh should be appropriately small to solve the differential equation accurately by choosing a large N. Algorithm 1 finds optimal values for the parameters a and b, defining the horizontal coverage, while Algorithm 2 is then used to get an appropriate value for the gain parameter g that fits the desired depth of the audience area.

For our simulation study we assume a listening area of 10 m × 10 m with the closest observation point at x_rV = 1.5 m. The simulations show the result of summing ideal point sources described by the Green’s function, that are positioned at the continuous source in 1mm intervals. This summation is performed for frequencies ranging from 50 Hz to 20 kHz with a logarithmic third-octave resolution of 27 points. To show broadband sound pressure maps over the listening area, the results are A-weighted [41] and summarized. Since the radiation is symmetrical, the maps are divided into two parts to show the influence of different values of b in a single map.

3.1 Continuously phased planar sources

The differential equation (18) was solved with uniform parameters z₀ = 1.35 m, inclination angle η = 0° and vertical source length V = 0.656 m which yields the number elements N = 657 for a step size Δv = 1 mm. The vertical source length was chosen to fit the length of eight mini line-array elements [42] which are used later in the measurements. For H, we choose 0.572 m so that the total source width corresponds to the length of seven mini line-array elements. The farthest point of observation in y direction is chosen to be at y_r = 5 m while in x direction it is x_r = 10 m. To show the influence of the parameter b, simulations are performed using b = 0 (digonalized Hessian matrix) on the one hand and on the other hand, by using an optimized value of b to reach an outermost observation point at x_r = [3.5 5 0]^⊺. Applying Algorithm 1 yields b = 0.2531.

Figures 8–10 show the delay length profiles w = c τ using different decay parameters β. Comparing the delay contours for different decay parameters, it is noticable that higher values of the decay parameter have rounder contours as the outcome of the source has to become more curved. Therefore, higher values of β also require higher delays.

Figure 8 Delay lengths for a continuously phased planar source (η = 0°), b = 0 (left) and b = 0.2531 (right) for β = 0 (0 dB/dod).

Figure 9 Delay lengths for a continuously phased planar source (η = 0°), b = 0 (left) and b = 0.2531 (right) for β = 0.25 (−1.5 dB/dod).

Figure 10 Delay lengths for a continuously phased planar source (η = 0°), b = 0 (left) and b = 0.2531 (right) for β = 0.5 (−3 dB/dod).

Figures 11–13 show the equalized A-weighted sound pressure maps of the continuously phased planar sources with different decay parameters β. For equalization, we use a filter whoose magnitude rises by 6 dB/octave to compensate the red frequency response of a planar source composed by ideal point sources, cf. Appendix B. As the integral of equation (1) and the stationary-phase approximation assume infinite integral bounds, the sound pressure level rolls off by −6 dB/dod for points of observation which lie outside the desired coverage area.

Figure 11 Map and contours of the simulated A-weighted sum of a continuous phased planar source (η = 0°), b = 0 (left) and b = 0.2531 (right) for β = 0 (0 dB/dod).

Figure 12 Map and contours of the simulated A-weighted sum of a continuous phased planar source (η = 0°), b = 0 (left) and b = 0.2531 (right) for β = 0.25 (−1.5 dB/dod).

Figure 13 Map and contours of the simulated A-weighted sum of a continuous phased planar source (η = 0°), b = 0 (left) and b = 0.2531 (right) for β = 0.5 (−3 dB/dod).

Considering the effect of the parameter b, higher values results in a wider coverage. As mentioned above, the observation points for b = 0 build a conical section on the listening area. This section gets warped for higher values of b, which can also be seen in the contours of the maps. Furthermore, higher values of b increase the amount of delay that is necessary.

Ideal phased line sources exhibit a rotationally symmetrical radiation pattern on the listening area without horizontal directivity if they are mounted straight at the wall (η = 0°). Tilting these optimum phased line sources causes the radiation to become horizontally directional, limiting the horizontal coverage. By contrast, the horizontal range for planar arrays can be set using the parameter a, even for planar arrays mounted on the wall, η = 0°, cf. Figures 14–16. Unlike phased line sources, the horizontal range is not a direct consequence of the source inclination; rather, it may be adjusted independently of it. In practice, line sources are employed with double bass arrangements and waveguides for higher frequencies to control the horizontal directivity.

Figure 14 Map and contours of the simulated A-weighted sum of a discrete phased planar source (η = 0°) with continuous delay profile, b = 0 (left) and b = 0.2531 (right) for β = 0 (0 dB/dod).

Figure 15 Map and contours of the simulated A-weighted sum of a discrete phased planar source (η = 0°) with continuous delay profile, b = 0 (left) and b = 0.2531 (right) for β = 0.25 (−1.5 dB/dod).

Figure 16 Map and contours of the simulated A-weighted sum of a discrete phased planar source (η = 0°) with continuous delay profile, b = 0 (left) and b = 0.2531 (right) for β = 0.5 (−3 dB/dod).

3.2 Discretization of geometry

Later, we use the miniature line array [42] with 8.2 cm × 10 cm × 10 cm 3D-printed enclosures and 2.5 inch broadband drivers to perform a sample measurement. Therefore, we choose the dimensions of the drivers for discretizing the geometry of the planar source. We assume that the planar source is composed by an array of 7 × 8 elements and simulate a single driver by a circular area with 3 cm radius leaving gaps in between the elements, cf. Figure 3. To comply with our simulations done in [11, 42, 43], the radius of these discrete elements gets shorter by the fifth power of the simulated frequency. For 50 Hz the radius is 3 cm, for 20 kHz it is 1.8 cm.

Figures 17–19 show the equalized A-weighted sound pressure maps over the desired area. For equalization, we use the filter that compensates the red frequency response, cf. Appendix B, but keep the magnitude constant above the spatial aliasing frequency for which endfire-radiation occurs whenever λ < 8.2 cm, i.e. f_spat.al. = 4.2 kHz. In comparison to the results of the continuous source, the main difference lies in the vicinity of the source where spatial aliasing has a high influence on the sound pressure. We denote that the contours are little narrower compared to those of the continuous source. Due to spatial aliasing, higher sound pressures are also to be expected for observation points which lie outside the desired coverage area. Figure 4 compares the radiation of a planar source with discretized geometry in 0 dB/dod configuration for 1 kHz and 6 kHz. Considering the maps for 1 kHz, it is noticeable that also for the continuous source grating lobes are observed. However, the maps of the continuous source and the discretized source are comparable. For 6 kHz, where spatial aliasing occurs for the discrete source, it is seen that the grating lobes, that lie outside the desired coverage area, are more pronounced and therefore higher sound pressure levels for close observation points x_r < 0.5 m are denoted.

Figure 17 Map and contours of the simulated A-weighted sum of a discrete phased planar source (η = 0°) with discretized delay profile, b = 0 (left) and b = 0.2531 (right) for β = 0 (0 dB/dod).

Figure 18 Map and contours of the simulated A-weighted sum of a discrete phased planar source (η = 0°) with discretized delay profile, b = 0 (left) and b = 0.2531 (right) for β = 0.25 (−1.5 dB/dod).

Figure 19 Map and contours of the simulated A-weighted sum of a discrete phased planar source (η = 0°) with discretized delay profile, b = 0 (left) and b = 0.2531 (right) for β = 0.5 (−3 dB/dod).

3.3 Discretization of phase

Both, the phase load and the source geometry are typically discrete in practical applications. Just like a smoothing filter in image processing [44], we apply a Gaussian filter as discretization kernel, to get an appropriate delay length w for a single element,

$$K[m,n]=e^{-\frac{m^2+n^2}{2 {\sigma }^2}}.$$(28)

σ² denotes the standard deviation and determines the width of the Gaussian kernel. An average value for the delay length that is employed as a discrete value is obtained by summing the continuous delay length values w_d belonging to a discrete loudspeaker element, weighted by the discretization kernel H[h, v],

$${\tilde{w}}_{\mathrm{d}}=\mathrm{rd}\left\{\frac{\sum _m\mathrm{}\sum _n\mathrm{}{w}_d[m,n] K[m,n]}{\sum _m\mathrm{}\sum _n\mathrm{} K[m,n]}\frac{f_{\mathrm{s}}}{c}\right\} \frac{c}{f_{\mathrm{s}}},$$(29)

whereby c = 343 m/s denotes the speed of sound. For our simulations, the discrete delay length is rounded to be realized with integer samples at f_s = 48 kHz. Normalization is necessary to prevent the delay length from being distorted whenever ∑_m ∑_n K[m, n] ≠ 1. In our case the diameter of a discrete element is 0.082 m with a resolution of 1mm which requires a discretization kernel of size 82 × 82. σ = 0.025 is found to be a useful value that lets the kernel fade out sufficiently towards the outside of a discrete element.

Figures 17–19 show the A-weighted measured level maps of the planar source with discrete geometry and phase load. Compared to the simulations before, the results for b = 0 differ in that lower sound pressure values are denoted for more distant and more outlying observation points, as well as that the contours are narrower. The results for the optimized value of b indicate that the curvature of the delay length cannot be achieved with the coarse discretization. The continuous delay lengths of the outermost element on the bottom of the source lie in a range between 0.156 m and 0.24 m for β = 0. The discrete delay length, 0.2 m, represents the curvature in this area very poorly, resulting in highly pronounced discretization effects. These effects could be minimized by using smaller drivers and smaller driver spacing correspondingly which would require a two-way system in practice. Figures 20–22 show the results of a simulation using half driver spacing, i.e. 0.041 m. These are comparable to the results of the simulation with discrete geometry only. For the majority of the listening area, only minor differences to the simulation with discrete geometry can be noted.

Figure 20 Map and contours of the simulated A-weighted sum of a discrete phased planar source (η = 0°) with discretized delay profile, β = 0 (0 dB/dod), b = 0 (left) and b = 0.2531 (right) using a smaller discretization.

Figure 21 Map and contours of the simulated A-weighted sum of a discrete phased planar source (η = 0°) with discretized delay profile, β = 0.25 (−1.5 dB/dod), b = 0 (left) and b = 0.2531 (right) using a smaller discretization.

Figure 22 Map and contours of the simulated A-weighted sum of a discrete phased planar source (η = 0°) with discretized delay profile, β = 0.5 (−3 dB/dod), b = 0 (left) and b = 0.2531 (right) using a smaller discretization.

A simpler solution is achieved by tilting the source which also results in smaller delay lengths to be realised, as performed later for the sample measurement. Furthermore, a tilt is preferable in practice because better coverage can be expected in high frequencies if the source is tilted towards the audience area, due to the directivity of the drivers.

3.4 Directivity

Another parameter to be considered beside the coverage, is the frequency-dependent directivity ratio Q(f) and the frequency-dependent directivity index DI(f) = 10 lgQ(f),

$$Q\left(f\right)=\frac{\max \left\{\right|p(f,\theta ,\varphi ){|}^2\}}{\frac{1}{4 \pi }{\int }_0^{2\pi }\mathrm{}{\int }_0^{\pi }\mathrm{}|p\left(f,\theta ,\varphi \right){|}^2\sin \theta \mathrm{d}\theta \mathrm{d}\varphi }.$$(30)

max{p(f, θ, ϕ)²} denotes the maximum squared sound pressure for the observed frequency evaluated at azimuth ϕ and elevation θ. Figure 5 shows the directivity index of phased planar sources for both values of b. As before, the simulation is based on summing the sound pressure of ideal point sources placed in 1 mm intervals on the continuous planar source. To comply with previous work [11], the simulation was done in 10 m distance to the source. It is noticable that for frequencies below 500 Hz the radiation of the source is independent of any settings of β and b. For frequencies below 250 Hz the radiation gets omnidirectional as the dimensions are too small. Above 500 Hz, the curves for different values of b run differently. Higher values of b result in lower values of the directivity index as the beam gets broader to cover a larger area.

4 Sample measurement

4.1 Experimental setup

To verify the theory in practice, impulse response measurements of a planar array consisting of 7 × 8 elements were taken using the miniature line-array enclosures presented in [42], cf. Figure 7. By contrast to the simulations before, the array elements cannot be placed as dense horizontally as they are vertically because the enclosures are wider and there are additional screws on the sides that are used for adjusting the line-array splay angles. In order to reduce the amount of required delay, an inclination of η = 6.5° was applied. Fiugre 6 shows the comparison of the continous delay lengths between a source that is mounted flat on the wall and a source that is inclined by η = 6.5°. The inclination significantly decreases the delays in vertical direction. At the lower end of the source, this results in a 6.85 cm shorter delay length, which corresponds to approximately 10 samples at a sampling rate of f_s = 48 kHz.

The measurements were conducted at the IEM CUBE, a 10.2 m × 12 m studio with T₃₀ ≈ 0.5s. To focus on the direct sound behavior, the impulse responses were truncated to the first 300 samples (6.25 ms at f_s = 48 kHz).

We recorded the responses of a 2 s sweep over a grid of 9 m × 6 m with a resolution of 0.5 m in order to capture the position-dependent direct sound radiation of the planar array. Each element was measured individually to allow the possibility of adding delays later in the analysis. To avoid floor reflections influencing the measurements, we used pressure zone microphones, AKG PZM 30D, which were positioned on the floor. For each element, we applied the filter proposed in [42], in order to make the frequency response of a single element flat. Furthermore, we used the filter that compensates the red frequency response of a planar source composed by point sources until the frequency for which endfire radiation occurs, f_spat.al. = 3.1 kHz. The results are averaged in third octaves and the curves are plotted as A-weighted sum.

4.2 Results

Figures 23–25 show the measured equalized A-weighted sound pressure level maps for the phased planar sources with different decay parameters β. As in the simulation, the results for b = 0 are shown on the left while the results for the optimized value of b to reach the outermost observation point at x_r = [3.5 5 0]^⊺ are shown on the right (b = 0.3324).

Figure 23 Map and contours of the measured A-weighted sum of a discrete phased planar source (η = 6.5°) with β = 0 (0 dB/dod) using b = 0 (left) and b = 0.3324 (right).

Figure 24 Map and contours of the measured A-weighted sum of a discrete phased planar source (η = 6.5°) with β = 0.25 (−1.5 dB/dod) using b = 0 (left) and b = 0.3324 (right).

Figure 25 Map and contours of the measured A-weighted sum of a discrete phased planar (η = 6.5°) source with β = 0.5 (−3 dB/dod) using b = 0 (left) and b = 0.3324 (right).

Considering the plots for b = 0, the width of the contours are similar to those of the simulations of the planar sources with discrete geometry only. The main difference lies in farther points of observation where jagged contours are visible. These are an effect of the discretization of geometry and phase as well as the directivity of the used drivers. Nevertheless, it can be seen that the desired level distribution is achieved with the coarse discretization when applying an inclination.

For b = 0.3324, similar discretization effects can be seen. Comparing to the simulation results of discrete geometry and phase load, it should be noted that the measurements results come closer to the simulation results of discrete geometry only. Also here, the advantage of the tilted source becomes apparent, as it requires less delay to achieve the desired level supply. However, further improvements can be expected by using a multi-way system, where smaller tweeters can be used and placed closer together.

5 Conclusion

With the aim of achieving sound pressure levels on the listening plane rolling-of by 6 ⋅ β dB/dod, this paper proposed a design method for the delay load of phased planar sources. In contrast to a vertical phased line, another objective was to be able to variably adjust the horizontal coverage of the source. Starting from the integral of a Green’s function along two dimensions, stationary-phase approximation was applied to obtain a second-order non-linear differential equation of the delay length profile. This equation is under-determined in terms of the three curvature parameters. To solve it, a manual separation into three design equations is proposed.

Simulations of the continuous source based on the summation of point sources showed the optimal −6 ⋅ β dB/dod maps. Furthermore, the simulation shows the effects caused by discretization of the geometry and also by discretization of the phase load that is necessary in practice but also outlined the limits of a coarse discretization. The simulations therefore indicate that acceptable broadband results can be achieved with smaller drivers and smaller driver spacing. It has been shown that an inclined source will reduce the amount of delay for easier practical applicability. Measurements of a prototypical setup showed the practical feasibility.

Future work should consider (virtual) listening experiments with different β profiles that compare the phased planar arrays with line arrays for medium-sized surround sound reinforcement for 50–250 listeners. Moreover, waveguide solutions should be investigated to improve radiation with reduced spatial aliasing. In order to reduce the spatial aliasing even further, multi-way systems could be investigated in which the smaller tweeters are spaced more closely. Furthermore, the impact of equalizing strategies for individual elements have to be considered.

Funding

Our research was funded by the Austrian Science Fund (FWF) under project number P 35254-N, Envelopment in Immersive Sound Reinforcement (EnImSo).

Conflicts of interest

The authors declared no conflicts of interests.

Data availability statement

The research data associated with this article are available in KUG Phaidra data repository, under the reference [45].

Appendix A

Deriving the z-component of the stationary-phase direction

From the derivatives of r in equation (7) and the definition of the direction u = [u_x u_y u_z]^⊺, we know that $\frac{\partial r}{\partial v}=-\sin \eta u_x-\cos \eta u_z$. The stationary-phase direction u must be a unit vector. We may replace u_x by $\sqrt{1-u_y^2-u_z^2}$ to get,

$$\frac{\partial r}{\partial v}+\cos \eta u_z=-\sin \eta \sqrt{1-u_y^2-u_z^2}.$$(A1)

Taking the square of both sides yields a quadratic equation for u_z,

$$u_z^2+2\frac{\partial r}{\partial v}\cos \eta u_z+{\left(\frac{\partial r}{\partial v}\right)}^2+{\sin }^2\eta (u_y^2-1)=0,$$(A2)

and by applying sin²η + cos²η = 1, the z component of the direction gets

$$u_z=-\frac{\partial r}{\partial v}\cos \eta \underline{(+)}\sin \eta \sqrt{1-{\left(\frac{\partial r}{\partial v}\right)}^2-u_y^2}.$$(A3)

Because we are only interested in a stationary-phase direction pointing downwards on the listening area, the negative sign must be selected.

Appendix B

Equalization of a planar source

We consider a phased planar source of infinite length and width that is expanded over the vertical axis z and the horizontal axis y. Integration the Green’s function over y and z yields its sound pressure,

$$p={\int }_{-\infty }^{\infty }\mathrm{}{\int }_{-\infty }^{\infty }\mathrm{}\frac{e^{-\mathrm{i} k (r+w)}}{4 \pi r}\mathrm{d}y \mathrm{d}z.$$(B1)

i denotes the imaginary unit, k the wave number, r the distance between source and receiver and w the delay expressed as delay length. By applying the stationary-phase approximation, we get

$$p\approx \frac{e^{-\mathrm{i} k (r+w{)}_{\mathrm{*}}-\mathrm{i} \frac{\pi }{4} \mathrm{sig} \{\mathit{H}(r+w{)}_{\mathrm{*}}\}}}{2k{r}_{\mathrm{*}}\sqrt{\left|\mathit{H}\right(r+w{)}_{\mathrm{*}}|}}{|}_{w \equiv 0}=\frac{e^{-\mathrm{i} k \left|x\right|}}{2 \mathrm{i} k},$$(B2)

where |H(·)| is the determinant of the Hessian and sig{H(·)} its signature. * is used here as explicit indicator for the evaluation at the stationary-phase point. For zero phasing w ≡ 0, the stationary-phase distance is the normal distance r_* = |x| to the y-z plane and the exact result of equation (B1) is obtained: the one-dimensional Green’s function $\frac{e^{-i k \left|x\right|}}{2 i k}$. In any case, the resulting sound pressure is proportional to $\frac{1}{k}$. So, its frequency response can be said to be colored red, in analogy to the colors of noise (red, pink, white). For a neutral frequency response, a filter proportional to $k=\frac{\omega }{c}$ is required, i.e. introducing a 6 dB/octave increase.

References

All Tables

All Figures

	Figure 1 Geometry of a planar source (gray rectangle) with horizontal extent H and vertical extent V showing the unit vectors eh in horizontal direction and ev in vertical direction.
In the text

Figure 2 Side view of a planar source (solid gray) inclined by η showing the stationary-phase direction u, cf. Section 2.3 and the beamforming angles ${\vartheta }_0=-\arcsin \frac{\partial w}{\partial v}{|}_{v =\mathrm{ }0, h = 0}$ on top and ${\vartheta }_V=-\arcsin \frac{\partial w}{\partial v}{|}_{v = -V, h = 0}$ on the bottom of the source for a desired depth coverage between xr0 and xrV, cf. Section 2.3.

In the text

	Figure 3 Discretization of geometry of a phased planar source.
In the text

	Figure 4 Side view of simulated sound pressure level maps at y = 0 for phased planar arrays in 0 dB/dod configuration comparing the continuous source to a source with discretized geometry for f = 1 kHz and f = 6 kHz.
In the text

	Figure 5 Directivity index of phased planar sources with b = 0 (solid lines) and b = 0.2531 (dashed lines) for different decay parameters β.
In the text

	Figure 6 Delay lengths of a continuous phased planar source mounted straight on the wall (left, η = 0°) compared to the delay lengths of a continuous phase planar source inclined by η = 6.5° with β = 0 and b = 0.3324.
In the text

	Figure 7 Experimental setup using line array enclosures of [42] mounted as 7 × 8 element planar array.
In the text

	Figure 8 Delay lengths for a continuously phased planar source (η = 0°), b = 0 (left) and b = 0.2531 (right) for β = 0 (0 dB/dod).
In the text

	Figure 9 Delay lengths for a continuously phased planar source (η = 0°), b = 0 (left) and b = 0.2531 (right) for β = 0.25 (−1.5 dB/dod).
In the text

	Figure 10 Delay lengths for a continuously phased planar source (η = 0°), b = 0 (left) and b = 0.2531 (right) for β = 0.5 (−3 dB/dod).
In the text

	Figure 11 Map and contours of the simulated A-weighted sum of a continuous phased planar source (η = 0°), b = 0 (left) and b = 0.2531 (right) for β = 0 (0 dB/dod).
In the text

	Figure 12 Map and contours of the simulated A-weighted sum of a continuous phased planar source (η = 0°), b = 0 (left) and b = 0.2531 (right) for β = 0.25 (−1.5 dB/dod).
In the text

	Figure 13 Map and contours of the simulated A-weighted sum of a continuous phased planar source (η = 0°), b = 0 (left) and b = 0.2531 (right) for β = 0.5 (−3 dB/dod).
In the text

	Figure 14 Map and contours of the simulated A-weighted sum of a discrete phased planar source (η = 0°) with continuous delay profile, b = 0 (left) and b = 0.2531 (right) for β = 0 (0 dB/dod).
In the text

	Figure 15 Map and contours of the simulated A-weighted sum of a discrete phased planar source (η = 0°) with continuous delay profile, b = 0 (left) and b = 0.2531 (right) for β = 0.25 (−1.5 dB/dod).
In the text

	Figure 16 Map and contours of the simulated A-weighted sum of a discrete phased planar source (η = 0°) with continuous delay profile, b = 0 (left) and b = 0.2531 (right) for β = 0.5 (−3 dB/dod).
In the text

	Figure 17 Map and contours of the simulated A-weighted sum of a discrete phased planar source (η = 0°) with discretized delay profile, b = 0 (left) and b = 0.2531 (right) for β = 0 (0 dB/dod).
In the text

	Figure 18 Map and contours of the simulated A-weighted sum of a discrete phased planar source (η = 0°) with discretized delay profile, b = 0 (left) and b = 0.2531 (right) for β = 0.25 (−1.5 dB/dod).
In the text

	Figure 19 Map and contours of the simulated A-weighted sum of a discrete phased planar source (η = 0°) with discretized delay profile, b = 0 (left) and b = 0.2531 (right) for β = 0.5 (−3 dB/dod).
In the text

	Figure 20 Map and contours of the simulated A-weighted sum of a discrete phased planar source (η = 0°) with discretized delay profile, β = 0 (0 dB/dod), b = 0 (left) and b = 0.2531 (right) using a smaller discretization.
In the text

	Figure 21 Map and contours of the simulated A-weighted sum of a discrete phased planar source (η = 0°) with discretized delay profile, β = 0.25 (−1.5 dB/dod), b = 0 (left) and b = 0.2531 (right) using a smaller discretization.
In the text

	Figure 22 Map and contours of the simulated A-weighted sum of a discrete phased planar source (η = 0°) with discretized delay profile, β = 0.5 (−3 dB/dod), b = 0 (left) and b = 0.2531 (right) using a smaller discretization.
In the text

	Figure 23 Map and contours of the measured A-weighted sum of a discrete phased planar source (η = 6.5°) with β = 0 (0 dB/dod) using b = 0 (left) and b = 0.3324 (right).
In the text

	Figure 24 Map and contours of the measured A-weighted sum of a discrete phased planar source (η = 6.5°) with β = 0.25 (−1.5 dB/dod) using b = 0 (left) and b = 0.3324 (right).
In the text

	Figure 25 Map and contours of the measured A-weighted sum of a discrete phased planar (η = 6.5°) source with β = 0.5 (−3 dB/dod) using b = 0 (left) and b = 0.3324 (right).
In the text