Research Article Generalized Broadband Beamforming Using a Modal Subspace Decomposition

We propose a new broadband beamformer design technique which produces an optimal receiver beam pattern for any set of field measurements in space and time. The modal subspace decomposition (MSD) technique is based on projecting a desired pattern into the subspace of patterns achievable by a particular set of space-time sampling positions. This projection is the optimal achievable pattern in the sense that it minimizes the mean-squared error (MSE) between the desired and actual patterns. The main advantage of the technique is versatility as it can be applied to both sparse and dense arrays, nonuniform and asynchronous time sampling, and dynamic arrays where sensors can move throughout space. It can also be applied to any beam pattern type, including frequency-invariant and spot pattern designs. A simple extension to the technique is presented for oversampled arrays, which allows high-resolution beamforming whilst carefully controlling input energy and error sensitivity.


Motivation and background
Beamforming is an important area of research with applications in acoustics, wireless communications, sonar, and radar [1,2]. A broadband beamformer at the receiver is based on a time-sampled sensor array-a finite set of wavefield samples taken throughout space and time. Weighted linear combinations of these space-time samples are used to filter far-field sources based on their angle of arrival and frequency content.
More formally, the response of a broadband beamformer can be expressed as a beam pattern-a 2D function of angle and frequency/wave number. A perfect beam pattern is designed to reject noise and interfering sources, but depending on sensor positions and time sampling, this desired pattern may not be achievable by a particular array. The beamformer design problem is to find an achievable pattern as close as possible to the desired pattern.
One popular technique for beam pattern design is the frequency decomposition method [3,4]. In this method, the time sequences measured at each sensor are projected into narrowband frequency bins using the discrete Fourier transform (DFT). A broadband beam pattern can then be created by using established narrowband beamforming techniques within each frequency bin.
One limitation of this technique is that sensors must be fixed, and closely spaced. To avoid spatial aliasing in the highest frequency bin, sensors must be placed no further than half of the corresponding wavelength apart [5][6][7].
Another limitation is that perfect frequency decomposition would require an infinite-length time sequence sampled at the Nyquist rate (twice the highest design frequency). Although infinite sequences are not possible in a real-world application, the decomposition can be well approximated using a windowed DFT if sufficiently long time sequences are available [8].
In many applications, however, getting a sequence of reasonable length may be impossible. Consider an environment where the source field is changing rapidly-past samples will quickly become useless, and short time sequences must be used to maintain relevance. Also, sequences may need to be shortened to minimize the computational complexity of filtering.
On top of these limitations, traditional broadband beamforming techniques suffer from a lack of versatility. Most implementations focus on the design of frequency-invariant beam patterns, where the angular response is constant across all frequencies, and there is little scope for more complex patterns. They are also generally based on the simple model of a static sensor array in space, sampled synchronously and uniformly throughout time. This can cause problems in real-world applications, such as sensor networks, which may require asynchronous sampling (where different sensors are sampled at different times), or dynamic arrays (where sensors are able to move throughout space).
One further interesting property of beamformers is the super-resolution effect, where infinite levels of directivity can be achieved using closely spaced sensor elements [9,10]. Super-resolution beamformers tend to be unsuitable for realworld applications as they require huge amounts of energy in weighting coefficients, and become extremely sensitive to small measurement errors [11,12].

Problem statement and contributions
The challenge is to develop a beamformer design technique which addresses the weaknesses and limitations identified in the previous section. This "ideal" technique would be applicable to any desired beam pattern (including both frequencyinvariant and more complex types), any sampling scheme (including both sparse and dense samplings, nonuniform and asynchronous time sampling, and moving sensors), and would always produce the optimal achievable pattern (where optimality is defined as minimizing the mean-square error, or MSE, between the desired and actual patterns).
In this paper, we introduce a method for "ideal" beamformer design based on a modal subspace decomposition (MSD). Although previous authors [13] have used projection into a modal subspace to optimize performance given a set of linear constraints, we are solving a slightly different problem. Rather than projecting valid beamformer weightings into an estimation space, we optimize the beamformer by projecting a desired (but perhaps unachievable) beam pattern into the subspace of achievable patterns.
(1) In Section 3, we define the subspace of achievable beam patterns for any set of space-time sampling positions, and derive a modal basis for this space. This basis is independent of the wavefield and desired beam pattern, and depends only on the sampling geometry in space and time. (2) In Section 4, we show that the optimal (in terms of minimizing the MSE) beamformer can be designed by orthogonally projecting any desired beam pattern into the subspace of achievable patterns. This projection can be performed by projection into the individual modal basis functions. (4) Optimal beam pattern design for closely spaced arrays will introduce the super-resolution effect, and the associated problems with error sensitivity and high input power. In Section 6, we show how a modified subspace projection based on a reduced set of modes can be used to find suitable tradeoff between resolution and these negative effects. (5) In Section 7, we discuss the simple extension of this method to more complex beamforming problems such as near-field and azimuth/elevation.

INTRODUCTION TO BEAMFORMING
In broadband far-field beamforming, the source distribution is a 2D function F(k, φ), where φ is the angle of arrival, k = 2π f /c is the wave number, c is the speed of wave propagation, and f is frequency. Consider wavefields f (x, t) in space and time. For notational simplicity, we consider only 2D space (denoted by the position vector x = (x, θ) in polar coordinates, where x = |x| and θ = ∠x). The extension to 3D space is considered in Section 7. A wavefield is generated by integrating the far-field distribution over all azimuth angles (φ = [−π, π]), and some bandlimited range of wave-numbers (k = [k 1 , k 2 ]) [14], This wavefield must now be sampled in space and time. An example of a traditional broadband beamforming array is depicted in Figure 1(a), where an array of Q sensors in space (x q , for q = 0, . . . , Q − 1) is sampled uniformly and synchronously at N relative snapshots throughout time (t n , for n = 0, . . . , N − 1). Thus, the space-time sampling can be completely specified by a uniform N × Q element grid-the Cartesian product of the time and space sampling.
This simple grid framework is inadequate to deal with more complex space-time sampling schemes, some examples of which are depicted in Figure 1(b). Consider asynchronous sampling where different sensors are sampled at different times, or consider dynamic arrays where the relative spatial position of sensors changes throughout time. To deal with these complexities, a more general model for space-time sampling is needed, which subsumes the traditional Cartesian product framework. An obvious choice is to treat each space-time sample independently. In this generalized model, the traditional array with Q sensors and N snapshots is represented as an M = N × Q dimensional space-time sampler, with m = qN + n. That is, each space-time sample is treated as an individual sensor element.
Michael I. Y. Williams et al. Following on from the traditional view of sensor arrays, the traditional view of a broadband beamformer is a set of length-N FIR filters attached to each of the Q sensors. We generalize this notation to match the notation of Definition 1. Any possible linear beamformer can be represented as a linear combination of the M field observations. ( It is worth noting that this more general model subsumes the traditional FIR filter model. Every traditional broadband beamformer can be equivalently expressed in the form of Definition 2 by simply treating each filter coefficient as an independent weighting coefficient. Combining (1) and (2), the beamformer output can be expressed as where we have defined an achievable beampattern [2], As long as every space-time sample is unique (i.e., no point in space-time is sampled more than once), every achievable beam pattern W ach (k, φ) will correspond uniquely with a particular vector of weighting coefficients w. Consider some desirable, but perhaps unachievable, beam pattern denoted by W des (k, φ). For a particular spacetime sampler, the beamforming problem is to find a set of weightings w which produce a beam pattern W ach (k, φ) as close as possible to the desired pattern. One measure of closeness is the mean-squared error (MSE) between the desired and achievable patterns:

Operators and spaces
In this section, we derive a modal basis for the subspace of achievable beam patterns. The first step is to formally define some vector and function spaces. Define S, the space of all finite-energy weight vectors, and therefore the space of all attainable beamformers based on the inner product where * denotes the complex conjugate, and the associated norm S is an M-dimensional complex vector space.
As mentioned in the previous section, the beam pattern design process is often based on some desirable beam pattern denoted by W des (k, φ). To formalize this concept, we define F , the space of desired patterns as those patterns with finite energy over the design ranges of k and φ,  based on the inner product and associated norm F is an infinite-dimensional, separable, Hilbert space 1 . Given a finite-dimensional space-time sampler, not all desired beam patterns will be achievable. With this in mind, we can partition F into two orthogonal subspaces-W is the space of achievable beampatterns and ⊥W is the space of unachievable patterns. Now, any desired patterm W des (k, φ) ∈ F can be expressed as where W ach (k, φ) ∈ W and W unach (k, φ) ∈ ⊥W . As mentioned earlier, each achievable pattern W ach (k, φ) ∈ W has a unique mapping with a weighting vector w ∈ S. Thus, since S is an M-dimensional space, W must be an M-dimensional proper subspace of F . The mapping between the two is defined by an invertible linear operator A : S → W that projects a set of weighting coefficients to its corresponding achievable beam pattern. From (4), the operator is defined by The relationships between these spaces and operators are summarized in Figure 2.

Modal bases
Given that the adjoint operator A * exists (see the appendix), operators A * A and AA * are known to have some particularly useful properties [15]. Specifically, the M eigenvectors of A * A denoted by u n form a complete orthonormal basis for S, and the M eigenfunctions of AA * denoted by U n (k, φ) form a complete orthonormal basis for W . Although infinitely many different sets of basis functions can be found for these spaces, these particular sets are uniquely significant due to their relation to the operator A, and so are known as modes.
The vector modes u n are found in the appendix to be the solutions to a matrix eigenvector equation of the form where the elements of the M × M matrix Z are given by and the real, nonnegative eigenvalues are ordered to form a monotonically decreasing series λ 0 ≥ λ 1 ≥ · · · ≥ λ M−1 . In general, this integral has no closed-form solution, but can be calculated numerically to any required degree of accuracy. Any weighting vector can be expressed as a linear combination of these vector modes, The continuous modes U n (k, φ) are found in the appendix to be given by where u n,m denotes the mth element of the nth vector mode. Any achievable beam pattern can be expressed as a linear combination of these continuous modes, From (13) and (17), the vector and continuous modes are related by the relationship where λ n are the matrix eigenvalues from (14). To summarize, given M space-time sampling locations, a set of M vector modes u n can be derived which form a basis for the space of allowable weight vectors. Each vector mode can be used to derive a continuous mode U n (k, φ) and the set of M continuous modes forms a basis for the space of achievable beam patterns. These modes are independent of the wavefield and desired beam pattern, and depend only on the geometry of the sampling. Thus, the modal bases need only be derived once for any particular array.

MSD BEAMFORMING
Given some desired pattern W des (k, φ) ∈ F , we need to find the achievable pattern which minimizes the MSE as defined in (5). From functional analysis, this optimal pattern will be the orthogonal projection of W des (k, φ) onto the subspace W [15]. This projection is denoted by B : F → W , and is shown in Figure 2. As the continuous modes U n (k, φ) form an orthonormal basis for W , this projection can be performed by projecting onto the modes, For all but the simplest patterns, this projection will need to be calculated numerically. Quadrature techniques can be used to calculate the 2D integral to any desired degree of accuracy with low computational complexity [16]. The modal coefficients b n can now be used to find the beamformer weighting coefficients. Combining (19) and (20), Equating both sides, the weighting coefficients are These coefficients can then be used in Definition 2 to create a beamformer with pattern W ach (k, φ).

The MSD beamforming algorithm
In summary, the steps of the algorithm are the following.
(1) Given a set of M space-time sampling positions (x m , t m ), build the matrix Z using (15). (2) Find the eigenvectors u n and eigenvalues λ n of Z, and use (17) to build the continuous modes U n (k, φ). (3) Given some desired beam pattern W des (k, φ), the optimal beamformer weighting coefficients will be Note that whilst the eigendecomposition in the second step is the most computationally intensive part of the process, the first two steps depend only on the space-time sampling positions, and need only be recalculated if the rare event that the space-time sampling geometry is changed. Thus, for an existing array, only the third step must be performed when the desired pattern changes. This is particularly important for adaptive beamforming, where the desired beam pattern may change over time. In this situation, the third step must be constantly recalculated to ensure that the designed beam pattern remains optimal.

DESIGN EXAMPLES
In this section, we apply the MSD design technique described in the previous section to an acoustic beamforming problem with c = 350 m/s, and a design frequency range between f 1 = 400 Hz and f 2 = 4 kHz.

Frequency-invariant beamforming
One of the more common problems in broadband beamforming is frequency-invariant beamforming, where the aim is to design a beam with identical angular response across all design frequencies [5][6][7]. To use the MSD technique, some desired pattern is needed. Thus, a frequency-invariant pattern is used with a beam centre at π/2, and a beam width of π/5. Consider a 5-element uniform circular array (UCA) designed for frequency f 2 , with a radius such that the elements are a half-wavelength apart, and 32 time samples are taken at the Nyquist rate. The frequency-invariant design pattern is projected into the modal basis, and the resulting achievable beam pattern is shown in Figure 3(a).
As mentioned earlier, one major strength of the MSD technique is that it is equally applicable to sparse arrays. To demonstrate this, a second array with the same spatial radius, but with only 3 elements, is used. Also, the sampling rate is halved so that only 16 samples are taken over the same time period. Figure 3(b) shows the resulting optimal achievable pattern. Whilst the pattern is clearly inferior to that for the denser array, it is still reasonably directional and frequency invariant.
A comparison with traditional frequency-invariant design techniques is difficult, as these techniques tend to be based on long time sequences, and perform poorly on the short sparse sequences used in this example. Note that since the MSD algorithm is optimized over all possible beamformers, it is impossible for any other design technique to produce a better pattern (in terms of the MSE).

Spot beamforming
Whereas most of the work in broadband beamforming has focused on frequency-invariant design, the MSD design technique is far more versatile. Consider the pattern shown in Figure 4(a), which filters sources based on both frequency and direction-effectively applying different bandpass filters to signals depending on angle of arrival. Projecting this pattern into the modal basis for the 5-element, 32-sample array used in the previous example results in the achievable pattern shown in Figure 4(b).

IMPLEMENTATION ISSUES: OVERSAMPLING AND SUPER-RESOLUTION
Traditional sensor arrays are designed around a particular frequency-sensors are placed a half-wavelength apart, and time samples are taken at the Nyquist rate. In broadband beamforming, the array is usually designed around the maximum design frequency [6]. Space-time sampling denser than this benchmark is known as oversampling.
The main advantage of oversampling is the large number of available space-time samples. The space of achievable beam patterns grows in dimension, and a wider range of patterns becomes achievable. In theory, we can get infinite levels of resolution by continually increasing sampling density.
This result is a due to the super-resolution (or supergain) effect where infinite levels of directivity can be achieved using closely spaced sensor elements [9,10]. Superresolution beamformers tend to be unsuitable for real-world applications, as they require huge amounts of energy in weighting coefficients, and become extremely sensitive to small measurement errors [11,12]. In this section, we will examine how these problems arise in the MSD technique, and how to manage the negative effects.
Consider a beam pattern defined by a set of modal coefficients b n in (20). The energy required to implement this pattern as a beamformer can be found by taking the norm of (23), The error sensitivity can be found by adding an error term, denoted by e m , to each field observation in (2). Consider an error modeled as spatially and temporally uncorrelated, with variance E{|e m | 2 } = σ 2 . The perturbed beamformer output will be Thus, the input energy w 2 is effectively the white noise gain [1]. Equations (25) and (28) show that the white noise gain of a modal beamformer depends directly on the inverse of the eigenvalues λ n . Figure 5 shows that as sampling becomes denser, the eigenvalues tend to converge towards zero-the sparse 3-element array has very few small eigenvalues, while more than half of the eigenvalues for the oversampled 6element array are close to zero. So, as space-time sampling is made denser, the eigenvalues approach zero, and the input energy and error sensitivity will approach infinity. To avoid this problem, but still retain the ability to exploit dense arrays and super-resolution, we must make some changes to the modal subspace projection. Consider the reduced set of eigenvalues λ n ≥ ζ, where ζ is some threshold. This set will have a reduced number of elements M ≤ M. Effectively, we are throwing away the modes which are causing the problems-those with eigenvalues closest to zero. We can now alter (23) to build beamformer weights from a reduced set of modes, The effect of threshold is shown in Figure 6 for the 5-element, 32-time-sample UCA considered earlier (and assuming a beam pattern with uniform modal coefficients). It is clear that a balance must be struck between a high threshold which will use more modes and produce the best beam pattern, and a low threshold which will use fewer modes and minimize the error sensitivity and input energy. In this way, superresolution can be achieved for high ζ if high input power and error sensitivity can be tolerated. Typically, a good choice of threshold is somewhere on the "knee" of the curve in Figure 5, where the eigenvalues converge rapidly to zero.

EXTENSIONS AND FURTHER WORK
Although the MSD technique has been presented in the context of the broadband far-field beamforming problem, the same basic concept can easily be applied to more complex beamforming problems. Equation (1), which maps a far-field source distribution to a 2D wavefield, can be altered to model a more complex mapping. Alternative mappings might allow near-field sources, sources distributed in both azimuth and elevation angles, or wavefields in 3D space. Although this altered mapping will change the form of the modal basis, the same basic methods and derivations can be applied. These extensions will be explored in a future paper.
We showed in Section 4.1 that one strength of this algorithm is that much of the modal decomposition processing can be performed offline. One important area for further work, however, is an investigation of how the computational complexity of the MSD technique compares with other comparable design techniques.

CONCLUSION
A new technique for broadband far-field beamforming has been proposed which produces an optimal beamformer for any set of space-time sampling positions. The MSD technique is based on projecting a desired pattern into the subspace of achievable patterns, and is optimal in the sense 8 EURASIP Journal on Advances in Signal Processing that it minimizes the mean-squared error (MSE). The main advantage of the technique is its versatility as it can be applied to both sparse and dense arrays, nonuniform and asynchronous time sampling, and dynamic arrays where sensors can move throughout space. Design examples were presented to show that the technique is applicable to both frequencyinvariant and spot beamformings. For oversampling arrays, a restricted modal expansion allows control over both superresolution and its associated negative effects. Future work will extend this basic technique to more complex situations including near-field sources, and azimuth/elevation beamforming. Since A is a bounded linear operator, the operators A * A and AA * are positive, selfadjoint linear operators.
The M orthonormal eigenvectors of the operator A * A, denoted by u n , form an orthonormal basis for the space S of weighting coefficients [15], λ n u n = A * Au n , ( A . 5 ) As A * A is positive and selfadjoint, the eigenvalues λ n are real and nonnegative. Thus, it makes sense to order the eigenvalues so that λ 0 ≥ λ 1 ≥ · · · ≥ λ M− In general, this integral has no closed-form solution, but can be calculated numerically to any required degree of accuracy.
From the same principles, the M orthonormal eigenfunctions of the operator AA * , denoted by U n (k, φ), form an orthonormal basis for the space W of achievable beampatterns, where u n,m denotes the mth element of the nth vector mode.