Improvement on EVESPA for Beamforming and Direction of Arrival Estimation

This paper presents an alternative estimation procedure of the generalized steering matrix of the sources in EVESPA, suitable for both beamforming and direction of arrival estimation. It is shown how the estimation of such a matrix can be restricted to that of its corresponding coe ﬃ cient matrix in the signal subspace, providing both performance enhancement and computational complexity reduction. Performance comparison through numerical simulations is presented to conﬁrm the e ﬀ ectiveness of the proposed procedure.


Introduction
The EVESPA algorithm was introduced by Gönen et al. in [1] as an adapted version of VESPA to handle the case of coherent signals. Using fourth-order cumulants, this algorithm properly generates an unbiased estimation of the generalized steering matrix (GSM) B of the sources from which an analysis of the signal parameters can be performed for each coherent group in an independent fashion. The same estimation procedure is also undertaken in [2] as a means of computing the weights of an optimal beamformer capable of maximizing the signal-to-interference plus noise ratio (SINR) in a coherent environment. EVESPA has also been considered in [3] under a slight variation for the problem of direction of arrival (DOA) estimation in mobile communications, and in [4] where it was adapted for a twodimensional scenario. Steps 1 to 7 in [2] describe the EVESPA algorithm where the M × G matrix B is considered in its whole. However, since the latter lies in the signal subspace, only the search of its corresponding G × G coefficient matrix proves sufficient. In this paper, we present a new estimation procedure of B based on this principle which is shown to bring significant performance enhancement both in terms of computational complexity and quality estimation.
Throughout the paper, " * and " † are, respectively, used as the conjugate and Hermitian transpose operators, and nonscalar quantities such as vectors and matrices are labeled in bold.

Background Theory
In this section we recall the main guidelines of the EVESPA algorithm [1,2], as it will be assumed that the reader has a sufficient knowledge on the matter. The signal model used in this paper is the same as that used by the original authors, namely, where x k , B, u k , and n k are, respectively, the M × 1 vector of the received signals, the M × G generalized steering matrix of the sources, the G × 1 elementary sources vector, and the M × 1 additive white Gaussian noise (AWGN) vector, which is assumed symmetric. Elements u g (t k ), g ∈ {1, 2, . . . , G}, of u k are modeled as uncorrelated zero-mean random processes, where G denotes the number of coherent groups 2 EURASIP Journal on Advances in Signal Processing of signals impinging on the M-element array. Without loss of generality, matrix B may be expressed as where A g and α g are the M × p g steering matrix and the p g ×1 coefficients vector of the g-th coherent group. The total number of sources impinging on the array is thus G g=1 p g . EVESPA may be suitable for any type of signals provided that a complex envelope of interest u g (t k ) admits a nonzero fourth-order cumulant, that is, The algorithm proceeds to the evaluation of the two cumulant matrices: from which an estimation B of B is obtained by following steps 1 through 7 in [2]. Upon calculation of B, whose columns match those of B to within scale and permutation, one can freely proceed to beamforming or DOA estimation as described in [1,2]. In the case of DOA estimation, the covariance matrix of the g-th coherent group of signals is first estimated as R xxg = b g b † g , where b g is the g-th column of B. Spatial smoothing is then applied to R xxg in order to estimate the DOAs of that group. Any other DOA estimation algorithm could also be applied once B has been obtained from EVESPA.
In the case of beamforming, an optimum weight vector is initially computed as w g,opt = c g R −1 xx b g , where R xx is the covariance matrix of the received signals, and where is a constant ensuring a unit response in the "look" direction. This beamforming vector is then used to recover the elementary signal u g (t k ) of that particular group. Here again, any other beamforming scheme can be implemented after estimation of B. One interesting example is that of the null steering beamformer [5], where the array response can be made null for all signals of a particular group instead of synthesizing nulls for each of these signals independently.
Note finally that although the EVESPA algorithm has been developed in a context of narrowband multipath propagation, its application remains valid in any environment provided that the received signals can be modeled as in (1), and that the appropriate requirements on B, u k , and n k be met. Those are essentially given by assumptions A1 through A6 in [1], but can be summarized in a more general way as follows.
(ii) B has G ≤ M columns linearly independent from one another.
(iii) The fourth-order cumulant of the noise vector n k is zero.
One interesting example of environment where the EVESPA algorithm may also be applied is that of narrowband nearfield sources [6], where the signal model also finds its correspondence to (1). In this paper, though, we will focus on the description of an enhanced estimation procedure of B with the aim of improving performance for both beamforming and direction of arrival estimation in a coherent narrowband scenario, since each of these subjects were covered in [1,2]. Note however that this estimation procedure is general, and may in fact be applied regardless of the type of subsequent processing.

Proposed Estimation Procedure
Consider the covariance matrix of the received signals: where E{·} denotes the expected value operator. Note that R uu is always diagonal since u k is a vector of zero-mean and uncorrelated elements. Expressing R xx in terms of its eigenvalue decomposition yields where E s represents a set of orthonormal basis vectors of the signal subspace. From (5) and (6), it follows that B can be expressed in terms of E s such that where Q is a G × G coefficient matrix ensured to be full rank under assumptions A1 to A4 in [2]. In the same issue, it is also shown that matrices C 1 and C 2 of (4) evaluate to where Λ and D are both full-rank diagonal matrices. Our alternative estimation procedure begins by forming the two G×G matrices (recall from (7) that E † s B = Q, since E † s E s = I) where E s is obtained from the eigenvalue decomposition (EVD) of R xx . We now apply an estimation procedure similar to steps 2 through 7 in [2], but in the aim of identifying Q. Consider for this the single value decomposition (SVD) of a 2G × G matrix C such that Step : Additional steps required for the covariance based improvement method in [1]. This improvement is applied by default in steps 1 to 4 of the proposed procedure.
where matrix U may be partitioned into and where U n1 and U n2 are both G × G. Since the G last rows of Σ are zeros, it follows that which implies that Hence, the eigenvalues of H must match the diagonal elements of D * . If these elements are distinct, there exists a unique mapping between E H , the eigenvectors matrix of H, and Q − † such that where Z is a scale permutation matrix containing only one non-zero element per line and column. Therefore, the columns of E − † H = QZ match those of Q to within scale and permutation and a straightforward estimation of B becomes Table 1 presents a summary of the main computational steps involved in both the original EVESPA and the proposed algorithm using their respective notations. It can be seen that the proposed procedure requires only one SVD, whereas two are required in the original EVESPA. Moreover, the only M-dependant decomposition involved in the proposed method is that of the initial EVD. Hence, as the number of sensors M increases for a fixed number of coherent groups G, there will obviously come a point where the proposed procedure outperforms the original EVESPA in terms of computational complexity. However, note that the final step of our proposed method involves the inverse of the nondiagonal matrix  E H , which for a minimum value of G represents a higher complexity operation than (25) in [2]. Thus, the gain in computational complexity of our proposed method does not appear evident for low values of M. In order to appreciate the computational complexity of each algorithm, Figure 1 displays their normalized average execution time curves for M ∈ {2, 3, . . . , 20} and M ≥ G ∈ {2, 4, 8, 12}, where each point was obtained from 20000 trials and randomly generated matrices C 1 , C 2 , and R xx . Note that even though the original EVESPA (without improvement) does not make use of R xx , its computation does not increase the complexity of the proposed method since it already represents an intermediate step in the evaluation of both C 1 and C 2 .

Computational Complexity.
No improvement was considered for the original EVESPA. As expected, the performance of the proposed procedure is at its worse for low values of M and G where it is slightly outperformed by the original EVESPA. However, this situation quickly changes as M and G increase where the gain in computational complexity obtained with the proposed procedure becomes obvious. This could also have been predicted from Table 1. All in all, the proposed estimation procedure thus constitutes an improvement in terms of computational complexity over the original one.

Statistical Performance.
Since E s can be estimated from second-order statistics, we now show that the proposed estimation procedure achieves a better statistical performance than the original EVESPA. Consider the complex angle β g (K) between b g , the g-th column of B, and b g (K), the g-th column of B(K) and corresponding estimation of b g Hence, under assumptions A1 to A4 in [2] and for G ≤ M, the MSE of this latter quantity can be computed as follow: Ideally, | cos(β g (K))| = 1 meaning that b g and b g (K) are collinear. Using such a performance criterion, the best algorithm is thus the one that maximizes E{| cos(β g (K))|} for all g and a given K < ∞. Taking the average of (17) over G, a global RMSE criterion can thus be defined as Figure 2 displays the RMSE curves of both the original EVESPA and the proposed algorithm obtained from 10000 runs of 50 snapshots for SNR values ranging from −10 dB to 12 dB. BPSK signals are considered using parameters of Table 2. The receiver consists of a ten-element uniform linear array (ULA) with equal power and independent AWGN on all elements. It can be seen that the proposed estimation procedure achieves a better performance than the original EVESPA for all SNRs, namely, because E s is estimated from second-order statistics which possess a lower variance than fourth-order statistics. Note however that the use of second-order statistics as a means of improving the quality estimate of B was also considered in Section 4 of [1], corresponding to steps 5 and 6 of Table 1. Upon a first estimation B of B, a new estimation B was computed such that B = E s E † s B. The performance of such an estimator would have been similar to that of the proposed estimation procedure in Figure 2. However, the computational complexity involved in a first evaluation of B from the original procedure followed by an additional EVD of R xx in order to compute E s would clearly become higher than that of the proposed algorithm. Hence, the latter does still represent an advantageous alternative in this context.

Conclusion
In this paper, we have shown how the original EVESPA algorithm could be improved both in terms of computational complexity and statistical performance by restricting the estimation of B to that of its corresponding coefficient matrix in the signal subspace. The use of E s as estimated from second-order statistics ensures a gain in statistical performance while the reduced dimensions of C through the use of Q accounts for a gain in computational complexity in a majority of scenarios.