An Efficient Broadband Adaptive Beamformer without Presteering Delays

Broadband adaptive beamformers have been widely used in many areas due to their ability of filtering signals in space domain as well as in frequency domain. However, the space-time array employed in broadband beamformers requires presteering delays to align the signals coming from a specific direction. Because the presteering delays are direction dependent, it is difficult to make precise delays in practice. A common way to eliminate the presteering delays is imposing constraints on the weight vector of the space-time array. However, the structure of the constraint matrix is not taken into account in the existing methods, leading to a computational complexity of O(N2) when updating the weight vector. In this paper, we describe a new kind of constraint method in time domain that preserves the block diagonal structure of the constraint matrix. Based on this structure, we design an efficient weight vector update algorithm that has a computational complexity of O(N). In addition, the proposed algorithm does not contain matrix operations (only scalar and vector operations are involved), making it easy to be implemented in chips such as FPGA. Moreover, the constraint accuracy of the proposed method is as high as the frequency constraint method when the fractional bandwidth of the signal is smaller than 10%. Numerical experiments show that our method achieves the same performance of the state-of-the-art methods while keeping a simpler algorithm structure and a lower computational cost.


Introduction
Beamforming has become a fundamental technique for sensor arrays and received considerable attention in recent decades [1][2][3][4][5]. Besides providing the ability of performing spatial filtering, the broadband beamformers can also control the frequency response in a given direction [6], which makes it a power tool in various areas such as radar, sonar, mobile communications, satellite navigations, and radio astronomy. The space-time adaptive processing is one of the widely used approaches to implement broadband beamforming. An adaptive broadband beamformer based on space-time array was proposed in [7], which uses the least mean square (LMS) algorithm to compute the weight vector of the array. This method has been used in a wide range of applications due to its low computational complexity and good numerical stability. Another commonly used technique for broadband beamforming is the space-frequency adaptive processing [8,9], which splits the frequency band of the incident signal into several narrow subbands and then uses narrowband beamformers to deal with the narrowband signals. The space-frequency processing has a faster convergence rate than space-time processing when the number of taps is large (e.g., hundreds or thousands). Nevertheless, space-frequency beamformer has a larger processing delay than space-time beamformer [10,11]. Therefore, space-time adaptive processing is still the first choice for delay sensitive applications such as satellite navigations [12,13].
To obtain a desired frequency response in the direction of the signal of interest (SOI), constraints should be imposed on the weight vector of the beamformer. Most of these constraints require that the array has been presteered toward the direction of SOI. However, it is difficult to implement precise delays in practice [14,15], especially for the situations that sensor array has to form multiple beams [16,17]. In addition, the presteering errors may cause signal cancellation in the direction of SOI, leading to performance degradation of the beamformer [18,19]. Therefore, it is very important to eliminate the presteering delays in broadband beamformers. A type of convolution constraint method was proposed in [20] to remove the presteering delays. This method has a simple procedure to implement; however, it is computationally expensive. The authors in [21] described another technique to remove the presteering delays, which multiplies the received signal by a matrix whose elements are the inverse Fourier transform of the steering vector. Thus, it is computationally expensive too. By introducing a set of frequency domain constraints (FDC), the authors in [22] developed a simple technique to eliminate the presteering delays, which requires less tapped delay-lines (TDLs) than the convolution constraint method. Moreover, it can incorporate the sensor patterns into the beamforming algorithm and thus can be applied to conformal arrays [23]. The FDC method based on generalized sidelobe canceller (FDC-GSC) was presented in [24], which is mathematically equivalent to FDC but with fewer computations. To further reduce the number of TDLs, infinite impulse response filter can be used [25,26]. However, all of these methods destroy the block diagonal structure of the constraint matrix, making the computational complexity of weight vector update increased from O(N) to O(N 2 ).
To reduce the computational complexity of weight vector update, we design a new kind of time domain approximate constraint (TDAC) method, which does not require the presteering delays while preserving the block diagonal structure of the constraint matrix. The constraint accuracy of TDAC is as high as the FDC method if the fractional bandwidth is smaller than 10%. This is a reasonable assumption because the fractional bandwidth of most radar and communication systems is smaller than 10%. Moreover, by exploiting the block diagonal structure of the constraint matrix, a new efficient weight vector update algorithm with a complexity of O(N) is also developed. In fact, the computations involved in each iteration of the proposed algorithm is only half of that of the conventional Frost algorithm (CFA). Both LMS and normalized LMS (NLMS) algorithms can be used to update the weight vector of the space-time array. We shall study the NLMS-based algorithm in this paper because the one based on LMS can be easily obtained from NLMS.
The rest of this paper is organized as follows. Section 2 reviews the signal model of space-time array and the conventional Frost algorithm. Section 3 describes the details of the proposed method, including algorithm design, geometrical interpretation, and complexity analysis. Numerical simulations are presented in Section 4, followed by conclusions in Section 5.

Space-Time Adaptive Beamformer
This section reviews the signal model of space-time array and the conventional Frost algorithm, including space-time steering vector, frequency-wavenumber response (FWR), and weight vector update methods of Frost algorithm based on LMS and NLMS. Figure 1 shows a space-time array with M sensors and each sensor is followed by a time delay unit and a finite impulse response (FIR) filter of L taps (the radio frequency chain and the analog-to-digital converter are omitted for simplicity). The aim of presteering delays τ m (θ, φ) is to align the signal coming from the direction of SOI so that there is no phase difference among the SOI received by different sensors. The total degrees of freedom of the array, i.e., the number of free weights [6], is N = ML, and the N × 1 baseband signal vector at the kth time instant is given by
If there are no presteering delays, the FWR of the array for the signal coming from (θ, φ) at frequency ω c + ω can be expressed as [27] where k = [(ω c + ω)/c]a is the wavevector of the impinging signal, ω c is the carrier frequency, ω is the baseband frequency, c is the speed of light, a = −[sin(θ) cos(φ), sin(θ) sin(φ), cos(θ)] T is the propagation direction, p m is the position vector of the mth sensor, and T s is the sampling period of the baseband signal. The steering vector v(ω, k) of the space-time array can be expressed as [28] v(ω, k) where is the temporal steering vector and is the spatial steering vector.

The Conventional Frost Algorithm
The linearly constrained minimum variance beamformer is a widely used technique to suppress the interferences while keeping the SOI, which minimizes the output power of the array subject to J linear constraints as follows where C is the N × J constraint matrix whose columns are linearly independent and g is the J × 1 gain vector. The solution of (7) is given by [27] It is not advisable to compute the weight vector by solving the above equation for real-time processing applications, because one has to estimate the covariance matrix R x and compute the matrix inversion (or solve linear systems of equations). Frost developed an adaptive approach to compute the weight vector based on the LMS algorithm as follows [7]: where µ is the step size parameter, is the quiescent weight vector in the column space of C, denoted by R(C), and is the projection matrix onto the orthogonal complement of R(C), denoted by R ⊥ (C). The normalized version of (9) is given by [29] w(k

The Proposed Method
Although the CFA (9) (or (12)) is simple, it requires the presteering delays. In this section, we first introduce a new kind of constraint for the space-time array, which eliminates the presteering delays and enables the weight vector to be updated efficiently. Then, inspired by [7], we give a geometric interpretation of the proposed algorithm. Finally, the comparison of computational complexity of our method and the existing methods is provided.

The Approximate Constraints in Time Domain
Suppose that the number of constraints J in (7) is equal to the number of taps L of the array, and let w H l v s (k) = g * l , then the FWR (3) of the space-time array can be expressed as where Ω = ωT s ∈ [−π, π] is the normalized frequency [30]. To keep the symbols consistent with (7), we define w H l v s (k) = g * l , i.e., v H s (k)w l = g l . If g l is independent of Ω, then Γ(ω, k) is equal to the discrete Fourier transform (DFT) of an FIR filter whose lth coefficient is equal to the inner product of w l and v s (k).
Equation (13) shows that the broadband beamformers can perform frequency filtering as well as spatial filtering. If the frequency response in the direction of SOI is determined by an FIR filter with coefficients g = [g 0 , · · · , g L−1 ] T , we can impose constraints on w l in the form of v H s (k)w l = g l . However, because the wavevector k depends on baseband frequency ω, the coefficients g l in (13) depend on ω too. Nevertheless, for systems that operate with a small fractional bandwidth, the difference between k and k 0 = (ω c /c)a is small and k can be approximated by k 0 . Thus, we obtain the following approximate constraints for the weight vector The above approximation is reasonable for many practical broadband systems that operate with a fractional bandwidth smaller than 10%, such as radar [28], satellite navigations [12], and wireless communications [31]. The benefit of this approximation is that it allows an efficient implementation of the weight vector update.
The constraints given by (14) can be rewritten as C H w = g, which takes the same form of the constraint given in (7). The constraint matrix C has a block diagonal structure as follows where I L is the L × L identity matrix. The main computations of (12) come from the projection operation. Thus, we should derive a simple form for the projection matrix P. and Let where s(k 0 ) = v s (k 0 )/ √ M is the normalized spatial steering vector. Then the projection matrix P can be written as Expressing the N × N projection matrix P by an M × M projection matrix Q is the key point of the efficient algorithm, which will be described in the next subsection.

An Efficient Weight Vector Update Algorithm
Because w q is independent of time instant k, only the second term on the right hand side of (12) needs updating. Let where Then w(k + 1) can be decomposed into two terms as follows We call w a (k + 1) the adaptive weight vector because it is updated adaptively according to the input signal vectors x(k).
Since P is the projection matrix onto R ⊥ (C) and w q ∈ R(C), we have Pw q = 0. Thus, By using the idempotent property P 2 = P of the projection matrix [33], we have Hence, Next we show how to compute z(k) = Px(k) and p(k) = x H (k)Px(k) efficiently. Because, as shown in Figure 1, x l (k) has a time delay structure x l (k) = x 0 (k − l), has a similar time delay structure. It follows from (1) and (19) that where Therefore, only z 0 (k) needs to be computed at each iteration, which involves 2M − 1 complex additions and 2M complex multiplications. Although p(k) can be computed by x H (k)z(k) with N − 1 complex additions and N complex multiplications. However, there exists a more efficient method to compute p(k) as follows where the time delay properties x l−1 (k − 1) = x l (k) and z l−1 (k − 1) = z l (k) are used. Since z 0 (k) has already been calculated and each term on the right hand side of the above equation is real, there are only M complex additions and M complex multiplications in computing p(k) if we store the quantities x H 0 (k)z 0 (k). A circular array [34] is employed in our method to store the latest L + 1 quantities of x H 0 (i)z 0 (i) for i = k, k − 1, · · · , k − L. The circular array, denoted by q, is shown in Figure 2, where the position front points at the current quantity q(k) and the position back points at the Lth previous quantity q(k − L). When new data q(k) arrives, the circular array overwrites q(k − L) with q(k).
The final description of our algorithm is summarized in Algorithm 1. Please note that there are no matrix operations in our algorithm, which means that the algorithm can be implemented at the level of scalar and vector operations. This feature is very important when the algorithm is implemented in chips such as field programmable gate array (FPGA) and digital signal processor (DSP).

Algorithm 1:
The TDAC algorithm based on NLMS.

Geometrical Interpretation
The new weight vector update algorithm has a simple geometrical interpretation as shown in Figure 3, where C H w = 0 and C H w = g define the constraint subspace and the constraint plane, respectively. It follows from (16) that w q ∈ R(C). Thus, it is perpendicular to the constraint subspace R ⊥ (C). Meanwhile, because w q satisfies the constraint equation C H w q = g, it terminates on the constraint plane. In addition, as defined by (20), w a (k) is the projection ofw(k) onto R ⊥ (C). Thus, it lies in the constraint subspace.
From Algorithm 1 we know that only w a (k) is updated during the iterations. To minimize the output power at the kth time instant, an increment −[µy * (k)/p(k)]x(k), based on the NLMS criterion, is added to the adaptive vector w a (k). However, this change may move w a (k) off the constraint subspace. Thus, the projection of the increment −[µy * (k)/p(k)]Px(k) is used to update w a (k) to ensure that w a (k) lies in the constraint subspace.

Computational Complexity Analysis
The computations involved in each iteration of CFA, FDC, FDC-GSC, and TDAC, in terms of complex additions and complex multiplications, are shown in Table 1. We see from the table that our method (TDAC) has the least computations, which is even more efficient than CFA. In contrast, FDC and FDC-GSC increase the weight vector update complexity from O(N) to O(N 2 ).
In addition, computing the projection matrix P at the initial stage of FDC requires a complexity of O(N J 2 ), where J is the number of constraint points in frequency band. Similarly, computing the blocking matrix B at the initial stage of FDC-GSC requires a complexity of O(N J 2 ). The convolution constraint method has the same update complexity of FDC. However, it needs DFT operations to form the constraint matrix, leading to a higher complexity than FDC. In contrast to FDC and FDC-GSC, the proposed method only requires N additions and 2N multiplications at the initial stage.

Simulation Results
In this section, we present some numerical experiments to compare the performance of the proposed algorithm with other methods. A 10 × 10 space-time array is employed in the first experiment. More specifically, a uniform linear array (ULA) located along the z-axis is considered. The array consists of 10 isotropic antennas whose interelement spacing is 0.5λ where λ is the wavelength corresponding to the carrier frequency. In baseband, each sensor is connected to an FIR filter with 10 TDLs. The distortionless gain vector is set to g = [0, 0, 0, 0, 1, 0, 0, 0, 0, 0] T . One desired signal and two uncorrelated interferences impinge on the array from 60 • , 80 • and 120 • , respectively. Noises are assumed to be spatially and temporally uncorrelated zero-mean Gaussian random processes with equal power. The signal-to-noise ratio (SNR) is 10 dB and the interference-to-noise ratios (INR) are 40 dB and 30 dB respectively. Both the signal and interferences occupy a bandwidth of 100 MHz around the carrier frequency of 1000 MHz, i.e., the fractional bandwidth is 10%. We compare our method (TDAC) with the methods of CFA [7], FDC [22], and FDC-GSC [24]. The NLMS-based update Equation (12) with µ = 0.02 is used for all of the tested algorithms. For the FDC method, 10 constraint points are uniformly selected in the frequency band.
The performance of broadband beamformer is measured in terms of array output power and signal-to-interference-plus-noise ratio (SINR). The output SINR is defined as where R s is the correlation matrix of the desired signal and R i+n is the correlation matrix of interferences plus noise. We assume that the power spectral densities of the desired signal and the interferences are flat in the bandwidth of considered. We also assume that the signal, interferences, and noise are statistically independent. Under these assumptions, the correlation between the lth tap after the mth sensor and the pth tap after the qth sensor for the dth impinging wave can be expressed as [35] where σ 2 d is the power of the dth signal, τ (d) m is the propagation delay of the dth signal at the mth sensor with respect to the origin, and B d is the bandwidth of the dth signal. For the simulation scenario described above, R s = R 1 , R i+n = R 2 + R 3 + σ 2 n I N , where σ 2 n is the noise power. Figure 4 shows the FWRs of CFA, FDC and TDAC when the adaptive algorithms have converged (for ULA located along the z-axis, FWR depends on frequency and polar angle). Because the FWR of FDC-GSC is the same as FDC, it is not shown here. It can be seen that all methods have constant magnitude responses at the constrained direction while placing nulls at the interference directions in the whole frequency band. Because the array is presteered toward the direction of SOI by CFA, the equivalent directions of interferences are also changed. In this example, the equivalent directions of the desired signal and interferences are 90 • , 109 • and 180 • respectively, which can be verified in Figure 4a. In addition, we see from Figure 4b,c that the FWR of TDAC is smoother than that of FDC. Hence TDAC may have better SINR performance than FDC in this experiment due to the better sidelobe structure.
The magnitude and phase errors in the constrained direction of CFA, FDC, FDC-GSC and TDAC are plotted in Figure 5. It can be seen that (i) the constrained response is equal to the distortionless response when the array is precisely presteered by CFA; (ii) FDC and FDC-GSC have the same frequency response since they are equivalent; (iii) there exist ripples in the frequency response of FDC, this is because only 10 frequency points are constrained by FDC and the responses between these constrained points are not guaranteed; (iv) both FDC and TDAC provide good approximations to the desired frequency response; and (v) TDAC has a smaller phase error than FDC.    Figure 6 shows the learning curves of the output power and SINR of different methods averaged over 100 independent trials. As shown in Figure 6, the four tested methods have the same convergence rate, but our method has the best steady-state SINR performance. The reason CFA has a smaller SINR than other methods is that when the array is presteered toward the direction of SOI, the directions of interferences are also changed. Thus, the array faces a different interference environment, leading to different SINR performance. Although FDC can impose accurate constraints in frequency domain, the frequency responses between the constrained points are not guaranteed, leading to ripples in frequency band as shown in Figure 5. Therefore, similar to the proposed method that approximates the FWR in time domain, FDC is also a type of approximate method that approximates the FWR in frequency domain. In the second experiment we change the bandwidth from 100 MHz to 200 MHz while keeping the carrier frequency fixed at 1000 MHz, i.e., the fractional bandwidth is 20%. We also change the length of FIR filter from 10 to 20, i.e., we use a 10 × 20 space-time array. The proposed algorithm is compared with the FDC methods with 20 constraint points (denoted by FDC1) and 30 constraint points (denoted by FDC2) respectively. The simulation results of FWRs, constraint errors, and learning curves are shown in Figures 7-9, respectively.

Conclusions
An efficient implementation of the broadband adaptive beamformer without presteering delays was studied in this paper. A new kind of approximate constraint in time domain has been proposed to eliminate the presteering delays of the space-time array. In addition, a new weight vector update algorithm was developed by using the block diagonal structure of the constraint matrix, leading to a computational complexity of O(N) in each iteration. In contrast, the computational complexity of the frequency constraint methods is O(N 2 ) in each iteration. Moreover, the new algorithm does not contain matrix operations and can be implemented at the level of scalar and vector operations. This feature is very important for real-time applications, in which the algorithms should be implemented in chips such as FPGA and DSP. Numerical experiments shown that the approximate accuracy of the proposed method is as high as the frequency constraint method for systems with a fractional bandwidth smaller than 10%, while our method has a simpler algorithm structure and a lower computational cost than the existing methods.