# Constant-Beamwidth Kronecker Product Beamforming With Nonuniform Planar Arrays Ariel Frank, et al.

May 13, 2022

## 1 Introduction

Speech, wireless communications, and radar are among the fields that use broadband beamforming (Van Veen and Buckley, 1988; Van Trees, 2004; Cohen et al., 2009). A broadband array is desirable since audible speech covers a wide range of frequencies (Snow, 1931). A constant beamwidth prevents distortion of the speech signal when the speaker moves away from the center of the mainlobe. If the beamformer weights are consistent versus frequency, the beamwidth narrows as the frequency increases (Ward et al., 1995). Therefore, the beamformer weights should be frequency-dependent to enable control of the beamwidth. Nevertheless, the sensor positions set bounds on the maximum and minimum frequency that can attain the desired beamwidth. A beamformer’s directivity is an indicator of its performance in various scenarios, such as reverberant environments. The purpose of our research is to design, with low computational complexity and a minimal number of sensors, a beamformer whose weights and positions are chosen to maximize the directivity while maintaining a constant beamwidth.

To calculate the weights for a frequency-invariant beamformer, several methods utilized the inverse Fourier transform (Ward et al., 1995, 2001; Liu and Weiss, 2004; Liu et al., 2007; Liu and Weiss, 2008; Pal and Vaidyanathan, 2010; Weiss et al., 2017), while others opted to defining an optimization problem. For simple objectives and constraints, closed form solutions were found (Parra, 2006; Zhao et al., 2009; Crocco and Trucco, 2011; Tourbabin et al., 2012; Zhao and Liu, 2012; Huang et al., 2020b), while others resorted to using optimization toolboxes (Gazor and Grenier, 1995; Yan and Ma, 2005; Chan and Chen, 2006; Yan, 2006; Zhao et al., 2008; Yan et al., 2010; Zhu and Wu, 2011; Yin, 2012; Liu et al., 2015; Li et al., 2017; Buchris et al., 2018; Liu et al., 2018; Yang X. et al., 2019; Liu et al., 2021; Son, 2021) and even deep learning methods (Aroudi and Braun, 2021; Ramezanpour et al., 2022). Recently, the Kronecker product has been utilized for beamforming to decompose the beamformer design problem into smaller problems (Yang W. et al., 2019; Benesty et al., 2019; Cohen et al., 2019; Itzhak et al., 2019; Wang et al., 2019; Huang et al., 2020a; Huang et al., 2020c; Sharma et al., 2020; Wang et al., 2020; Yang et al., 2020; Itzhak et al., 2021; Kuhn et al., 2021; Wang et al., 2021), yet it has not been applied for constant-beamwidth beamforming.

Most frequency-invariant beamformers were designed for uniform linear arrays (ULA) (Liu and Weiss, 2004; Liu et al., 2007; Liu and Weiss, 2008; Zhao et al., 2008, 2009; Liu and Weiss, 2010; Haupt and Moosbrugger, 2012; Zhao and Liu, 2012; Li et al., 2017; Rosen et al., 2017; Weiss et al., 2017; Yang X. et al., 2019; Long et al., 2019; Huang et al., 2020b; Erokhin et al., 2020; Liu et al., 2021), yet some were also developed for linear arrays (Ward et al., 1998; Parra, 2006; Pal and Vaidyanathan, 2010; Crocco and Trucco, 2011; Liu et al., 2015, 2018; Son, 2021), logarithmically spaced arrays (Ward et al., 2001), circular arrays (Chan and Chen, 2006; Yan et al., 2010; Buchris et al., 2020; Sharma et al., 2021a; Sharma et al., 2021b; Kleiman et al., 2021), and arbitrary geometry arrays (Ward et al., 1995; Yan and Ma, 2005; Yan, 2006; Zhu and Wu, 2011; Tourbabin et al., 2012; Yin, 2012). Different geometries yield varying performances. Therefore, choosing where to position the sensors has been addressed in (Gazor and Grenier, 1995; Liu et al., 2015; Buchris et al., 2018; Liu et al., 2018; Son, 2021) to maximize performance measures such as white noise gain (WNG). However, these methods require intensive computations. To alleviate the complexity of designing a two-dimensional (2D) array, we show how a constant beamwidth can be attained based on designing one-dimensional (1D) arrays. We introduce two different techniques. The first one utilizes Kronecker product properties, whereas the second technique uses a matrix transform to convert the linear array weights into planar array weights. With the second technique, we can tradeoff between WNG and directivity even when removing some sensors.

The beamwidth in the XZ plane may differ from the beamwidth in the YZ plane. This can be accomplished by designing two linear arrays, each with a different beamwidth. We position the planar array sensors based on the positions of the linear array sensors. We transform the weights used for the linear arrays to get weights for the planar beamformer. The transform ensures that the beamwidth remains constant in the XZ and YZ planes. The proposed techniques enable constant-beamwidth beamforming using planar arrays with many sensors by designing simpler linear arrays. Our design generalizes the window-based technique introduced by Long et al. (Long et al., 2019) for symmetric nonuniform linear arrays. Allowing nonuniform spacing between sensors generates additional degrees of freedom, which enables maximizing the directivity while maintaining a constant beamwidth over a broader range of frequencies. The beamformer weights are designed to take full advantage of the nonuniform interelement spacing of the symmetric array. Compared to existing methods, our generalization improves the directivity and attains a constant beamwidth over a broader range of frequencies, even for ULAs. We present two different procedures for choosing the symmetric linear array sensor positions. The first procedure is computationally intensive and involves iteratively updating the sensor positions. The purpose of this procedure is to find an upper bound on the achievable performance. Afterward, we present the second procedure that is computationally simple while maintaining on-par performance. Simulations demonstrate the properties, low design complexity, and superior performance of our design. To facilitate reproducibility, our code is available online1.

The paper is structured as follows: Section 2 describes the signal model, performance measures, and the problem formulation. Section 3 presents beamformers for planar array geometries based on beamformers designed for linear arrays. In Section 4, we develop algorithms for finding optimal placements of the sensors for a symmetric linear array geometry and optimal design of the beamforming weights. Section 5 contains simulations and experimental results, and Section 6 concludes.

## 2 Signal Model and Problem Formulation

### 2.1 Signal Model

Consider a nonuniform planar array whose sensors are positioned on an M × N grid where the x-coordinates are

${xm}m=1M$

and the y-coordinates are

${yn}n=1N$

. The sensor at (m, n) = (1, 1) is referred to as the reference sensor. We assume a far-field scenario, where a plane wave propagates in an anechoic environment at the speed of sound in the air, v = 343 m/s. Its direction of arrival is parameterized by the azimuth angle φ, 0° ≤ φ ≤ 360°, and elevation angle θ, 0° ≤ θ ≤ 180°. In the frequency domain, the observed signal at frequency f can be written in vector form as

where X(f) is the source signal arriving from the direction (φd, θd), q(f) is the additive noise vector of length MN, and a (f, φ, θ) is the steering vector of length MN whose entries are given by

$e−J2πfvxm⁡cosφsinθ+yn⁡sinφsinθ;1≤m≤M,1≤n≤N;(2)$

where

$J=‐1$

.

The beamformer output is given by

where

$hf$

is the beamforming weight vector of length MN and H denotes conjugate transpose. A frequency-invariant gain is desired in the direction (φd, θd). Therefore, we impose the common distortionless constraint:

The beampattern expresses the response of the beamformer to a plane wave and is defined as

$Bf,φ,θ,hf=aHf,φ,θhf.(5)$

To improve readability, we sometimes use a more concise notation for the steering vector, a, beamforming weight vector, h, and beampattern,

$B(h)$

, by omitting the variables f, φ, and θ.

### 2.2 Performance Measures

We now develop wideband performance measures to quantify the WNG and directivity over a frequency band. The narrowband and wideband input signal-to-noise ratios (SNR) are

and

$iSNR=∫fϕXfdf∫fϕQ0fdf,(7)$

respectively, where

$ϕXf$

is the variance of

$Xf$

and

$ϕQ0f$

is the variance of the additive noise at the reference sensor. From (1) and (3), the narrowband and wideband output SNRs are

$oSNRf=ϕXfhHfaf,φd,θd2hHfΦqfhf=4ϕXfhHfΦqfhf(8)$

and

$oSNR=∫fϕXfhHfaf,φd,θd2df∫fhHfΦqfhfdf=4∫fϕXfdf∫fhHfΦqfhfdf,(9)$

respectively, where

$Φqf$

is the correlation matrix of

$qf$

. The narrowband and wideband array gains are given by

and

respectively.

When

$qf$

is spherically isotropic noise,

$Φqf$

is given by

where the elements of

$Γdf$

are given by

$Γdfr,r′=sinc2πfvdr,r′;1≤r,r′≤MN;(13)$

where dr,r is the Euclidean distance between the corresponding sensors.

The narrowband directivity factor (DF) is defined as the narrowband array gain in the presence of spherically isotropic noise (Benesty et al., 2018) and is given by

We define the wideband directivity factor over the frequency band

$fL,fH$

as the wideband array gain in the presence of spherically isotropic noise. For scenarios where

$ϕQ0f$

is frequency-invariant, i.e.,

$ϕQ0f≡ϕQ0$

, a compact expression for the wideband directivity factor is given by:

$DfL,fH=11fH−fL∫f=fLfHhHfΓdfhfdf=1411fH−fL∫f=fLfH1Dfdf.(15)$

From here, we define the wideband directivity index (DI) over the frequency band

$fL,fH$

as

$DIfL,fH=10⁡log10DfL,fH.(16)$

WNG is a commonly used performance measurement that reflects the robustness of an array to physical imperfections (Benesty et al., 2018). It is defined as the array gain for white spatial noise, i.e.,

$Φq(f)=ϕQ0(f)IMN$

, where IMN is the MN × MN identity matrix. For scenarios where

$ϕQ0f$

is frequency-invariant, the narrowband and wideband WNGs are given by the compact expressions:

and

$WfL,fH=11fH−fL∫f=fLfH1Wfdf,(18)$

respectively.

### 2.3 Problem Formulation

We assume that the number of sensors, MN, is given. Our goal is to find optimal placements for the sensors and optimal beamformer weights that maximize the DI while maintaining a constant beamwidth and acceptable WNG.

Let bX and bY denote the half-power beamwidths in the XZ and YZ planes, respectively. Then, our optimization problem is defined as

$maximizexm,1≤m≤Myn,1≤n≤Nhf,fL≤f≤fHDIfL,fHsubjecttohHfaf,φd,θd=1,fL≤f≤fHbXf=θX,fL≤f≤fHbYf=θY,fL≤f≤fHWfL,fH=ℵ,(19)$

where θX and θY are the desired half-power beamwidths in the XZ and YZ planes, respectively, and ≥1 for acceptable wideband WNG.

## 3 Planar Array Constant-Beamwidth Beampattern Design

This section shows how a constant beamwidth can be developed for a planar array based on designs of linear arrays. The beamwidth in the XZ plane may differ from the beamwidth in the YZ plane. This can be accomplished by designing two linear arrays, each with a different beamwidth.

A linear array beampattern is designed such that when the linear array lays along the x-axis, the beampattern maintains the distortionless constraint in the broadside direction, i.e., (φd, θd) = (90°, 90°). Its beamwidth, bφ, is measured in the XY plain. Our goal is to design a beampattern with a half-power beamwidth of θX in the XZ plane and θY in the YZ plane. Assume that the two linear arrays are given (in the next section, we suggest a method for designing linear arrays, yet the designs in this section are not limited to that method). The first linear array has M sensors with a constant beamwidth of θX. Denote the positions by

${xm}m=1M$

and the frequency-dependent weights by hX(f). The second linear array has N sensors with a constant beamwidth of θY. Denote the positions by

${yn}n=1N$

and the frequency-dependent weights by hY(f). To simplify the analysis, we assume that xm = 0 for some m and that yn = 0 for some n. We position the planar array sensors based on the positions of the linear array sensors. We transform the weights used for the linear arrays to get weights for the planar beamformer. The transform ensures that the beamwidth remains constant in the XZ and YZ planes. The transform does not depend on the frequency, only on the physical locations of the sensors. The transform enables constant-beamwidth beamforming using planar arrays with many sensors by designing simpler linear arrays.

First, we show how to use Kronecker product beamforming (Benesty et al., 2019) for the nonuniform planar array. Then we offer how the beamwidth can remain constant using fewer sensors; the proposed method enables a tradeoff between WNG and directivity.

### 3.1 Kronecker Product Beamforming

From the two given linear arrays, we construct an M × N grid where the x-coordinates are

${xm}m=1M$

and the y-coordinates are

${yn}n=1N$

—creating a planar array with a total of MN sensors. Let aP (f, φ, θ) denote the steering vector of length MN whose entries are given in (2). Denote by aX (f, φ, θ) and aY (f, φ, θ) the steering vectors (with respective lengths M and N) of the subarrays that lay on the x and y axes whose entries are given by

$e−J2πfvxm⁡cos(φ)sin(θ),1≤m≤M,(20)$

and

$e−J2πfvyn⁡sin(φ)sin(θ),1≤n≤N,(21)$

respectively. Notice that

where ⊗ denotes the Kronecker product.

The beampatterns of the subarrays are given by

and

For a linear array, because it exists along a 1D line, every 2D slice (that contains the 1D line) of its beampattern is identical. For the subarray lying on the x-axis (y-axis), every 2D slice (containing the x-axis (y-axis)) of its beampattern is identical; For example, the XZ (YZ) and XY planes. Therefore, by the constant-beamwidth design of the linear arrays:

$BhX=1,φ=90°orθ=0°,12,φ,θ=0°,12θX,BhY=1,φ=0°,12,φ,θ=90°,12θY.(25)$

Let h3 be a 2D window constructed as the outer product of two 1D windows:

$h3m,n=h1mh2n;1≤m≤M,1≤n≤N;(26)$

where h1 and h2 are 1D windows with respective lengths M and N. For example, each 1D window can be a Kaiser window—and then the 2D window is a 2D Kaiser window. By the construction, h3 is separable in the two variables. Therefore, if we arrange its values into a vector of length MN, the resulting vector can be formulated as the Kronecker product of the two 1D windows. We design the planar array beamforming weight vector, hP, as

Following the methodology in (Benesty et al., 2019) (while generalizing for a nonuniformly spaced rectangular array), the beampattern of the planar array is given by

$BhP=aPHf,φ,θhPf=aY⊗aXHhY⊗hX=aXHhXaYHhY=BhXBhY,(28)$

where for the third equality we used the Kronecker product identity

$u1⊗u2Hv1⊗v2=u2Hv2u1Hv1.(29)$

The beamwidth in the XZ and YZ planes is attained by substituting (25) into (28):

$BhP=1,φ,θ=0°,0°,12,φ,θ=0°,12θX,12,φ,θ=90°,12θY.(30)$

The resulting beampattern maintains the desired beamwidth in each plane and satisfies the distortionless constrain. This shows that we can decompose the 2D problem of designing a constant-beamwidth beamformer with a planar array to two 1D problems of designing constant-beamwidth beamformers with linear arrays. The Kronecker product is valid only when sensors are placed on all the MN grid positions. Next, we develop a method for attaining a constant-beamwidth beamformer even when only a subset of the grid positions have sensors.

Computational Complexity of the Kronecker Beamformer:

In this paper, we consider the computational complexity as the number of multiplications needed for the method. Assuming that the two linear arrays are given, the planar array beamforming weights are calculated with (27), requiring MN multiplications per frequency. Let F denote the number of frequency bins. Therefore, the computational complexity is

$O(MNF)$

.

### 3.2 Linear to Planar Weight Transform

This section presents a method that does not require a sensor at each grid position and can tradeoff between WNG and directivity. From the two linear arrays, we construct the M × N grid, but this time choose only S out of the possible MN positions to place a total of S sensors. Conditions regarding where to place the sensors are discussed later.

The beampattern is given by

$Bf,φ,θ,hPf=∑m∑naP*m,nhPm,n,(31)$

where the summation is only over S values (only the index pairs where there is a sensor) and * denotes complex conjugate. Here, hP [m, n] is the beamformer weight at the sensor positioned at (xm, yn). The other expressions in this section with square braces are similarly defined.

Let us examine the beampattern in the XZ plane. At elevation

$θ=12θX$

:

$Bf,φ=0,θ=12θX,hPf=∑m∑neJ2πfvxm⁡sin12θXhPm,n=∑meJ2πfvxm⁡sin12θX∑nhPm,n.(32)$

We desire that the absolute value of (32) will be equal to

$1/2$

, which is precisely the amplitude of the linear array beampattern at

$12θX$

(25):

$BhX=∑meJ2πfvxm⁡sin12θXhXm=12.(33)$

Therefore, we require:

$∑meJ2πfvxm⁡sin12θX∑nhPm,n=∑meJ2πfvxm⁡sin12θXhXm.(34)$

To attain (34), it is enough to require:

$∑nhPm,n=hXm,1≤m≤M.(35)$

Similarly, to attain the desired beamwidth in the YZ plane, we require:

$∑mhPm,n=hYn,1≤n≤N.(36)$

Altogether, (35) holds M constraints and (36) holds N constraints. Therefore, we have to solve a linear system of M + N equations with S parameters (the weights for the sensors, hP). In matrix notation:

where C is an (M + N) × S matrix,

$hL≐hXT,hYTT$

, and T denotes transpose. We only need to calculate C once; notice that it is not a function of frequency—all the values of C are 0 or 1.

A solution to (37) exists when the rank of

$C,hL$

equals rank{C}. When S = MN, then rank{C} = M + N − 1 (the rows of C are linearly dependent because the sum of the first M rows equals the sum of the last N rows). The sum of the elements in hX equals the sum of the elements in hY and equals 1 to fulfill the distortionless constraint. It follows that

$rank{C,hL}=M+N−1$

. Therefore, when considering removing a sensor, we must ensure that without the corresponding column, the rank of C remains M + N−1. The beamformer weights are found following the methodology in (Benesty et al., 2018).

The following optimization problem maximizes the WNG under constraint (37):

$hPW=argminhhHhs.t.Ch=hL.(38)$

Its solution is given by:

which is also the least-squares solution to (37). The following optimization problem maximizes the directivity under constraint (37):

$hPD=argminhhHΓhs.t.Ch=hL,(40)$

where the elements of Γ are given by

$Γfs,s′=sinc2πfvds,s′;1≤s,s′≤S;(41)$

where ds,s is the Euclidean distance between the sensors corresponding to the columns s and s′ of C. Its solution is given by:

$hPD=Γ−1CHCΓ−1CH−1hL.(42)$

A tradeoff between maximizing the WNG and directivity can be attained with weights given by

$hPα=Γα−1CHCΓα−1CH−1hL,(43)$

where Γα = (1 − α)Γ + αIS for 0 ≤ α ≤ 1. Selecting α controls the tradeoff between WNG and directivity: a bigger α increases the WNG while a smaller α increases the directivity. We can set α to be frequency-dependent.

Computational Complexity of the Tradeoff Beamformers:

We assume that the two linear arrays are given. First, we analyze the computational complexity of the maximum WNG beamformer (39). Recall that C is not a function of frequency. Therefore,

$CH(CCH)−1$

can be calculated once with

$OSM2+N2+M3+N3$

operations. As SM + N − 1 for a viable solution, the complexity can be written in a more compact form:

$OSM2+N2$

. Afterward, multiplying by hL requires

$SM+N$

multiplications per frequency. Therefore, the total computational complexity is

$OSM+NM+N+F$

.

Next, we analyze the computational complexity of the maximum directivity beamformer (42) and tradeoff beamformer (43). The matrix multiplications and inversions require

$OS3$

computations per frequency. Therefore, the total computational complexity is

$OS3F$

.

## 4 Linear Array Constant-Beamwidth Beampattern Design

This section presents a method for obtaining beamforming weights for a symmetric nonuniform linear array, assuming that the sensor positions are given. Then we show how to choose the sensor positions optimally.

Consider a symmetric linear array with M = 2L + 1 sensors. We use a symmetric array so that the array phase center is frequency invariant (Ward et al., 2001) and so that we can simply generalize the window-based technique. The sensor positions are denoted by xl, − LlL. Due to the symmetry of the array,

$xl=0,l=0,x−l=−xl,1≤l≤L.(44)$

### 4.1 Beamforming Weight Vector Design

In this section, we assume that the sensor positions are given. Our goal is to design the beamformer weights for maximum DI while maintaining a constant beamwidth of φB. Formally,

$maximizehf,fL≤f≤fHDIfL,fHsubjecttohHfaf,φd,θd=1,fL≤f≤fHbφ=φB,fL≤f≤fH.(45)$

A simple method for constant-beamwidth beamforming is to constrain the beamformer weights to equal the taps of a window. The window has a parameter that is tuned so that the beamwidth is as desired. This method was introduced by Long et al. in (Long et al., 2019) for ULAs. This method compromises between approximating the desired beampattern and computational complexity. We generalize this method for symmetric linear arrays while using the Kaiser window. Compared to the existing form, our generalization improves the DI. It enables maintaining the desired beamwidth over a broader range of frequencies, even for ULAs, as is shown in the results.

The Kaiser window depends on one parameter, the window shape factor, β. For each frequency,

$βf$

is chosen so that the beamwidth is as desired. Due to the inverse relationship between the window width and the beamformer width (Ward et al., 1995), increasing β widens the beampattern. This inverse relationship can be intuitively understood from (5) by noticing that the beampattern is the (nonuniform) discrete Fourier transform of the beamformer weights. If the weights are frequency independent, the beampattern narrows as the frequency increases. Therefore, higher values of β are used at higher frequencies to widen the beampattern so that the beamwidth remains constant.

We do not limit the array to be uniform; therefore, we sample the continuous Kaiser window at the sensor positions. The continuous Kaiser window (Kaiser and Schafer, 1980) is given by

$wKx=I0β1−xW/22I0β,x≤W/2,0,x>W/2,(46)$

where

$I0x$

is the zeroth-order modified Bessel function of the first kind, and W is the window support. The Kaiser window (of length R + 1) used for digital signal processing is given by

$wULARl=wKx=WRl=I0β1−lR/22I0β,l≤R/2,0,l>R/2,(47)$

where R + 1 = M. We suggest sampling the continuous Kaiser window at the sensor positions, i.e.,

$wLAl=wKx=xl=I0β1−xlW/22I0β,xl≤W/2,0,xl>W/2.(48)$

Algorithm 1. Algorithm for calculating the beamformer weights.

By choosing the window support to equal W = 2xL (as done in (Long et al., 2019)), we are left with a single parameter to tune, β. We increase the degrees of freedom of the problem by optimizing over the parameter W. We consider multiple values of W (per frequency), and for each W, tune β so that the beamwidth is as desired. Then, per frequency, we select the W that achieves the highest DF. To reduce the number of possible values of W we limit the window support to the set

$Wii=1L$

of size L where Wi = 2xi. Here, Wi corresponds to having only the 2i + 1 center sensors having nonzero weights. We refer to these sensors as active sensors. In the results, we show that the performance improvement of not limiting the values of W to this set is negligible. Therefore, we limit the values of W to this set to reduce the computational complexity of finding the weights. Notice that W is optimized per frequency, i.e., W = W(f).

To reflect that some sensors are closely spaced, we use the trapezoidal integration technique (TIT), i.e., we multiply each weight by

$Δl=xl+1−xl−12$

as was done in (Ward et al., 2001). We multiply the weight of the outermost sensors by ΔL = ΔL = xLxL−1. Altogether,

$Δl=xl+1−xl−12,l

By using the TIT at frequencies where all of the sensors are used (W = 2xL), the center sensors that are close together have less influence. This enables maintaining the desired beamwidth at low frequencies.

Last, we normalize the weights to fulfill (4). To conclude, the beamformer weights are:

$hl=ΔlwKxl∑l′=−LLΔl′wKxl′,−L≤l≤L.(50)$

After the weights are found for each frequency, the DI can be evaluated. Algorithm 1 summarizes the process. The resulting beamformer weights,

$ĥ(f)$

, satisfy the constraints in (45). So, we proceed with determining where to place the sensors.

Computational Complexity of Algorithm 1:

Calculating the weights with (50) requires M multiplications for computing the products

$ΔlwKxl$

, − LlL. Because we begin with the largest β, the beamwidth is initially wider than desired, and thus the absolute value of the beampattern at azimuth

$φ=φd+12φB$

is greater than

$1/2$

. Therefore, to check whether bφ > φB, it is sufficient to evaluate the absolute value of the beampattern at azimuth

$φ=φd+12φB$

and check whether it is greater than

$1/2$

, requiring M multiplications. Let B denote the number of window shape factor options. The beampattern is evaluated up to F + B times per window support option, thus evaluated for up to a total of

$F+BL$

times. Altogether, the first part of the algorithm requires

$OM2F+B$

computations.

In the second part of the algorithm, calculating the DF requires M multiplications for calculating the weights and M2 + M multiplications for the matrix multiplications, per frequency and window support option. Therefore, requires

$OM3F$

computations.

In total, Algorithm 1 requires

$OM2MF+B$

computations.

### 4.2 Symmetric Linear Array Sensor Positioning

Our goal is to find where to position the sensors such that the DI is maximized while using the weight design proposed in the previous subsection. Formally,

$maximizexl,−L≤l≤LDIfL,fH,(51)$

where the weights are attained using Algorithm 1; these weights fulfill the constraints in (45). A brute force search over all possible positions is infeasible. So, we first present an iterative algorithm for determining the sensor positions. The resulting performance is an upper bound on the DI (using the weight design proposed in the previous section). Afterward, we present a simpler algorithm that achieves competitive performance.

#### 4.2.1 Iterative Algorithm

The sensor positions are iteratively updated to maximize the DI using the Nelder-Mead simplex algorithm (Lagarias et al., 1998). A set of L variables represent the sensor positions

$xll=1L$

. The

$DIfL,fH$

is evaluated in each iteration by constructing the symmetric array by (44), finding the beamformer’s weights according to Algorithm 1, and finally calculating the DI with (16). The algorithm is initialized multiple times with different starting points to ensure that a global maximum is found. Therefore, the resulting DI is an upper bound on the DI (using the weight design proposed in the previous section).

#### 4.2.2 Non-Iterative Algorithm

In this section, we outline steps to achieve a constant beamwidth over a maximum range of frequencies. Afterward, we modify one of the steps to improve the DI at the expense of reducing the frequency range.

The algorithm is initialized with a symmetric linear array with M0 = 2L0 + 1 (<M) sensors. When these sensors are active, the lowest frequency that can maintain the desired beamwidth is attained using the widest window (β = 0). Denote this frequency by

$fmin(L0)$

. We add a pair of sensors with symmetric positions. Denote by

$d(L0+1)$

the interelement spacing between the new sensor to the previously outermost sensor. It is desirable for

$d(L0+1)$

to be large so that the array can support low frequencies. However, if it is too large, the beamwidth at higher frequencies (specifically, at

$fmin(L0)$

) is too narrow. So,

$d(L0+1)$

is chosen such that above

$fmin(L0)$

the new pair of sensors are inactive and below

$fmin(L0)$

the new pair of sensors are active. Therefore, we increase

$d(L0+1)$

so long as the array can maintain the desired beamwidth at

$fmin(L0)$

. Because

$fmin(L0)$

is the highest frequency that the new sensors are active at, the desired beamwidth is attained there with the narrowest window: β = βmax (= 10). The process of adding a pair of sensors is repeated until the array consists of M sensors.

As displayed in the results, the sidelobes are high for low values of β: decreasing the DI. Therefore, we limit the widest window to β = βmin > 0. This modification improves the DI but reduces the frequency range because the beampattern cannot be as narrow as when β = 0. We need to initialize with M0 = 5 sensors because three cannot attain a narrow enough beamwidth at fH = 8 kHz without exhibiting high sidelobes. An option how to choose the initial M0 sensors: a ULA with M0 sensors whose interelement spacing is determined to maximize the DF at fH. The steps are detailed in Algorithm 2.

Algorithm 2. Algorithm for computing the optimal sensor positions.

The algorithm can be used to determine how many sensors are necessary to maintain a constant beamwidth over the frequency band [fL, fH]. Once ffL, we can finish adding new pairs of sensors. Once the positions are found with Algorithm 2, the beamformer weights are obtained with Algorithm 1. Computational Complexity of Algorithm 2:

As detailed in the complexity analysis of Algorithm 1,

$OM$

operations suffice for checking whether bφ > φB. Let D denote the number of position options. D is derived from the maximum permissible array aperture and the resolution in which d(l) is increased in each iteration. The beampattern is evaluated up to D times in the part Find where to place the next pair of sensors, and up to F times in the part Find lowest frequency that maintains the desired beamwidth. Altogether, the beampattern is evaluated up to D + F times. Therefore, Algorithm 2 requires

$OMD+F$

computations.

## 5 Experimental Results

In this section, we examine the performance of our proposed approach and compare it to existing methods. We first examine linear arrays and compare our proposed iterative and non-iterative algorithms to existing methods. Next, we design planar arrays using the Kronecker and tradeoff beamformers.

### 5.1 Comparison of Linear Array With Previous Methods

In this section, we examine the performance of our proposed approach and compare it to previous approaches using M = 11 sensors. The desired half-power beamwidth is arbitrarily chosen to be φB = 15°. We compare our design method to three other methods: the sparse array design of Son (Son, 2021), the ULA design of Long et al. (Long et al., 2019), and the logarithmically spaced array design of Ward et al. (Ward et al., 2001). The array layouts are displayed in Figure 1.

FIGURE 1. Configurations of the arrays: proposed (triangles), Son (circles), Long (squares), and Ward (diamonds).

Using the iterative algorithm proposed in the previous section, the resulting sensor positions of the symmetric linear array (for l ≥ 1) were [3.8, 7.9, 14.3, 29.2, 74.8] cm. To reproduce the array of Son (Son, 2021), we set the reference beampattern to equal

which has a half-power beamwidth of 15°. We also checked other reference beampatterns, yet they resulted in similar performances as using (52). An initial sensor with the highest spanning energy was chosen from a ULA with 2001 sensors and an interelement spacing of 0.1 cm. The remaining 2000 sensors were divided into five groups. From each group, two sensors closest to the group centroid were chosen. Finally, the sparse array was created with the 11 chosen sensors. Next, the weights were found. The mainlobe error tolerances were 0.006E1 for φ ∈ [82.5°, 97.5°]. The sidelobe error tolerances were 0.028E2 for φ ∈ [0°, 82.5°) ∪ (97.5°, 180°]. Here, E1 = 31 and E2 = 330 were the number of azimuth angles sampled (the resolution was 0.5°). To ensure that the beamwidth was as desired, we added two constraints: that the error tolerances were 0.0001 for φ ∈ {82.5°, 97.5°}. The resulting sparse array was a symmetric linear array with real weights. The ULA weights used by Long et al. (Long et al., 2019) were reproduced using Algorithm 1, where the window support was limited to only equal the support of the array, i.e., W ≡ 2xL. For the comparison, we used an interelement spacing of 2.8 cm because it resulted in the largest

$DI0,8kHz$

while attaining a beamwidth of at least 15° at all frequencies. To reproduce the array of Ward et al. (Ward et al., 2001), we set the reference beampattern to equal (52). We used a logarithmic ratio of 3 (denoted by Q in (Ward et al., 2001)) because it resulted in the highest DI.

The optimal window parameters (β and W for our array and β for Long) were calculated using Algorithm 1 and are plotted in Figure 2 as a function of frequency. As expected, for frequencies with the same number of active sensors, β increases as the frequency increases. This is necessary to counteract the narrowing of the beampattern as the frequency increases. Increasing β narrows the Kaiser window, which in turn widens the beamwidth due to the inverse relationship between the window width and the beamformer width. We also observe that the number of active sensors (2i + 1) decreases as the frequency increases. Reducing the number of active sensors narrows the window, yielding an effect similar to increasing β.

FIGURE 2. Continuous Kaiser window parameters, (A) the window shape factor, and (B) the index of the window support from the set

$2xii=15$

, as a function of frequency: proposed (solid line) and Long (dashed line).

The beamformer weights (computed by (50) for our array and by (47) followed by normalization to fulfill (4) for Long) are illustrated in Figure 3. For the proposed and Ward arrays, the inner sensors are effectively inactive at low frequencies due to the TIT (the TIT does not affect ULAs due to their uniform interelement spacing). The outer sensors are given greater weight due to their greater distance between adjacent sensors (49). This enables maintaining the desired beamwidth at low frequencies by effectively widening the window. With these weights, the beampatterns were constructed according to (5) and are illustrated in Figure 3. The beamwidth, WNG, and DF are plotted in Figure 4 as a function of frequency. Son’s method performance is only reported at frequencies where the WNG is at least 0 dB (f ≥ 350 Hz). Our array maintained the desired beamwidth down to 620 Hz, while Long could not hold the desired beamwidth below 3.77 kHz. Ward’s array does not exhibit a constant beamwidth; Ward needed more sensors to succeed. Son’s array maintained the desired beamwidth down to 350 Hz but had lower WNG and DI than ours. Son’s and ours arrays’ WNG was lower than Long’s and Ward’s due to fewer active sensors at each frequency, as shown in Figure 3A.

FIGURE 3. The left subplots (A) illustrate the weight values of the arrays: (1) proposed, (2) Son, (3) Long, and (4) Ward. For each array, the weights are given per sensor as a function of frequency. The right subplots (B) illustrate the beampatterns. The overlayed dashed line is the half-power contour of the mainlobe.

FIGURE 4. (A) Beamwidth, (B) WNG, and (C) narrowband directivity factor as a function of frequency: proposed (solid line), Son (dashed line), Long (dash-dot line), and Ward (dotted line). The obtained

$W350Hz,8kHz$

are respectively 7.8, 7.6, 9.5, and 9.2 dB. The obtained

$DI350Hz,8kHz$

are respectively 8.1, 6.7, 6.5, and 6.3 dB.

From the DF, the DI was calculated using (16). The

$DI0,8kHz$

of the proposed, Long, and Ward arrays are 7.6, 5.9, and 5.7 dB, respectively; the

$DI350Hz,8kHz$

of the proposed array is 1.4 dB better than Son’s. This demonstrates the advantage of our proposed design.

### 5.2 Beamformer Weights Design

In this section, we discuss the properties of the proposed weight design. We report the degradation in performance of the proposed configuration when one of the design steps is omitted.

The TIT improved the

$DI0,8kHz$

by 0.6 dB and enabled maintaining the desired beamwidth at lower frequencies (an extra 410 Hz), proving the technique’s added value. The TIT gives larger weight to the outer sensors because of their greater distance between adjacent sensors (49). This enables maintaining the desired beamwidth at low frequencies by effectively widening the window.

Sampling the continuous Kaiser window at the element positions according to (48) led to a modest improvement of 0.05 dB over uniformly sampling the window according to (47). The

$DI0,8kHz$

improvement when using continuous values of W without limiting W to the set

$2xii=1L$

was negligible (0.01 dB), therefore, it is sufficient to limit the window support to this set. Although, if we were to limit the window support to only equal the support of the array, i.e., W ≡ 2xL, the beamwidth would be narrower than 15° at frequencies above 1.63 kHz. This shows the merit of setting a subset of the sensors to be active to widen the beampattern at higher frequencies effectively.

Last, we examine the performance when we limited the array configuration to be an M = 11 sensor ULA. To find the ULA’s optimal configuration, the iterative algorithm proposed in the previous section optimized one variable that represented the interelement spacing. The resulting interelement spacing of the ULA was 4.1 cm and its

$DI350Hz,8kHz$

was 7.2 dB—still better than Son, Long, and Ward.

### 5.3 Non-iterative Algorithm’s Performance

Here we display the properties of our proposed non-iterative algorithm for determining the sensor positions. We compare the iterative and non-iterative algorithms. For the iterative algorithm, we used the same array design as in the previous section. For the non-iterative algorithm, we initialized with a ULA of M0 = 5 sensors with an interelement spacing of 3.4 cm. We used this spacing because it maximized the DF at fH = 8 kHz. We initialized with M0 = 5 sensors because three cannot attain the desired beamwidth at 8 kHz without exhibiting high sidelobes. We set fH = 8 kHz, M = 11, βmax = 10, and φB = 15°.

We constructed two arrays using the non-iterative algorithm (detailed in Algorithm 2). For one array, we set βmin = 0, and for the other array, we choose βmin = 1.36 so that the minimum frequency for the desired beamwidth is the same as that obtained with the iterative algorithm (620 Hz). The resulting sensor positions are displayed in Figure 5. After the positions were chosen, the optimal window parameters were found using Algorithm 1 and are plotted in Figure 6. The beamwidth, WNG, and DF are plotted in Figure 7 as a function of frequency.

FIGURE 5. Configurations of the arrays: iterative (triangles), non-iterative with βmin = 0 (asterisks), and non-iterative with βmin = 1.36 (crosses).

FIGURE 6. Continuous Kaiser window parameters, (A) the window shape factor, and (B) the index of the window support from the set

$2xii=15$

, as a function of frequency: iterative (solid line), non-iterative with βmin =0 (dashed line), and non-iterative with βmin =1.36 (dash-dot line).

FIGURE 7. (A) Beamwidth, (B) WNG, and (C) narrowband directivity factor as a function of frequency: iterative (solid line), non-iterative with βmin = 0 (dashed line), and non-iterative with βmin = 1.36 (dash-dot line). The obtained

$W0,8kHz$

are respectively 7.8, 7.5, and 7.8 dB. The obtained

$DI0,8kHz$

are respectively 7.6, 7.3, and 7.4 dB.

For the array designed with βmin = 0, we see that as the frequency decreases, right before another sensor is activated, the window is as wide as possible (β = 0). For the array designed with βmin = 1.36, we see that the optimal window parameter, β, is sometimes less than 1.36. This is plausible because βmin was used only for choosing where to place the sensors. Afterward, Algorithm 1 may use any value of β.

We see that setting βmin > 0 resulted in improved directivity compared to βmin = 0. As expected, the iterative algorithm produced an array with a larger DI than the non-iterative algorithm. Nevertheless, the

$DI0,8kHz$

obtained with the non-iterative algorithm was only 0.2 dB smaller. This is beneficial considering the computational simplicity of the non-iterative algorithm compared to the iterative algorithm: no optimization toolbox is needed for the non-iterative algorithm.

### 5.4 Planar Array

This section compares our proposed techniques for constant-beamwidth beamforming for planar arrays to Son’s design (Son, 2021) using 99 sensors. The array layouts are displayed in Figure 8. The desired half-power beamwidths were arbitrarily chosen to be θX = 15° in the XZ plane and θY = 30° in the YZ plane.

FIGURE 8. (A) Configurations of the arrays: proposed planar (triangles) and proposed star (squares). (B) Configuration of Son’s array (circles).

To construct the planar array, we first design two linear arrays. To find the positions of the sensors, we use the non-iterative algorithm proposed in the previous section: we design a symmetric linear array with M = 11 sensors with a constant beamwidth of θX = 15° and another symmetric linear array with N = 9 sensors with a constant beamwidth of θY = 30°. We initialize with a ULA of M0 = 5 sensors with an interelement spacing of 3.4 cm. For the first array, we set βmin = 1.36, and for the second array, we set βmin = 1.00 so that the minimum frequencies that the beampatterns still maintain the desired beamwidths are both 620 Hz. Fewer sensors are needed in the second array because the desired beamwidth is wider. We use the non-iterative algorithm instead of the iterative algorithm because of its low computational complexity while attaining on-par performance. The beamformer weights for the two symmetric linear arrays are computed using Algorithm 1.

Based on the linear arrays’ sensor positions, we construct a symmetric planar array with a total of MN = 99 sensors. For this proposed layout, four different beamformer weight vectors are compared. For Kronecker product beamforming, the weights are given by the Kronecker product of the linear arrays’ weights (27). The beamformer weights for maximum WNG, tradeoff, and maximum directivity are given by (43), where we set α ∈ {1, 0.5, 0.01}, respectively.

To reproduce the array of Son (Son, 2021), we set the reference beampattern to equal

$sincθ−0.1477⁡cos2φ+0.4432,(53)$

which has a half-power beamwidth of 15° in the XZ plane and 30° in the YZ plane. From a uniform planar array with 201 × 201 sensors and an interelement spacing of 1 cm, an initial 30 sensors with the highest spanning energy were chosen. The remaining sensors were divided into 23 groups. From each group, three sensors closest to the group centroid were chosen. Finally, the sparse array was created with the 99 chosen sensors. Next, the weights were found. The mainlobe error tolerances were 0.006E1 for θ ∈ [0°, 15°]. The sidelobe error tolerances were 0.028E2 for θ ∈ (15°, 90°]. Here, E1 = 31 × 181 and E2 = 150, ×, 181 were the number of angles sampled (the azimuth and elevation resolutions were 2° and 0.5°, respectively). To ensure that the beamwidths were as desired, we added four constraints: that the error tolerances were 0.0001 for (φ, θ) ∈ {(0°, 7.5°), (90°, 15°), (180°, 7.5°), (270°, 15°)}. The resulting sparse array was not symmetric and had complex-valued weights.

With these weights, the beampatterns were constructed according to (5) and illustrated for a few frequencies in Figure 9. The beamwidths in the XZ and YZ planes of the proposed and Son’s beamformers are plotted in Figure 10 as a function of frequency. The beamwidths of the proposed beamformers were the same as the beamwidths of the linear arrays upon which they were built. For the Kronecker beamformer, this is due to property (28), and for the other beamformers, this is due to the constraint (37). Therefore, the proposed techniques maintained the desired beamwidth down to 620 Hz. Son’s array maintained the desired beamwidths at even lower frequencies. We see in Figure 10 that Son’s beamwidth was narrower than desired at some frequencies. This was because the beampattern dipped below

$1/2$

and then peaked back to

$1/2$

at

$12θY$

as constrained.

FIGURE 9. Planar arrrays’ beampatterns using the beamformers: (A) proposed Kronecker, (B) proposed tradeoff with α = 0.5, and (C) Son, for several frequencies: (1) 2 kHz, (2) 5 kHz, and (3) 8 kHz.

FIGURE 10. Beamwidths in the XZ and YZ planes as a function of frequency of Son’s and our proposed planar beamformers.

The WNG and DF are plotted in Figure 11 as a function of frequency. Son’s performance is only reported at frequencies where the WNG is at least 0 dB (f ≥ 890 Hz). For most frequencies, the WNGs of our design techniques were lower than Son’s. This was expected because Son’s objective was to maximize the WNG. We also see that the larger α was for the proposed techniques, the larger the WNG.

FIGURE 11. (A) WNG, and (B) narrowband directivity factor as a function of frequency: Kronecker (solid line), max WNG (dashed line), tradeoff (dash-dot line), max directivity (dotted line), and Son (thick solid line). The obtained

$W890Hz,8kHz$

are respectively 13.5, 15.0, 14.9, 13.0, and 13.1 dB. The obtained

$DI890Hz,8kHz$

are respectively 14.4, 15.2, 15.5, 15.8, and 6.1 dB.

From the DF, the DI was calculated using (16). The

$DI890Hz,8kHz$

of the proposed techniques were higher than achieved by Son. This demonstrates the advantage of the proposed techniques. We also see that the smaller α was for the proposed techniques, the larger the DI. It is important to note that the tradeoff beamformer has greater WNG and DI than the Kronecker beamformer. The Kronecker beamformer’s performance is inferior because it obtains the desired beamwidths without solving an optimization problem.

A major advantage of our technique over Son’s is the computational complexity. We report the running times on an Intel Core™ i7-8550U CPU with 8 GB RAM with the code implemented in MATLAB. For the tradeoff beamformer, it took 3 s to find the sensor positions for the two linear arrays using the non-iterative algorithm, 17 s to compute the weights with Algorithm 1 (the resolution of β was 0.01), and 33 s to calculate the beamformer weights (for all 801 frequencies). In total: less than 1 min. On the other hand, for Son, it took 279 s to choose the positions of the sensors and 286 s to calculate the beamformer weights for one frequency. To calculate for all 801 frequencies, Son’s technique took more than 1,000 times longer than ours.

Last, we examine the performance when limiting the array configuration to an 11 × 11 planar array with uniform spacing. We construct two ULAs (each with M = 11 sensors) with an interelement spacing of 4.2 cm to find the sensor positions. We use this spacing because it results in the highest

$DI0,8kHz$

while maintaining a beamwidth at least as wide as desired for all frequencies. We use Algorithm 1 to find the weights for both ULAs: for one the desired beamwidth was θX = 15° and for the other the desired beamwidth was θY = 30°. From these two ULAs, we construct a planar array (with an interelement spacing of 4.2 cm) where the weights are given by the tradeoff beamformer (43) with α = 0.01 to maximize the directivity. Its

$DI0,8kHz$

was 11.2 dB—lower by at least 1.3 dB than our proposed nonuniformly spaced planar array beamformers, even though it had 121 sensors. This is similar to what we have seen for the linear arrays: the nonuniformly spaced array achieves greater DI.

### 5.5 Planar Array With Missing Sensors

In this section, we simulate a planar array with sensors only at some of the grid positions. We use the same planar array design from the previous section. However, we place sensors at only S = 39 out of the MN = 99 grid positions. In addition, we choose to place the sensors in a way that resembles a star. The selected grid positions are displayed in Figure 8 and called the star array. Because only a subset of the grid positions is used, the Kronecker product beamformer is invalid here. Therefore, we only compare three different beamformers: the beamformers given by (43) with α ∈ {1, 0.5, 0.01}, for maximum WNG, tradeoff, and maximum directivity, respectively. With the beamformers’ weights, the beampatterns were constructed according to (5) and are illustrated for a few frequencies in Figure 12. The beamwidths are shown in Figure 10. Even though we only placed sensors at a subset of the grid positions, the beamformers still maintain the desired beamwidth down to 620 Hz. The WNG and DF are plotted in Figure 13 as a function of frequency. As expected, the larger is α, the higher is the WNG and the lower is the DI. The

$DI890Hz,8kHz$

of the S = 39 sensor star array is still higher than achieved by Son’s sparse 99 sensor array. This demonstrates the advantage of the proposed technique.

FIGURE 12. Star array’s beampatterns using the beamformers: (A) max WNG, (B) tradeoff, and (C) max directivity, for several frequencies: (1) 2 kHz, (2) 5 kHz, and (3) 8 kHz.

FIGURE 13. (A) WNG, and (B) narrowband directivity factor as a function of frequency for the star array beamformers: max WNG (dashed line), tradeoff (dash-dot line), and max directivity (dotted line). The obtained

$W890Hz,8kHz$

are respectively 10.4, 10.3, and 9.2 dB. The obtained

$DI890Hz,8kHz$

are respectively 10.5, 10.8, and 11.0 dB.

We compare the performance of the star array to a 7 × 7 planar array with uniform spacing. To find the positions of the sensors, we construct two ULAs (each with M = 7 sensors) with an interelement spacing of 4.3 cm. We use this spacing because it results in the highest

$DI0,8kHz$

while maintaining beamwidths that are at least as wide as desired for all frequencies. We use Algorithm 1 to find the weights for both ULAs: for one the desired beamwidth is θX = 15°, and for the other the desired beamwidth is θY = 30°. From these two ULAs, we construct a uniform planar array (with an interelement spacing of 4.3 cm) where the weights are given by the tradeoff beamformer (43) with α = 0.01 to maximize the directivity. Its

$DI0,8kHz$

is 8.7 dB—lower by at least 0.4 dB than our nonuniformly spaced star array beamformers, even though it has 49 sensors. This once again shows the value of using a nonuniformly spaced design.

## 6 Conclusion

We have proposed planar beamformers that yield different beamwidths in the XZ and YZ planes by designing constant-beamwidth linear arrays. The planar beamwidths equal the linear array beamwidths for all frequencies. The first technique involves utilizing Kronecker product properties, whereby no matrix inversion is needed to find the weights. Our second technique provides a closed-form solution that allows a tradeoff between WNG and directivity. The second technique is applicable even when only a subset of the sensors is used. Further research is needed to choose the optimal subset of the sensors.

As a building block for the planar beamformer, we proposed a method to find optimal sensors placements in a symmetric nonuniform linear array and optimal beamformer to maximize the DI. The design involves sampling the continuous Kaiser window, varying the window support, and using the TIT. Sampling the continuous Kaiser window at the sensor positions is appropriate given the nonuniform interelement spacing of the symmetric array. Varying the window support and the window shape factor enables more control over the effective width of the window. Using the TIT allows for maintaining the desired beamwidth at low frequencies. Experimental results have shown that the proposed techniques yield higher directivity and involve lower computational complexity than the state-of-the-art. Additionally, the resulting symmetric nonuniform arrays outperform uniformly-spaced arrays. Future work will focus on generalizing the proposed methods to three-dimensional arrays and addressing the problem of constant-beamwidth steering with nonuniform arrays.

## Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

## Author Contributions

AF contributed to methodology, investigation, formal analysis, software, validation, and writing—original draft. IC contributed to conceptualization, methodology, supervision, and writing—review and editing.

## Funding

This work was supported by the Pazy Research Foundation.

## Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

## Acknowledgments

The authors would like to thank Stem Audio for providing equipment and technical guidance.