Coarrays, MUSIC, and the Cramér

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 65, NO. 4, FEBRUARY 15, 2017 933

Coarrays, MUSIC, and the Cram

er–Rao Bound

Mianzhi Wang, Student Member, IEEE, and Arye Nehorai, Life Fellow, IEEE

Abstract—Sparse linear arrays, such as coprime arrays and

nested arrays, have the attractive capability of providing enhanced

degrees of freedom. By exploiting the coarray structure, an

augmented sample covariance matrix can be constructed and

MUtiple SIgnal Classiﬁcation (MUSIC) can be applied to identify

more sources than the number of sensors. While such a MUSIC

algorithm works quite well, its performance has not been theoret-

ically analyzed. In this paper, we derive a simpliﬁed asymptotic

mean square error (MSE) expression for the MUSIC algorithm

applied to the coarray model, which is applicable even if the source

number exceeds the sensor number. We show that the directly

augmented sample covariance matrix and the spatial smoothed

sample covariance matrix yield the same asymptotic MSE for

MUSIC. We also show that when there are more sources than the

number of sensors, the MSE converges to a positive value instead

of zero when the signal-to-noise ratio (SNR) goes to inﬁnity. This

ﬁnding explains the “saturation” behavior of the coarray-based

MUSIC algorithms in the high-SNR region observed in previous

studies. Finally, we derive the Cram

er–Rao bound for sparse

linear arrays, and conduct a numerical study of the statistical

efﬁciency of the coarray-based estimator. Experimental results

verify theoretical derivations and reveal the complex efﬁciency

pattern of coarray-based MUSIC algorithms.

Index Terms—MUSIC, Cram

er–Rao bound, coarray, sparse

linear arrays, statistical efﬁciency.

I. INTRODUCTION

STIMATING source directions of arrivals (DOAs) using

sensors arrays plays an important role in the ﬁeld of ar-

ray signal processing. For uniform linear arrays (ULA), it is

widely known that traditional subspace-based methods, such as

MUSIC, can resolve up to N −1 uncorrelated sources with N

sensors [1]–[3]. However, for sparse linear arrays, such as mini-

mal redundancy arrays (MRA) [4], it is possible to construct an

augmented covariance matrix by exploiting the coarray struc-

ture. We can then apply MUSIC to the augmented covariance

matrix, and up to O(N

) sources can be resolved with only N

sensors [4].

Recently, the development of co-prime arrays [5]–[8] and

nested arrays [9]–[11], has generated renewed interest in sparse

linear arrays, and it remains to investigate the performance

of these arrays. The performance of the MUSIC estimator

and its variants (e.g., root-MUSIC [12], [13]) was thoroughly

analyzed by Stoica et al. in [2], [14] and [15]. The same authors

also derived the asymptotic MSE expression of the MUSIC

Manuscript received March 31, 2016; revised September 16, 2016; accepted

October 19, 2016. Date of publication November 8, 2016; date of current ver-

sion December 5, 2016. The associate editor coordinating the review of this

manuscript and approving it for publication was Prof. Wee Peng Tay. This work

was supported by the AFOSR under Grant FA9550-11-1-0210 and the ONR

under Grant N000141310050.

The authors are with the Preston M. Green Department of Electrical and

Systems Engineering, Washington University in St. Louis, St. Louis, MO 63130

USA (e-mail: [email protected]; [email protected]).

Color versions of one or more of the ﬁgures in this paper are available online

at http://ieeexplore.ieee.org.

Digital Object Identiﬁer 10.1109/TSP.2016.2626255

estimator, and rigorously studied its statistical efﬁciency.

In [16], Li et al. derived a uniﬁed MSE expression for common

subspace-based estimators (e.g., MUSIC and ESPRIT [17])

via ﬁrst-order perturbation analysis. However, these results

are based on the physical array model and make use of the

statistical properties of the original sample covariance matrix,

which cannot be applied when the coarray model is utilized.

In [18], Gorokhov et al. ﬁrst derived a general MSE expression

for the MUSIC algorithm applied to matrix-valued transforms

of the sample covariance matrix. While this expression is

applicable to coarray-based MUSIC, its explicit form is

rather complicated, making it difﬁcult to conduct analytical

performance studies. Therefore, a simpler and more revealing

MSE expression is desired.

In this paper, we ﬁrst review the coarray signal model

commonly used for sparse linear arrays. We investigate two

common approaches to constructing the augmented sample

covariance matrix, namely, the direct augmentation approach

(DAA) [19], [20] and the spatial smoothing approach [9]. We

show that MUSIC yields the same asymptotic estimation er-

ror for both approaches. We are then able to derive an explicit

MSE expression that is applicable to both approaches. Our MSE

expression has a simpler form, which may facilitate the perfor-

mance analysis of coarray-based MUSIC algorithms. We ob-

serve that the MSE of coarray-based MUSIC depends on both

the physical array geometry and the coarray geometry. We show

that, when there are more sources than the number of sensors,

the MSE does not drop to zero even if the SNR approaches

inﬁnity, which agrees with the experimental results in previous

studies. Next, we derive the CRB of DOAs that is applicable to

sparse linear arrays. We notice that when there are more sources

than the number of sensors, the CRB is strictly nonzero as the

SNR goes to inﬁnity, which is consistent with our observation

on the MSE expression. It should be mentioned that during the

review process of this paper, Liu et al. and Koochakzadeh et al.

also independently derived the CRB for sparse linear arrays in

[21], [22]. In this paper, we provide a more rigorous proof the

CRB’s limiting properties in high SNR regions. We also include

various statistical efﬁciency analysis by utilizing our results on

MSE, which is not present in [21], [22]. Finally, we verify our

analytical MSE expression and analyze the statistical efﬁciency

of different sparse linear arrays via numerical simulations. We

we observe good agreement between the empirical MSE and

the analytical MSE, as well as complex efﬁciency patterns of

coarray-based MUSIC.

Throughout this paper, we make use of the following nota-

tions. Given a matrix A,weuseA

, A

, and A

∗

to denote

the transpose, the Hermitian transpose, and the conjugate of

A, respectively. We use A

to denote the (i, j)-th element of A,

and a

to denote the i-th column of A.IfA is full column rank,

we deﬁne its pseudo inverse as A

†

=(A

−1

.Wealso

deﬁne the projection matrix onto the null space of A as Π

⊥

I − AA

†

.LetA =[a

... a

] ∈ C

M × N

, and we deﬁne

See http://www.ieee.org/publications

standards/publications/rights/index.html for more information.

934 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 65, NO. 4, FEBRUARY 15, 2017

the vectorization operation as vec(A)=[a

... a

]

, and

mat

M,N

(·) as its inverse operation. We use ⊗ and  to denote

the Kronecker product and the Khatri-Rao product (i.e., the

column-wise Kronecker product), respectively. We denote by

R(A) and I(A) the real and the imaginary parts of A.IfA is

a square matrix, we denote its trace by tr(A). In addition, we

use T

to denote a M × M permutation matrix whose anti-

diagonal elements are one, and whose remaining elements are

zero. We say a complex vector z ∈ C

is conjugate symmetric

if T

z = z

∗

.Wealsousee

to denote the i-th natural base vec-

tor in Euclidean space. For instance, Ae

yields the i-th column

of A, and e

A yields the i-th row of A.

II. C

OARRAY SIGNAL MODEL

We consider a linear sparse array consisting of M sensors

whose locations are given by D = {d

,...,d

}. Each sen-

sor location d

is chosen to be the integer multiple of the smallest

distance between any two sensors, denoted by d

. Therefore we

can also represent the sensor locations using the integer set

D = {

,...,

}, where

= d

for i =1, 2,...,M.

Without loss of generality, we assume that the ﬁrst sensor

is placed at the origin. We consider K narrow-band sources

,θ

,...,θ

impinging on the array from the far ﬁeld. Denot-

ing λ as the wavelength of the carrier frequency, we can express

the steering vector for the k-th source as

a(θ



1 e

··· e



, (1)

where φ

=(2πd

sin θ

)/λ. Hence the received signal vectors

are given by

y(t)=A(θ)x(t)+n(t),t =1, 2,...,N, (2)

where A =[a(θ

) a(θ

) ... a(θ

)] denotes the array steering

matrix, x(t) denotes the source signal vector, n(t) denotes ad-

ditive noise, and N denotes the number of snapshots. In the

following discussion, we make the following assumptions:

A1 The source signals follow the unconditional model

[15] and are uncorrelated white circularly-symmetric

Gaussian.

A2 The source DOAs are distinct (i.e., θ

= θ

∀k = l).

A3 The additive noise is white circularly-symmetric Gaus-

sian and uncorrelated from the sources.

A4 The is no temporal correlation between each snapshot.

Under A1–A4, the sample covariance matrix is given by

R = AP A

+ σ

I, (3)

where P = diag(p

,...,p

) denotes the source covariance

matrix, and σ

denotes the variance of the additive noise. By

vectorizing R, we can obtain the following coarray model:

r = A

p + σ

i, (4)

where A

= A

∗

 A, p =[p

,...,p

]

, and i = vec(I).

It has been shown in [9] that A

corresponds to the steer-

ing matrix of the coarray whose sensor locations are given by

= {d

− d

|1 ≤ m, n ≤ M}. By carefully selecting rows

of (A

∗

 A), we can construct a new steering matrix repre-

senting a virtual ULA with enhanced degrees of freedom. Be-

cause D

is symmetric, this virtual ULA is centered at the

origin. The sensor locations of the virtual ULA are given by

Fig. 1. A coprime array with sensors located at [0, 2, 3, 4, 6, 9]λ/2 and its

coarray: (a) physical array; (b) coarray; (c) central ULA part of the coarray.

[−M

+1, −M

+2,...,0,...,M

− 1]d

, where M

is de-

ﬁned such that 2M

− 1 is the size of the virtual ULA. Fig. 1

provides an illustrative example of the relationship between the

physical array and the corresponding virtual ULA. The obser-

vation vector of the virtual ULA is given by

z = Fr = A

p + σ

Fi, (5)

where F is the coarray selection matrix, whose detailed deﬁni-

tion is provided in Appendix A, and A

represents the steering

matrix of the virtual array. The virtual ULA can be divided

into M

overlapping uniform subarrays of size M

. The output

of the i-th subarray is given by z

= Γ

z for i =1, 2,...,M

where Γ

=[0

× (i−1)

× M

× (M

−i)

] represents the

selection matrix for the i-th subarray.

Given the outputs of the M

subarrays, the augmented covari-

ance matrix of the virtual array R

is commonly constructed

via one of the following methods [9], [20]:

=[z

−1

··· z

], (6a)



i=1

, (6b)

where method (6a) corresponds to DAA, while method (6b)

corresponds to the spatial smoothing approach.

Following the results in [9] and [20], R

and R

are related

via the following equality:

+ σ

, (7)

where A

corresponds to the steering matrix of a ULA whose

sensors are located at [0, 1,...,M

− 1]d

. If we design a sparse

linear array such that M

>M, we immediately gain enhanced

degrees of freedom by applying MUSIC to either R

or R

instead of R in (3). For example, in Fig. 1, we have a co-prime

array with M

=8> 6=M. Because MUSIC is applicable

only when the number of sources is less than the number of

sensors, we assume that K<M

throughout the paper. This

assumption, combined with A2, ensures that A

is full column

rank.

It should be noted that the elements in (6a) are obtained

via linear operations on the elements in R, and those in (6b)

are obtained via quadratic operations. Therefore the statistical

properties of R

and R

are different from that of R.

Consequently, traditional performance analysis for the MUSIC

algorithm based on R cannot be applied to the coarray-based

MUSIC. For brevity, we use the term direct augmentation

based MUSIC (DA-MUSIC), and the term spatial smoothing

based MUSIC (SS-MUSIC) to denote the MUSIC algorithm

WANG AND NEHORAI: COARRAYS, MUSIC, AND THE CRAM

ER–RAO BOUND 935

applied to R

and R

, respectively. In the following section,

we will derive a uniﬁed analytical MSE expression for both

DA-MUSIC and SS-MUSIC.

III. MSE

OF COARRAY-BASED MUSIC

In practice, the real sample covariance matrix R is

unobtainable, and its maximum-likelihood estimate

R =

1/N



t=1

x(t)x(t)

is used. Therefore z, R

, and R

are

also replaced with their estimated versions

, and

Due to the estimation error ΔR =

R − R, the estimated noise

eigenvectors will deviate from the true one, leading to DOA

estimation errors.

In general, the eigenvectors of a perturbed matrix are not

well-determined [23]. For instance, in the very low SNR sce-

nario, ΔR may cause a subspace swap, and the estimated noise

eigenvectors will deviate drastically from the true ones [24].

Nevertheless, as shownin [16], [18] and [25], given enough sam-

ples and sufﬁcient SNR, it is possible to obtain the closed-form

expressions for DOA estimation errors via ﬁrst-order analysis.

Following similar ideas, we are able to derive the closed-form

error expression for DA-MUSIC and SS-MUSIC, as stated in

Theorem 1.

Theorem 1: Let

(1)

and

(2)

denote the estimated values

of the k-th DOA by DA-MUSIC and SS-MUSIC, respectively.

Let Δr = vec(

R − R). Assume the signal subspace and the

noise subspace are well-separated, so that Δr does not cause a

subspace swap. Then

(1)

− θ

(2)

− θ

= −(γ

)

−1

R(ξ

Δr), (8)

where

= denotes asymptotic equality, and

= F

(β

⊗ α

), (9a)

= −e

†

, (9b)

= Π

⊥

(θ

), (9c)

(θ

)Π

⊥

(θ

), (9d)

Γ =



−1

···Γ



, (9e)

(θ

∂a

(θ

)

∂θ

. (9f)

Proof: See Appendix B. 

Theorem 1 can be reinforced by Proposition 1. β

=0en-

sures that γ

−1

exists and (8) is well-deﬁned, while ξ

=0

ensures that (8) depends on Δr and cannot be trivially zero.

Proposition 1: β

, ξ

= 0 for k =1, 2,...,K.

Proof: We ﬁrst show that β

= 0 by contradiction.

Assume β

= 0. Then Π

⊥

(θ

)=0, where D =

diag(0, 1,...,M

− 1). This implies that Da

(θ

) lies in the

column space of A

.Leth = e

−jφ

(θ

). We immedi-

ately obtain that [A

h] is not full column rank. We now add

− K − 1 distinct DOAs in (−π/2,π/2) that are different

from θ

,...,θ

, and construct an extended steering matrix

of the M

− 1 distinct DOAs, θ

,...,θ

−1

.LetB =[

h].

It follows that B is also not full column rank. Because B is

a square matrix, it is also not full row rank. Therefore there

exists some non-zero c ∈ C

such that c

B = 0.Lett

jφ

for l =1, 2,...,M

, where φ

=(2πd

sin θ

)/λ. We can

express B as

⎡

⎢

⎣

11··· 10

··· t

−1

··· t

−1

··· t

−1

− 1)t

−2

⎤

⎥

⎦

We deﬁne the complex polynomial f (x)=



l=1

l−1

.It

can be observed that c

B = 0 is equivalent to f(t

)=0for

l =1, 2,...,M

− 1, and f



)=0. By construction, θ

are

distinct, so t

are M

− 1 different roots of f(x). Because c = 0,

f(x) is not a constant-zero polynomial, and has at most M

− 1

roots. Therefore each root t

has a multiplicity of at most one.

However, f



)=0implies that t

has a multiplicity of at least

two, which contradicts the previous conclusion and completes

the proof of β

= 0.

We now show that ξ

= 0. By the deﬁnition of F in

Appendix A, each row of F has at least one non-zero ele-

ment, and each column of F has at most one non-zero element.

Hence F

x = 0 for some x ∈ C

−1

if and only of x = 0.

It sufﬁces to show that Γ

(β

⊗ α

) = 0. By the deﬁnition of

Γ, we can rewrite Γ

(β

⊗ α

) as

, where

⎡

⎢

⎣

0 ··· 0

k(M

−1)

··· 0

··· β

0 β

··· β

k(M

−1)

00··· β

⎤

⎥

⎦

and β

is the l-th element of β

. Because β

= 0

and K<

is full column rank. By the deﬁnition of pseudo inverse,

we know that α

= 0. Therefore

= 0, which completes

the proof of ξ

= 0. 

One important implication of Theorem 1 is that DA-MUSIC

and SS-MUSIC share the same ﬁrst-order error expression, de-

spite the fact that R

is constructed from the second-order

statistics, while R

is constructed from the fourth-order statis-

tics. Theorem 1 enables a uniﬁed analysis of the MSEs of

DA-MUSIC and SS-MUSIC, which we present in Theorem 2.

Theorem 2: Under the same assumptions as in Theorem 1,

the asymptotic second-order statistics of the DOA estimation

errors by DA-MUSIC and SS-MUSIC share the same form:

E[(

− θ

)(

− θ

)]

R[ξ

(R ⊗ R

)ξ

]

. (10)

Proof: See Appendix C. 

By Theorem 2, it is straightforwardto write the uniﬁed asymp-

totic MSE expression as

(θ

(R ⊗ R

)ξ

. (11)

936 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 65, NO. 4, FEBRUARY 15, 2017

Therefore the MSE

depends on both the physical array ge-

ometry and the coarray geometry. The physical array geometry

is captured by A, which appears in R ⊗ R

. The coarray geom-

etry is captured by A

, which appears in ξ

and γ

. Therefore,

even if two arrays share the same coarray geometry, they may

not share the same MSE because their physical array geometry

maybedifferent.

It can be easily observed from (11) that (θ

)→0 as N →∞.

However, because p

appears in both the denominator and

numerator in (11), it is not obvious how the MSE varies

with respect to the source power p

and noise power σ

.Let

¯p

= p

/σ

denote the signal-to-noise ratio of the k-th source.

Let

P = diag(¯p

, ¯p

,...,¯p

), and

R = A

+ I. We can

then rewrite (θ

) as

(θ

(

R ⊗

)ξ

N ¯p

. (12)

Hence the MSE depends on the SNRs instead of the absolute

values of p

or σ

. To provide an intuitive understanding how

SNR affects the MSE, we consider the case when all sources

have the same power. In this case, we show in Corollary 1 that

the MSE asymptotically decreases as the SNR increases.

Corollary 1: Assume all sources have the same power p.Let

¯p = p/σ

denote the common SNR. Given sufﬁciently large N,

the MSE (θ

) decreases monotonically as ¯p increases, and

lim

¯p→∞

(θ

Nγ

ξ

(A ⊗ A

∗

)

. (13)

Proof: The limiting expression can be derived straightfor-

wardly from (12). For monotonicity, without loss of generality,

let p =1,so¯p =1/σ

. Because f(x)=1/x is monotonically

decreasing on (0, ∞), it sufﬁces to show that (θ

) increases

monotonically as σ

increases. Assume 0 <s

, and we

have

(θ

− (θ

Nγ

Qξ

where Q =(s

− s

)[(AA

) ⊗ I + I ⊗(AA

)+(s

)I]. Because AA

is positive semideﬁnite, both (AA

) ⊗ I

and I ⊗(AA

) are positive semideﬁnite. Combined with our

assumption that 0 <s

, we conclude that Q is positive

deﬁnite. By Proposition 1 we know that ξ

= 0. Therefore

Qξ

is strictly greater than zero, which implies the MSE

monotonically increases as σ

increases. 

Because both DA-MUSIC and SS-MUSIC work also in cases

when the number of sources exceeds the number of sensors, we

are particularly interested in their limiting performance in such

cases. As shown in Corollary 2, when K ≥ M, the correspond-

ing MSE is strictly greater than zero, even though the SNR

approaches inﬁnity. This corollary explains the “saturation” be-

havior of SS-MUSIC in the high SNR region as observed in [8]

and [9]. Another interesting implication of Corollary 2 is that

when 2 ≤ K<M, the limiting MSE is not necessarily zero.

Recall that in [2], it was shown that the MSE of the traditional

MUSIC algorithm will converge to zero as SNR approaches

inﬁnity. We know that both DA-MUSIC and SS-MUSIC will

For brevity, when we use the acronym “MSE” in the following discussion,

we refer to the asymptotic MSE, (θ

), unless explicitly stated.

be outperformed by traditional MUSIC in high SNR regions

when 2 ≤ K<M. Therefore, we suggest using DA-MUSIC

or SS-MUSIC only when K ≥ M.

Corollary 2: Following the same assumptions in Corollary 1,

1) When K =1, lim

¯p→∞

(θ

)=0;

2) When 2 ≤ K<M, lim

¯p→∞

(θ

) ≥ 0;

3) When K ≥ M, lim

¯p→∞

(θ

) > 0.

Proof: The right-hand side of (13) can be expanded into

Nγ



m =1



n=1

ξ

[a(θ

) ⊗ a

∗

(θ

)]

By the deﬁnition of F , F [a(θ

) ⊗ a

∗

(θ

)] becomes

j(M

−1)φ

j(M

−2)φ

,...,e

−j(M

−1)φ

Hence ΓF [a(θ

) ⊗ a

∗

(θ

)] = a

(θ

) ⊗ a

∗

(θ

). Observe

that

[a(θ

) ⊗ a

∗

(θ

)] = (β

⊗ α

)

(θ

) ⊗ a

∗

(θ

))

=(β

(θ

))(α

∗

(θ

))

(θ

)Π

⊥

(θ

))(α

∗

(θ

))

=0.

We can reduce the right-hand side of (13) into

Nγ



1≤m,n≤K

m =n

ξ

[a(θ

) ⊗ a

∗

(θ

)]

Therefore when K =1, the limiting expression is exactly zero.

When 2 ≤ K<M, the limiting expression is not necessary

zero because when m = n, ξ

[a(θ

) ⊗ a

∗

(θ

)] is not neces-

sarily zero.

When K ≥ M, A is full row rank. Hence A ⊗A

∗

is also

full row rank. By Proposition 1 we know that ξ

=0, which

implies that (θ

) is strictly greater than zero. 

IV. C

RAM

ER-RAO BOUND

The CRB for the unconditional model (2) has been well stud-

ied in [15], but only when the number of sources is less than

the number of sensors and no prior knowledge of P is given.

For the coarray model, the number of sources can exceed the

number of sensors, and P is assumed to be diagonal. Therefore,

the CRB derived in [15] cannot be directly applied. Based on

[26, Appendix 15C], we provide an alternative CRB based on

the signal model (2), under assumptions A1–A4.

For the signal model (2), the parameter vector is deﬁned by

η =[θ

,...,θ

,...,p

,σ

]

, (14)

and the (m, n)-th element of the Fisher information matrix

(FIM) is given by [15], [26]

FIM

= N tr



∂R

∂η

−1

∂R

∂η

−1



. (15)

WANG AND NEHORAI: COARRAYS, MUSIC, AND THE CRAM

ER–RAO BOUND 937

Observe that tr(AB)=vec(A

)

vec(B), and that vec

(AXB)=(B

⊗ A) vec(X). We can rewrite (15) as

FIM

= N



∂r

∂η



⊗ R)

−1

∂r

∂η

Denote the derivatives of r with respect to η as

∂r

∂η



∂r

∂θ

···

∂r

∂θ

∂r

∂p

···

∂r

∂p

∂r

∂σ



. (16)

The FIM can be compactly expressed by

FIM =



∂r

∂η



⊗ R)

−1

∂r

∂η

. (17)

According to (4), we can compute the derivatives in (16) and

obtain

∂r

∂η





, (18)

where

∗

 A + A

∗



A, A

and i follow the same def-

initions as in (4), and

A =



∂a(θ

)

∂θ

∂a(θ

)

∂θ

···

∂a(θ

)

∂θ



Note that (18) can be partitioned into two parts, speciﬁcally,

the part corresponding to DOAs and the part corresponding to

the source and noise powers. We can also partition the FIM.

Because R is positive deﬁnite, (R

⊗ R)

−1

is also positive

deﬁnite, and its square root (R

⊗ R)

−1/2

also exists. Let

=(R

⊗ R)

−1/2

P ,

=(R

⊗ R)

−1/2





We can write the partitioned FIM as

FIM = N





The CRB matrix for the DOAs is then obtained by block-wise

inversion:

CRB

⊥

)

−1

, (19)

where Π

⊥

= I −M

)

−1

. It is worth noting

that, unlike the classical CRB for the unconditional model in-

troduced in [15, Remark 1], expression (19) is applicable even

if the number of sources exceeds the number of sensors.

Remark 1: Similar to (11), CRB

depends on the SNRs in-

stead of the absolute values of p

or σ

.Let¯p

= p

/σ

, and

P = diag(¯p

, ¯p

,...,¯p

).Wehave

⊗

−1/2

P , (20)

= σ

−2

(

⊗

−1/2





. (21)

Substituting (20) and (21) into (19), the term σ

gets canceled,

and the resulting CRB

depends on the ratios ¯p

instead of

absolute values of p

or σ

Remark 2: The invertibility of the FIM depends on the coar-

ray structure. In the noisy case, (R

⊗ R)

−1

is always full

rank, so the FIM is invertible if and only if ∂r/∂η is full col-

umn rank. By (18) we know that the rank of ∂r/∂η is closely

related to A

, the coarray steering matrix. Therefore CRB

not valid for an arbitrary number of sources, because A

may

not be full column rank when too many sources are present.

Proposition 2: Assume all sources have the same power p,

and ∂r/∂η is full column rank. Let ¯p = p/σ

1IfK<M, and lim

¯p→∞

CRB

exists, it is zero under mild

conditions.

2IfK ≥ M, and lim

¯p→∞

CRB

exists,it is positivedeﬁnite.

Proof: See Appendix D. 

While inﬁnite SNR is unachievable from a practical stand-

point, Proposition 2 gives some useful theoretical implications.

When K<M, the limiting MSE (13) in Corollary 1 is not

necessarily zero. However, Proposition 2 reveals that the CRB

may approach zero when SNR goes to inﬁnity. This observa-

tion implies that both DA-MUSIC and SS-MUSIC may have

poor statistical efﬁciency in high SNR regions. When K ≥ M,

Proposition 2 implies that the CRB of each DOA will converge

to a positive constant, which is consistent with Corollary 2.

V. N UMERICAL ANALYSIS

In this section, we numerically analyze of DA-MUSIC and

SS-MUSIC by utilizing (11) and (19). We ﬁrst verify the MSE

expression (10) introduced in Theorem 2 through Monte Carlo

simulations. We then examine the application of (8) in predicting

the resolvability of two closely placed sources, and analyze the

asymptotic efﬁciency of both estimators from various aspects.

Finally, we investigate how the number of sensors affect the

asymptotic MSE.

In all experiments, we deﬁne the signal-to-noise ratio (SNR)

SNR = 10 log

min

k=1,2,...,K

where K is the number of sources.

Throughout Section V-A, V-B and V-C, we consider the fol-

lowing three different types of linear arrays with the following

sensor conﬁgurations:

Co-prime Array [5]: [0, 3, 5, 6, 9, 10, 12, 15, 20, 25]λ/2

Nested Array [9]: [1, 2, 3, 4, 5, 10, 15, 20, 25, 30]λ/2

MRA [27]: [0, 1, 4, 10, 16, 22, 28, 30, 33, 35]λ/2

All three arrays share the same number of sensors, but differ-

ence apertures.

A. Numerical Veriﬁcation

We ﬁrst verify (11) via numerical simulations. We consider 11

sources with equal power, evenly placed between −67.50

◦

and

56.25

◦

, which is more than the number of sensors. We compare

the difference between the analytical MSE and the empirical

MSE under different combinations of SNR and snapshot num-

bers. The analytical MSE is deﬁned by

MSE



k=1

(θ

938 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 65, NO. 4, FEBRUARY 15, 2017

Fig. 2. |MSE

− MSE

|/MSE

for different types of arrays under

different numbers of snapshots and different SNRs.

and the empirical MSE is deﬁned by

MSE



l=1



k=1



(l)

− θ

(l)



where θ

(l)

is the k-th DOA in the l-th trial, and

(l)

is the

corresponding estimate.

Fig. 2 illustrates the relative errors between MSE

and

MSE

obtained from 10,000 trials under various scenarios.

It can be observed that MSE

and MSE

agree very well

given enough snapshots and a sufﬁciently high SNR. It should

be noted that at 0 dB SNR, (8) is quite accurate when 250

snapshots are available. In addition. there is no signiﬁcant dif-

ference between the relative errors obtained from DA-MUSIC

and those from SS-MUSIC. These observations are consistent

with our assumptions, and verify Theorem 1 and Theorem 2.

We observe that in some of the low SNR regions, |MSE

−

MSE

|/MSE

appears to be smaller even if the number of

snapshots is limited. In such regions, MSE

actually “satu-

rates”, and MSE

happens to be close to the saturated value.

Therefore, this observation does not imply that (11) is valid in

such regions.

B. Prediction of Resolvability

One direct application of Theorem 2 is predicting the re-

solvability of two closely located sources. We consider two

sources with equal power, located at θ

=30

◦

− Δθ/2, and

=30

◦

+Δθ/2, where Δθ varies from 0.3

◦

to 3.0

◦

.Wesay

the two sources are correctly resolved if the MUSIC algorithm

is able to identify two sources, and the two estimated DOAs

satisfy |

− θ

| < Δθ/2,fori ∈{1, 2}. The probability of res-

olution is computed from 500 trials. For all trials, the number of

Fig. 3. Probability of resolution vs. source separation, obtained from 500

trials. The number of snapshots is ﬁxed at 500, and the SNR is set to 0 dB.

snapshots is ﬁxed at 500, the SNR is set to 0 dB, and SS-MUSIC

is used.

For illustration purpose, we analytically predict the resolv-

ability of the two sources via the following simple criterion:

(θ

)+(θ

)

Unresovalble



Resolvable

Δθ. (22)

Readers are directed to [28] for a more comprehensive criterion.

Fig. 3 illustrates the resolution performance of the three ar-

rays under different Δθ, as well as the thresholds predicted by

(22). The MRA shows best resolution performance of the three

arrays, which can be explained by the fact that the MRA has the

largest aperture. The co-prime array, with the smallest aperture,

shows the worst resolution performance. Despite the differences

in resolution performance, the probability of resolution of each

array drops to nearly zero at the predicted thresholds. This con-

ﬁrms that (11) provides a convenient way of predicting the

resolvability of two close sources.

C. Asymptotic Efﬁciency Study

In this section, we utilize (11) and (19) to study the asymp-

totic statistical efﬁciency of DA-MUSIC and SS-MUSIC under

different array geometries and parameter settings. We deﬁne

their average efﬁciency as

κ =

tr CRB



k=1

(θ

)

. (23)

For efﬁcient estimators we expect κ =1, while for inefﬁcient

estimators we expect 0 ≤ κ<1.

We ﬁrst compare the κ value under different SNRs for

the three different arrays. We consider three cases: K =1,

K =6, and K =12.TheK sources are located at {−60

◦

[120(k −1)/(K −1)]

◦

|k =1, 2,...,K}, and all sources have

the same power. As shown in Fig. 4(a), when only one source is

present, κ increases as the SNR increases for all three arrays.

However, none of the arrays leads to efﬁcient DOA estimation.

Interestingly, despite being the least efﬁcient geometry in the

low SNR region, the co-prime array achieves higher efﬁciency

than the nested array in the high SNR region. When K =6,

we can observe in Fig. 4(b) that κ decreases to zero as SNR

increases. This rather surprising behavior suggests that both

DA-MUSIC and SS-MUSIC are not statistically efﬁcient meth-

ods for DOA estimation when the number of sources is greater

than one and less than the number of sensors. It is consistent

WANG AND NEHORAI: COARRAYS, MUSIC, AND THE CRAM

ER–RAO BOUND 939

Fig. 4. Average efﬁciency vs. SNR: (a) K =1,(b)K =6,(c)K =12.

with the implication of Proposition 2 when K<M. When

K =12, the number of sources exceeds the number of sensors.

We can observe in Fig. 4(c) that κ also decreases as SNR in-

creases. However, unlike the case when K =6, κ converges to

a positive value instead of zero.

The above observations imply that DA-MUSIC and

SS-MUSIC achieve higher degrees of freedom at the cost of

decreased statistical efﬁciency. When statistical efﬁciency is

concerned and the number of sources is less than the number

of sensors, one might consider applying MUSIC directly to the

original sample covariance R deﬁned in (3) [29].

Next, we then analyze how κ is affected by angular separa-

tion. Two sources located at −Δθ and Δθ are considered. We

compute the κ values under different choices of Δθ for all three

arrays. For reference, we also include the empirical results ob-

tained from 1000 trials. To satisfy the asymptotic assumption,

the number of snapshots is ﬁxed at 1000 for each trial. As shown

in Fig. 5(a)–(c), the overall statistical efﬁciency decreases as the

SNR increases from 0 dB to 10 dB for all three arrays, which

is consistent with our previous observation in Fig. 4(b). We can

also observe that the relationship between κ and the normalized

angular separation Δθ/π is rather complex, as opposed to the

traditional MUSIC algorithm (c.f. [2]). The statistical efﬁciency

Fig. 5. Average efﬁciency vs. angular separation for the co-prime array:

(a) MRA, (b) nested array, (c) co-prime array. The solid lines and dashed lines

are analytical values obtained from (23). The circles and crosses are emprical

results averaged from 1000 trials.

of DA-MUSIC and SS-MUSIC is highly dependent on array

geometry and angular separation.

D. MSE vs. Number of Sensors

In this section, we investigate how the number of sensors

affect the asymptotic MSE, (θ

). We consider three types of

sparse linear arrays: co-prime arrays, nested arrays, and MRAs.

In this experiment, the co-prime arrays are generated by co-

prime pairs (q,q +1) for q =2, 3,...,12. The nested arrays

are generated by parameter pairs (q +1,q) for q =2, 3,...,12.

The MRAs are constructed according to [27]. We consider two

cases: the one source case where K =1, and the under deter-

mined case where K = M. For the former case, we placed the

only source at the 0

◦

. For the later case, we placed the sources

uniformly between −60

◦

and 60

◦

.WesetSNR = 0 dB and

N = 1000. The empirical MSEs were obtained from 500 trials.

SS-MUSIC was used in all the trials.

In Fig. 6(a), we observe that when K =1, the MSE de-

creases at a rate of approximately O(M

−4.5

) for all three arrays.

940 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 65, NO. 4, FEBRUARY 15, 2017

Fig. 6. MSE vs. M :(a)K =1,(b)K = M . The solid lines are analytical

results. The “+”, “◦”, and “” denote empirical results obtains from 500 trials.

The dashed lines are trend lines used for comparison.

In Fig. 6(b), we observe that when K = M, the MSE only de-

creases at a rate of approximately O(M

−3.5

). In both cases, the

MRAs and the nested arrays achieve lower MSE than the co-

prime arrays. Another interesting observation is that for all three

arrays, the MSE decreases faster than O(M

−3

). Recall that for a

M-sensor ULA, the asymptotic MSE of traditional MUSIC de-

creases at a rate of O(M

−3

) as M →∞[2]. This observation

suggests that given the same number of sensors, these sparse

linear arrays can achieve higher estimation accuracy than ULAs

when the number of sensors is large.

VI. C

ONCLUSION

In this paper, we reviewed the coarray signal model and de-

rived the asymptotic MSE expression for two coarray-based

MUSIC algorithms, namely DA-MUSIC and SS-MUSIC. We

theoretically proved that the two MUSIC algorithms share the

same asymptotic MSE error expression. Our analytical MSE

expression is more revealing and can be applied to various types

of sparse linear arrays, such as co-prime arrays, nested arrays,

and MRAs. In addition, our MSE expression is also valid when

the number of sources exceeds the number of sensors. We also

derived the CRB for sparse linear arrays, and analyzed the sta-

tistically efﬁciency of typical sparse linear arrays. Our results

will beneﬁt to future research on performance analysis and opti-

mal design of sparse linear arrays. Throughout our derivations,

we assume the array is perfectly calibrated. In the future, it will

be interesting to extend the results in this paper to cases when

model errors are present. Additionally, we will further investi-

gate how the number of sensors affect the MSE and the CRB for

sparse linear arrays, as well as the possibility of deriving closed

form expressions in the case of large number of sensors.

PPENDIX A

EFINITION AND PROPERTIES OF THE COARRAY

SELECTION MATR IX

According to (3),



k=1

exp[j(

−

)φ

]+δ

where δ

denotes Kronecker’s delta. This equation implies

that the (m, n)-th element of R is associated with the difference

(

−

). To capture this property, we introduce the difference

matrix Δ such that Δ

−

. We also deﬁne the weight

function ω(n):Z → Z as (see [9] for details)

ω(l)=|{(m, n)|Δ

= l}|,

where |A| denotes the cardinality of the set A . Intuitively, ω(l)

counts the number of all possible pairs of (

) such that

−

= l. Clearly, ω(l)=ω(−l).

Deﬁnition 1: The coarray selection matrix F is a (2M

−

1) × M

matrix satisfying

m,p+(q−1)M

⎧

⎨

⎩

ω(m − M

)

, Δ

= m − M

0, otherwise,

(24)

for m =1, 2,...,2M

− 1,p=1, 2,...,M,q =1, 2,...,M.

To better illustrate the construction of F , we consider a toy

array whose sensor locations are given by {0,d

, 4d

}.The

corresponding difference matrix of this array is

Δ =

⎡

⎣

0 −1 −4

10−3

43 0

⎤

⎦

The ULA part of the difference coarray consists of three sensors

located at −d

,0,andd

. The weight function satisﬁes ω(−1) =

ω(1) = 1, and ω(0) = 3,soM

=2. We can write the coarray

selection matrix as

F =

⎡

⎢

⎣

000100000

000

010000000

⎤

⎥

⎦

If we pre-multiply the vectorized sample covariance matrix r by

F , we obtain the observation vector of the virtual ULA (deﬁned

in (5)):

z =

⎡

⎣

⎤

⎦

⎡

⎢

⎣

+ R

)

⎤

⎥

⎦

It can be seen that z

is obtained by averaging all the el-

ements in R that correspond to the difference m − M

,for

m =1, 2,...,2M

− 1.

Based on Deﬁnition 1, we nowderive several useful properties

of F .

Lemma 1: F

m,p+(q−1)M

= F

−m,q+(p−1)M

for m =

1, 2,...,2M

− 1,p=1, 2,...,M,q =1, 2,...,M.

Proof: If F

m,p+(q−1)M

=0, then Δ

= m − M

.Be-

cause Δ

= −Δ

, Δ

= −(m − M

). Hence (2M

− m) − M

= −(m − M

) =Δ

, which implies that

−m,q+(p−1)M

is also zero.

If F

m,p+(q−1)M

=0, then Δ

=m−M

and F

m,p+(q−1)

M =

1/ω(m − M

). Note that (2M

− m) − M

= −(m−M

−Δ

=Δ

. We thus have F

−m,q+(p−1)M

=1/ω(−(m −

)) = 1/ω(m − M

)=F

m,p+(q−1)M

. 

Lemma 2: Let R ∈ C

be Hermitian symmetric. Then z =

F vec(R) is conjugate symmetric.

WANG AND NEHORAI: COARRAYS, MUSIC, AND THE CRAM

ER–RAO BOUND 941

Proof: By Lemma 1 and R = R



p=1



q=1

m,p+(q−1)M



q=1



p=1

−m,q+(p−1)M

∗

= z

∗

−m

. 

Lemma 3: Let z ∈ C

−1

be conjugate symmetric. Then

mat

M,M

z) is Hermitian symmetric.

Proof: Let H = mat

M,M

z). Then

−1



m =1

m,p+(q−1)M

. (25)

We know that z is conjugate symmetric, so z

= z

∗

−m

Therefore, by Lemma 1

−1



m =1

∗

−m

−m,q+(p−1)M



−1





,q+(p−1)M



∗

= H

∗

. 

(26)

PPENDIX B

ROOF OF THEOREM 1

We ﬁrst derive the ﬁrst-order expression of DA-MUSIC. De-

note the eigendecomposition of R

= E

+ E

where E

and E

are eigenvectors of the signal subspace and

noise subspace, respectively, and Λ

, Λ

are the correspond-

ing eigenvalues. Speciﬁcally, we have Λ

= σ

Let

= R

+ΔR

= E

+ΔE

, and

+ΔΛ

be the perturbed versions of R

, E

, and Λ

The following equality holds:

+ΔR

)(E

+ΔE

)=(E

+ΔE

)(Λ

+ΔΛ

If the perturbation is small, we can omit high-order terms and

obtain [16], [23], [25]

ΔE

= −P

−1

†

ΔR

. (27)

Because P is diagonal, for a speciﬁc θ

,wehave

(θ

)ΔE

= −p

−1

†

ΔR

, (28)

where e

is the k-th column of the identity matrix I

K × K

. Based

on the conclusion in Appendix B of [2], under sufﬁciently small

perturbations, the error expression of DA-MUSIC for the k-th

DOA is given by

(1)

− θ

= −

R[a

(θ

)ΔE

(θ

))]

(θ

)

(29)

where

(θ

)=∂a

(θ

)/∂θ

Substituting (28) into (29) gives

(1)

− θ

= −

R[e

†

ΔR

(θ

)]

(θ

)

. (30)

Because vec(AXB)=(B

⊗ A) vec(X) and E

⊥

, we can use the notations introduced in (9b)–(9d) to ex-

press (30) as

(1)

− θ

= −(γ

)

−1

R[(β

⊗ α

)

Δr

], (31)

where Δr

=vec(ΔR

Note that

is constructed from

R. It follows that ΔR

actually depends on ΔR, which is the perturbation part of the

covariance matrix R. By the deﬁnition of R

Δr

=vec



Δz ··· Γ

Δz Γ

Δz



= ΓF Δr,

where Γ =[Γ

−1

···Γ

]

and Δr = vec(ΔR).

Let ξ

= F

(β

⊗ α

). We can now express (31) in

terms of Δr as

(1)

− θ

= −(γ

)

−1

R(ξ

Δr), (32)

which completes the ﬁrst part of the proof.

We next consider the ﬁrst-order error expression of SS-

MUSIC. From (7) we know that R

shares the same eigen-

vectors as R

. Hence the eigendecomposition of R

can be

expressed by

= E

+ E

where Λ

and Λ

are the eigenvalues of the signal subspace

and noise subspace. Speciﬁcally, we have Λ

= σ

I.Note

that R

=(A

+ σ

. Following a similar ap-

proach to the one we used to obtain (27), we get

ΔE

= −M

−1

(PA

+2σ

−1

†

ΔR

where ΔE

is the perturbation of the noise eigenvectors pro-

duced by ΔR

. After omitting high-order terms, ΔR

given by

ΔR



k=1

Δz

+Δz

According to [9], each subarray observation vector z

can be

expressed by

= A

−k

p + σ

−k+1

, (33)

for k =1, 2,...,M

, where i

is a vector of length M

whose

elements are zero except for the l-th element being one, and

Ψ = diag(e

−jφ

,...,e

−jφ

Observe that



k=1

−k+1

Δz

= σ

ΔR

942 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 65, NO. 4, FEBRUARY 15, 2017

and



k=1

−k

pΔz

= A

⎡

⎢

⎣

−j(M

−1)φ

−j(M

−2)φ

··· 1

−j(M

−1)φ

−j(M

−2)φ

··· 1

−j(M

−1)φ

−j(M

−2)φ

··· 1

⎤

⎥

⎦

⎡

⎢

⎣

Δz

⎤

⎥

⎦

= A

P (T

)

ΔR

= A

ΔR

where T

is a M

× M

permutation matrix whose anti-

diagonal elements are one, and whose remaining elements are

zero. Because ΔR =ΔR

, by Lemma 2 we know that Δz is

conjugate symmetric. According to the deﬁnition of R

,itis

straightforward to show that ΔR

=ΔR

also holds. Hence

ΔR



+2σ

I)ΔR

+ΔR



Substituting ΔR

into the expression of A

ΔE

, and uti-

lizing the property that A

= 0, we can express A

ΔE

−P

−1

(PA

+2σ

−1

†

+2σ

I)ΔR

Observe that

†

+2σ

I)=(A

)

−1

+2σ

=[PA

+2σ

)

−1

]

=(PA

+2σ

I)A

†

Hence the term (PA

+2σ

I) gets canceled and we obtain

ΔE

= −P

−1

†

ΔR

, (34)

which coincides with the ﬁrst-order error expression of

ΔE

PPENDIX C

ROOF OF THEOREM 2

Before proceeding to the main proof, we introduce the fol-

lowing deﬁnition.

Deﬁnition 2: Let A =[a

...a

] ∈ R

N × N

, and B =

...b

] ∈ R

N × N

. The structured matrix C

∈

× N

is deﬁned as

⎡

⎢

⎣

... a

⎤

⎥

⎦

We now start deriving the explicit MSE expression. Accord-

ingto(32),

E[(

− θ

)(

− θ

)]

=(γ

)

−1

(γ

)

−1

E[R(ξ

Δr) R(ξ

Δr)]

=(γ

)

−1

(γ

)

−1



R(ξ

)

E[R(Δr) R(Δr)

]

× R(ξ

)

+ I(ξ

)

E[I(Δr) I(Δr)

] I(ξ

)

− R(ξ

)

E[R(Δr) I(Δr)

] I(ξ

)

− R(ξ

)

E[R(Δr) I(Δr)

] I(ξ

)



, (35)

where we used the property that R(AB)=R(A) R(B) −

I(A) I(B) for two complex matrices A and B with proper

dimensions.

To obtain the closed-form expression for (35), we need to

compute the four expectations. It should be noted that in the case

of ﬁnite snapshots, Δr does not follow a circularly-symmetric

complex Gaussian distribution. Therefore we cannot directly

use the properties of the circularly-symmetric complex Gaus-

sian distribution to evaluate the expectations. For brevity, we

demonstrate the computation of only the ﬁrst expectation in

(35). The computation of the remaining three expectations fol-

lows the same idea.

Let r

denote the i-th column of R in (3). Its estimate,

is given by



t=1

y(t)y

∗

(t), where y

(t) is the i-th element of

y(t). Because E[

]=r

E[R(Δr

) R(Δr

)

]

= E[R(

) R(

)

] − R(r

) R(r

)

(36)

The second term in (36) is deterministic, and the ﬁrst term in

(36) can be expanded into







s=1

y(s)y

∗

(s)







t=1

y(t)y

∗

(t)









s=1



t=1

R(y(s)y

∗

(s)) R(y(t)y

∗

(t))





s=1



t=1





R(y(s)) R(y

∗

(s)) − I(y(s)) I(y

∗

(s))





R(y(t))

R(y

∗

(t)) − I(y(t))

I(y

∗

(t))







s=1



t=1



E[R(y(s)) R(y

(s)) R(y(t))

R(y

(t))]

+ E[R(y(s)) R(y

(s)) I(y(t))

I(y

(t))]

+ E[I(y(s)) I(y

(s)) R(y(t))

R(y

(t))]

+ E[I(y(s)) I(y

(s)) I(y(t))

I(y

(t))]



. (37)

We ﬁrst consider the partial sum of the cases when s = t.By

A4, y(s) and y(t) are uncorrelated Gaussians. Recall that for

WANG AND NEHORAI: COARRAYS, MUSIC, AND THE CRAM

ER–RAO BOUND 943

x ∼CN(0, Σ),

E[R(x) R(x)

R(Σ), E[R(x) I(x)

]=−

I(Σ)

E[I(x) R(x)

I(Σ), E[I(x) I(x)

R(Σ).

We have

E[R(y(s)) R(y

(s)) R(y(t))

R(y

(t))]

= E[R(y(s)) R(y

(s))]E[R(y(t))

R(y

(t))]

R(r

) R(r

)

Similarly, we can obtain that when s = t,

E[R(y(s)) R(y

(s)) I(y(t))

I(y

(t))] =

R(r

) R(r

)

E[I(y(s)) I(y

(s)) R(y(t))

R(y

(t))] =

R(r

) R(r

)

E[I(y(s)) I(y

(s)) I(y(t))

I(y

(t))] =

R(r

) R(r

)

(38)

Therefore the partial sum of the cases when s = t is given by

(1 − 1/N ) R(r

) R(r

)

We now consider the partial sum of the cases when s = t.We

ﬁrst consider the ﬁrst expectation inside the double summation

in (37). Recall that for x ∼N(0, Σ), E[x

]=σ

+ σ

. We can express the (m, n)-th element of the

matrix E[R(y(t)) R(y

(t)) R(y(t))

R(y

(t))] as

E[R(y

(t)) R(y

(t))]

= E[R(y

(t)) R(y

(t))]

= E[R(y

(t)) R(y

(t))]E[R(y

(t)) R(y

(t))]

+ E[R(y

(t)) R(y

(t))]E[R(y

(t)) R(y

(t))]

+ E[R(y

(t)) R(y

(t))]E[R(y

(t)) R(y

(t))]

[R(R

) R(R

)+ R(R

) R(R

)+R(R

) R(R

)].

Hence

E[R(y(t)) R(y

(t)) R(y(t))

R(y

(t))]

[R(r

) R(r

)

+ R(r

) R(r

)

+ R(R) R(R

)].

Similarly, we obtain that

E[I(y(t)) I(y

(t)) I(y(t))

I(y

(t))]

[R(r

) R(r

)

+ R(r

) R(r

)

+ R(R) R(R

)],

E[R(y(t)) R(y

(t)) I(y(t))

I(y

(t))]

= E[I(y(t)) I(y

(t)) R(y(t))

R(y

(t))]

[R(r

) R(r

)

− I(r

) I(r

)

+ I(R) I(R

)].

Therefore the partial sum of the cases when s = t

is given by (1/N ) R(r

) R(r

)

+(1/2N)[R(R) R(R

)+ I

(R) I(R

)+R(r

) R(r

)

− I(r

) I(r

)

]. Combined with

the previous partial sum of the cases when s = t, we obtain

that

E[R(Δr

) R(Δr

)

]

[R(R) R(R

)+I(R) I(R

)

+ R(r

) R(r

)

− I(r

) I(r

)

(39)

Therefore

E[R(Δr) R(Δr)

]

[R(R) ⊗ R(R)+I(R) ⊗ I(R)

+ C

R(R) R(R)

− C

I(R) I(R)

(40)

which completes the computation of ﬁrst expectation in (35).

Utilizing the same technique, we obtain that

E[I(Δr) I(Δr)

]

[R(R) ⊗ R(R)+I(R) ⊗ I(R)

+ C

I(R) I(R)

− C

R(R) R(R)

], (41)

and

E[R(Δr) I(Δr)

]

[I(R) ⊗ R(R) − R(R) ⊗I(R)

+ C

R(R) I(R)

+ C

I(R) R(R)

]. (42)

Substituting (40)–(42) into (35) gives a closed-form MSE ex-

pression. However, this expression is too complicated for analyt-

ical study. In the following steps, we make use of the properties

of ξ

to simply the MSE expression.

Lemma 4: Let X, Y , A, B ∈ R

N × N

satisfying X

(−1)

X, A

=(−1)

A, and B

=(−1)

B, where

∈{0, 1}. Then

vec(X)

(A ⊗ B) vec(Y )

=(−1)

vec(X)

vec(Y ),

vec(X)

(B ⊗ A) vec(Y )

=(−1)

vec(X)

vec(Y ).

Proof: By Deﬁnition 2,

vec(X)

vec(Y )



m =1



n=1



m =1



n=1





p=1





p=1





m =1



n=1



p=1



q=1

944 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 65, NO. 4, FEBRUARY 15, 2017

=(−1)



p=1



n=1



m =1



q=1

=(−1)



p=1



n=1

=(−1)

vec(X)

(A ⊗ B) vec(Y ).

The proof of the second equality follows the same idea. 

Lemma 5: T

⊥

=(Π

⊥

)

∗

Proof: Since Π

⊥

= I −A

)

−1

, it sufﬁces

to show that T

)

−1

=(A

)

−1

)

∗

. Because A

is the steering matrix of a ULA with

sensors, it is straightforward to show that T

Φ)

∗

, where Φ = diag(e

−j(M

−1)φ

−j(M

−1)φ

,...,

−j(M

−1)φ

Because T

= I, T

= T

)

−1

= T

)

−1

=(A

Φ)

∗

((A

Φ)

∗

)

−1

Φ)

=(A

)

−1

)

∗



Lemma 6: Let Ξ

= mat

M,M

(ξ

). Then Ξ

= Ξ

for k =

1, 2,...,K.

Proof: Note that ξ

= F

(β

⊗ α

). We ﬁrst prove that

⊗ α

is conjugate symmetric, or that (T

⊗ T

)(β

⊗

)=(β

⊗ α

)

∗

. Similar to the proof of Lemma 5, we utilize

the properties that T

=(A

Φ)

∗

and that T

(θ

−j(M

−1)φ

)

∗

to show that

†

)

(θ

=[(A

†

)

(θ

)]

∗

. (43)

Observe that

(θ

)=j

(θ

), where

=(2πd

cos θ

)/λ and D = diag(0, 1,...,M

− 1).Wehave

⊗ T

)(β

⊗ α

)=(β

⊗ α

)

∗

⇐⇒ T

=(α

)

∗

⇐⇒ T

[(A

†

)

(θ

)DΠ

⊥

]

∗

= −(A

†

)

(θ

)DΠ

⊥

Since D = T

, combining with Lemma 5

and (43), it sufﬁces to show that

†

)

(θ

⊥

= −(A

†

)

(θ

)DΠ

⊥

. (44)

Observe that T

+ D =(M

− 1)I.Wehave

⊥

+ D)a

(θ

)=0,

or equivalently

(θ

⊥

= −a

(θ

)DΠ

⊥

. (45)

Pre-multiplying both sides of (45) with (A

†

)

leads to

(44), which completes the proof that β

⊗ α

is conjugate

symmetric. According to the deﬁnition of Γ in(9e),itis

straightforward to show that Γ

(β

⊗ α

) is also conjugate

symmetric. Combined with Lemma 3 in Appendix A, we con-

clude that mat

M,M

(β

⊗ α

)) is Hermitian symmet-

ric, or that Ξ

= Ξ

. 

Given Lemma 4–6, we are able continue the simpliﬁcation.

We ﬁrst consider the term R(ξ

)

E[R(Δr) R(Δr)

] R(ξ

)

in (35). Let Ξ

= mat

M,M

(ξ

), and Ξ

= mat

M,M

(ξ

By Lemma 6, we have Ξ

= Ξ

, and Ξ

= Ξ

. Observe

that R(R)

= R(R), and that I(R)

= I(R). By Lemma 4

we immediately obtain the following equalities:

R(ξ

)

(R(R) ⊗ R(R)) R(ξ

)

= R(ξ

)

R(R) R(R)

R(ξ

)

(I(R) ⊗ I(R)) R(ξ

)

= −R(ξ

)

I(R) I(R)

R(ξ

Therefore R(ξ

)

E[R(Δr) R(Δr)

] R(ξ

) can be com-

pactly expressed as

R(ξ

)

E[R(Δr) R(Δr)

] R(ξ

)

R(ξ

)

[R(R) ⊗ R(R)+I(R) ⊗ I(R)] R(ξ

)

R(ξ

)

R(R

⊗ R) R(ξ

), (46)

where we make use of the properties that R

= R

∗

, and

R(R

∗

⊗ R)=R(R) ⊗ R(R)+I(R) ⊗ I(R). Similarly, we

can obtain that

I(ξ

)

E[I(Δr) I(Δr)

] I(ξ

)

I(ξ

)

R(R

⊗ R) I(ξ

), (47)

R(ξ

)

E[R(Δr) I(Δr)

] I(ξ

)

= −

R(ξ

)

I(R

⊗ R) I(ξ

), (48)

R(ξ

)

E[R(Δr) I(Δr)

] I(ξ

)

= −

R(ξ

)

I(R

⊗ R) I(ξ

). (49)

Substituting (46)–(49) into (35) completes the proof.

PPENDIX D

ROOF OF PROPOSITION 2

Without loss of generality, let p =1and σ

→ 0. For brevity,

we denote R

⊗ R by W . We ﬁrst consider the case when K<

M. Denote the eigendecomposition of R

−1

by E

−1

−2

.Wehave

−1

= σ

−4

+ σ

−2

+ K

WANG AND NEHORAI: COARRAYS, MUSIC, AND THE CRAM

ER–RAO BOUND 945

where

= E

∗

⊗ E

= E

∗

−1

⊗ E

+ E

∗

⊗ E

−1

= E

∗

−1

⊗ E

−1

Recall that A

= 0.Wehave

=(E

∗

⊗ E

)(

∗

 A + A

∗



= E

∗

 E

A + E

∗

 E

= 0. (50)

Therefore

−1

= σ

−2

+ σ

)

. (51)

Similar to W

−1

, we denote W

−

= σ

−2

+ σ

−1

+ K

where

= E

∗

−

⊗ E

+ E

∗

⊗ E

−

= E

∗

−

⊗ E

−

Therefore

−

= σ

−2

(σ

+ K

)Π

(σ

+ K

)

where Π

= M

†

. We can then express the CRB as

CRB

= σ

+ σ

)

−1

, (52)

where

− K

)

= −

+ K

)

− K

)

When σ

=0, R reduces to AA

. Observe that the eigende-

composition of R always exists for σ

≥ 0.WeuseK



–K



denote the corresponding K

–K

when σ

→ 0.

Lemma 7: Let K<M. Assume ∂r/∂η is full column rank.

Then lim

→0

exists.

Proof: Because A

= 0,

=(E

∗

−1

⊗ E

)(A

∗

 A)

+(E

∗

⊗ E

−1

)(A

∗

 A)

= E

∗

−1

∗

 E

+ E

∗

 E

−1

= 0

Similarly, we can show that K

= 0, i

i = i

i =

0, and i

i =rank(E

)=M − K. Hence



−1



Because ∂r/∂η is full column rank, M

is full rank and

positive deﬁnite. Therefore the Schur complements exist, and

we can inverse M

block-wisely. Let V = A

and v = i

−1

i. After tedious but straightforward computa-

tion, we obtain

= K

−1

− s

−1

+ σ

−2

)

− v

−1

+ σ

−2

)ii

−1

+ s

−1

+ σ

−2

)ii

+ σ

−2

where S and s are Schur complements given by

S = V − v

−1

s = v −i

−1

Observe that

v = i

−1

i = σ

−4

(M − K)+i

We know that both v

−1

and s

−1

decrease at the rate of σ

.As

→ 0,wehave

S → A



−1

+ σ

−2

) → 0,

−1

+ σ

−2

) → 0,

−1

+ σ

−2

)ii

+ σ

−2

) →



M − K

We now show that A



is nonsingular. Denote the

eigendecomposition of AA

by E



)

. Recall that

for matrices with proper dimensions, (A  B)

(C  D)=

C) ◦ (B

D), where ◦ denotes the Hadamard product.

We can expand A



into





(Λ



)

−1



)



∗

◦





(Λ



)

−1



)



Note that AA



(Λ



)

−1



)

A = E



)

A = A, and

that A is full column rank when K<M. We thus

have A



(Λ



)

−1



)

A = I. Therefore A



= I,

which is nonsingular.

Combining the above results, we obtain that when σ

→ 0,

→ K



M − K



For sufﬁciently small σ

> 0, it is easy to show that K

–K

are bounded in the sense of Frobenius norm (i.e., K



≤ C

for some C>0,fori ∈{1, 2, 3, 4, 5}). Because ∂r/∂η is

full rank, M

is also full rank for any σ

> 0, which implies

that Π

is well-deﬁned for any σ

> 0. Observe that Π

is positive semideﬁnite, and that tr(Π

) = rank(M

).We

know that Π

is bounded for any σ

> 0. Therefore Q

and

are also bounded for sufﬁciently small σ

, which implies

that σ

+ σ

→ 0 as σ

→ 0.

By Lemma 7, we know that Q

→ Q



as σ

→ 0, where



− K



)

946 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 65, NO. 4, FEBRUARY 15, 2017

and Π



= lim

→0

as derived in Lemma 7. As-

sume Q



is nonsingular

. By (52) we immediately obtain that

CRB

→ 0 as σ

→ 0.

When K ≥ M, R is full rank regardless of the choice of

. Hence (R

⊗ R)

−1

is always full rank. Because ∂r/∂η is

full column rank, the FIM is positive deﬁnite, which implies its

Schur complements are also positive deﬁnite. Therefore CRB

is positive deﬁnite.

EFERENCES

[1] R. Schmidt, “Multiple emitter location and signal parameter estimation,”

IEEE Trans. Antennas Propag., vol. 34, no. 3, pp. 276–280, Mar. 1986.

[2] P. Stoica and A. Nehorai, “MUSIC, maximum likelihood, and Cramer-

Rao bound,” IEEE Trans. Acoust., Speech, Signal Process., vol. 37, no. 5,

pp. 720–741, May 1989.

[3] B. Ottersten, M. Viberg, P. Stoica, and A. Nehorai, “Exact and large sample

maximum likelihood techniques for parameter estimation and detection in

array processing,” in Radar Array Processing, T. S. Huang, T. Kohonen,

M. R. Schroeder, H. K. V. Lotsch, S. Haykin, J. Litva, and T. J. Shepherd,

Eds. Berlin, Germany: Springer, 1993, vol. 25, pp. 99–151.

[4] A. Moffet, “Minimum-redundancy linear arrays,” IEEE Trans. Antennas

Propag., vol. 16, no. 2, pp. 172–175, Mar. 1968.

[5] P. Pal and P. P. Vaidyanathan, “Coprime sampling and the MUSIC algo-

rithm,” in Proc. 2011 IEEE Digital Signal Process. Workshop IEEE Signal

Process. Education Workshop, Jan. 2011, pp. 289–294.

[6] Z. Tan and A. Nehorai, “Sparse direction of arrival estimation using co-

prime arrays with off-grid targets,” IEEE Signal Process. Lett., vol. 21,

no. 1, pp. 26–29, Jan. 2014.

[7] Z. Tan, Y. C. Eldar, and A. Nehorai, “Direction of arrival estimation

using co-prime arrays: A super resolution viewpoint,” IEEE Trans. Signal

Process., vol. 62, no. 21, pp. 5565–5576, Nov. 2014.

[8] S. Qin, Y. Zhang, and M. Amin, “Generalized coprime array conﬁgurations

for direction-of-arrival estimation,” IEEE Trans. Signal Process., vol. 63,

no. 6, pp. 1377–1390, Mar. 2015.

[9] P. Pal and P. Vaidyanathan, “Nested arrays: A novel approach to array pro-

cessing with enhanced degrees of freedom,” IEEE Trans. Signal Process.,

vol. 58, no. 8, pp. 4167–4181, Aug. 2010.

[10] K. Han and A. Nehorai, “Wideband gaussian source processing using a

linear nested array,” IEEE Signal Process. Lett., vol. 20, no. 11, pp. 1110–

1113, Nov. 2013.

[11] K. Han and A. Nehorai, “Nested vector-sensor array processing via tensor

modeling,” IEEE Trans. Signal Process., vol. 62, no. 10, pp. 2542–2553,

May 2014.

[12] B. Friedlander, “The root-MUSIC algorithm for direction ﬁnding with

interpolated arrays,” Signal Process., vol. 30, no. 1, pp. 15–29, Jan. 1993.

[13] M. Pesavento, A. Gershman, and M. Haardt, “Unitary root-MUSIC with

a real-valued eigendecomposition: A theoretical and experimental perfor-

mance study,” IEEE Trans. Signal Process., vol. 48, no. 5, pp. 1306–1314,

May 2000.

[14] P. Stoica and A. Nehorai, “MUSIC, maximum likelihood, and Cramer–Rao

bound: Further results and comparisons,” IEEE Trans. Acoust., Speech,

Signal Process., vol. 38, no. 12, pp. 2140–2150, Dec. 1990.

[15] P. Stoica and A. Nehorai, “Performance study of conditional and un-

conditional direction-of-arrival estimation,” IEEE Trans. Acoust., Speech,

Signal Process., vol. 38, no. 10, pp. 1783–1795, Oct. 1990.

[16] F. Li, H. Liu, and R. Vaccaro, “Performance analysis for DOA estimation

algorithms: uniﬁcation, simpliﬁcation, and observations,” IEEE Trans.

Aerosp. Electron. Syst., vol. 29, no. 4, pp. 1170–1184, Oct. 1993.

[17] R. Roy and T. Kailath, “ESPRIT-estimation of signal parameters via ro-

tational invariance techniques,” IEEE Trans. Acoust., Speech, Signal Pro-

cess., vol. 37, no. 7, pp. 984–995, Jul. 1989.

[18] A. Gorokhov, Y. Abramovich, and J. F. Bohme, “Uniﬁed analysis of DOA

estimation algorithms for covariance matrix transforms,” Signal Process.,

vol. 55, no. 1, pp. 107–115, Nov. 1996.

[19] Y. Abramovich, N. Spencer, and A. Gorokhov, “Detection-estimation of

more uncorrelated Gaussian sources than sensors in nonuniform linear

antenna arrays .I. Fully augmentable arrays,” IEEE Trans. Signal Process.,

vol. 49, no. 5, pp. 959–971, May 2001.

The condition when Q



is singular is difﬁcult to obtain analytically. In

numerical simulations, we have veriﬁed that it remains nonsingular for various

parameter settings.

[20] C.-L. Liu and P. Vaidyanathan, “Remarks on the spatial smoothing

step in coarray MUSIC,” IEEE Signal Process. Lett., vol. 22, no. 9,

pp. 1438–1442, Sep. 2015.

[21] A. Koochakzadeh and P. Pal, “Cram

er–Rao bounds for underdeter-

mined source localization,” IEEE Signal Process. Lett., vol. 23, no. 7,

pp. 919–923, Jul. 2016.

[22] C.-L. Liu and P. Vaidyanathan, “Cram

er–Rao bounds for coprime and

other sparse arrays, which ﬁnd more sources than sensors,” Digital

Signal Process., 2016. [Online]. Available: http://www.sciencedirect.

com/science/article/pii/S1051200416300264

[23] G. W. Stewart, “Error and perturbation bounds for subspaces associated

with certain eigenvalue problems,” SIAM Rev., vol. 15, no. 4, pp. 727–764,

Oct. 1973.

[24] M. Hawkes, A. Nehorai, and P. Stoica, “Performance breakdown of

subspace-based methods: prediction and cure,” in Proc. IEEE Int. Conf.

Acoustics, Speech, Signal Process., 2001, vol. 6, pp. 4005–4008.

[25] A. Swindlehurst and T. Kailath, “A performance analysis of subspace-

based methods in the presence of model errors. I. The MUSIC algorithm,”

IEEE Trans. Signal Process., vol. 40, no. 7, pp. 1758–1774, Jul. 1992.

[26] S. M. Kay, Fundamentals of Statistical Signal Processing. Englewood

Cliffs, NJ, USA: Prentice-Hall, 1993.

[27] M. Ishiguro, “Minimum redundancy linear arrays for a large number of

antennas,” Radio Sci., vol. 15, no. 6, pp. 1163–1170, Nov. 1980.

[28] Z. Liu and A. Nehorai, “Statistical angular resolution limit for point

sources,” IEEE Trans. Signal Process., vol. 55, no. 11, pp. 5521–5527,

Nov. 2007.

[29] P. P. Vaidyanathan and P. Pal, “Direct-MUSIC on sparse arrays,” in Proc.

2012 Int. Conf. Signal Process. Commun., Jul. 2012, pp. 1–5.

Mianzhi Wang (S’15) received the B.Sc. degree

in electronic engineering from Fudan University,

Shanghai, China, in 2013. He is currently working

toward the Ph.D. degree in the Preston M. Green

Department of Electrical and Systems Engineering,

Washington University in St. Louis, St. Louis, MO,

USA under the guidance of Dr. A. Nehorai. His re-

search interests include the areas of statistical signal

processing for sensor arrays, optimization, and ma-

chine learning.

Arye Nehorai (S’80–M’83–SM’90–F’94–LF’17)

received the B.Sc. and M.Sc. degrees from the Tech-

nion, Haifa, Israel, and the Ph.D. degree from Stan-

ford University, Stanford, CA, USA. He is currently

the Eugene and Martha Lohman Professor of electri-

cal engineering in the Preston M. Green Department

of Electrical and Systems Engineering, Washington

University in St. Louis (WUSTL), St. Louis, MO,

USA. He served as the Chair of this department from

2006 to 2016. Under his Department Chair Lead-

ership, the undergraduate enrollment has more than

tripled in four years and the master’s enrollment grew seven-fold in the same

time period. He is also a Professor in the Department of Biomedical Engineering

(by courtesy), a Professor in the Division of Biology and Biomedical Studies,

and the Director of the Center for Sensor Signal and Information Processing

at WUSTL. Prior to serving at WUSTL, he was a Faculty Member at the Yale

University and the University of Illinois at Chicago.

Dr. Nehorai served as an Editor-in-Chief of the IEEE T

RANSACTIONS ON

SIGNAL PROCESSING from 2000 to 2002. From 2003 to 2005, he was the Vice

President (Publications) of the IEEE Signal Processing Society, the Chair of the

Publications Board, and a Member of the Executive Committee of this Society.

He was the Founding Editor of the special columns on Leadership Reﬂections in

IEEE S

IGNAL PROCESSING MAGAZINE from 2003 to 2006. He received the 2006

IEEE SPS Technical Achievement Award and the 2010 IEEE SPS Meritorious

Service Award. He was elected Distinguished Lecturer of the IEEE SPS for a

term lasting from 2004 to 2005. He received several Best Paper Awards in the

IEEE journals and conferences. In 2001, he was named University Scholar of the

University of Illinois. He was the Principal Investigator of the Multidisciplinary

University Research Initiative project titled Adaptive Waveform Diversity for

Full Spectral Dominance from 2005 to 2010. He has been a Fellow of the Royal

Statistical Society since 1996 and Fellow of the AAAS since 2012.