Design of Differential Loudspeaker Line Array for Steerable Frequency-Invariant Beamforming

1. Introduction

Beamforming is a fundamental technique in array signal processing and has garnered significant attention in research [1,2,3,4,5,6]. In microphone arrays, beamforming suppresses spatial interferences and improves the output signal-to-noise ratio [7,8,9,10,11,12]. In loudspeaker arrays, beamforming creates highly directional beam patterns, reducing room reflections and enhancing the desired sound's delivery [13,14,15,16,17,18]. Array systems can be classified into two categories: additive arrays and differential arrays. Additive arrays generate directional beam patterns through the synchronize-and-add principle [19]. However, they cannot achieve high directivity at low frequencies with limited array aperture. In contrast, differential arrays, which respond to the spatial derivatives of the sound field, can generate narrow beam patterns and achieve high directional gains even with small apertures [20].

Differential arrays offer advantages in compact size, frequency-invariant beampatterns, and high spatial directivity, and have been extensively studied in microphone array applications over the last two decades. Differential microphone arrays can be implemented in various geometric shapes, including line-shaped [20,21,22,23,24,25], planar [26,27,28,29,30,31], and volumetric designs [32,33,34,35]. Among these, the line-shaped differential microphone array (LDMA) is particularly well-studied due to its ease of integration with electronic devices. While most research focuses on directing the main lobe in the endfire direction, this is often unsuitable for applications like smart TVs or tablets, where beam steering is needed, e.g., to accurately pick up voice commands from users speaking from different positions.

Efforts have been made to design first- and second-order steerable LDMAs [36,37], including the development of general conditions for steerable differential beamformers [38]. However, the ideal beam pattern can become ineffective due to improper null position selection. To address this, extremum and interference suppression constraints have been used in designing steerable LDMAs [39]. Nonetheless, selecting the optimal trade-off parameter between the directivity factor (DF) and interference suppression energy remains challenging. Recently, a novel method was proposed using both omnidirectional and directional microphones to design fully steerable LDMAs, achieving steering- and frequency-invariant beampatterns with a linear super array (LSA) [40,41]. However, mismatches between omnidirectional and directional microphones can impact beamformer performance in practical applications.

Differential beamformers have also been applied to linear loudspeaker arrays to produce highly directional patterns [42,43,44,45,46,47,48,49]. For example, a three-element line array with a second-order differential broadside beam has been studied to create a near-field sound zone in car cabins [44]. For higher-order broadside differential beampatterns, a null-constrained method has been proposed as a convenient approach for designing differential beamformers [45]. However, its frequency-invariant beampattern is limited, which may be undesirable for certain applications. Our recent research extends the null-constrained method using a series expansion approach for broadside linear loudspeaker arrays, enhancing the beamformer's robustness and better preserving the frequency-invariant beam pattern [48,49]. In the study of beam steering for differential loudspeaker arrays (DLAs), A linear loudspeaker array composed of miniature omnidirectional speakers and dipole speakers is used to design the steerable first-order DLAs [47]. Our recent study develops steerable frequency-invariant differential beampatterns using a circular loudspeaker array with a rigid baffle [46]. However, miniature loudspeakers and circular arrays are not widely applicable, such as in public broadcasting and outdoor sound reinforcement scenarios. Therefore, there is a need to explore beam steering with loudspeaker line arrays.

In this paper, we propose a method for designing steerable differential beamformers using a loudspeaker line array. The key contributions include: 1) calculating the target steerable differential beampattern based on the desired direction, main lobe width, and beampattern order; 2) formulating the desired differential beampattern in the modal domain for representation; and 3) developing a multi-constraint optimization problem to design a steerable frequency-invariant beamformer using the Jacobi-Anger expansion, while minimizing the mean square error (MSE) between the synthesized and ideal beampatterns. Both simulation and experimental studies are provided. Although this work focuses on differential loudspeaker line arrays, the proposed method is also applicable to differential microphone arrays, making it suitable for designing LDMAs with steerable differential beampatterns for high-quality acoustic signal acquisition.

The paper is organized as follows: Section 2 introduces the signal model, problem formulation, and key definitions. Section 3 details the proposed method for designing steerable frequency-invariant beampatterns with a linear array. Section 4 presents simulations with design examples and discusses the impacts of constraint parameters and element mismatch, as well as a performance comparison with existing methods. Section 5 provides experimental results that align with the simulations, validating the proposed method. Finally, Section 6 concludes the paper.

2. Problem Formulation

Figure 1 illustrates a loudspeaker array designed to radiate a steerable differential beampattern towards a listener. The line array consisting of L loudspeakers with an interelement spacing $\delta$ lies on the x axis. The loudspeaker array is centers at the origin of the coordinate system. Each loudspeaker is modelled as omnidirectional point source, and the far-field sound pressure at a listener position $(r,\theta )$ generated by the loudspeaker array can be expressed as

$$p(k,r,\theta )\approx {\displaystyle \frac{e^{ikr}}{4\pi r}}{\displaystyle \sum _{l=1}^L{w}_l^{*}(k)e^{-ik{x}_l\cos \theta }}$$

(1)

where $i$ is the imaginary unit, with $i^2=-1$, $k=2\pi f/c$ is the wave number, f denotes the frequency, c is the speed of sound in air, the superscript ${\left(·\right)}^{*}$ is the complex-conjugate operator,$w_l(k)$ denotes the weight of the l-th loudspeaker at $(x_l,0)$, where$x_l=-(L+1)\delta /2+l\delta$, $l$ is the loudspeaker index, l = 1, 2, …, L, r is the distance from the origin of the coordinate system to the listener, the angle $\theta$ is defined with the respect to the positive x axis.

In the far-field, the directivity of the loudspeaker array does not change with the distance r. To evaluate the directivity of radiation, the normalized far-field radiation pattern is

$$B(k,\theta )=p(k,r,\theta )/(e^{ikr}/4\pi r)={\displaystyle \sum _{l=1}^L{w}_l^{*}(k)e^{-ik{x}_l\cos \theta }},$$

(2)

which can be written in a vector form

$$B(k,\theta )={\mathrm{w}}^{\mathrm{H}}(k)\mathrm{g}(k,\theta ),$$

(3)

where

$$\mathrm{w}(k)={[w_1(k),\mathrm{...},w_L(k)]}^{\mathrm{T}},$$

(4)

$$\mathrm{g}(k,\theta )={[e^{-ik{x}_1\cos \theta },\cdots ,e^{-ik{x}_L\cos \theta }]}^{\mathrm{T}},$$

(5)

where the superscripts ${\left(·\right)}^{\mathrm{H}}$ and ${\left(·\right)}^{\mathrm{T}}$ are the conjugate-transpose operator and the transpose operator respectively, $\mathrm{w}$ is the beamforming filter to be designed to achieve a steerable differential radiation pattern.

Apart from the normalized far-field radiation pattern, there are two commonly used metrics to evaluate the performance of a beamformer: white noise gain (WNG) and directivity factor (DF). The WNG quantifies the robustness of the beamformer against the white noise that

$$\mathrm{WNG}(\mathrm{k})={\displaystyle \frac{{\left|B(k,{\theta }_s)\right|}^2}{{\mathrm{w}}^H(k)\mathrm{w}(k)}},$$

(6)

where ${\theta }_s$ is the desired radiation direction. For the loudspeaker array, the WNG also represents the radiation efficiency of the array.

The DF evaluates the directional characteristics of the beamformer. The 2-dimensional DF is defined as the ratio of the sound power radiated in the desired direction to the average sound power across the half-plane in front of the loudspeaker array, and is expressed as

$$DF(\mathrm{k})={\displaystyle \frac{\pi {\left|\mathrm{B}(\mathrm{k},{\theta }_s)\right|}^2}{{\displaystyle {\int }_0^{\pi }{\left|\mathrm{B}(\mathrm{k},\theta )\right|}^{2}d\theta }}}={\displaystyle \frac{{\left|\mathrm{B}(\mathrm{k},{\theta }_s)\right|}^2}{{\mathrm{w}}^H(\mathrm{k})\mathsf{\Gamma }(\mathrm{k})\mathrm{w}(\mathrm{k})}},$$

(7)

where

$$\mathsf{\Gamma }(k)={\displaystyle \frac{1}{\pi }}{\displaystyle {\int }_0^{\pi }{\mathrm{g}}^{\mathrm{H}}(k,\theta )}\mathrm{g}(k,\theta )d\theta$$

(8)

is a $L\times L$ square matrix, whose elements are ${\left[\mathsf{\Gamma }\right]}_{m,n}(k)=\begin{array}{cc}{J}_0[k(i-j)\delta ],& m,n\in \{1,2,\mathrm{...},L\}\end{array}$. Here, $J_0(\cdot )$ is the zero order Bessel function of the first kind. For brevity, we will omit the dependence on k in the following text.

3. Methods

The objective of the proposed method is to form a directional steerable loudspeaker array such that the energy radiated by the array is mainly concentrated in the main lobe area while minimizing the output energy of the sidelobe area. It is achieved in three steps: (1) calculating the target radiation pattern according to the radiation steering angle, main lobe area and sidelobe area; (2) formulating the target radiation pattern in the modal domain; (3) designing the steerable differential beamformer using the modal matching with Jacobi-Anger expansion.

3.1. Target Radiation Pattern with the Mainlobe Steering

The ideal N-th order differential beampattern for a differential microphone array is

$${\tilde{B}}^{(N)}(\theta )={\displaystyle \sum _{n=0}^N{\alpha }_{N,n}{\cos }^n\theta }$$

(9)

• which assumes the microphone spacing is much smaller than the wavelength, and
• ${\alpha }_{N,n},n=0,1,\mathrm{...},N$ are real coefficients. To allow mainlobe steering to the desired direction
• ${\theta }_s$, two conditions should be satisfied, that ${\tilde{B}}^{(N)}({\theta }_s)=1$ and $\partial {\tilde{B}}^{(N)}({\theta }_s)/\partial {\theta }_s=0$. In the vector form,

$${\mathbf{A}}^{\mathrm{T}}{\mathbf{\alpha }}_N=\mathbf{b}$$

(10)

Where

$$\mathbf{A}={\left[\begin{array}{cccc}1& \cos {\theta }_s& \dots & {\cos }^N{\theta }_s\\ 0& \sin {\theta }_s& \dots & N{\cos }^{N-1}{\theta }_s\sin {\theta }_s\end{array}\right]}^{\mathrm{T}},$$

(11)

$${\alpha }_N={\left[\begin{array}{cccc}{\alpha }_{N,0}& {\alpha }_{N,1}& \dots & {\alpha }_{N,N}\end{array}\right]}^{\mathrm{T}},\mathrm{b}={\left[\begin{array}{cc}1& 0\end{array}\right]}^{\mathrm{T}}.$$

(12)

Also, (9) in the vector form is

$${\tilde{B}}^{(N)}(\theta )={\alpha }_N^{\mathrm{T}}{c}_N,$$

(13)

where

$$c_N={[\begin{array}{cccc}1& \cos \theta & \dots & {\cos }^N\theta \end{array}]}^{\mathrm{T}}.$$

(14)

The corresponding beam power can be expressed as

$${\left|{\tilde{B}}^{(N)}(\theta )\right|}^2={\alpha }_N^{\mathrm{T}}{c}_N{c}_N^{\mathrm{T}}{\alpha }_N={\alpha }_N^{\mathrm{T}}\mathrm{C}{\alpha }_N,$$

(15)

where

$$C=c_N{c}_N^{\mathrm{T}}.$$

(16)

Due to the acoustic reciprocity principle, our proposed method applies equation (13) as the ideal N-th order differential radiation pattern for a differential loudspeaker array, enabling precise steering of sound propagation towards the desired direction ${\theta }_s$. Additionally, it is necessary to minimize the radiation of sound energy in the sidelobe area, defined as

$${\mathrm{\Theta }}_{SL}\in \left[\begin{array}{cc}{0}^{\circ }，& {\theta }_s-\Delta \end{array}/2\right]\cup \left[\begin{array}{cc}{\theta }_s+\Delta /2,& {180}^{\circ }\end{array}\right],$$

(17)

where $\Delta$ is the mainlobe width. Therefore, once the desired direction and mainlobe width are determined, the optimal solution ${\alpha }_N$ can be obtained by solving the following optimization problem

$$\begin{array}{cc}\underset{{\mathsf{\alpha }}_N}{\min }{\left|{\tilde{B}}^{(N)}({\Theta }_{SL})\right|}^2& s.t.\end{array}\begin{array}{c}{\mathbf{A}}^{\mathrm{T}}{\mathbf{\alpha }}_N=\mathbf{b}\end{array}.$$

(18)

This can be solved using the CVX toolbox [50]. Moreover, the selection of the mainlobe width and the desired direction ${\theta }_s$ must satisfy

$$\Delta \le 2\cdot \min ({\theta }_s,{180}^{\circ }-{\theta }_s)$$

(19)

The optimal target radiation pattern can be obtained by substituting the resultant into (13). Figure 2 provides examples of target radiation patterns designed using the above method at different desired angles. This indicates that the proposed method can design effective radiation patterns for different orders in various desired directions. Figure 2 also presents the target radiation patterns of the same order and desired directions designed using the method stated in Ref [39]. Yu's method finds the optimal vector ${\mathbf{\alpha }}_N$ by maximizing the directivity factor of the beampattern under the constraints of the distortionless constraint and the extremum constraint in the desired direction [as (10)]. This method needs a regularization parameter to trade off the radiation power in the end-fire direction and the spatial average of the radiated power over the entire space. The regularization parameter is set to 0.15, as given in Ref [39] and applied in the following simulation.

Figure 2a–Figure 2d present the comparison results of the fourth-order target differential radiation patterns at desired angles of 30°, 60°, 90°, and 135°, using Yu's method and the proposed method, respectively. Though both methods can design effective fourth-order differential radiation patterns at the desired angles, the proposed method has smaller sidelobes (excluding that in the end-fire direction). This is because Yu's method designs the target radiation pattern by adjusting the regularization parameter to balance the array's energy radiated in the end-fire direction and the average energy radiated into the entire space. Since this method does not account for the energy distribution of sidelobes other than the end-fire direction, it can result in excessively large sidelobes. Thus, the regularization parameter selection is required for different desired directions. However, selecting the appropriate regularization parameter is quite challenging, which hinders the practical application of this method. In contrast, the proposed method designs the target radiation patterns by minimizing the total energy in the sidelobe area, which includes the end-fire direction. As a result, the target beampattern radiates less energy into the sidelobe area. Additionally, once the desired direction and main lobe width are determined, we can design the optimal radiation pattern without adjusting other parameters, facilitating practical applications.

Figure 2e-Figure 2f further compare the fifth-order target radiation patterns by Yu’s and the proposed methods, with the desired directions of ${\theta }_s={45}^{\circ }$ and ${\theta }_s={90}^{\circ }$, respectively. The above findings hold. In addition, when the desired direction is in the broadside direction, ${\theta }_s={90}^{\circ }$, the 5th-order target radiation pattern (Figure 2f) is the same as the 4th-order target pattern (Figure 2c). This occurs because the target radiation pattern in the broadside direction is symmetric about the y axis. As a result, the odd terms of the optimal vector ${\mathbf{\alpha }}_N$ are zero, leading to pairs of nulls in the target beampattern that are symmetrically distributed about the y axis. Therefore, when designing N-th order broadside radiation patterns (where N is odd), the process automatically reduces to designing an (N−1)-th order pattern.

3.2. Formulating the Target Radiation Pattern in the Modal Domain

Using the N-th order target radiation pattern obtained in Section 3.1, we will transform the ideal radiation pattern from equation (9) into the modal domain, depending on the parity of N.

3.2.1.N. is Even

Let $N=2\tilde{N}$ and ${\alpha }_{N,-1}=0$, (9) will be

$${\tilde{B}}^{(N)}(\theta )={\displaystyle \sum _{n=0}^N{\alpha }_{N,n}{\cos }^n\theta }={\displaystyle \sum _{n=0}^{\tilde{N}}{\alpha }_{N,2n}{\cos }^{2n}\theta +}{\displaystyle \sum _{n=0}^{\tilde{N}}{\alpha }_{N,2n-1}{\cos }^{2n-1}\theta }.$$

(20)

According to the following relations [51]:

$${\cos }^{2n}\theta ={\displaystyle \frac{1}{2^{2n}}}\left\{{\displaystyle \sum _{t=0}^{n-1}2\left(\begin{array}{c}2n\\ t\end{array}\right)\cos 2\left(n-t\right)\theta +}\right.\left.\left(\begin{array}{c}2n\\ n\end{array}\right)\right\},$$

(21)

$${\cos }^{2n-1}\theta ={\displaystyle \frac{1}{2^{2n-2}}}\left\{{\displaystyle \sum _{t=0}^{n-1}\left(\begin{array}{c}2n-1\\ t\end{array}\right)\cos \left[2\left(n-t\right)-1\right]\theta }\right\},$$

(22)

where $\left(\begin{array}{c}·\\ ·\end{array}\right)$ is combinations. Combining the Euler’s formula, i.e.,

$$\cos (n\theta )={\displaystyle \frac{e^{in\theta }+e^{-in\theta }}{2}},$$

(23)

then (20) can be expressed as

$${\tilde{B}}^{(N)}(\theta )={\displaystyle \sum _{n=-N}^N{\gamma }_n{e}^{in\theta }}.$$

(24)

where

$${\gamma }_n=\left\{\begin{array}{c}\begin{array}{cc}{\displaystyle \sum _{p=\frac{\left|n\right|+1}{2}}^{\tilde{N}}{\alpha }_{N,2p-1}{\tilde{\eta }}_{\frac{n+\mathrm{sgn}(n)}{2}}(p)}& n=\pm 1,\pm 3,\mathrm{...},\pm (2\tilde{N}-1)\begin{array}{c},\end{array}\end{array}\\ \begin{array}{cc}{\displaystyle \sum _{p=\left|n\right|/2}^{\tilde{N}}{\alpha }_{N,2p}{\eta }_{n/2}(p)}& n=0,\pm 2,\pm 4\end{array},\mathrm{...},\pm 2\tilde{N}\begin{array}{c},\end{array}\end{array}\right.$$

(25)

$${\tilde{\eta }}_t(n)=\left\{\begin{array}{c}\begin{array}{cc}0& t=0\begin{array}{c},\end{array}\end{array}\\ \begin{array}{cc}{\displaystyle \frac{1}{2^{2n-1}}}\left(\begin{array}{c}2n-1\\ n-\left|t\right|\end{array}\right)& t=\pm 1,\pm 2,\mathrm{...},\pm n\begin{array}{c},\end{array}\end{array}\end{array}\right.$$

(26)

$$\mathrm{sgn}(t)=\left\{\begin{array}{c}\begin{array}{cc}1& t>0\begin{array}{c},\end{array}\end{array}\\ \begin{array}{cc}0& t=0\begin{array}{c},\end{array}\end{array}\\ \begin{array}{cc}-1& t<0\begin{array}{c},\end{array}\end{array}\end{array}\right.$$

(27)

$${\eta }_t(n)={\displaystyle \frac{1}{2^{2n}}}\left(\begin{array}{c}2n\\ n-\left|t\right|\end{array}\right)\begin{array}{c}t=0,\pm 1,\mathrm{...},\pm n\end{array}.$$

(28)

where $\left|·\right|$ represents absolute value.

3.2.2.N. is odd

Let $N\text{'}=(N+1)/2$, and

$${\alpha }_{2N\text{'},n}=\left\{\begin{array}{c}\begin{array}{cc}{\alpha }_{N,n}& n=0,1,\mathrm{...},N\begin{array}{c},\end{array}\end{array}\\ 0\begin{array}{c}n=N+1\begin{array}{c},\end{array}\end{array}\end{array}\right.$$

(29)

then we can express the odd (N-th) order target radiation pattern as a form of the even ($2N\text{'}$) order target beampattern:

$${\tilde{B}}^{(N)}(\theta )={\tilde{B}}^{(2N\text{'})}(\theta )={\displaystyle \sum _{n=0}^{2N\text{'}}{\alpha }_{2N\text{'},n}{\cos }^n\theta }={\displaystyle \sum _{n=0}^{N\text{'}}{\alpha }_{2N\text{'},2n}{\cos }^{2n}\theta +}{\displaystyle \sum _{n=0}^{N\text{'}}{\alpha }_{2N\text{'},2n-1}{\cos }^{2n-1}\theta },$$

(30)

Similar to the derivation steps in Section 3.2.1, we can formulate (30) in the modal domain with a symmetric form

$${\tilde{B}}^{(N)}(\theta )={\tilde{B}}^{(2N\text{'})}(\theta )={\displaystyle \sum _{n=-2N\text{'}}^{2N\text{'}}{\gamma }_n{e}^{in\theta }}$$

(31)

where

$${\gamma }_n=\left\{\begin{array}{c}\begin{array}{cc}{\displaystyle \sum _{p=\frac{\left|n\right|+1}{2}}^{N\text{'}}{\alpha }_{2N\text{'},2p-1}{\tilde{\eta }}_{\frac{n+\mathrm{sgn}(n)}{2}}(p)}& n=\pm 1,\pm 3,\mathrm{...},\pm (2N\text{'}-1)\end{array}\\ \begin{array}{cc}{\displaystyle \sum _{p=\left|n\right|/2}^{N\text{'}}{\alpha }_{2N\text{'},2p}{\eta }_{n/2}(p)}& n=0,\pm 2,\pm 4\end{array},\mathrm{...},\pm 2N\text{'}\end{array}\right.,$$

(32)

where ${\tilde{\eta }}_t(n)$, $\mathrm{sgn}(t)$ and ${\eta }_t(n)$ are defined in (26), (27) and (28), respectively.

3.3. Beamformer Design

The objective is to determine the loudspeaker weighting so that the array’s radiation pattern closely approximates the target radiation pattern. Based on (2) and (24), we use the Jacobi-Anger series expansion method to accomplish this task.

3.3.1. Modal Matching Method with Maximum WNG

In the series expansion method, the Jacobi–Anger expansion is used to decompose the resulting radiation pattern into a linear combination of circular harmonics. The Jacobi-Anger expansion is

$$e^{-i\sigma \cos \theta }={\displaystyle \sum _{n=-\infty }^{+\infty }{\beta }_n(\sigma )e^{in\theta }},$$

(33)

where ${\beta }_n(\sigma )={(-i)}^n{J}_n(\sigma )$, and $J_n(\cdot )$ is the nth-order Bessel function of the first kind. Substituting (33) into (2), we obtain

$$B(\theta )={\displaystyle \sum _{l=1}^L{w}_l^{*}{\displaystyle \sum _{n=-\infty }^{+\infty }{\beta }_n(k{x}_l)e^{in\theta }}}.$$

(34)

In order to obtain a Nth-order target radiation pattern, the infinite summation in (34) is truncated to the order N, that

$$B(\theta )\approx {\displaystyle \sum _{n=-N}^N{e}^{in\theta }}{\displaystyle \sum _{l=1}^L{\beta }_n(k{x}_l)w_l^{*}}.$$

(35)

When the normalized far-field radiation pattern in (35) is in consistent with the target radiation pattern in (24), there is

$${\displaystyle \sum _{l=1}^L{\beta }_n(k{x}_l)w_l^{*}}={\gamma }_n\begin{array}{c}n=0,\pm 1,\mathrm{...},\pm N.\end{array}.$$

(36)

For the design of a steerable differential beamformer, the distortionless constraint in the desired direction is needed

$${\mathbf{w}}^{\mathrm{H}}\mathbf{g}({\theta }_s)=1.$$

(37)

Combing (36) and (37) yields

$$\mathbf{\Phi }\mathbf{w}=\mathsf{\nu }.$$

(38)

where

$$\mathbf{\Phi }={\left[\begin{array}{ccc}{\mathbf{\beta }}_{-N}& \dots & \begin{array}{cc}{\mathbf{\beta }}_N& \mathbf{g}({\theta }_s)\end{array}\end{array}\right]}^{\mathrm{H}},$$

(39)

$${\mathbf{\beta }}_n={[\begin{array}{ccc}{\beta }_n(k{x}_1)& \cdots & {\beta }_n(k{x}_L)\end{array}]}^{\mathrm{T}},$$

(40)

and

$$\mathsf{\nu }={\left[\begin{array}{cc}\begin{array}{ccc}{\gamma }_{-N}& \mathrm{...}& \begin{array}{cc}{\gamma }_0& \begin{array}{cc}\mathrm{...}& {\gamma }_N\end{array}\end{array}\end{array}& 1\end{array}\right]}^{\mathrm{T}}.$$

(41)

Using the symmetric property of the Bessel function, $i^n{J}_n(\cdot )=i^{-n}{J}_{-n}(\cdot )$, (38) can be simplified to

$$\tilde{\mathbf{\Phi }}\mathbf{w}=\tilde{\mathbf{\nu }},$$

(42)

where

$$\tilde{\mathbf{\Phi }}={\left[\begin{array}{ccc}{\mathbf{\beta }}_0& \cdots & \begin{array}{cc}{\mathbf{\beta }}_N& \mathbf{g}({\theta }_s)\end{array}\end{array}\right]}^{\mathrm{H}}$$

(43)

is $\left(N+2\right)\times L$ full matrix and

$$\tilde{\mathbf{\nu }}={[\begin{array}{cc}\begin{array}{ccc}{\gamma }_0& \mathrm{...}& {\gamma }_N\end{array}& 1\end{array}]}^{\mathrm{T}}.$$

(44)

With $L>\left(N+2\right)$ loudspeakers to generate Nth-order steerable target radiation pattern, one approach to design the beamformer is to maximize WNG with the constraints of (42). Due to (37) holds, The WNG, defined in (6), can be written as $\mathrm{WNG}=1/{\mathbf{w}}^H\mathbf{w}$.

The optimization problem can be formulated as

$$\begin{array}{cc}\underset{\mathbf{w}}{\min }{\mathbf{w}}^{\mathrm{H}}\mathbf{w}& \begin{array}{cc}\mathrm{s}.\mathrm{t}.& \tilde{\mathbf{\Phi }}\mathbf{w}=\tilde{\mathbf{\nu }}\end{array}\end{array}.$$

(45)

The solution is

$$\mathbf{w}={\tilde{\mathbf{\Phi }}}^H{(\tilde{\mathbf{\Phi }}{\tilde{\mathbf{\Phi }}}^H)}^{-1}\tilde{\mathbf{\nu }}.$$

(46)

3.3.2. Modal Matching Method with WNG Constraint

Since the resulting beamformer (46) is obtained with some approximations, the mean square error (MSE) is adopted to evaluate the accuracy of these approximations. The MSE of the approximations of the beamformer to the Nth-order target differential radiation pattern is defined as

$$d(\mathbf{w})={\displaystyle \frac{1}{\pi }}{\displaystyle {\int }_0^{\pi }{\left|{\mathbf{w}}^{\mathrm{H}}\mathbf{g}(\theta )-{\tilde{B}}^{(N)}(\theta )\right|}^2}d\theta .$$

(47)

Substituting (9) into (47), the MSE can be written as a quadratic function,

$$d(\mathbf{w})={\mathbf{w}}^{\mathrm{H}}\mathsf{\Gamma }\mathbf{w}-{\mathbf{w}}^{\mathrm{H}}\mathrm{q}-{\mathrm{q}}^{\mathrm{H}}\mathbf{w}+\xi$$

(48)

where

$$\mathbf{q}={\displaystyle \frac{1}{\pi }}{\displaystyle {\int }_0^{\pi }\mathbf{g}(\theta ){\tilde{B}}^{(N)}(\theta )d\theta }=\mathrm{Q}{\mathbf{\alpha }}_N,$$

(49)

$$\mathbf{Q}={\displaystyle \frac{1}{\pi }}{\displaystyle {\int }_0^{\pi }\mathbf{g}(\theta ){\mathrm{c}}_N^{\mathrm{T}}}d\theta ,$$

(50)

$$\xi ={\displaystyle {\int }_0^{\pi }{\left|{\tilde{B}}^{(N)}(\theta )\right|}^2}d\theta ={\alpha }_N^{\mathrm{T}}\tilde{\mathrm{C}}{\alpha }_N,$$

(51)

$$\tilde{\mathrm{C}}={\displaystyle \frac{1}{\pi }}{\displaystyle {\int }_0^{\pi }{c}_N{c}_N^{T}d\theta }.$$

(52)

The matrix $\mathbf{\Gamma }$ is the array spatial correlation matrix, which is defined in (8). $\mathbf{Q}$ is an complex matrix. Its (m, n + 1)th (m = 1, 2, …, L; n = 0, 1, …, N) element is

$$Q_{m,n+1}={\displaystyle \frac{1}{2^{n-1}}}{\displaystyle \sum _{p=0}^{⌊n/2⌋}\left(\begin{array}{c}n\\ p\end{array}\right)i^{n-2p}{J}_{n-2p}(-k{x}_m)}-{\displaystyle \frac{1}{2^n}}{\delta }_{0,n-2⌊n/2⌋}\left(\begin{array}{c}n\\ ⌊n/2⌋\end{array}\right)i^{n-2⌊n/2⌋}{J}_{n-2⌊n/2⌋}(-k{x}_m),$$

(53)

where $⌊\cdot ⌋$ represents the floor function applied to a real number, $J_n(\cdot )$ is the nth-order Bessel function of the first kind and ${\delta }_{m,n}$ represents the Kronecker delta. $\tilde{\mathbf{C}}$ is an $(L+1)\times (L+1)$ real matrix, whose (m + 1, n + 1)th element is

$${\tilde{\mathrm{C}}}_{m+1,n+1}={\displaystyle \frac{1}{\pi }}{\displaystyle {\int }_0^{\pi }{\cos }^{m+n}\theta d\theta =\frac{(1+{(-1)}^{(m+n)})\Gamma \left((1+m+n)/2\right)}{2\sqrt{\pi }\Gamma \left((2+m+n)/2\right)}}\begin{array}{c}\left(m=0,1,\mathrm{...},N;n=0,1,\mathrm{...},N\right)\end{array}$$

(54)

where $\mathrm{\Gamma }(\cdot )$ is the Euler’s Gamma function.

After calculating the optimal coefficients of the target differential beampattern (9) according to the actual requirements, the objective for designing the beamformer are twofold. On the one hand, we aim to design the resulting beampattern approaches the target one as much as possible. On the other hand, we hope the system is robust enough to against unknown channel mismatch in real applications. The MSE, defined in (47), measures how closely the synthesized beampattern approximates the target beampattern. The WNG, defined in (6), is a metric for measuring the robustness of the designed beamformer. Therefore, attention should be paid to these two metrics when designing the beamformer.

The solution in (46) gives the maximum WNG. The higher the value of WNG, the more robust is the radiation. However, it does not consider the MSE in the design of the beamformer. In practice, the WNG greater than a suitable value should be ensured. It indicates the beamformer is robust enough for usage. In addition, the degrees of freedom provided by the weighting functions can be used to minimize MSE. The optimization problem can be described as follows

$$\begin{array}{l}\begin{array}{cc}\underset{\mathbf{w}}{\min }& {\mathbf{w}}^{\mathrm{H}}\mathbf{\Gamma }\mathbf{w}-{\mathbf{w}}^{\mathrm{H}}\mathbf{q}-{\mathbf{q}}^{\mathrm{H}}\mathbf{w}\end{array}\\ \begin{array}{c}s.t.\end{array}\tilde{\mathbf{\Phi }}\mathit{w}=\tilde{\mathbf{\nu }},\\ \begin{array}{c}{\mathbf{w}}^{\mathrm{H}}\mathbf{w}\le {\epsilon }_{\mathrm{WNG}}.\end{array}\end{array}$$

(55)

The parameter ${\epsilon }_{\mathrm{WNG}}$represents the lower limit of the WNG of the designed beamformer. It is noteworthy that ${\epsilon }_{\mathrm{WNG}}$ must not exceed ${\epsilon }_{\mathrm{WNG}}^{\max }$, which is the WNG obtained from the solution in (46). The optimal weighting vector $\mathbf{w}$ can be obtained by solving the optimization problem in (55) using the CVX toolbox [50].

4. Simulations

In this section, we will verify the effectiveness of the proposed method and investigate the performance of the beamformer through simulations. A linear array of 21 loudspeakers with a spacing of 0.04 m is used in the following simulation. The frequency range of interest is from 300 Hz to 4 kHz, which covers the frequency range of speech.

4.1. Performance Study and Comparison

We first conduct the simulations to verify the effectiveness of the proposed method. For ease of identification, we refer to the method proposed in Section 3.3.1 as the proposed method I and the method proposed in Section 3.3.2 as the proposed method II. The parameter ${\epsilon }_{\mathrm{WNG}}$ used in the proposed method II is set to 0 dB in this simulation.

The second-order and third-order target radiation pattern, with the desired direction is set as ${\theta }_s={30}^{°}$ and the main lobe width $\Delta ={60}^{\text{°}}$, are synthesized by the proposed method I and the proposed method II, respectively. Figure 3 presents comparisons between the beam patterns synthesized using the two proposed methods and the target beampattern at 2 kHz. The target beampatterns have maximum values of one, which located at the preset angle ${\theta }_s={30}^{°}$. The second-order target beampattern has a null at ${138}^{°}$. The third-order target beampattern has two nulls at ${100}^{°}$ and ${154}^{°}$. Figure 3a shows the comparison result of the second-order beampattern. It is observed that the beampattern synthesized using the proposed method I deviates from the target beampattern. Despite the value of the synthesized beampattern is one at the preset direction, the values of the resulting beampattern near ${20}^{°}$ and ${60}^{°}$ are both greater than one. This indicates that although the distortionless constraint is satisfied in the proposed method I, the finite orders of modal matching cannot ensure the maximum value of the generated beampattern in the preset direction. However, the beampattern of the proposed method II matches well the target beampattern shown with a black solid line. The same results can be seen in the comparison of the third-order beampattern shown in Figure 3b. It is noteworthy that, compared to the second-order synthesized beampattern, the third-order beampattern of the proposed method I is closer to the target beampattern. This is because the truncation error caused by (35) decreases as the order of the target beampattern increases.

The broadband synthesized radiation patterns are shown in Figure 4. Figure 4a-Figure 4b shows the broadband beampatterns of the second-order and third-order synthesized beampatterns of the proposed method I. There are two issues with the beampattern generated by the proposed method I: Firstly, the beampatterns cannot maintain consistent, as both the main lobe width and the null positions change with frequency. Secondly, the maximum values of the generated beampatterns deviate from the desired direction, which means that the generated beampattern does not closely approximate the desired beampattern. Figure 4c-4d plot the broadband beampatterns synthesized using the proposed method II. As seen, the proposed method II can design the different orders of target radiation patterns. The resulting beampatterns maintain almost frequency-invariant across the evaluated frequency range.

The performance comparisons of the two proposed methods for designing the second- and third order target differential beampatterns are shown in Figure 5. The performance measures are DF, WNG and MSE. Figure 5a shows the DF increases with the order of target beampattern increases. the DF of the proposed method I cannot remain constant values and exhibits some fluctuations at certain frequencies. However, the proposed method II can maintain constant DF values over the frequency range of interest (300 Hz – 4 kHz) for designing the different order target beampatterns. This indicates the proposed method II can achieve the frequency-invariant beampatterns in the whole evaluated frequency range. From Figure 5b, the WNGs of the proposed method II is smaller than that of the proposed method I. This is because the proposed method I uses all the remaining degrees of freedom of the array to optimize the WNG while satisfying the linear constraints. Therefore, the proposed method I has the highest WNG. The proposed method II, while meeting the linear constraints and ensuring the WNG exceeds ${\epsilon }_{\mathrm{WNG}}$, uses the remaining degrees of freedom of the array to optimize the MSE. The preset value ${\epsilon }_{\mathrm{WNG}}$ plays a role in a tradeoff between the WNG and the MSE. As can be seen in Figure 5c, The MSE of the proposed method II is significantly lower than that of the proposed method I. Especially in the range of 1-3.5 kHz, compared to the proposed method I, the proposed method II achieves a reduction of over 40 dB in MSE by only slightly decreasing the WNG value. It is interesting to note that the performance of the proposed method I in designing the third-order target beampattern is superior to its performance in designing the second-order target beampattern. It has a larger DF factor, higher WNG, and lower MSE. This means that the beamformer of the third-order ideal beampattern has stronger directivity, higher robustness, and lower mean square error. Nevertheless, the MSEs of the proposed method I are still only around -20 dB, which does not meet the requirements for practical applications. However, the proposed method II can significantly improve the mean square error by sacrificing a slight degree of robustness, which is highly valuable in real applications. Therefore, we will only study the performance of the proposed method II in the following sections. Any mentioned the proposed method hereafter refers to the beamformer obtained by solving the optimization problem (55).

4.2. Impact of the Parameter ${\epsilon }_{\mathrm{WNG}}$

${\epsilon }_{\mathrm{WNG}}$ is a key parameter to trade off the robustness and the approximation error of the proposed method. This simulation study the impact of the value of ${\epsilon }_{\mathrm{WNG}}$ on the performance of the proposed method with the third-order differential beampattern (see in Figure 3b) being the ideal radiation pattern. The level of ${\epsilon }_{\mathrm{WNG}}$ is set to four different values: 1) ${\epsilon }_{\mathrm{WNG}}={\epsilon }_{\mathrm{WNG}}^{\max }$dB; 2) ${\epsilon }_{\mathrm{WNG}}=({\epsilon }_{\mathrm{WNG}}^{\max }-2)$dB; 3) ${\epsilon }_{\mathrm{WNG}}=0$dB; and 4) ${\epsilon }_{\mathrm{WNG}}=-10$dB.

The corresponding broadband radiation patterns of the proposed method are shown in Figure 6. Figure 6a presents the synthesized beampattern when ${\epsilon }_{\mathrm{WNG}}$ equals ${\epsilon }_{\mathrm{WNG}}^{\max }$dB. It gives the same result shown in Figure 4b, indicating the solutions obtained by solving the optimized problem in (55) equal to the solution in (46) when ${\epsilon }_{\mathrm{WNG}}={\epsilon }_{\mathrm{WNG}}^{\max }$dB. Figure 6b displays the synthesized beampattern when ${\epsilon }_{\mathrm{WNG}}=({\epsilon }_{\mathrm{WNG}}^{\max }-2)$dB. With a slight decrease in the value of ${\epsilon }_{\mathrm{WNG}}$, the mean square error (MSE) between the synthesized beampattern and the ideal beampattern is reduced above 500 Hz (shown in Figure 6c), resulting in an improvement of the synthesized beampattern, although there is a slight deviation from the ideal beampattern at low frequencies. Furthermore, from Figure 6c and 6d, we can observe that as the value of ${\epsilon }_{\mathrm{WNG}}$ is further reduced, the synthesized beampattern can be closer to the ideal one across the entire evaluated frequency range.

Table 7. Figure 7a shows that the DF of the beamformer with ${\epsilon }_{\mathrm{WNG}}={\epsilon }_{\mathrm{WNG}}^{\max }$ fluctuates at certain frequencies and the DFs of the beamformers with other ${\epsilon }_{\mathrm{WNG}}$ values can maintain a constant over the whole evaluated frequencies range. This indicates that the frequency-invariant beampattern can be achieved by the proposed method with an appropriate ${\epsilon }_{\mathrm{WNG}}$ value. Figure 7b and 7c plot the WNG and the MSE of the proposed method with different ${\epsilon }_{\mathrm{WNG}}$ values. It can be observed that at low frequencies, especially below 1.4 kHz, the WNG and MSE of the beamformer decrease as ${\epsilon }_{\mathrm{WNG}}$ decreases. This means that, with the decrease in the value of ${\epsilon }_{\mathrm{WNG}}$, the synthesized beampattern more closely approaches the target beampattern, but this improvement comes at the expense of a decrease in the WNG. Above 1.4 kHz, the WNG of these beamformers dose not differ significantly. However, the MSE of the beamformers when ${\epsilon }_{\mathrm{WNG}}={\epsilon }_{\mathrm{WNG}}^{\max }-2$, ${\epsilon }_{\mathrm{WNG}}=0$, and ${\epsilon }_{\mathrm{WNG}}=-10$ is much smaller than the MSE of the beamformer when ${\epsilon }_{\mathrm{WNG}}={\epsilon }_{\mathrm{WNG}}^{\max }$. Hence, we can tune the tradeoff between the WNG and MSE of the beamformer by setting ${\epsilon }_{\mathrm{WNG}}$ to slightly less than ${\epsilon }_{\mathrm{WNG}}^{\max }$.

4.3. Impact of Loudspeaker Mismatch

It was assumed in the previous section that the loudspeakers used for designing the beamformer are ideal, without any mismatch issues among the loudspeaker units. However, in practice, there exists the uncertainty in the loudspeaker characteristics (magnitude, phase and position) due to the variations in the response of the drivers. This simulation investigates the impact of the driver mismatch on the proposed method with different ${\epsilon }_{\mathrm{WNG}}$ values. The target beampattern and ${\epsilon }_{\mathrm{WNG}}$ values are set the same as in Section 4.2. The perturbations are added to the spatial responses of the loudspeakers, for which the error has a multiplicative form with uniform distribution between −3 and +3 dB in magnitude and uniform distribution between −10° and +10° in phase. The performance measures are averaged over 1000 Monte Carlo trails.

Broadband beampatterns are plotted in Figure 8. By comparing Figure 8a and Figure 8b with Figure 6a and Figure 6b, we can observe that, aside from the difference in amplitude range, the basic shape of the beampattern remains unchanged. This indicates that when ${\epsilon }_{\mathrm{WNG}}={\epsilon }_{\mathrm{WNG}}^{\max }$ and ${\epsilon }_{\mathrm{WNG}}={\epsilon }_{\mathrm{WNG}}^{\max }-2$, the beamformer is robust enough to against the perturbations in the frequency response of the loudspeaker. However, Figure 8c and Figure 8d show significant differences from Figure 6c and Figure 6d at low frequencies, suggesting that when ${\epsilon }_{\mathrm{WNG}}=0$ and ${\epsilon }_{\mathrm{WNG}}=-10$, the beamformer’s robustness is reduced, causing the synthesized beampattern to deviate from the target beampattern due to the errors added to the loudspeaker.

Figure 9 shows the averaged DF, WNG and MSE of the beamformers. Figure 9a shows that the DF of the beamformers with ${\epsilon }_{\mathrm{WNG}}={\epsilon }_{\mathrm{WNG}}^{\max }$ and ${\epsilon }_{\mathrm{WNG}}={\epsilon }_{\mathrm{WNG}}^{\max }-2$ are almost identical over the whole evaluated frequency range. When ${\epsilon }_{\mathrm{WNG}}=0$ and ${\epsilon }_{\mathrm{WNG}}=-10$, the DFs decrease with the ${\epsilon }_{\mathrm{WNG}}$ value decrease below 1 kHz. Figure 9b and 9c plot the WNG and MSE, respectively. Unlike the conclusion drawn in Section 4.2 that a smaller WNG of the beamformer results in a smaller MSE, when there exists a mismatch in the drivers, a smaller WNG value leads to a larger MSE which means large deviations from the target beampattern at low frequencies. Hence, the appropriate selection of the ${\epsilon }_{\mathrm{WNG}}$ value is crucial for the proposed method in practical applications.

4.4. Validation of the Steering Flexibility

In this simulation, we investigate the steering flexibility of the proposed method. For the broadside array illustrated in Figure 1, we consider four different desired directions: ${\theta }_s={30}^{\circ }$, ${\theta }_s={45}^{\circ }$, ${\theta }_s={60}^{\circ }$ and ${\theta }_s={90}^{\circ }$. The coefficients of the third-order differential beampattern with mainlobe width equals 60° and the different desired directions need to be determined by solving (18) first, and then the synthesized beampattern is obtained using the proposed method. The parameter ${\epsilon }_{\mathrm{WNG}}$ used in the simulation is set to $({\epsilon }_{\mathrm{WNG}}^{\max }-2)$dB.

Figure 10 shows the synthesized broadband beampatterns. As seen, the mainlobe of the synthesized beampatterns can be steered to the desired directions. This indicates that the proposed method can design steerable different beamformers with a linear array. The synthesized beampatterns slightly deviate from the ideal beampatterns at low frequencies for ${\theta }_s={30}^{\circ }$, ${\theta }_s={45}^{\circ }$, ${\theta }_s={60}^{\circ }$ and ${\theta }_s={90}^{\circ }$. When ${\theta }_s={90}^{\circ }$, the proposed method can achieve frequency-invariant beampattern over the whole evaluated frequency range.

Figure 11 shows the DF, WNG and MSE of the beamformers with different desired directions. As seen from Figure 11a, when ${\theta }_s$ are 30°, 45°, and 60° respectively, the DF of the beamformer increases as ${\theta }_s$ increases. However, when the desired direction ${\theta }_s$ is 90°, the DF does not increase as expected. That is because, in this case, the ideal beampattern is symmetric with respect to the direction , the coefficient of the ideal beampattern satisfies ${\alpha }_{N,m}=0$ when $m$ is odd, leading to a reduction in the DF value. From Figure 11b, the WNG decrease as ${\theta }_s$ increases at low frequencies. Figure 11c shows that when ${\theta }_s={90}^{\circ }$, the MSE is below -40 dB over the whole evaluated frequency range, which is considered the ideal beampattern has been well approximated.

4.5. Comparison with Other Steerable Beamforming Methods

To demonstrate the advantages of the proposed method in designing the steerable broadband differential beampattern with linear array, we compare the performance of the proposed method with three other steerable beamforming methods: (i) the Delay and Sum (DS) method, which increases the acoustic energy in the desired direction by phase-aligning [43]. (ii) the Minimum Variance Distortionless Response (MVDR) method, which aims to minimize the output power of the array while maintaining a distortionless response in the desired direction [52]. (iii) the Null Constrained (NC) method, which utilizes the null position for the ideal beampattern to design the differential beamformer [38]. The MVDR method requires a regularization parameter during the solving process. We set the regularization parameter to 0.001 in the simulation. The steerable differential beamformer we designed aims to approach the fourth-order differential beampattern with the mainlobe width of 60° and the desired direction being steered to 120°.

Figure 12 plots the broadband beampatterns synthesized by these four methods. As seen, the mainlobes of the DS and MVDR methods become narrower as the frequencies increases. The NC method can only maintain a constant main lobe width at low frequencies. Above 600 Hz, like the DS method, its mainlobe width also narrows as the frequency increases. In contrast, the proposed method maintains the frequency-invariant beampattern over the whole frequency range.

The DF and WNG of the different beamforming methods are plotted in Figure 13. Figure 13a shows the DF of DS and MVDR increase with frequency. DF of NC maintains almost the same below 600 Hz but increases with frequencies above 600 Hz. The DF of the proposed method can maintain a constant across the whole frequency range, indicating the frequency-invariant pattern can be achieved by the proposed method. In Figure 13b, DS has the maximum WNG and remains constant across the entire frequency range. The WNG of the NC method increases with frequency at low frequencies, and at high frequencies, its value matches that of the DS. In the 500-3k Hz frequency range, the MVDR method has the lowest WNG, indicating the worst anti-perturbation ability within this frequency band. Above 1 kHz，the WNG of the proposed method is lower than that of the DS and NC methods, but it still greater than 5 dB, which is considered a proper level of robustness for practical usage.

5. Experiment and Discussion

Experiments were carried out in an anechoic chamber using a linear loudspeaker array shown in Figure 14a consisting of 31 moving-coil speaker units (B1S of HiVi Inc.) spaced 3.8 cm apart. The experimental setup is illustrated in Figure 14b. The array was mounted on a turntable and positioned in the front when the turntable angle was 0°. An omni-directional microphone (BSWA MPA201) was placed at 3 m from the geometrical center of the array. Both the microphone and the linear array were placed at a height of 2 m from the floor.

A linear swept sine signal in the frequency range of 300-4k Hz with a resolution of 5 Hz was reproduced by each loudspeaker to measure its transfer function. The signals recorded by the microphone were processed by an AP2720 audio analyzer to obtain the transfer function at a specific angle. The turntable was then rotated by 5° in the counter-clockwise direction, and the same measurement procedure was repeated. This process was repeated when the array was rotated from 0° to 180°, giving a total of 37 positions with 31 measurements each. Finally, the transfer functions were used to calculate the measured beampattern to evaluate the performance of the proposed method.

Figure 15 illustrates comparisons between the measured beampatterns and the ideal beampattern at frequencies: 500 Hz, 1 kHz, 2 kHz and 4 kHz. As can be seen, the measured beampatterns closely resemble the ideal beampattern, although there exist some minor differences due to array imperfections and measurement errors.

To further demonstrate the good performance of the proposed method, Figure 16 gives the measured broadband beampatterns of the four different methods: DAS, MVDR, NC and the proposed method. The regularization parameter used in the MVDR method is 0.001 and the WNG constraint value of the proposed method is set to ${\epsilon }_{\mathrm{WNG}}^{\max }-2$ over the whole evaluated frequency range. It can be observed that the main lobe of the beampattern generated by the DS method became narrower as the frequency increases. The MVDR method cannot generate an effective beampattern due to its poor robustness. The main lobe of the NC method is well-maintained at low frequencies. However, at high frequencies, its beampattern becomes similar to that of the DS method, with the main lobe narrowing as the frequency increases. In contrast, the main lobe of the proposed method kept almost the same in the evaluated frequency range of up to 4 kHz, which demonstrates that the frequency-invariant beampattern can be synthesized by the proposed method in a broadband of frequencies.

6. Conclusions

In this paper, we propose a method for designing the steerable frequency-invariant beamformer using a differential loudspeaker line array. The design process begins by determining the target differential beampatterns based on the desired direction, mainlobe width, and beampattern order. These target beampatterns are then represented in the modal domain. The Jacobi-Anger series expansion is employed to design the beamformer, ensuring that the resulting beampattern closely aligns with the target differential beampattern. To further enhance beampattern matching and robustness, a multi-constraint optimization problem is formulated. This approach introduces a white noise gain constraint value, enabling a trade-off between the white noise gain performance and the mean square error of the synthesized beampatterns. We also derive the upper limit for the proposed constraint value. Simulation results demonstrate that setting the proposed constraint value slightly below this upper limit can significantly improve the mean square error performance. Both simulations and experimental results demonstrate that the proposed method outperforms existing steerable beamforming techniques, achieving steerable frequency-invariant beamforming across the frequency range of 300 Hz to 4 kHz.