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-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-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 ∈ ×1 at time t can be expressed as: where (t) ∈ ×1 is the vector of incident signals at time t, (t) ∈ ×1 is the vector of additive noises, the steering matrix (θ) = [ (θ1) (θ2) ⋯ (θ) ] is the ideal array response and: Here, operator (·) stands for transpose, θ1, θ2, ⋯ , θ 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 = 1, ⋯ , M denotes the deterministic unknown gain and phase error of sensor i. The vector γ = [γ1,·⋯ , γ] is constructed out of the diagonal elements of matrix Γ(γ). For notational convenience, we omit the arguments of Γ(γ) and (θ) sometimes.

Two assumptions need to be made. Firstly, signal (t) is a temporally-complex white Gaussian random vector with mean zero, and its covariance matrix has full rank p. Secondly, noise (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 and noise subspace can be easily obtained from the eigenvalue decomposition of array output covariance matrix: where E{·} and (·) denote the statistical expectation and the complex conjugate transpose, respectively. The eigenvalues are ordered, such as λ1 ≥ λ2 ≥ ⋯ ≥ λ > λ+1 = ⋯ = λ. The corresponding signal subspace and noise subspace are given by = [ 1 ⋯ ] and = [ +1 ⋯ ].

The focus of this paper is the estimation of spatial signature matrix = Γ(γ)(θ) from the N snapshots of the array output. One ambiguity for this problem can be observed between the unknown signal vector (t) and (i.e., 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 in this section. The Vandermonde structure of the ideal array steering matrix (θ) 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 and elements of the array, extracting a subarray from the array is equivalent to picking up rows of matrix . 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 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 − 1), l denotes that the i–th element of the array is used l 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 can be partitioned into two subarrays. where:

Γ, Γ and 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 is p, i.e., Rank(Γ) = p. It is clear that the matrix Γ spans the same space as the signal subspace . The relationship implies that there exists a nonsingular p × p matrix satisfying: Combining Equations (5) and (11) leads to: where diagonal matrix Φ̄ = Φ. Eliminating in Equation (12) yields to: where the q × q diagonal matrix and Ξ = −1Φ̄. 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 ∥ · ∥ 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: where , 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 . 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 Γ and Γ can be expressed as: and: Recall ; 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.

Estimate the signal subspace from the eigendecomposition of the sample covariance matrix , where N is the finite number of snapshots.

Calculate and according to Equation (11), then determine an estimate of vector γ̄ from Equation (16).

Calculate the matrix Ξ according to Equation (15), then extract the rotational DOA θ̂, i = 1, ⋯ , p.

Construct matrix Γ(γ̂), and (θ̂), then the spatial signature matrix is estimated as ^ = Γ(γ̂)(θ̂).

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 by a m × M selection matrix . The full row rank matrix 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 1 is a nonsingular p × p matrix. Considering n subarrays, we have: where ε denotes the distance between the (i + 1)–th subarray and the reference subarray; the matrix ¯ is formed by the selection matrix and , , and Γ and 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 is a Hermitian positive definite weighting matrix and γ̄ and γ͂, respectively, denote the real part and image part of γ. Since the signal subspace is usually substituted for its estimated value , 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 Φ = Φ and Γ = Γ are used.

This is a nonlinear least squares problem for unknown vector μ (see [29]). Solving by minimizing the cost function Ω, we can obtain: Substituting Equation (28) back to Equation (27) yields: 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 (ΨΨ) is nonsingular.

The method considered herein requires a multidimensional search over parameter vector μ ∈ 2(+−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 μ denotes the estimate at k—th iteration, ξ is an iterative step length and the gradient vector and the Hessian matrix are evaluated at μ.

We define a vector 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 μ, the i-th element of the gradient vector is given by: where . 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 ) 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 −−1 to be a decent direction. The other is that we only need to calculate the first derivative of instead of the second one. Furthermore, a useful modification to ill-conditioned cases is to usually use ( + ζ) in lieu of , where the diagonal loading factor ζ is usually small.

Next, we focus on the first derivative of vector in Equation (32). By the aid of the differentiation of the projection matrix [29], can be expressed as: where ψ† = ψ(ψψ)−1 and . The expression of ψ 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 if γ is not included in matrix Γ, otherwise, and , where is a matrix with one in the kk–th position and zeros elsewhere, and γ is the k–th element at the main diagonal of matrix Γ. Moreover, and:

From the above analysis, the vector can be conveniently calculated.

This multidimensional search procedure is briefly described as follows: where 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.

Estimate the signal subspace .

Initialize vector μ Equation (23) utilizing our first method or ESPRIT-like [22].

Update μ according to Equation (31) fromthe estimates of vector γ̂ and θ̂, i = 1, ⋯ , p.

Calculate matrix Γ(γ̂) and (θ̂); estimate spatial signature matrix = Γ(γ̂)(θ̂).

Evaluate:

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 . The identifiability of from signal subspace can be found in [22].

Remark B: It is an intractable problem how to choose weighting matrix in Equation (22). Refer to the weighted subspace fitting (WSF) method [31]; we define: where the diagonal matrix Λ̂ contains the so-called signal eigenvalues estimated from the array output sample covariance matrix and denotes the consistent estimate of noise power. Simulation results using the proposed algorithm show that = 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 matrix , which takes NM2 flops. The process of obtaining signal subspace costs (M2p) flops.

For the first method, the calculation of in Equation (16) requires p(q − p/3) operations if QR-decomposition using Householder transformations is performed, while the calculation of the minimum eigenvalue of takes (q2) flops. Another main computational load is primarily due to the eigendecomposition of Ξ̂ in Equation (15), and it requires about p2q + (p3) flops. The above scheme takes about q2p + 2p2q − p3/3 + (p3) flops. When q ≫ p, our first method requires NM2 + q2p + (M2p) flops.

The main computational load of our second method is in the Gauss–Newton iterative process Equation (31) and the first derivative of vector Equation (36). The computation of needs about m3n3 + 3m2n2p operations, and the inverse matrix of in Equation (31) costs approximately (8(M + p − 1)3) flops. If M = mn and M ≫ p, the above two steps take M4 + 8M3p operations per iteration. Then, our second method requires about NM2 + (M4 + 8M3p) flops.

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.

Algorithm	Computational Complexity
ESPRIT-like	NM2 + (2M2p)
Method 1	NM2 + q2p + (M2p)
Method 2	NM2 + (M4 + 8M3p)

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.10e10°, 0.90e−5°, 1.25e20°, 0.80e−9°, 0.96e15°, 1.18e−23°, 0.88e−2° and 0.85e4°. 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 , where and , 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 and and 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.

Example 1. Value of the Stopping Criterion Versus Iteration Number

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.

Figure 1

Value of the stopping criterion versus iteration number with different initial values. The number of snapshots N = 500, and SNR = 15 dB.

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 = 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.

Figure 2

RMSE of the spatial signature matrix estimation versus SNR. The number of snapshots N = 500.

Example 3. Performance Versus Different

In the third example, we give the simulation results of the multidimensional search procedure (Method 2) with different weighting matrix = opt, and = . 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 opt is clearly demonstrated by Figure 3. The introduction of weighting matrix 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.

Figure 3

RMSE of the spatial signature matrix estimation versus SNR with different . The number of snapshots N = 500.

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 = opt. As depicted in Figure 4, our proposed two algorithms perform excellently for correlated signals. However, because the weighting matrix is used in Method 2, the multidimensional search algorithm is superior to Method 1.

Figure 4

RMSE of the spatial signature matrix estimation versus correlation factor ρ. The number of snapshots N = 1000 and SNR = 15 dB.

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 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].