A recursive Bayesian beamforming for steering vector uncertainties

A recursive Bayesian approach to narrowband beamforming for an uncertain steering vector of interest signal is presented. In this paper, the interference-plus-noise covariance matrix and signal power are assumed to be known. The steering vector is modeled as a complex Gaussian random vector that characterizes the level of steering vector uncertainty. Applying the Bayesian model, a recursive algorithm for minimum mean square error (MMSE) estimation is developed. It can be viewed as a mixture of conditional MMSE estimates weighted by the posterior probability density function of the random steering vector given the observed data. The proposed recursive Bayesian beamformer can make use of the information about the steering vector brought by all the observed data until the current short-term integration window and can estimate the mean and covariance of the steering vector recursively. Numerical simulations show that the proposed beamformer with the known signal power and interference-plus-noise covariance matrix outperforms the linearly constrained minimum variance, subspace projection, and other three Bayesian beamformers. After convergence, it has similar performance to the optimal Max-SINR beamformer with the true steering vector.


Introduction
Digital beamforming is widely used in array signal processing for enhancing a desired signal while suppressing interference and noise at the output of an array of sensors [1,2]. It has applications in fields, such as radar, sonar, radio astronomy, speech processing, and wireless communications. [1][2][3][4][5]. It is well known, however, that the digital beamformers are sensitive to error in the estimated signal steering vector. Any errors in the steering vector will lead to signal distortions and degrade the beamforming performance severely. The causes of steering vector error in practical applications include improper array modeling, pointing error, miscalibration, and source motion as well as other effects [6][7][8][9][10]. Robustness with respect to steering vector uncertainties in digital beamformers is desirable.
Several approaches are known to partly overcome the problem of arbitrary steering vector error [11]. The most popular of them are the diagonal loading approaches [12] *Correspondence: hanyb@mail.njust.edu.cn 1 School of Electronic and Optical Engineering, Nanjing University of Science & Technology, Nanjing, 210094, China Full list of author information is available at the end of the article and constrained minimum variance approaches [13][14][15][16]. The diagonal loading can be viewed as a means either to equalize the least significant eigenvalues of the sample covariance matrix or to constrain the array gain. The constrained minimum variance beamforming includes directional constraints [13], derivative constraints [14], quadratic constraints [15], and soft constraints on the norm of the weight vector [16]. In these techniques, robustness to steering vector uncertainty is increased at the expense of a reduction in noise and interference suppression. Recently, some developments based on worst-case performance optimization [17,18] and subspace projection [19][20][21][22] were proposed. The worst-case approaches ensure that the response of the beamformer is above a given level for all steering vectors whose distance to the presumed steering vector is less than a certain distance. The subspace methods estimate the signal-plusinterference subspace to reduce mismatch.
Stochastic methods have also been proposed to tackle the uncertainty of the steering vector [23][24][25][26][27][28]. The most popular of them is Bayesian beamforming which derives from minimum mean square error (MMSE) estimation http://asp.eurasipjournals.com/content/2013/1/108 [25][26][27][28]. In this method, the uncertain steering vector or direction-of-arrival (DOA) is assumed to be a random vector or random variable with a prior distribution that describes the level of uncertainty. The corresponding MMSE estimator can be viewed as a mixture of conditional MMSE estimators combined according to the data-driven posterior distribution function of the steering vector or DOA given the data. More recently, a Bayesian beamforming with order recursive implementation for steering vector uncertainties was proposed in [28], which has the form of a Kalman filter that is recursive in order instead of time.
In this paper, we develop a narrowband beamformer using the Bayesian approach based on [25,26,28]. In this approach, the interference-plus-noise covariance matrix and signal power are assumed to be known, and the steering vector is assumed to be a complex Gaussian random vector that characterizes the level of steering vector uncertainty. Different from the Bayesian beamformer in [28], the proposed method is a time-recursive Bayesian beamformer which can make use of all previous observed data rather than just a recent short-term integration (STI) window and can estimate the steering vector recursively to approach the optimum performance. In comparison with the beamformers of [25,26], the steering vector uncertainty is addressed and modeled as a complex Gaussian random vector in this paper, while the DOA uncertainty is considered in [25,26] and modeled as a random variable. Unlike DOA, random modeling for the steering vector can address not only the uncertainties due to pointing error but also the uncertainties due to scattering around the source, miscalibration, array deformation, different gain, and phase responses of sensors in the array, etc. These three methods all apply the maximum a posteriori estimation for the uncertain DOA or steering vector and have the assumption that the interferers are located far away from the main lobe of the expected beamformer.
The rest of the paper is organized as follows: In Section 2, we present the signal model used in the paper. Section 3 presents the derivation of MMSE estimation and Bayesian beamformer. The recursive implementation of the Bayesian beamformer is developed in Section 4. Numerical simulations are reported in Section 5, and conclusions are given in Section 6.

Signal model
Let us consider an array of N sensors. Assuming narrowband processing, the N×1 complex receiving sensor signal at a snapshot k can be given by where a is the N × 1 complex steering vector of interest, s k is the desired signal with known power E{s k s * k } = σ 2 s .
i k and n k are the N × 1 interference and noise components with known covariance matrix R i+n = E{(i k + n k )(i k +n k ) H }. s k , i k , and n k are assumed to be zero-mean, complex Gaussian random processes that are wide-sense stationary and mutually independent, and their successive snapshots are also statistically independent. In practice, the true steering vector often deviates from its presumed value for various reasons. It is often reasonable to model these errors collectively as a random error vector associated some prior information that is often available in statistical form. Similar to [27,28], we assume that the steering vector a has a complex Gaussian priori probability density function (PDF) with mean a 0 and covariance matrix C 0 , that is, The array output signal is processed in STI windows with K time samples. In each STI, the steering vector a in (1) is assumed to be time-invariant, and the received samples are X j = (x jK , · · · , x (j+1)K−1 ), where j is the index of STI. The goal is to design a beamformer to estimate the signal of interest s j = (s jK , · · · , s (j+1)K−1 ), which is a row vector with length of K.

Bayesian beamformer
In this section, we consider a recursive Bayesian method to estimate the desired signal vector s j with the optimization criterion of MMSE. Given the received samples X 0:j = (X 0 , · · · , X j ), the MMSE estimate of the signal vector is described by the conditional mean of s j given X 0:j , which can be expressed as [25,26] s j = E{s j |X 0:j } = s j p(s j |X 0:j )ds j = s j p(s j |X 0:j , a)p(a|X 0:j )ds j da = p(a|X 0:j )E{s j |X 0:j , a}da.
Based on the assumptions of s k , i k , and n k , we can derive that s k is independent of x l for a different k and l when the steering vector a is known. Then, s j and X 0:j−1 are independent given a. In addition, all elements of s j are independent when X j and a are given. So, we have E{s j |X 0:j , a}= E{s j |X j , a} = E{(s jK , · · · , s (j+1)K−1) )|X j , a} = (E{s jK |X j , a}, · · · , E{s (j+1)K−1 |X j , a}) = (E{s jK |x jK , a}, · · · , E{s (j+1)K−1 |x (j+1)K−1 ,a}).
(4) http://asp.eurasipjournals.com/content/2013/1/108 For arbitrary jK ≤ k < (j + 1)K, the conditional mean E{s k |x k , a} is an optimal MMSE estimator when it is constrained to be linear [25]. That is where E{x k s H k |a} = aσ 2 s is owing to the assumptions of zero mean and the mutual independence of s k , i k , and n k . R x (a) = σ 2 s aa H + R i+n is the data covariance matrix given a. Obviously, it is a spatial Wiener filter with weights [25]. So, we have and the MMSE estimate of the signal vector iŝ It is a beamformer-like data processing with weights and known as the Bayesian beamformer. The Bayesian beamformer is a weighted average of spatial Wiener filters with different steering vectors, which are combined according to the value of the a posteriori probability for each steering vector.

Recursive implementation of the Bayesian beamformer
In this section, we compute the weights of the Bayesian beamformer. According to the Bayesian principle, the posteriori PDF p(a|X 0:j ) can be given by Because successive snapshots of s k , i k , and n k are all statistically independent, it can be seen from (1) that x k is sample independent at different snapshots when a is given. Then, X j and X 0:j−1 are independent given a since they are in the different STI windows. So, we have p(X j |a, X 0:j−1 ) = p(X j |a) and where p(X j |X 0:j−1 ) = p(a|X 0:j−1 )p(X j |a)da is the regularization probability. Assume that the posteriori PDF p(a|X 0:j−1 ) follows a complex Gaussian distribution with mean a j−1 and covariance matrix C j−1 , that is, As stated in [27,28], this Gaussian random model for steering vector is general because it can address the uncertainties due to many reasons, including DOA pointing error, array calibration error, scattering around the source, propagation through an inhomogeneous medium, and other systematic problems. Another reason of this assumption is that the Gaussian distribution can be expressed analytically, and we can obtain the recursive expressions for a j−1 and C j−1 in the following discussions. The likelihood function p(X j |a) can be given by [25] p(X j |a) = The determinant |R x (a)| has the form Expanding R −1 x (a) using the matrix inversion lemma, it yields We note the sample autocorrelation matrix So, the likelihood function is given by where According to the derivation in [26], (16) can be rewritten as where a r is the ideal steering vector, which can be approximated by a j−1 . Calculating this likelihood PDF presents a bigger difficulty. As in [25,26], we assume that the expected projection of the steering vector into the interference subspace is small, or equivalently, the interferers are located far away from the main lobe of the beamformer using the steering vector a r . Therefore, the quadratic functional β(a) can be approximated by the constant β(a j−1 ) [25,26,28], which is defined by The likelihood function can be alternatively approximated as where and γ = α(1 + σ 2 s β(a j−1 )) −K is a normalization factor that ensures the function integrates to one. Substituting (11) and (19) into (10), the posterior PDF p(a|X 0:j ) can be given by It can be seen that the posterior PDF p(a|X 0:j ) is a complex Gaussian with mean a j and covariance matrix C j , which can be recursively expressed by Because β(a) can be approximated by β(a j−1 ), we use (14) and have Substituting (24) into (8), the weights of the Bayesian beamformer are approximated by where a j = ap(a|X 0:j )da.
In summary, we assume that p(a|X 0:j−1 ) is a complex Gaussian PDF with mean a j−1 and covariance matrix C j−1 . The proposed Bayesian beamformer in STI j can be described as follows: 1) ComputeR x,j using (15) and obtainR −1 x,j ; 2) Compute β(a j−1 ) and μ(a j−1 ) using (18) and (20); 3) Compute using (22) and update a j and C j using (23); 4) Compute w MMSE j using (25). http://asp.eurasipjournals.com/content/2013/1/108 At the end of this section, we state again that the signal power σ 2 s and interference-plus-noise covariance matrix R i+n are assumed to be known in the proposed beamformer. However, in practice, the signal power and interference-plus-noise covariance are unknown. In each STI, a simple method for estimating σ 2 s is the minimum variance spatial spectral estimation using a j−1 [25], that is, For the interference-plus-noise covariance matrix, we can consider the output of the array in the absence of the signal of interest and collect a long-term sample covariance estimate for R i+n [3].

Numerical simulation
In this section, we provide numerical illustrations of the performances of the proposed Bayesian beamformer. A uniform linear array with N = 10 omnidirectional sensors spaced half a wavelength apart is considered. Even though arbitrary noise is applicable to the proposed method, for simplicity, we assume that the noise is a white complex Gaussian with power σ 2 n = 1 in the simulations. A desired source signal is from the broadside and generated as a random complex Gaussian process with power σ 2 s . There are  two interferers with the interference-to-noise ratio (INR) 10 dB and 30°and 60°away from the desired signal. Similar to [27], we define the uncertainty ratio (UR) and the signal-to-noise ratio (SNR), respectively, as UR = 10 log 10 tr{C 0 } a H r a r where tr{·} denotes the matrix trace. For algorithm initialization, we set C 0 = 0.001I, and a 0 is drawn from a complex Gaussian distribution with mean a r and covariance matrix C 0 . The performance evaluation criterion is the output signal-to-interference-plus-noise ratio (SINR) SINR(w) = 10 log 10 In all experiments, the STI block length is set to be K = 512, and the optimum Max-SINR beamformer with the true steering vector a r is provided as reference whose weight is w MaxSINR = R −1 i+n a r . In the first experiment, we investigate the convergence of the proposed method. Method of [28] Method of [25] Method of [26] Proposed method MaxSINR Method of [28] Method of [25] Method of [26] Proposed method MaxSINR Method of [28] Method of [25] Method of [26] Proposed method MaxSINR Method of [28] Method of [25] Method of [26] Proposed method MaxSINR Method of [28] Method of [25] Method of [26] Proposed method MaxSINR (e) Method of [28] Method of [25] Method of [26] Proposed method MaxSINR  Figure 1 shows the output SINR versus STI index for different URs under different SNRs. It can be seen that the proposed recursive Bayesian beamformer has good convergence performance. With the increase of time, the output SINR converges to that of the Max-SINR beamformer for different URs and SNRs. In the case of high SNR and small UR, the convergence is very fast; otherwise, the proposed Bayesian beamformer has slow convergence speed.
For comparison purposes, we also display the performances of the linearly constrained minimum variance (LCMV) beamformer [29], subspace projection beamformer [22], and other three Bayesian beamformers proposed in [25,26,28], where the beamformers of [22,28,29] are non-recursive STI block-based methods and the beamformers of [25,26] are recursive. For the LCMV beamformer, the weight is x,j a 0 , where we constrain the weight to satisfy a H 0 w LCMV = 1 and do not use the diagonal loading. For the subspace projection beamformer, the time-varying (on the STI scale) weight vector is computed asP j a 0 , whereP j = I − U d,j U H d,j is the perpendicular projection matrix for the interference subspace, and U d,j contains normalized eigenvectors corresponding to the two largest eigenvalues in the decom-positionR x,j = U j U H j such that U j =[ U d,j |U s+n,j ]. The Bayesian beamforming weight of [27] is σ 2 s (R x,j + Kσ 2 s C 0 ) −1 a 0 . For the Bayesian methods in [25,26], they cannot address the uncertainties due to some systematic problems such as array calibration error and drift in sensor gains or phases, so the DOA uncertainty is considered here. The source DOA is modeled as a random variable with prior statistics, and the a priori PDF of DOA is fixedly set to be uniform at 81 evenly spaced points over [ −10 • : 10 • ] in all simulations. Five hundred Monte Carlo trials were run for different SNRs and URs. Figure 2 shows the output SINR versus SNR for different beamformers under different URs, where the results are obtained after 80,000 STI frames. The beampatterns of different beamformers from one trial is shown in Figure 3, where SNR = 10 dB, UR = −10 dB, and the STI index is 100. The DOAs of the desired signal and interferers are indicated by the vertical solid and dashed lines, respectively.
From Figures 2 and 3, we can observe that the proposed recursive Bayesian beamformer outperforms the other methods and has similar performance to the optimal Max-SINR beamformer. Compared to the LCMV, subspace projection, and Bayesian method of [28], the proposed method produces higher output SINR and better beampattern shape. This is because the proposed beamformer can exploit information about the steering vector contained in past STI windows, while the LCMV, subspace projection, and method of [28] only simply use the information in the current STI window. Compared to the recursive Bayesian beamformers of [25,26], our method shows some improvements in terms of output SINR and beampattern shape, especially in the case of high SNR. The reason is that we use (25) to calculate the beamformer weights where R i+n is assumed to be known, while this assumption is not necessary for the beamformers of [25,26]. In addition, from Figure 2, we can see that the differences between our proposed beamformer and the beamformers of [25,26] are not really significant until some higher SNRs are encountered. The reason is that there is less information about the steering vector of interest contained in the received data when the SNR is much lower. These three Bayesian beamformers have similar performance after convergence, and the SINR improvement of our method is insignificant in the case of lower SNRs.
Finally, we assess the performance of the proposed beamformer when the signal power and the interferenceplus-noise covariance matrix are inaccurate. As stated in Section 4, we use (26) to estimate σ 2 s and collect a longterm sample covariance estimate for R i+n considering the output of the array in the absence of the signal of interest. Besides the output SINR, the vectorial angle error between w MMSE j and w MaxSINR is introduced to evaluate the accuracy of beamforming weight, which is defined by where abs{·} denotes the absolute value of a complex number. When this value is close to 0, the proposed beamforming weight in (25) approaches to the Max-SINR beamformer weight. Figure 4 is shown to assess the performance of the proposed beamformer when the signal power is estimated according to (26), where the interferenceplus-noise covariance matrix is known exactly and SNR = 0 dB, UR = 20 dB. Figure 4a is the output SINR versus STI index for the Max-SINR beamformer, the proposed beamformer with true signal power, and estimated signal power. Figure 4b compares the vectorial angle error between w MMSE j and w MaxSINR when the signal power is known and estimated by using (26).
In Figure 5, some results are presented to assess the performance of the proposed beamformer when http://asp.eurasipjournals.com/content/2013/1/108 the interference-plus-noise covariance matrix is inaccurate, where the signal power is known exactly and SNR = 0 dB, UR = 20 dB. In simulations, some longterm samples of the array in the absence of the signal of interest are collected offline to estimate the interference-plus-noise covariance matrix. Here, the length of samples are L = 1 × K, 10 × K, and 100 × K, respectively, where K = 512 is the length of STI. Figure 5a is the output SINR versus STI index for the Max-SINR beamformer and the proposed beamformer with true INCM and different estimated INCMs, where the INCM is the abbreviation of 'interference-plus-noise covariance matrix' . Figure 5b shows the vectorial angle error between w MMSE j and w MaxSINR .
From Figures 4 and 5, it can be seen that the proposed Bayesian beamformer with estimated signal power using (26) has similar performance to the beamformer with the true signal power. For the effect of the interferenceplus-noise covariance matrix, when the sample length is not long enough, the estimation of R i+n is inaccurate, which will result in the inaccuracy of w MMSE j and performance degradation of the proposed algorithm. With the increase of sample length, the estimation of R i+n is more accurate, and the proposed beamformer has better SINR evaluation and smaller vectorial angle error. Compared to the inaccuracy of interference-plus-noise covariance matrix, the signal power estimation error affects the performance of the proposed beamformer more slightly, which means that (26) can be used to estimate the signal power in practice. In other words, the accuracy of R i+n is predominant in the evaluation of the proposed beamformer. Figure 6 shows the overall performance of the proposed beamformer when signal power is estimated using (26) and 100 × K samples without the signal of interest are adopted to estimate R i+n . From Figure 6, we can see that there is a small performance degradation compared to the proposed beamformer with the true signal power and interference-plus-noise covariance matrix.

Conclusions
In this paper, a recursive Bayesian approach is proposed to mitigate uncertainty in the steering vector for narrowband beamforming. By assuming that the steering vector is a complex Gaussian random vector, the beamformer can be viewed as a mixture of conditional MMSE estimates weighted by the posterior PDF of the random steering vector. To make use of the information about a steering vector contained in past STI windows, a recursive algorithm is developed to estimate the posterior PDF of the steering vector. Simulation results show better performance of the proposed beamformer compared with the LCMV, subspace projection, and other three Bayesian beamformers. After convergence, it has similar performance to the optimal Max-SINR beamformer with true steering vector. A future direction of research consists of generalizing this algorithm to consider the unknown signal power and unknown interference-plus-noise covariance matrix.