Spatial Signature Estimation with an Uncalibrated Uniform Linear Array

In this paper, the problem of spatial signature estimation using a uniform linear array (ULA) with unknown sensor gain and phase errors is considered. As is well known, the directions-of-arrival (DOAs) can only be determined within an unknown rotational angle in this array model. However, the phase ambiguity has no impact on the identification of the spatial signature. Two auto-calibration methods are presented for spatial signature estimation. In our methods, the rotational DOAs and model error parameters are firstly obtained, and the spatial signature is subsequently calculated. The first method extracts two subarrays from the ULA to construct an estimator, and the elements of the array can be used several times in one subarray. The other fully exploits multiple invariances in the interior of the sensor array, and a multidimensional nonlinear problem is formulated. A Gauss–Newton iterative algorithm is applied for solving it. The first method can provide excellent initial inputs for the second one. The effectiveness of the proposed algorithms is demonstrated by several simulation results.


Introduction
In wireless communication, a multiple antenna array at base stations is usually used for combating fading, suppressing co-channel interference and increasing capacity [1,2]. If the spatial signatures are exactly known, signals from different users can be separated at the base station array. However, the user spatial signatures are usually unknown in practice and, therefore, must be estimated.
A natural approach to spatial signature estimation makes use of training sequences transmitted by each user. This is a time-consuming method, and the problems of information transmission rate and synchronization are difficult to overcome. As a result, researchers recently have been paying attention to blind spatial signature estimation techniques. In [3], a blind approach based on parallel factor (PARAFAC) analysis is proposed for spatial signature estimation. This method resorts to the assumption of the diagonal covariance matrix of the signals (i.e., uncorrelated signals) to formulate the PARAFAC model. The time-varying user power loading is exploited to identify this model. The strict requirements for signal and user power can limit practical applications of this method.
Unlike the technique proposed in [3], the most common one utilizes the parametric model of spatial signature. This method is involved in determining the directions-of-arrival (DOAs) of incident signals. For an ideal array output model, many high-resolution methods, such as multiple signal classification (MUSIC) [4] and estimation of signal parameters via rotational invariance techniques (ESPRIT) [5], can achieve high estimating accuracy [6]. However, sensor perturbations (e.g., mutual coupling, gain/phase errors and sensor position uncertainty) are not easily eliminated or compensated [7]. The performance of the subspace-based methods will be severely degraded by these model errors [8,9], and we know from [10][11][12] that model errors also affect the Cramér-Rao lower bound. Therefore, sensor perturbations must be taken into account for renewing the capability of these approaches. In other words, once the incident signal DOAs and the unknown model errors are obtained, the problem of spatial signature or channel estimation is solved.
Two parametric errors (i.e., statistical and deterministic) [13,14] are usually used to extend the ideal array model. This paper only focuses on the DOA-based spatial signature estimation method under deterministic unknown gain and phase errors. The term deterministic here implies that the unknown error parameter at each element of the array is a complex stable constant during the period of observation.
A way to compensate for the unknown model errors is to make use of signal sources in known directions. The methods for calibrating sensor gain and phase errors are proposed in [15,16]. A combined calibration of gain, phase and mutual coupling, as well as sensor position perturbations using the maximum likelihood approach has been studied in [7]. All of the above calibration methods require knowledge of exact source localizations or the covariance matrix of the ideal array response. This is a practical and valid approach when the model errors keep constant with time. Meanwhile, several so-called auto-calibration approaches have also been developed in the literature. Here, the auto-calibration indicates that array calibration can be accomplished without employing any dummy elements or transmitters at known directions. When the first and second order statistic model of sensor errors is given, an optimal maximum a posteriori (MAP) estimator is presented [14,17,18]. The problem of DOA estimation is taken apart from optimization function, and alternative iteration is avoided. Despite such techniques based on MAP being called auto-calibration methods, they need the prior distribution of array perturbation parameters.
In [19], an iterative eigenstructure-based technique of estimating DOAs in the presence of unknown gain and phase errors is presented, and this method can apply to arbitrary array geometries, except a uniform linear array (ULA). For the ULA, a phase ambiguity exists between the diagonal error matrix and the ideal array steering matrix (see [20][21][22]). Thus, it is impossible to estimate DOA and gain and phase errors simultaneously. A Hermitian Toeplitz structure of the ideal ULA output covariance matrix is exploited in [23,24]. This method takes advantage of the elements' equivalence at every diagonal line of covariance matrix to form two equations, and the sensor gain and phase error vectors can be separately obtained from these two equations. The main drawback of the approach is that the Toeplitz matrix assumption is only established under infinite sampling conditions. This means that if the number of snapshots is fixed, the performance of the algorithm does not improve after increasing the signal-to-noise ratio (SNR) to a certain point.
In addition, ESPRIT can also be extended to this case. The idea firstly emerged in [15] and is applied for estimating spatial signatures in [22,25]. In [22], a closed-form ESPRIT-like method is proposed, and the maximum overlapping two subarrays are selected from the array to construct an estimator. This method is also used for DOA estimation with a partly-calibrated ULA [26]. However, it has been pointed out that any two subarrays' configuration may be inherently suboptimal [27]. This provides the motivation for the estimation strategy that fully exploits the multiple invariances in the array. In fact, a multiple invariance (MI) ESPRIT has been proposed for DOA estimation under an ideal array output model [28]. Because model errors exist, the MI-ESPRIT method will fail in our problem.
In this paper, we consider the problem of spatial signature estimation with the ULA in the presence of unknown gain and phase errors. As is well known, the incident signal DOAs can only be determined within an arbitrary rotation angle in this case. However, this rotational ambiguity has no impact on the identification of spatial signature [22]. The sensor gain and phase errors are different from mutual coupling among elements of the array, and it cannot be affected by other sensors. In addition, the ideal steering matrix of the ULA has a Vandermonde structure. The two properties provide the possibility of selecting subarrays from the uncalibrated array.
In practice, two eigenstructure-based algorithms are presented in this paper. The first one extracts two subarrays from the ULA to construct an estimator, and a closed-form solution is obtained. The difference from the ESPRIT-like method is that one element of the ULA can be used several times in one subarray. The second method fully exploits the multiple invariances in the ULA. It is natural to expect that an estimator formulated by multiple invariances in the array performs better than the one involving only one invariance (e.g., the ESPRIT-like method). In this method, a multidimensional procedure of minimizing a cost function is presented. The Gauss-Newton algorithm [29] is suggested to solve the multidimensional search problem. The choice of the initial values is essential for the global convergence of the iterative method. Fortunately, the excellent initial inputs can be provided by our first method. Our algorithms firstly estimate the rotational DOAs (not absolute DOAs) and error parameters, then obtain the estimation of the spatial signature. The effectiveness of the proposed algorithms is compared with the ESPRIT-like method [22] by several computer simulations. Furthermore, unlike the method proposed in [3], the subsequent use of these two methods, as done here, does not depend on the assumption of uncorrelated incident signals and time-varying user power.

Data Model
A uniform linear array with M sensors receiving narrowband signals from p far-field sources and the array output vector y ∈ C M ×1 at time t can be expressed as: where s(t) ∈ C p×1 is the vector of incident signals at time t, n(t) ∈ C M ×1 is the vector of additive noises, the steering matrix A(θ) = a(θ 1 ) a(θ 2 ) · · · a(θ p ) is the ideal array response and: Here, operator (·) T stands for transpose, θ 1 , θ 2 , · · · , θ p are the directions-of-arrival of incident signals and d and λ represent the distance between two consecutive sensors and the identical wavelength for all signals, respectively. The diagonal matrix Γ(γ) is given by: and γ i , i = 1, · · · , M denotes the deterministic unknown gain and phase error of sensor i. The vector γ = [γ 1 , · · · , γ M ] T is constructed out of the diagonal elements of matrix Γ(γ). For notational convenience, we omit the arguments of Γ(γ) and A(θ) sometimes. Two assumptions need to be made. Firstly, signal s(t) is a temporally-complex white Gaussian random vector with mean zero, and its covariance matrix R ss has full rank p. Secondly, noise n(t) is a temporally and spatially complex white Gaussian random vector with mean zero and uncorrelated with incident signals. Then, the so-called signal subspace E s and noise subspace E n can be easily obtained from the eigenvalue decomposition of array output covariance matrix: where E{·} and (·) H denote the statistical expectation and the complex conjugate transpose, respectively. The eigenvalues are ordered, such as λ 1 ≥ λ 2 ≥ · · · ≥ λ p > λ p+1 = · · · = λ M . The corresponding signal subspace E s and noise subspace E n are given by E s = e 1 · · · e p and E n = e p+1 · · · e M .
The focus of this paper is the estimation of spatial signature matrix V = Γ(γ)A(θ) from the N snapshots of the array output. One ambiguity for this problem can be observed between the unknown signal vector s(t) and V (i.e., V s(t) = αV · ( 1 α s(t)) for an unknown non-zero scaling α). A reasonable constraint for solving this scaling ambiguity is to let the first element of diagonal matrix Γ(γ) be equal to one.

Estimation Algorithms
Two new subspace-based approaches are presented for estimating spatial signature matrix V in this section. The Vandermonde structure of the ideal array steering matrix A(θ) provides the opportunity for optimally exploiting one invariance or multiple invariances in the ULA, even though the gain and phase errors exist. Since there is a one-to-one relationship between rows in matrix V and elements of the array, extracting a subarray from the array is equivalent to picking up rows of matrix V . We make use of these particular properties to construct our estimators. Here, we have to point out that the one-to-one relationship in rows of V and the elements of the array is not held under mutual coupling errors. The reason is that the response of one element is affected by other elements of the array.

A Novel ESPRIT-Like Method
Define a q × (M − ε) selection matrixĴ consisting of zeros and ones, where q = M − ε + (l i − 1), l i denotes that the i−th element of the array is used l i times in one subarray, and the distance between the subarrays is ε times the element spacing. Only one element at every row is equal to one, and it corresponds to the selected sensor of the array. However, the number one at each column inĴ indicates the number of repetitions of selected element in one subarray.
Then, the spatial signature matrix V can be partitioned into two subarrays. where: Γ x , Γ y and A x can be calculated by the following equations: Here, the operator diag{·} denotes forming a diagonal matrix by taking the given vector as the main diagonal. The displacement diagonal matrix Φ is given by: Here, we give an example to illustrate how to choose the selection matrix. Consider a ULA with six elements; a possible selection matrixĴ 9×5 can be expressed as: In this example, the third element of the array is used three times, while the second and fourth one are selected twice in one subarray.
The determinant of the matrix Γ(γ) is given by: then the rank of the matrix V is p, i.e., Rank(ΓA) = p. It is clear that the matrix ΓA spans the same space as the signal subspace E s . The relationship implies that there exists a nonsingular p × p matrix T satisfying: Combining Equations (5) and (11) leads to: where diagonal matrixΦ = Φ ε . Eliminating A x in Equation (12) yields to: where the q × q diagonal matrixΓ(γ) = Γ −1 y Γ x and Ξ = T −1Φ T . The vectorγ can be constructed from the diagonal elements of matrixΓ(γ).
Similar to the ESPRIT-like method in [22],Γ(γ) and Ξ can be estimated from the least squares problem as follows:Γ where · F denotes the Frobenius norm. The solution for Ξ is given by: Substituting Equation (15) back to Equation (14) and after some manipulation, the least squares problem becomes:γ = arg min γγ and stands for the element-wise product of two matrices. From the above equation, an estimate of vectorγ may be obtained from the eigenvector associated with the smallest eigenvalue of the matrix Π ⊥ Ex E y E H y T . Here, we note that the constraint may be added to matrix Ξ, and some discussions about this question can be found in [22]. Consider the example in Equation (9), the two diagonal matrices Γ x and Γ y can be expressed as: and: RecallΓ(γ) = Γ −1 y Γ x ; we can see that once an estimate of vectorγ is obtained from Equation (16), it is convenient to obtain the elements of vector γ. Furthermore, the matrix Ξ can be obtained by substituting the matrixΓ(γ) reconstructed fromγ back to Equation (15).
The proposed algorithm now is briefly outlined below.
1. Estimate the signal subspaceÊ s from the eigendecomposition of the sample covariance matrix where N is the finite number of snapshots. 2. Calculate E x and E y according to Equation (11), then determine an estimate of vectorγ from Equation (16). 3. Calculate the matrix Ξ according to Equation (15), then extract the rotational DOAθ i , i = 1, · · · , p. 4. Construct matrix Γ(γ) and A(θ), then the spatial signature matrix is estimated aŝ V = Γ(γ)A(θ).

An MI-ESPRIT-Like Method
Assume that the array comprises n identical subarrays of m sensors. Obviously, overlapping subarrays make M ≤ mn. In addition, it has been pointed out that extracting a subarray from the array is equivalent to picking up m rows of matrix V by a m × M selection matrix J i . The full row rank matrix J i consists of zeros and ones. Only one element at every row is equal to one, and it corresponds to the selected sensor of the array.
Similar to Equation (11), for the i-th subarray, we have: where T 1 is a nonsingular p × p matrix. Considering n subarrays, we have: where ε i d denotes the distance between the (i + 1)−th subarray and the reference subarray; the matrix J is formed by the selection matrix J i andJ = J T 1 J T 2 · · · J T n T , E si , and Γ i and A 1 can be calculated by the following equations: Analogous to the MI-ESPRIT method in [28], the unknown parameter vector µ can be obtained from the following least squares problem: and: where W is a Hermitian positive definite weighting matrix andγ i andγ i , respectively, denote the real part and image part of γ i . Since the signal subspace E s is usually substituted for its estimated valueÊ s , the unitary matrix Φ in Equation (8) becomes: . . .
Next, for ease of notation, we define: and: The expression in Equation (22) is reformulated as: Here, the equations Φ ε i T = Φ ε i and Γ i T = Γ i are used.
This is a nonlinear least squares problem for unknown vector µ (see [29]). Solving T T 1 by minimizing the cost function Ω, we can obtain:T Substituting Equation (28) where: Note that it is reasonable to involve all of the elements of the array to estimate unknown parameters. The limiting case is that there are no overlapping subarrays (i.e., M = mn). Thus, the matrix Ψ has full row rank, except that some elements of the array are not always selected (i.e., M < mn). Therefore, the matrix (ΨΨ H ) is nonsingular.
The method considered herein requires a multidimensional search over parameter vector µ ∈ R 2(M +p−1)×1 . It is well known that the Gauss-Newton algorithm [29,30] is a classical way to obtain the solution of this problem. Considering the cost function Ω in Equation (29), the unknown parameters can be iteratively obtained by: where µ k denotes the estimate at k−th iteration, ξ k is an iterative step length and the gradient vector Q = ∂Ω ∂µ and the Hessian matrix H = ∂ 2 Ω ∂µ∂µ T are evaluated at µ k . We define a vector r formed by columns scanning matrix ΥΠ ⊥ Ψ , i.e., Then, the cost function Ω becomes: In the following, based on the above equation, the expressions for the gradient vector and an approximate Hessian matrix are presented. Firstly, considering the gradient of Ω w.r.t µ i , the i-th element of the gradient vector Q is given by: where r i = ∂r ∂µ i . The ij-th component of the Hessian matrix may be expressed as (see [30,31]): The reasons for this approximation (i.e., ignoring the second derivative of r) can be found in [28,30], and it brings two benefits for the Gauss-Newton iterative algorithm. One is that the approximate Hessian matrix is positive semidefinite and it guarantees −H −1 Q to be a decent direction. The other is that we only need to calculate the first derivative of r instead of the second one. Furthermore, a useful modification to ill-conditioned cases is to usually use (H + ζI) in lieu of H, where the diagonal loading factor ζ is usually small.
Next, we focus on the first derivative of vector r in Equation (32). By the aid of the differentiation of the projection matrix [29], r i can be expressed as: The expression of Ψ i can be obtained from Equations (21), (23), (24) and (26). This is a tiresome calculation process. However, it can be simplified with the observations that where I kk is a matrix with one in the kk−th position and zeros elsewhere, and γ j is the k−th element at the main diagonal of matrix Γ j . Moreover, and: From the above analysis, the vector r i can be conveniently calculated. This multidimensional search procedure is briefly described as follows: 1. Estimate the signal subspaceÊ s . 2. Initialize vector µ Equation (23) utilizing our first method or ESPRIT-like [22]. 3. Update µ according to Equation (31) fromthe estimates of vectorγ andθ i , i = 1, · · · , p. 4. Calculate matrix Γ(γ) and A(θ); estimate spatial signature matrixV = Γ(γ)A(θ).

Evaluate:
whereV k denotes the estimated spatial signature matrix at the k−th step. If it is smaller than a pre-set tolerance, terminate the iterations; if not, repeat Steps 3, 4 and 5.
Remark A: Although our first method and the ESPRIT-like proposed in [22] can only provide rotational values over true DOA θ and true gain and phase error matrix Γ, it makes no difference to the estimation of spatial signature V . The identifiability of V from signal subspace E s can be found in [22].
Remark B: It is an intractable problem how to choose weighting matrix W in Equation (22). Refer to the weighted subspace fitting (WSF) method [31]; we define: where the diagonal matrixΛ s contains the so-called signal eigenvalues estimated from the array output sample covariance matrix andσ 2 n denotes the consistent estimate of noise power. Simulation results using the proposed algorithm show that W = W opt results in a better performance.
Remark C: The Gauss-Newton algorithm is introduced for solving our nonlinear least squares problem Equation (29). The initial estimate is essential for the global convergence of this iterative method. Fortunately, our first method proposed in Section 3.1 or the ESPRIT-like [22] can provide excellent initial values. Simulation results indicate that only two or three iterative steps are required to achieve a solution close to the global minimum.
Remark D: The computational load of the proposed two methods in this paper is usually dominated by forming the sample covariance matrixR, which takes N M 2 flops. The process of obtaining signal subspaceÊ s costs O(M 2 p) flops.
For the first method, the calculation of Π ⊥ Ex in Equation (16)  The main computational complexities of the proposed two algorithms are compared in Table 1 with that of the ESPRIT-like method. In this table, Method 1 denotes the method proposed in Section 3.1, and Method 2 is the approach proposed in Section 3.2. As can be seen from Table 1, the complexity of Method 2 is higher than that of the ESPRIT-like and Method 1. This is because more unknown parameters need to be estimated in Method 2. However, more accurate results can be obtained from this multidimensional search process. Table 1. Computational complexity of the ESPRIT-like and the proposed two techniques.

Experiments
In this section, computer simulations are conducted for evaluating the performance of the proposed algorithms. The approach proposed in Section 3.1 will be referred to as Method 1, while we refer to the method proposed in Section 3.2 as Method 2. In all scenarios, a ULA of 9 elements with half of wavelength element spacing is used. The deterministic gain and phase errors are given by 1, 1.10e j10 • , 0.90e −j5 • , 1.25e j20 • , 0.80e −j9 • , 0.96e j15 • , 1.18e −j23 • , 0.88e −j2 • and 0.85e j4 • . The ULA is divided into 5 subarrays in Method 2, and the i-th subarray is selected by matrix: However, the two subarrays of Method 1 are selected by matrices: and:Ĵ Two equal-power signals located at 25 • and 30 • are considered, and the SNR per element for each source is defined by SNR = 10 log 10 (σ 2 s /σ 2 n ), where σ 2 s and σ 2 n , respectively, denote the power of the incident signal and that of additive noise at each sensor. RMSE is used as the performance measure and defined by: where K is the number of trials andV i and V i are the estimated spatial signature and the true one at the i−th experiment, respectively. The result of the ESPRIT-like method in [22] is also included. A total of 200 trials are performed for each simulation scenario. Method 2 terminates after 10 iterations in all simulations.  Since our proposed Method 2 requires an iterative optimization, it is of interest to evaluate how the estimates improve with each iteration. In this example, Method 2 is initialized separately by Method 1 and ESPRIT-like, and the two incident signals are uncorrelated. Figure 1 shows the value of the stopping criterion Equation (39) as a function of iteration number. It is obvious that all improvement comes in the first two or three iterations. Furthermore, the initial values provided by Method 1 are better than those of ESPRIT-like. This is because more accurate parameter values are estimated by Method 1, and we will show it in the next example.

Example 2. Performance Versus SNR
In the second example, the proposed two methods and the ESPRIT-like in [22] are compared against the SNR. We assume that the incident two signals are uncorrelated, and the number of snapshots is N = 500. The SNR is varied from −10 to 10 dB. The weighting matrix W = W opt in Method 2 is calculated by Equation (40). Moreover, the initial value is given separately by ESPRIT-like and the Method 1.  The simulation result in Figure 2 shows that our proposed two algorithms perform better than ESPRIT-like, even at low SNR. This means that the maximum overlapping subarray is not the optimal choice in the presence of unknown sensor gain and phase response. In other words, fully exploiting invariance in the array should be taken into account to construct an estimator. In addition, the behavior of the multidimensional search procedure (Method 2) has been improved when the initial parameter values are provided by our Method 1. Simulation results also indicate that our method only requires 2 to 4 Gauss-Newton iterations, and the main performance improvement comes from the first two iterations.

Example 3. Performance Versus Different W
In the third example, we give the simulation results of the multidimensional search procedure (Method 2) with different weighting matrix W = W opt , and W = I. The two sources are also uncorrelated, and the snapshots are fixed at N = 500. Method 2 is only initialized by the ESPRIT-like method. The advantage of using the weighting matrix W opt is clearly demonstrated by Figure 3. The introduction of weighting matrix W can improve the performance of the proposed algorithm especially under low SNR. Meanwhile, this example also shows that an algorithm that fully exploits the physical structure of the array can give better performance.

Example 4. Performance Versus Correlation between Signals
In the last one, we test the estimation performance in terms of the correlation between the two incident signals, where the correlation factor is denoted by ρ. The number of snapshots is N = 1000 and SNR is fixed at 15 dB. The generations of the two correlated signals can be found in [32]. The inputs of the multidimensional search process are given by Method 1 and the weighting matrix W = W opt . As depicted in Figure 4, our proposed two algorithms perform excellently for correlated signals. However, because the weighting matrix W is used in Method 2, the multidimensional search algorithm is superior to Method 1.

Conclusions
The ESPRIT-like method proposed in [22] just makes use of the maximum overlapping subarray of ULA with unknown gain and phase errors to formulate an estimator. However, there exist many invariances in this kind of ULA. Thus, we firstly presented a simple eigenstructure-based algorithm for spatial signature estimation. One element of the ULA can be selected several times in one subarray, while the same element can only be used once in one subarray of the ESPRIT-like method. Another multidimensional search algorithm was also proposed. This method fully excavates the multiple invariances of the ULA and is implemented with the help of the Gauss-Newton iterative algorithm. Because excellent initial parameter values can be provided by the ESPRIT-like or our first method, this algorithm can converge rapidly to an appropriate solution. On the other hand, the introduction of a weighting matrix W further improves the performance of the multidimensional search algorithm. However, how to select subarrays from ULA, making the algorithm lead to optimal parameter estimation, is not addressed in this paper. Some pragmatic discussion about this question can be found in [27].

Author Contributions
Xiang Cao and Jingmin Xin provided the idea of this work. Xiang Cao and Yoshifumi Nishio conceived of and designed the experiments. Xiang Cao performed the experiments and provided all of the figures and data for the paper. Xiang Cao, Jingmin Xin, and Nanning Zheng prepared the literature and analyzed the data. Xiang Cao wrote the paper. Correspondence and requests for the paper should be addressed to Jingmin Xin.

Conflicts of Interest
The authors declare no conflict of interest.