Sensors 2015, 15, 13899-13915; doi:10.3390/s150613899

OPEN ACCESS

sensors
ISSN 1424-8220
www.mdpi.com/journal/sensors
ArticleSpatial Signature Estimation with an Uncalibrated Uniform
Linear Array †Xiang Cao 1 , Jingmin Xin 1 , Yoshifumi Nishio 2, * and Nanning Zheng 1
1

Institute of Artificial Intelligence and Robotics, Xi’an Jiaotong University, Xi’an 710049, China;
E-Mails: cao.xiang@stu.xjtu.edu.cn (X.C.); jxin@mail.xjtu.edu.cn (J.X.);
nnzheng@mail.xjtu.edu.cn (N.Z.)
2
Department of Electrical and Electronic Engineering, Tokushima University,
Tokushima 770-8506, Japan†

This paper is an extended version of our paper published in Sensor Signal Processing for Defence(SSPD 2012), 25–27 September 2012, London.

* Author to whom correspondence should be addressed; E-Mail: nishio@ee.tokushima-u.ac.jp;
Tel.: +81-88-656-7470.
Academic Editor: Assefa M. Melesse
Received: 13 March 2015 / Accepted: 10 June 2015 / Published: 12 June 2015

Abstract: 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 identificationof 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 byseveral simulation results.
Keywords: array signal processing; auto-calibration; parameter estimation

Sensors 2015, 15

13900

1. 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 areexactly 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 byeach user. This is a time-consuming method, and the problems of information transmission rate andsynchronization 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 assumptionof 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], canachieve high estimating accuracy [6]. However, sensor perturbations (e.g., mutual coupling, gain/phaseerrors and sensor position uncertainty) are not easily eliminated or compensated [7]. The performanceof 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 theincident signal DOAs and the unknown model errors are obtained, the problem of spatial signature orchannel 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, theauto-calibration indicates that array calibration can be accomplished without employing any dummyelements 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Sensors 2015, 15

13901

such techniques based on MAP being called auto-calibration methods, they need the prior distributionof 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 andphase 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 ofcovariance 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 ofsnapshots 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, andthe maximum overlapping two subarrays are selected from the array to construct an estimator. Thismethod 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 presenceof unknown gain and phase errors. As is well known, the incident signal DOAs can only be determinedwithin 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 expectthat 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 ofminimizing a cost function is presented. The Gauss–Newton algorithm [29] is suggested to solve themultidimensional search problem. The choice of the initial values is essential for the global convergenceof the iterative method. Fortunately, the excellent initial inputs can be provided by our first method. Ouralgorithms 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

Sensors 2015, 15

13902

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.
2. Data Model
A uniform linear array with M sensors receiving narrowband signals from p far-field sources and thearray output vector y ∈ C M ×1 at time t can be expressed as:y(t) = Γ(γ)A(θ)s(t) + n(t)

(1)

M ×1
where s(t) ∈ C p×1 is the vector ofhincident signals at time t, n(t)is the vector of additivei ∈ Cnoises, the steering matrix A(θ) = a(θ1 ) a(θ2 ) · · · a(θp ) is the ideal array response and:a(θ) =h2π

2π

1 ej λ d sin(θ) · · · ej λ (M −1)d sin(θ)

iT

(2)

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:


γ1




γ2

Γ(γ) = 
(3)
..


.


γM
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 Rss has full rank p. Secondly, noise n(t)is a temporally and spatially complex white Gaussian random vector with mean zero and uncorrelatedwith 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:
H

R = E{y(t)y (t)} =MX

λi ei eH
i

(4)

i=1

where E{·} and (·)H denote the statistical expectation and the complex conjugate transpose, respectively.
The eigenvalues are ordered, such as λ1 ≥ λ2 ≥ · · · ≥ λph > λp+1 = · · · i= λM . The corresponding
signalhisubspace E s and noise subspace E n are given by E s = e1 · · · ep and E n = ep+1 · · · eM .
The focus of this paper is the estimation of spatial signature matrix V = Γ(γ)A(θ) from the Nsnapshots 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.

Sensors 2015, 1513903

3. Estimation Algorithms
Two new subspace-based approaches are presented for estimating spatial signature matrix V in thissection. 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 phaseerrors exist. Since there is a one-to-one relationship between rows in matrix V and elements of thearray, extracting a subarray from the array is equivalent to picking up rows of matrix V . We make useof 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.
3.1. A Novel ESPRIT-Like Method
P
Define a q × (M − ε) selection matrix Jˆ consisting of zeros and ones, where q = M − ε + (li − 1),li denotes that the i−th element of the array is used li 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 Jˆ indicatesthe number of repetitions of selected element in one subarray.
Then, the spatial signature matrix V can be partitioned into two subarrays.
"
#
Γx Ax
J ΓA =(5)
Γ y Ax Φ εwhere:
"
J=

Jˆ1
Jˆ2

#

"
=

Jˆq×(M −ε)
0q×ε
0q×ε
Jˆq×(M −ε)#(6)Γx , Γy and Ax can be calculated by the following equations:Γx = diag{Jˆ1 γ}, Γy = diag{Jˆ2 γ}, Ax = Jˆ1 A(7)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: 2π
ej λ d sin(θ1 )
2π




ej λ d sin(θ2 )


(8)
Φ=
..

.


2π

ej λ d sin(θp )Here, we give an example to illustrate how to choose the selection matrix. Consider a ULA with six
elements; a possible selection matrix Jˆ9×5 can be expressed as:

Sensors 2015, 15

13904



1
 0


 0

 0


ˆ
J 9×5 =  0

 0

 0


 0
0

0
1
0
1
0
0
0
0
0

0
0
1
0
1
0
1
0
0

0
0
0
0
0
1
0
1
0


0
0 


0 

0 


0 

0 

0 


0 
1

(9)

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:det(Γ(γ)) =MY

γi 6= 0(10)i=1then 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:"#Ex
J Es == J ΓAT(11)
Ey
Combining Equations (5) and (11) leads to:E x = Γx Ax T , E y = Γy Ax Φ̄T(12)where diagonal matrix Φ̄ = Φε . Eliminating Ax in Equation (12) yields to:Γ̄(γ̄)E y = E x Ξ

(13)

−1
Φ̄T . The vector γ̄ can be constructed
where the q × q diagonal matrix Γ̄(γ̄) = Γ−1y Γx and Ξ = T
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:ˆ , Ξ̂ = arg min Γ̄(γ̄)E − E Ξ 2
Γ̄
y
x
F
¯
Γ
,Ξ
where k · kF denotes the Frobenius norm. The solution for Ξ is given by:
−1 HΞ̂ = E H
E x Γ̄(γ̄)E y
x Ex

(14)

(15)

Substituting Equation (15) back to Equation (14) and after some manipulation, the least squares
problem becomes:
h
 i
H T
ˆ = arg min γ̄ H Π⊥
γ̄
E
E
γ̄
(16)
y y
Ex
γ̄

Sensors 2015, 15

13905

−1 H
H
where Π⊥
=
I
−
E
E
E
E x , and stands for the element-wise product of two matrices.
xxx
Ex
From the above equation, an estimate of vector γ̄ may be obtained from the eigenvector associated withT
the smallest eigenvalue of the matrix Π⊥
EyEH. Here, we note that the constraint may be addedyEx
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:
Γx = diag {1, γ 2 , γ 3 , γ 2 , γ 3 , γ 4 , γ 3 , γ 4 , γ 5 }(17)Γy = diag {γ 2 , γ 3 , γ 4 , γ 3 , γ 4 , γ 5 , γ 4 , γ 5 , γ 6 }

(18)

and:

Recall Γ̄(γ̄) = Γ−1y Γ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ˆ back to Equation (15).
substituting the matrix Γ̄(γ̄) reconstructed from γ̄
The proposed algorithm now is briefly outlined below.
1. Estimate the signal subspace Ê s from the eigendecomposition of the sample covariance matrix
P
HR̂ = N1 N
n=1 y(n)y (n), 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 asV̂ = Γ(γ̂)A(θ̂).
3.2. 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 consistsof 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:
J i E s = J i Γ(γ)A(θ)T 1
where T 1 is a nonsingular p × p matrix. Considering n subarrays, we have:

 

E s1
Γ1 A1

 

ε
 E s2   Γ2 A1 Φ 1   =  Γ3 A1 Φε2  T 1
EJ¯E s = 
s3

 

 ..   ..
 .   .
E sn

Γn A1 Φεn−1

(19)

(20)

Sensors 2015, 15

13906

where εi d denotes the distance between the (i + 1)−th subarray and the reference subarray; the matrixhiTJ¯ is formed by the selection matrix J i and J¯ = J T J T · · · J T, E si , and Γi and A1 can be12ncalculated by the following equations:
E si = J i E s , Γi = J i Γ(γ)J Ti , A1 = J 1 A(θ)(21)

Analogous to the MI-ESPRIT method in [28], the unknown parameter vector µ can be obtained fromthe following least squares problem:



2
E s1
Γ1 A 1ε E s2  Γ2 A 1 Φ 1 
 1/2 
 W −  Γ3 A1 Φε2  T 1
EΩ(µ) = s3 .. 
 ..

 . 
 .

εn−1
E sn
Γn A 1 Φ
F


(22)

and:
µ = [ γ̄2 · · · γ̄M γ̃2 · · · γ̃Mρ1 · · · ρp θ1 · · · θp ]T(23)

where W is a Hermitian positive definite weighting matrix and γ̄i and γ̃i , respectively, denote the realpart and image part of γi . Since the signal subspace E s is usually substituted for its estimated value Ê s ,the unitary matrix Φ in Equation (8) becomes:2π
ρ1 ej λ d sin(θ1 )
2π




ρ2 ej λ d sin(θ2 )


(24)
Φ=
..

.


2π

ρp ej λ d sin(θp )Next, for ease of notation, we define:
Υ=W

T /2

h

E Ts1

E Ts2

···

E Tsn

i

(25)

and:
Ψ=

h

AT1 Γ1 Φε1 AT1 Γ2 · · · Φεn−1 AT1 Γn

i

(26)

The expression in Equation (22) is reformulated as:2Ω(µ) = Υ − T T1 Ψ F(27)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 T1 by minimizing
the cost function Ω, we can obtain:TT̂ 1 = ΥΨH ΨΨH

−1

(28)

Sensors 2015, 15

13907

Substituting Equation (28) back to Equation (27) yields:2µ̂ = arg min Ω = arg min ΥΠ⊥
Ψ F
µ
µ

(29)

H
H −1
Π⊥
Ψ = I − Ψ (ΨΨ ) Ψ(30)

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 fullrow rank, except that some elements of the array are not always selected (i.e., M < mn). Therefore, thematrix (ΨΨH ) is nonsingular.
The method considered herein requires a multidimensional search over parameter vectorµ ∈ R2(M +p−1)×1 . It is well known that the Gauss–Newton algorithm [29,30] is a classical way toobtain the solution of this problem. Considering the cost function Ω in Equation (29), the unknown
parameters can be iteratively obtained by:µk+1 = µk − ξk H −1 Q

(31)

where µk denotes the estimate at k−th iteration, ξk is an iterative step length and the gradient vector∂2Ωk
Q = ∂∂Ωµ and the Hessian matrix H = ∂ µ∂ µT are evaluated at µ .
We define a vector r formed by columns scanning matrix ΥΠ⊥Ψ , i.e.,no⊥r = vec ΥΠΨ(32)Then, the cost function Ω becomes:Ω = rH r(33)

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 elementof the gradient vector Q is given by:
Qi =


∂Ω
= 2Re r H
i r ,
∂µi

(34)

where r i = ∂∂µr . The ij-th component of the Hessian matrix may be expressed as (see [30,31]):i

H ij =


∂ 2Ω' 2Re r Hi rj∂µi ∂µj(35)

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 usefulmodification to ill-conditioned cases is to usually use (H + ζI) in lieu of H, where the diagonal loading
factor ζ is usually small.

Sensors 2015, 15

13908

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:
 

H 
†
⊥
†
⊥
r i = −vec Υ Ψ Ψi ΠΨ + Ψ Ψi ΠΨ
(36)
Ψ . The expression of Ψ can be obtained from Equations (21),where Ψ† = ΨH (ΨΨH )−1 and Ψi = ∂∂µii(23), (24) and (26).
This is a tiresome calculation process. However, it can be simplified with the observations that
T
εi−1
ε
T
∂ Φ i−1 A1 Γi= 0 if γ is not included in matrix Γ , otherwise, ∂ Φ A1 Γi = Φεi−1 AT I and
j

∂γj

∂Φ

εi−1

i

∂γ̄j

1

kk

T
1

A Γi = jΦεi−1 AT I , where I is a matrix with one in the kk−th position and zeros elsewhere,kk
1 kk
and γj is the k−th element at the main diagonal of matrix Γj . Moreover, εi−1T∂Φεi−1 AT1 Γi∂ΦTεi−1 ∂A1=
A1 + Φ
Γi
(37)
∂θj
∂θj
∂θj
∂γ̃j

and:
∂Φεi−1 AT1 Γi∂Φεi−1 T=
A 1 Γi
∂ρj
∂ρj

(38)

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 Ê 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 matrix V̂ = Γ(γ̂)A(θ̂).
5. Evaluate:kV̂k+1 − V̂k kF(39)

kV̂k kF

where V̂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 befound in [22].
Remark B: It is an intractable problem how to choose weighting matrix W in Equation (22). Referto the weighted subspace fitting (WSF) method [31]; we define:−1

W opt = (Λ̂s − σ̂2n I)2 Λ̂s(40)

where the diagonal matrix Λ̂s contains the so-called signal eigenvalues estimated from the array output
sample covariance matrix and σ̂2n denotes the consistent estimate of noise power. Simulation results
using the proposed algorithm show that W = W opt results in a better performance.

Sensors 2015, 15

13909

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 toachieve a solution close to the global minimum.
Remark D: The computational load of the proposed two methods in this paper is usually dominatedby forming the sample covariance matrix R̂, which takes N M 2 flops. The process of obtaining signal
subspace Ê s costs O(M 2 p) flops.
2
For the first method, the calculation of Π⊥E x in Equation (16) requires p (q − p/3) operations
if QR-decomposition using Householder transformations is performed, while the calculation of theH T
minimum eigenvalue of (Π⊥EE) takes O(q 2 ) flops. Another main computational load isyyExprimarily due to the eigendecomposition of Ξ̂ in Equation (15), and it requires about p2 q + O(p3 ) flops.
The above scheme takes about q 2 p + 2p2 q − p3 /3 + O(p3 ) flops. When q  p, our first method requires
N M 2 + q 2 p + O(M 2 p) flops.
The main computational load of our second method is in the Gauss–Newton iterative process
Equation (31) and the first derivative of vector r Equation (36). The computation of r i needs
about m3 n3 + 3m2 n2 p operations, and the inverse matrix of H in Equation (31) costs approximately
O(8(M + p − 1)3 ) flops. If M = mn and M  p, the above two steps take M 4 + 8M 3 p operations periteration. Then, our second method requires about N M 2 + O(M 4 + 8M 3 p) flops.
The main computational complexities of the proposed two algorithms are compared in Table 1 withthat 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 complexityof 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 thismultidimensional search process.
Table 1. Computational complexity of the ESPRIT-like and the proposed two techniques.
Algorithm

Computational Complexity

ESPRIT-like
Method 1
Method 2

N M 2 + O(2M 2 p)
N M 2 + q 2 p + O(M 2 p)
N M 2 + O(M 4 + 8M 3 p)

4. 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 tothe 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.10ej10 ,
◦
◦
◦
◦
◦
◦
◦
0.90e−j5 , 1.25ej20 , 0.80e−j9 , 0.96ej15 , 1.18e−j23 , 0.88e−j2 and 0.85ej4 . The ULA is divided into
5 subarrays in Method 2, and the i-th subarray is selected by matrix:

Sensors 2015, 1513910

J i = [05×(i−1) , I 5×5 , 05×(5−i) ], i = 1, · · · , 5(41)

However, the two subarrays of Method 1 are selected by matrices:h
iTJˆ1 = J T1 J T2 J T3 J T4

(42)

and:
Jˆ2 =

h

J T2 J T3 J T4 J T5

iT

(43)

Two equal-power signals located at 25◦ and 30◦ are considered, and the SNR per element for eachsource is defined by SNR = 10 log10 (σ2s /σ2n ), where σ2s and σ2n , respectively, denote the power of theincident signal and that of additive noise at each sensor. RMSE is used as the performance measure anddefined by:
K r
2
1 X
RMSE =
Vˆ i − V i
(44)
K i=1
F
where K is the number of trials and V̂ 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 inall simulations.
−6

10

Method2 initialized by ESPRIT−like
Method2 initialized by method1−7

Value (dB)

10

−8

10

−9

10

1

2

3

4
5
6
Iteration Number

7

8

9

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

Sensors 2015, 15

13911

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 theestimates 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 thefirst 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 willshow 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 andthe Method 1.
0.8
ESPRIT-like
Method2 initialized by method1
Method1
Method2 initialized by ESPRIT-like

0.7

0.6

RMSE

0.5

0.4

0.3

0.2

0.1

0
-10

-8

-6

-4

-2

0
SNR

2

4

6

8

10

Figure 2. RMSE of the spatial signature matrix V estimation versus SNR. The number ofsnapshots N = 500.
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 behaviorof the multidimensional search procedure (Method 2) has been improved when the initial parametervalues are provided by our Method 1. Simulation results also indicate that our method only requires

Sensors 2015, 15

13912

2 to 4 Gauss–Newton iterations, and the main performance improvement comes from the firsttwo 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 physicalstructure of the array can give better performance.
2
ESPRIT-like
Method2 W=Wopt

1.8

Method2 W=I

1.6

1.4

RMSE

1.2

1

0.8

0.6

0.4

0.2

0
-10

-8

-6

-4

-2

0
SNR

2

4

6

8

10

Figure 3. RMSE of the spatial signature matrix V estimation versus SNR with different W .
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 ofthe 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.

Sensors 2015, 15

13913
0.4
method1
method2 W=Wopt
0.35

0.3

RMSE

0.25

0.2

0.15

0.1

0.05

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7



Figure 4. RMSE of the spatial signature matrix V estimation versus correlation factor ρ.
The number of snapshots N = 1000 and SNR = 15 dB.
5. Conclusions
The ESPRIT-like method proposed in [22] just makes use of the maximum overlapping subarrayof 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. Anothermultidimensional 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 introductionof 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].
Acknowledgments
This work was supported in part by the National Natural Science Foundation of China under Grants
61172162 and 61231018, the Programme of Introducing Talents of Discipline to Universities underGrant B13043, and the Specialized Research Fund for the Doctoral Program of Higher Education ofChina under Grant 20130201110016.

Sensors 2015, 15

13914

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.
References
1. Paulraj, A.J.; Papadias, C.B. Space-time processing for wireless communications. IEEE Signal
Process. Mag. 1997, 14, 49–83.
2. Winters, J.H. Smart antennas for wireless systems. IEEE Pers. Commun. Mag. 1998, 5, 23–27.
3. Rong, Y.; Vorobyov, S.A.; Gershman, A.B.; Sidiropoulos, N.D. Blind spatial signature estimation
via time-varying user power loading and parallel facotr analysis. IEEE Trans. Signal Process.
2005, 53, 1697–1710.
4. Schmidt, R.O. Multiple emitter location and signal parameters estimation. IEEE Trans.
Antennas Propag. 1986, 34, 267–280.
5. Roy, R.; Kailath, T. ESPRIT-estimation of signal parameters via rotational invariance techniques.
IEEE Trans. Acoust. Speech Signal Process. 1989, 37, 984–995.
6. Stoica, P.; Nehorai, A. MUSIC, maximum likelihood, and cramer-rao bound. IEEE Trans. Acoust.
Speech Signal Process. 1989, 37, 720–741.
7. Ng, B.C.; See, C.M.S. Sensor-array calibration using a maximum-likelihood approach. IEEE Trans.
Antennas Propag. 1996, 44, 827–835.
8. Swindlehurst, A.; Kailath, T. A performance analysis of subspace-based methods in the presenceof model errors, part I: The MUSIC algorithm. IEEE Trans. Signal Process. 1992, 40, 1758–1774.
9. Swindlehurst, A.; Kailath, T. A performance analysis of subspace-based methods in the presenceof model errors, part II: Multidimensional algorithms. IEEE Trans. Signal Process. 1993, 41,
2882–2890.
10. Zhu, J.X.; Wang, H. Effects of sensor position and pattern perturbations on CRLB for
direction finding of multiple narrow-band sources. In Proceedings of the 4th Acoustics Speech
Signal Processing Workshop on Spectral Estimation and Modeling, Minneapolis, MN, USA,
3–5 August 1988; pp. 98–102.
11. Rockah, Y.; Schultheiss, P. Array shape calibration using sources in unknown locations–part I:
Far-field sources. IEEE Trans. Acoust. Speech Signal Process. 1987, 35, 286–299.
12. Rockah, Y.; Messer, H.; Schultheiss, P. Localization performance of arrays subject to phase errors.
IEEE Trans. Aerosp. Electron. Syst. 1988, 24, 402–410.
13. Weiss, A.J.; Friedlander, B. Eigenstructure methods for direction finding with sensor gain and phase
uncertainties. Circuits Syst. Signal Process. 1990, 9, 272–300.

Sensors 2015, 1513915

14. Viberg, M.; Swindlehurst, A. A bayesian approach to auto-calibration for parametric array signal
processing. IEEE Trans. Signal Process. 1994, 42, 3495–3507.
15. Soon, V.C.; Tong, L.; Huang, Y.F.; Liu, R. A subspace method for estimating sensor gains andphases. IEEE Trans. Signal Process. 1994, 42, 973–976.
16. Fuhrmann, D.R. Estimation of sensor gain and phase. IEEE Trans. Signal Process. 1994, 42,
77–87.
17. Wahlberg, B.; Ottersten, B.; Viberg, M. Robust signal parameter estimation in the presence of array
perturbations. In Proceedings of the 1991 IEEE International Conference on Acoustics, Speech,
and Signal Processing, Toronto, ON, Canada, 14–17 April 1991; pp. 3277–3280.
18. Jansson, M.; Swindlehurst, A.; Ottersten, B. Weighted subspace fitting for general array error
models. IEEE Trans. Signal Process. 1998, 46, 2484–2498.
19. Friedlander, B.; Weiss, A.J. Direction finding in the presence of mutual coupling. IEEE Trans.
Antennas Propag. 1991, 39, 273–284.
20. Pierre, J.; Kaveh, M. Experimental performance of calibration and direction-finding algorithms.
In Proceedings of the 1991 IEEE International Conference on Acoustics, Speech, and Signal
Processing, Toronto, ON, Canada, 14–17 April 1991; pp. 1365–1368.
21. Hung, E.H.L. A critical study of a self-calibrating direction-finding method for arrays. IEEE Trans.
Signal Process. 1994, 42, 471–474.
22. Astély, D.; Swindlehurst, A.; Ottersten, B. Spatial signature estimation for uniform linear arrayswith unknown receive gains and phases. IEEE Trans. Signal Process. 1999, 47, 2128–2138.
23. Paulraj, A.; Kailath, T. Direction of arrival estimation by eigenstructure methods with unknown
sensor gain and phase. In Proceedings of the 1985 IEEE International Conference on Acoustics,
Speech, and Signal Processing, Tampa, FL, USA, 26–29 April 1985; pp. 640–643.
24. Li, Y.; Er, M.H. Theoretical analyses of gain and phase error calibration with optimal
implementation for linear equispaced array. IEEE Trans. Signal Process. 2006, 54, 712–723.
25. Gunther, J.; Swindlehurst, A. Algorithms for blind equalization with multiple antennas basedon frequency domain subspaces. In Proceedings of the 1996 IEEE International Conference on
Acoustics, Speech, and Signal Processing, Atlanta, GA, USA, 7–10 May 1996; pp. 2419–2422.
26. Liao, B.; Chan, S.C. Direction finding with partly calibrated uniform linear arrays. IEEE Trans.
Antennas Propag. 2012, 60, 922–929.
27. Ottersten, B.; Viberg, M.; Kailath, T. Performance analysis of the total least squares ESPRIT
algorithm. IEEE Trans. Signal Process. 1991, 39, 1122–1135.
28. Swindlehurst, A.; Ottersten, B.; Roy, R.; Kailath, T. Multiple invariance ESPRIT. IEEE Trans.
Signal Process. 1992, 40, 867–881.
29. Golub, G.H.; Pereyra, V. The differentiation of pseudo-inverses and nonlinear least squares
problems whose variables separate. SIAM J. Numer. Anal. 1973, 10, 413–432.
30. Viberg, M.; Ottersten, B.; Kailath, T. Detection and estimation in sensor arrays using weighted
subspace fitting. IEEE Trans. Signal Process. 1991, 39, 2436–2449.
31. Viberg, M.; Ottersten, B. Sensor array processing based on subspace fitting. IEEE Trans.
Signal Process. 1991, 39, 1110–1121.

Sensors 2015, 1513916

32. Xin, J.; Sano, A. Computationally Efficient Subspace-Based Method for Direction-of-Arrival
Estimation without Eigendecomposition. IEEE Trans. Signal Process. 2004, 52, 876–893.
c 2015 by the authors; licensee MDPI, Basel, Switzerland. This article is an open access article
distributed under the terms and conditions of the Creative Commons Attribution license(http://creativecommons.org/licenses/by/4.0/).