A Method for Spatial Signature Estimation with
Uncalibrated Uniform Linear Array
Xiang Cao #∗1 , Jingmin Xin # , Yoshifumi Nishio ∗2 , Nanning Zheng # , Akira Sano †3#

Institute of Artiﬁcial Intelligence and Robotics, Xi’an Jiaotong University
Xi’an 710049, China
1∗cao.xiang@stu.xjtu.edu.cn

Department of Electrical and Electronic Engineering, The University of Tokushima
Tokushima 770-8506, Japan2†nishio@ee.tokushima-u.ac.jp

Department of System Design Engineering, Keio University
Yokohama 223-8522, Japan3

sano@sd.keio.ac.jp

Abstract—In this paper, a multidimensional search procedure
is presented for estimating spatial signature with uniform linear
array in the presence of unknown gain and phase errors.
The proposed method fully exploits multiple invariances in theinterior of the sensor array and it is desirable that the quality of
estimate improves as more and more array invariances are used.
A Gauss-Newton iterative algorithm is applied for solving themultidimensional nonlinear problem. Simulation results indicatethat a solution close to the minimum is obtained after only twoor three iterations. Furthermore, the introduction of a weighting
matrix, similar to weighted subspace ﬁtting, can tone up the
performance of our method.

I. I NTRODUCTION
In sensor array signal processing, many high-resolution
methods such as ESPRIT [1] have been advanced for estimation of unknown parameters embedded in array output
model. It is a natural motivation of array calibration, since theperformance of these methods will degrade dramatically when
the ideal array model is damaged by some unknown errors(e.g., mutual coupling and sensor position uncertainties).
Two parametric errors (i.e., statistical and deterministic)
were usually used to extend the ideal array model. This paper
only focus on the deterministic unknown gain and phase errors.
The term deterministic here implies that the unknown error ateach of elements of array is a complex stable constant during
the period of observation.
Several so-called auto-calibration approaches have beendeveloped in literature. Here, the auto-calibration indicatesthat array calibration may be accomplished without employing
any dummy elements or transmitters at known direction. In[2], an iterative eigenstructure-based technique of estimating
Directions-of-Arrival (DOA’s) in the presence of unknown
gain and phase errors is presented and it can apply to arbitraryarray geometries except uniform linear array (ULA). For ULA,a phase ambiguity exists between the diagonal error matrix
and the ideal array steering matrix (see [3] and [4]). Thus,
it is impossible to estimate DOA and gain and phase errors
simultaneously for ULA. A Hermitian Toeplitz structure of

ULA output covariance matrix in the absence of gain andphase errors is exploited in [5] and [6]. This method takes
advantage of the elements equivalence at every diagonal line
to form two equations for estimation of gain and phase errors.
The main drawback of the approach is that the Toeplitz
matrix assumption is only established under inﬁnite sampling
condition. It means that if the number of snapshots is ﬁxed, theperformance of the algorithm does not improve after increasing
the signal-to-noise ratio (SNR) to a certain point. In addition,
ESPRIT method also can be extended to the case. The idea
ﬁrstly emerged from [7], and is applied for estimating spatial
signature in [8] and [4]. A closed-form ESPRIT-like method(see [4]) constructs an estimator based on the maximum
overlapping two subarrays. However, it has been pointed out
in [9] that any two subarrays conﬁguration may be inherently
suboptimal.
In this paper, we consider the spatial signature estimation
problem with ULA in the presence of unknown gain and phaseerrors. The practicality of spatial signature can be found in[8] and [4]. Unlike ESPRIT-like method, our algorithm fully
exploits the multiple invariances into ULA even though thearray model errors exist. The sensor gain and phase errors
are different from mutual coupling among elements of array
and it can not be affected by other sensors. Furthermore, thesteering matrix of ideal ULA has a Vandermonde structure.
Both of them provide the possibility of selecting subarraysfrom uncalibrated array. It is natural to expect that an estimator
formulated by multiple invariances in the array performs better
than the one involved only one invariance. In practice, amultidimensional procedure of minimizing a cost function ispresented. Our method ﬁrstly estimates the rotational DOA’s(not the absolute DOA’s) and error parameters, then obtains theestimation of spatial signature. The Gauss-Newton algorithm[10] is suggested for solving the multidimensional search
problem. The choice of the initial values is essential for the
global convergence of the iterative method. The excellent
initial inputs are provided by ESPRIT-like method. Computersimulations show that only two or three iterations are required
to achieve acceptable estimation values.

The determinant of the matrix Γ(γ) is given bydet(Γ(γ)) =II. DATA MODEL

y(t) = Γ(γ)A(θ)s(t) + n(t),(1)

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 idealarraysteering matrix A(θ) = a(θ1 ) a(θ2 ) · · · a(θp ) , anda(θ) =

2π

1 ej λ d sin(θ)···


j 2π (M −1)d sin(θ) Te λγi = 0,(3)

i=1

A uniform linear array with M sensors received narrowband
signals from p far-ﬁeld sources and the vector response y ∈C M ×1 of the array at time t can be expressed as



M


.

(2)

Here, θ1 , θ2 , · · · , θp are the Directions-of-Arrival of signals, dand λ represent the distance between two consecutive sensors
and the identical wavelength
for all signals, respectively. The
matrix Γ(γ) = diag γ1 γ2 · · · γM
and |γi | > 0
denotes the deterministic unknown gain and phase error of
sensor i.
Two assumptions need to be made. Firstly, s(t) is a temporally complex white Gaussian random vector with mean zeroand its covariance matrix Rss has full rank p (assuming nocorrelate signals). Secondly, n(t) is a temporally and spatiallycomplex white Gaussian random vector with mean zero anduncorrelated with incident signals.
The focus of this paper is the estimation of Spatial Signature matrix V = Γ(γ)A(θ) from the N snapshots of
the array outputs. One ambiguity for this problem may beobserved 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 ﬁrst element of diagonal matrix Γ(γ) be equal toone.

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 relation implies that there exists
a nonsingular p × p matrix T satisfying
J i E s = J i Γ(γ)A(θ)T .
Considering n subarrays, we have⎡⎤ ⎡
Γ 1 A1E s1⎢ E s2 ⎥ ⎢ Γ2 A1 Φε1
⎢
⎥ ⎢
ε2
⎢
⎥ ⎢
J E s = ⎢ E s3 ⎥ = ⎢ Γ3 A1 Φ⎢ .. ⎥ ⎢ ..
⎣ . ⎦ ⎣ .
E sn

⎤
⎥
⎥
⎥
⎥T,
⎥
⎦

(5)

Γn A1 Φεn−1

where the matrix J is formed by selection matrix J i and J = T
T, E si , Γi and A1 may be calculated
J 1 J T2 · · · J Tn
by the following equations:
E si = J i E s , Γi = J i Γ(γ)J Ti , A1 = J 1 A(θ).
The matrix Φ is given by 2π
Φ = diag ej λ d sin(θ1 )2πej λ d sin(θp )

···



.

(6)

(7)

εi d denotes the distance between the (i + 1)th subarray andthe reference subarray. Analogous to MI ESPRIT method in[11], the unknown parameter vector μ may now be obtained
from the following least squares problem
⎡
⎡
⎤
⎤
2
E s1
Γ 1 A1⎢ E s2 ⎥⎢ Γ 2 A 1 Φε 1
⎥
⎢
⎢
⎥
⎥
⎢ E s3 ⎥ 1/2 ⎢ Γ3 A1 Φε2⎥Ω(μ) = ⎢
−⎢
(8)
⎥W
⎥T
⎢ .. ⎥
⎢ ..
⎥
⎣ . ⎦
⎣ .
⎦
Γn A1 Φεn−1

E snIII. ESTIMATION ALGORITHM
A new subspace approach is presented for estimating spatial
signature matrix V in this section. The so-called signal
subspace E s and noise subspace E n can be obtained fromthe array outputs covariance matrix R = E{y(t)y H (t)}.
The Vandermonde structure of the ideal array steeringmatrix A(θ) permits optimally exploiting multiple invariancesfor ULA even though the gain and phase of the sensor has notbeen calibrated. It is expected that an algorithm, fully making
use of the particular properties, should obtain superior results
than ESPRIT-like estimator proposed in [4]. Assume that thearray comprises n identical subarrays of m sensors. Obviously,
overlapping subarrays make M ≤ mn. In addition, since the
one-to-one correspondence between rows in V and elementsof array, extracting a subarray from array is equivalent topicking up m rows of matrix V by a m × M selection
matrix J i . The full row rank matrix J i consists of zeros andones. Only one element at every row is equal to one and it
corresponds to the selected sensor of the array.

(4)

F

and
μ

=

[ γ̄2

···

γ̄M

ρ1

···

ρp

γ̃2
θ1

···
···

γ̃M
θp ] T ,(9)where W is a Hermitian positive deﬁnite weighting matrix,
γ̄i and γ̃i , respectively, denote the real part and image part ofscaling γi . Since the signal subspace E s is usually substituted
with Ê s obtained from the eigendecomposition of the sample
covariance matrix, the unitary matrix Φ in (7) becomes


2π
2π
Φ = diag ρ1 ej λ d sin(θ1 ) · · · ρp ej λ d sin(θp ) . (10)Next, for ease of notation, we deﬁneΥ = W T /2 E Ts1 E Ts2 · · · E Tsn
and
Ψ=



AT1 Γ1Φε1 AT1 Γ2···

Φεn−1 AT1 Γn

(11)


.

(12)

The expression in (8) is reformulated as2

Ω(μ) = Υ − T T Ψ

F

.

(13)

Here, the equations Φεi T = Φεi and Γi T = Γi are used.
This is a nonlinear least squares problem for unknown vector
μ (see, [10]). Solving T T by minimizing the cost function Ω,
we can obtainT̂T= ΥΨH ΨΨH

−1

.

(14)

Substituting this back to (13) yields to2μ̂ = arg min Ω = arg min ΥΠ⊥
Ψ F
μ
μ

(15)

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

(16)

where

Note that, it is reasonable to involve all of the elements of thearray for estimating unknown parameters. The limiting caseis that there are no overlapping subarrays (i.e., M = mn).
Thus, the matrix Ψ has full row rank except for deleting some
sensors from the array (i.e., M < mn). So the matrix (ΨΨH )is nonsingular.
IV. ALGORITHM IMPLEMENTATION
The method considered herein requires a multidimensional
search over parameter vector μ ∈ R2(M +p−1)×1 . It is wellknown that the Gauss-Newton algorithm [10], [12] is a classical way to obtain the solution of this problem. Considering
the cost function Ω in (15), the unknown parameters can be
iteratively obtained byμk+1= μ − ξk H
k

−1

Q,

(17)

where μk denotes the estimate at k iteration, ξk is a iterative
step length, the gradient vector Q = ∂∂Ωμ and the Hessian2kmatrix H = ∂ μ∂∂Ω
μT are evaluated at μ .
We deﬁne a vector r formed by columns scanning matrixΥΠ⊥Ψ , i.e.,
r = vec ΥΠ⊥
(18)
Ψ .

Then, the cost function Ω becomesΩ = r H r.

(19)

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 ith element of the gradient vector Q is given by


∂Ω
Qi =
= 2Re r H
(20)
i r ,
∂μi
where r i = ∂∂μr . The ijth component of the Hessian matriximay be expressed as (see [12], [13])H ij =


∂2Ω
 2Re r Hi rj .
∂μi ∂μj(21)

Reasons for this approximation (i.e., ignoring the second
derivative of r ) can be found in [11], [12] and it brings two
beneﬁts for the Gauss-Newton iterative algorithm. One is thatthe approximate Hessian matrix is positive semideﬁnite andit guarantees −H −1 Q to be decent direction. The other isthat we only need to calculate the ﬁrst derivative of r insteadof the second one. Furthermore, a useful modiﬁcation to illconditioned cases is to usually use (H + ζI) in lieu of H,
where ζ is a scaling.
Next, we focus on the ﬁrst derivative of vector r in (18).
By the aid of the differentiation of the projection matrix [10],r i may be expressed as

 
H
†
⊥
†
⊥
r i = −vec Υ Ψ Ψi ΠΨ + Ψ Ψi ΠΨ
(22)
Ψ . The expression
where Ψ† = ΨH (ΨΨH )−1 and Ψi = ∂∂μi
of Ψi can be obtained from (6), (9), (10) and (12). This isa tiresome calculation process. However, it can be simpliﬁed
εT∂ Φ i−1 A1 Γiwith the observation that
= 0 if γj can not be in∂γjεT∂ Φ i−1 A1 Γi= Φεi−1 AT1 I kk
cluded in matrix Γi , otherwise,
∂ γ̄j
εi−1
T
∂Φ
A1 Γ iand= jΦεi−1 AT1 I kk , where I kk is a matrix∂ γ̃jwith one in the kkth position and zeros elsewhere, γj is thekth element at the main diagonal of matrix Γj . Moreover,T∂Φεi−1 AT1 Γi∂Φεi−1 T
εi−1 ∂A1
Γi , (23)=
A1 + Φ
∂θj
∂θj
∂θj
and
∂Φεi−1 AT1 Γi∂Φεi−1 T=
A 1 Γi .
∂ρj
∂ρj

(24)

From the above analyses, the vector r i can be calculated
conveniently.
The proposed algorithm now is brieﬂy outlined below.
1) Estimate the signal subspace Ê s from the eigendecompositionof the sample covariance matrix R̂ =
N
1
H
y(n)y
(n), where N is the ﬁnite number ofn=1
N
snapshots.
2) Initialize vector μ utilizing ESPRIT-like method [4].
3) Update μ according to (17) until the stopping criterion
is satisﬁed, then form the estimates of γ̂ and θ̂.
4) Calculate matrix Γ(γ̂) and A(θ̂), then spatial signaturematrix is estimated as V̂ = Γ(γ̂)A(θ̂).
Note that, although the ESPRIT-like method can only provide rotational values over true DOA’s θ and true gain andphase error matrix Γ, it makes no difference to estimationof spatial signature V . The identiﬁability of V from signal
subspace E s can be found in [4]. In addition, it is an
intractable problem on how to choose weighting matrix W .
Refer to the weighted subspace ﬁtting (WSF) method, we
deﬁne
−1

W opt = (Λ̂s − σ̂n2 I)2 Λ̂s ,

(25)

1

10

−2.1

ESPRIT−like
Proposed method Wopt
Proposed method W=I

ESPRIT−like
Proposed method Wopt
Proposed method W=I

10

−2.2

10

−2.3

10

RMSE (dB)

RMSE (dB)

−2.4

0

10

10

−2.5

10

−2.6

10

−2.7

10

−2.8

10

−1

10

0
−10

−8

−6

−4

−2

0
SNR

2

4

6

8

10

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

where the diagonal matrix Λ̂s contains the so-called signal
eigenvalues estimated from the array outputs sample covariance matrix, σ̂n2 denotes the consistent estimate of noise power.
Simulation results using the proposed algorithm show that
W = W opt results in a better performance for parameter
estimates than W = I.
V. EXPERIMENTS
In this section, computer simulations were conducted forevaluating the performance of the proposed algorithm. In all
scenarios, a ULA of 9 elements with half of wavelength
element spacing is used. The ULA is divided into 5 subarrays and the ith subarray is selected by matrix J i =[05×(i−1) , I 5×5 , 05×(5−i) ], i = 1, · · · , 5. The deterministic◦◦
gain and phase errors are given by 1, 1.10ej10 , 0.90e−j5 ,
j20◦
−j9◦
j15◦
−j23◦
−j2◦
1.25e
, 0.80e, 0.96e, 1.18e, 0.88e
and◦
0.85ej4 . Two equal-power uncorrelated signals are considered
and the SNR per element for each source is deﬁned by
SNR = 10 log10 (σs2 /σn2 ), where σs2 and σn2 , respectively,
denote the power of incident signal and that of additive noiseat each sensor. RMSE is used as the performance measure anddeﬁned byK 
2
1 
RMSE =
Vˆ i − V i ,
(26)
K i=1
F
where K is the number of trials, V̂ i and V i are the estimatedspatial signature and the true one at the ith experiment,
respectively. A total of 200 trials were performed for eachsimulation scenario.
Example 1: Performance versus SNR
In this example, two closed sources located at 25◦ and 30◦
were used. For comparison, we give the results of ESPRIT-like
method [4] and the proposed method in this paper. Figure 1
shows that our method performs better than ESPRIT-like evenin the weighting matrix W = I. It means that the maximum
overlapping subarray is not the optimal choice in the presence

500

1000

1500

2000

2500
3000
Snapshots N

3500

4000

4500

5000

Fig. 2. RMSE of the spatial signature matrix H estimation versus the numberof snapshots N . SNR = 15dBof unknown sensor gain and phase responses. In other words,
fully exploiting the multiple invariances in the array should
be taken into account. In addition, adding a weighting matrix
can evidently improve the performance of our algorithm.
Example 2: Performance versus Number of Snapshots
In the second case, the performance of different algorithms
versus the number of snapshots N were compared. The DOA’s
assumed to be 25◦ and 40◦ . Similar conclusions to the Example 1 can be made from Figure 2. It is worth pointing out thatthe behavior of our method has been improved signiﬁcantly atlarge number of snapshots. The reason may be that excellent
initial input values are provided by the ESPRIT-like method.
Simulation results also show that our method only requires
2-4 Gauss-Newton iterations and the main performance improvement comes from the ﬁrst two iterations.
VI. CONCLUSIONSA multidimensional search algorithm is presented for spatial
signature estimation for ULA with unknown gain and phaseerrors. Our method fully excavates the multiple invariancesof the ULA and is implemented with the help of the GaussNewton iterative algorithm. Because of excellent initial parameter values provided by ESPRIT-like, our method canconverge rapidly to a appropriate solution. On the other hand,
the introduction of a weighting matrix W further improve
the performance of the proposed algorithm. How to selectsubarrays from ULA, making the algorithm lead to optimal
parameter estimation, is not addressed in this paper. Somepragmatic discussion about this question can be found in [9].
ACKNOWLEDGMENT
This work was supported in part by the National Nature
Science Foundation of China under Grant 61172162.
R EFERENCES[1] R. Roy and T. Kailath, “Esprit-estimation of signal parameters via
rotational invariance techniques,” IEEE Trans. Acoust.,Speech, Signal
Process., vol. 37, no. 7, pp. 984–995, 1989.

[2] A. J. Weiss and B. Friedlander, “Eigenstructure methods for direction
ﬁnding with sensor gain and phase uncertainties,” Circuits Syst. Signal
Process., vol. 9, no. 3, pp. 272–300, 1990.
[3] J. Pierre and M. Kaveh, “Experimental performance of calibration anddirection-ﬁnding algorithms,” in Proc. of IEEE ICASSP’91, Toronto,
Ont., Apr. 1991, pp. 1365–1368.
[4] D. Astély, A. Swindlehurst, and B. Ottersten, “Spatial signature estimation for uniform linear arrays with unknown receive gains and phases,”IEEE Trans. Signal Prosess., vol. 47, pp. 2128–2138, Aug. 1999.
[5] A. Paulraj and T. Kailath, “Direction of arrival estimation by eigenstructure methods with unknown sensor gain and phase,” in Proc. of IEEEICASSP’85, Tampa, FL, Apr. 1985, pp. 640–643.
[6] Y. Li and M. H. Er, “Theoretical analyses of gain and phase error
calibration with optimal implementation for linear equispaced array,”IEEE Trans. Signal Process., vol. 54, pp. 712–723, Feb. 2006.
[7] V. C. Soon, L. Tong, Y. F. Huang, and R. Liu, “A subspace method
for estimating sensor gains and phases,” IEEE Trans. Signal Process.,
vol. 42, pp. 973–976, Apr. 1994.
[8] J. Gunther and A. Swindlehurst, “Algorithms for blind equalization withmultiple antennas based on frequency domain subspaces,” in Proc. ofIEEE ICASSP’96, Atlanta, GA, May 1996, pp. 2419–2422.
[9] B. Ottersten, M. Viberg, and T. Kailath, “Performance analysis of
the total least squares esprit algorithm,” IEEE Trans. Signal Process.,
vol. 39, pp. 1122–1135, May 1991.
[10] G. H. Golub and V. Pereyra, “The differentiation of pseudo-inverses
and nonlinear least squares problems whose variables separate,” SIAM
J. Numer. Anal., vol. 10, no. 2, pp. 413–432, 1973.
[11] A. Swindlehurst, B. Ottersten, R. Roy, and T. Kailath, “Multiple invariance esprit,” IEEE Trans. Signal Prosess., vol. 40, pp. 867–881, Apr.
1992.
[12] M. Viberg, B. Ottersten, and T. Kailath, “Detection and estimationin sensor arrays using weighted subspace ﬁtting,” IEEE Trans. Signal
Process., vol. 39, pp. 2436–2449, Nov. 1991.
[13] M. Viberg and B. Ottersten, “Sensor array processing based on subspaceﬁtting,” IEEE Trans. Signal Process., vol. 39, pp. 1110–1121, May 1991.