Covariance Matching Estimation Techniques for Array Signal Processing Applications

DIGITAL SIGNAL PROCESSING ARTICLE NO.

8, 185–210 (1998)

SP980316

Covariance Matching Estimation Techniques for Array Signal Processing Applications B. Ottersten,*,† P. Stoica,‡ and R. Roy§ *Signal Processing/S3, Royal Institute of Technology, SE-100 44 Stockholm, Sweden; ‡Systems and Control, Uppsala University, Box 27, SE-751 03 Uppsala, Sweden; and §ArrayComm, Inc., 3141 Zanker Road, San Jose, California 95134

1. INTRODUCTION Ottersten, B., Stoica, P., and Roy, R., Covariance Matching Estimation Techniques for Array Signal Processing Applications, Digital Signal Processing 8 (1998), 185–210.

Signal parameter estimation from sensor array data is a problem that frequently arises in many applications. Examples include target localization with phased array radars, detection and classification of sources with underwater arrays, and, more recently, interference rejection and signal separation in wireless communications using antenna arrays. A large class of estimation and detection methods has been developed; most of these methods are based on a model of the sensor array response. This model is either derived analytically or obtained from calibration measurements. Thus, the sensor array receiver outputs are assumed to be a known function of the signal parameters which may include azimuth, elevation, polarization, etc. A large number of algorithms have been devised for the sensor array problem. The so-called subspacebased algorithms such as MUSIC, ESPRIT, MODE, and WSF [1–5], and variations thereof [6–10], are based on a vector space formulation of the signal parameter estimation problem. Maximum likelihood estimators, though often leading to comparatively difficult optimization problems, have also been proposed [11–14]. Herein, we will present a general framework for a novel class of covariance matching estimation techniques (COMET) which have gained some interest recently [15–24]. This approach may be traced back to the early 1970s in the statistical literature (see [25] and references therein) but it has only recently received attention in the signal processing community. This class of estimators may be applied to a wider range of signal processing problems than the subspace methods. The covariance matching approach is well suited to problems where the param-

A class of covariance matching estimation techniques (COMET) has recently attracted interest in the signal processing community. These techniques have their roots in the statistical literature where they are sometimes referred to as generalized least squares methods. Covariance matching is an alternative to maximum likelihood estimation, providing the same large sample properties often at a lower computational cost. Herein, we present a general framework for covariance matching techniques and show that they are well suited to solve several problems arising in array signal processing. A straightforward derivation of the COMET criterion from first principles is presented, which also establishes the large sample properties of the estimator. Closed form compact expressions for the asymptotic covariance of the estimates of the parameters of interest are also derived. Some detection schemes are reviewed and two COMET-based detection schemes are proposed. The main part of the paper treats three applications where the COMET approach proves interesting. First, we consider the localization of underwater sources using a hydro-acoustic array. The background noise is often spatially correlated in such an application and this must be taken into account in the estimation procedure. Second, the problem of channel estimation in wireless communications is treated. In digital communications, an estimate of the channel is often required to perform accurate demodulation as well as spatially selective transmission. Finally, a radar detection problem is formulated and the proposed detection schemes are evaluated. r 1998 Academic Press

† To whom correspondence should be addressed. E-mail: [email protected].

1051-2004/98 $25.00 Copyright r 1998 by Academic Press All rights of reproduction in any form reserved.

185

eter set can be separated into linear and non-linear parameters: the array covariance matrix depends linearly on the linear parameters and non-linearly on the non-linear parameters. This allows for a separation of the estimation problem where the linear parameters may be solved analytically in terms of the non-linear parameters, resulting in a lower dimensional optimization problem. The aforementioned type of covariance matrix structure arises, for example, in situations with linearly parameterized noise fields [23,26,27], distributed source signals [19,28,29], and uncorrelated source signals [30]. A straightforward derivation of the covariance matching approach is presented herein based on the extended invariance principle (EXIP) [31]. This derivation also establishes the asymptotic efficiency of the COMET method without any additional calculations. Furthermore, a compact expression for the Crame´r–Rao lower bound on the estimation error variance is derived from the COMET cost function. The problem of determining the correct number of parameters required to model the measured signal is also addressed. By deriving the limiting distribution of the minimum value of the COMET cost function, a hypothesis testing scheme is formulated for signal detection. Three examples are presented where the COMET estimation and detection method may be applied. The next paragraphs briefly describe these applications. We consider the problem of underwater source localization in unknown noise fields. The noise field in underwater applications is often spatially colored and unknown. If a linearly parameterized model of the noise covariance matrix is used, the COMET approach is well suited for this problem. A significant reduction in estimation bias and root-mean-squared error can be achieved as is demonstrated on experimental data collected from a full scale hydroacoustic array in the Baltic Sea. In wireless communications, the use of base station antenna arrays allows spatially selective reception and transmission. However, non-line-of-sight propagation and local scattering at the terminals is commonplace resulting in a wave field that can be seen as generated by distributed sources. We introduce a model for the channel distribution and estimate its parameters using COMET on an experimental data set collected in a typical mobile communication environment. We consider a radar application where a target is to be detected and estimated in the presence of strong interference. The performance of the COMET-

based detection schemes and a traditional likelihood ratio test is demonstrated on a simulation example as well as on data collected by an experimental radar system. It is shown that COMET-based detection can provide accurate detection and estimation capabilities in the presence of interference by exploiting the signal covariance structure. In the following section, the data model and some assumptions are introduced. Section 3 presents the COMET framework and Section 4 provides a straightforward statistical analysis of COMET based on first principles. In Section 5 the detection problem is treated. The emphasis of the paper is on the three applications presented in Section 6.

2. DATA MODEL Consider the signal model y(t) 5 Ax(t) 1 e(t),

(1)

where the complex observation vector, y(t), is of dimension m 3 1, Ax(t) models the signals of interest, e(t) is additive measurement noise, and t is a time index. In sensor array signal processing, the columns of the m 3 nu matrix A are often referred to as array response vectors and depend on, for example, the location parameters of the emitters. The nu emitter signals received at the array are collected in the vector x(t). The model in (1) requires a narrowband emitter signal. The propagation delays between sensors in the array may be modeled as phase shifts provided the bandwidth of the emitter signal is sufficiently small [32]. However, for wideband signals, an equation similar to (1) may be formed in the frequency domain [11]. It is normally assumed that the array response matrix, A, is deterministic and known as a function of the unknown signal parameters. Here, a more general model will be considered allowing for stochastic response vectors with a parameterized distribution function. The reason for this will become transparent when considering the application in Section 6.3. Let the observation vectors, 5y(t)6t51,2,..., be independent identically distributed circular Gaussian random variables with zero mean. The measurement noise is assumed to be uncorrelated with the signals x(t), which leads to the following model for the covariance matrix of 5y(t)6, R(u, m, s) 5 E5y(t)y*(t)6 5 Rs (u, m) 1 Q(s),

186

(2)

where E5·6 denotes the expected value, (·)* denotes the conjugate transpose, and

has been proposed [16,23,27,34] K

Rs(u, m) 5 E5Ax(t)x*(t)A*6

(3)

Q(s) 5 E5e(t)e*(t)6.

(4)

Q5

m

(6)

k

o s vec (Q ) k

k

k50

s0 · 5 [vec (Q0 ) · · · vec (QK )] · 5 Ss, · sK

34

(7)

where sk are unknown parameters and Qk are known Hermitian weighting matrices. This will guarantee that Q is Hermitian but in general positive (semi-) definiteness is not ensured. By choosing K 5 0, Q0 5 I, and s0 5 s2 we obtain the special case where the noise is spatially white. The more general parameterization in (6) allows for incorporation of a priori knowledge on the background noise and interference into the model. This model is used for direction estimation in Section 6.2 where an example of how to choose 5Qk 6 is also given. j In order to ensure that the parameterization used is ‘‘uniquely identifiable,’’ some restrictions are necessary. Let the dimensions of the parameter vectors u, m, and s be nu, nµ, and ns respectively. Note that we assume the case of one parameter per emitter signal, hence

r(u, m, s)5 vec (R(u, m, s)) 5 C(u)m1 Ss

3s4 0 F(u)a,

k

K

vec (Q) 5

With the assumption of independent vector samples, the maximum likelihood estimator is easily formulated; see Section 3. The estimation methods that are developed in the following may of course be applied successfully when the vector samples are temporally correlated. Also, the temporal correlation may be exploited by considering additional covariance matrix lags and not only the covariance matrix in (4) as will be the case here. It will further be assumed that the signal and noise covariance matrices are linearly parameterized by the real-valued vectors m and s, respectively, which is often the case in practical applications (as will be illustrated). This assumption is crucial for the development that follows. The elements of m are referred to as the emitter signal covariance parameters and those of s are the noise covariance parameters. In general, the signal parameters, u, enter in a highly non-linear fashion in the covariance matrix in (3). Under the assumptions previously made, (2) can be written as

5 [C(u) S]

osQ

k50

nu 5 dim [x(t)] 5 dim (u).

(5)

(8)

The extension to the case where dim [x(t)] Þ dim (u) is straightforward. Later it will be shown how closedform estimates of the linear parameters m and s can be obtained, hence reducing the problem to the estimation of u. Determining nu is termed the detection problem and is addressed in Section 5. To uniquely determine u, one obvious condition is that

where vec (·) is a vector obtained by stacking the columns of the argument on top of each other, a is obtained by concatenating the parameter vectors m and s, S is a known matrix, and C(u) is a given function of the unknown signal parameter vector u. An example of this parameterization is given below for the case of an unknown spatial noise covariance.

C(u8) 5 C(u9)

EXAMPLE 2.1 (Array processing with linear noise models). Consider an array processing situation where A(u) represents a deterministic array response matrix and x(t) contains the nu , m emitter signals with covariance matrix Rx 5 E5x(t)x*(t)6. The Hermitian emitter covariance matrix, Rx, may be linearly parameterized by the matrix entries (without imposing the condition of positive semidefinitedness). To ensure identifiability from the second-order statistics, a structured noise covariance matrix model must be introduced [33]. The following linear model

implies u8 5 u9.

(9)

This condition guarantees an unambiguous parameterization of C(u). To guarantee the uniqueness of a, we require that F(u) is full rank at the true values u. Another transparent condition for identifiability is that the number of unknowns is not larger than the number of estimating equations available, i.e., nu 1 n µ 1 ns # m2,

(10)

where m 2 equals the number of distinct real-valued

187

entries of R(u, m, s). This condition is necessary but it is not sufficient to ensure identifiability; see the example below. Sufficient conditions for identifiability will be application specific since they will depend on the structure of parameterization. The following example also illustrates that the parameterization of this paper may have weaker identifiability requirements than the commonly used (unstructured) parameterizations in most of the previous literature.

(ULA) with seven elements, i.e., the same structure as in Fig. 1 but without ‘‘missing’’ sensors, we would also be able to identify at most six emitters: even with uncorrelated emitters, a seven element ULA can only resolve six emitters owing to the invariance structure of the ULA. Indeed since the emitter signals are uncorrelated, the ULA covariance matrix is Toeplitz and there are only seven (m) covariance lags that carry independent information; furthermore the number of equations (Toeplitz covariance defined by 2m 2 1 real entries) must be larger than the number of unknowns (nu directions, nu powers, and the noise variance, have 2nu 1 1 unknowns) which results in m . nu. The sparse array structure chosen in the example is a so-called minimum redundancy array which in effect allows the estimation of all seven covariance lags with only four sensors [40]. The previous discussion also shows that the condition (10), while necessary for identifiability, is not sufficient. Indeed, with nu 5 7 emitters we would have 15 unknown parameters, which is less than m2 5 16. Yet as explained above, the array in Fig. 1 can j uniquely estimate the directions of only nu 5 6.

EXAMPLE 2.2 (Locating more sources with fewer sensors.) Most array signal processing methods require that the number of incoming wavefronts, nu, be less than the number of sensors, m [32,35]. This is the identifiability condition when the correlation among the wavefronts is completely unknown, i.e., E5x(t)x*(t)6 is unstructured. In some situations, the incoming emitter signals are known to be uncorrelated [30,36,37]. It is not possible to exploit this information in the traditional array processing techniques, unless they are modified significantly [38]. However, it is trivial to use this information in the current model by parameterizing E5x(t)x*(t)6 as a diagonal matrix. This also enables us to identify a number of emitter parameters (such as directions of arrival) that is larger than the number of sensors [39]. A necessary condition is of course that the number of equations be greater than or equal to the number of unknowns, i.e., m 2 $ 2nu 1 ns; see [9] and also [21]. For example, consider the four element sparse array structure in Fig. 1 with nominal half-wavelength element spacing. A number of uncorrelated emitter signals are received at the array from directions equally spaced between 0° and 180°. For example, with four signals the directions of arrival are 36°, 72°, 108°, and 144°. The signal-to-noise ratio for each signal is 20 dB and the noise is assumed to be spatially white. In Fig. 2 the average standard deviation of the direction estimates is displayed for increasing number of emitters assuming 100 vector samples. The average standard deviation is obtained as the trace of the Crame´r–Rao lower bound matrix derived in Section 4, divided by the number of emitters. Note that we can accurately estimate the directions to six emitters with only four sensors. With seven emitters, the problem becomes unidentifiable. To see this, note that if we had a uniform linear array

3. PARAMETER ESTIMATION VIA COVARIANCE MATCHING Based on the model presented in the previous section, we can formulate a maximum likelihood (ML) estimator for the unknown parameters. Under regularity conditions, the ML estimator is expected to give optimal estimates in the sense that the Crame´r–Rao bound (CRB) is achieved asymptotically (in the number of data samples). Under the data model introduced in the previous section, the ML parameter estimates are given by [7,25] ˆ ML, s ˆ ML ) 5 arg min (log 0R(u, m, s)0 (uˆ ML, m u,m,s

ˆ 6), 1 Tr 5R21 (u, m, s)R

(11)

where 0·0 denotes the determinant, Tr 5·6 denotes ˆ is the sample covariance matrix the trace, R ˆ 5 R

FIG. 1. A four element sparse linear array. A circle indicates a ‘‘missing’’ sensor.

1

N

o y(t)y*(t),

N t51

(12)

and N is the number of available snapshots. When

188

FIG. 2. The average standard deviation of the direction estimates for increasing number of emitters uniformly space between 0° and 180°. The emitters are equi-powered 20 dB stronger than the noise. A four element minimum redundancy array is used and 100 snapshots are assumed.

the covariance matrix is modeled as R 5 A(u)RxA*(u) 1 s2 I with Rx being a general Hermitian matrix, explicit expressions for the ML estimates of Rx and s2 can be obtained [41–43]. This allows a concentration of the likelihood function in (11) with respect to these parameters and the resulting function depends only on u. In general, however, it is not possible to concentrate m and s out from the negative log-likelihood function in (11) which must therefore be minimized with respect to all nu 1 nµ 1 ns parameters. This is not an attractive optimization problem and alternative approaches are of great interest.

belong to a compact domain. This condition rules out any model with rank constraints on its matrices. A most natural choice of a model set, which includes the original set and in which the MLE is consistent and readily derived, corresponds to the parameterization of R via its elements (subject to the condition of Hermitian symmetry, only). The ML estimate of the parameters in this ‘‘model’’ of R is ˆ in well known to be the sample covariance matrix, R (12) [25]. In the second step, the parameter estimates from the first step are viewed as input data from which Markov-like estimates of the original parameters are obtained. Let

3.1. Covariance Matching Estimation Technique ˆ ). rˆ 5 vec (R

The extended invariance principle [31] has been introduced with the purpose of simplifying complicated ML problems by means of large-sample approximations. According to EXIP, first we should determine the exact ML estimator of the unknown parameters in a model of R which is less detailed than the original one (2) (and hence more easily handled) and which is such that the ML estimate of its parameters is consistent. Note that, additionally, the parameters of the less detailed model should

(13)

Note that rˆ has several related entries due to the ˆ . Also, the off-diagonal enHermitian property of R ˆ are complex-valued. Form the real-valued tries of R vector gˆ which contains the real parts of entries on ˆ and the imaginary and below the main diagonal of R ˆ. parts of the entries below the main diagonal of R The relationship between the real-valued vector gˆ and the complex-valued vector rˆ can be described by

189

a simple linear transformation ˆ ) 5 Jrˆ, gˆ 5 J vec (R

E5rˆirˆ*j 6 5 (14)

N

N

o o E5y(t)y*(t)y*(s)y (s)6

N 2 t51

i

m2

ˆ 1R ˆ ) ˆ ij ) 5 1(R real (R ij ji 2 ˆ ij ) 5 imag (R

1 2j

ˆ ij 2 R ˆ ji ). (R

C 0 E5(rˆ 2 r)(rˆ 2 r)*6 5

(15)

(16)

(17)

(RT # R),

(21)

5 r˜*(R2T # R21 )r˜,

(22)

(23)

where we have used the fact that J is a nonsingular matrix. The cost function can be further simplified by using properties of the vec and # operators [44]:

(18)

˜) r˜*(R2T # R21 )r˜ 5 r˜(R2T # R21 ) vec (R ˜ R21 ) 5 r˜* vec (R21R

(24)

˜ R21 ) 5 Tr 5R ˜ R21R ˜ R21 6. ˜ )* vec (R21R 5 vec (R

(25)

This form of the cost function (for the real-valued ˆ is a case) was used for example in [25]. Since R consistent estimate of R, a large-sample ML estimator of u and a is obtained, according to EXIP, by the minimization of the function

As a first step toward deriving an expression for C, we determine the second-order moments of the covariance estimates. The ith (m 3 1) subvector of the (m 2 3 1) vector rˆ is given by

ˆ 21 )r˜(u, a) ˆ 2T # R r˜*(u, a)(R

N

No

N

N 21g˜ TC21g˜ 5 N 21r˜*J*(JCJ*)21Jr˜ 5 N21r˜* C21r˜

3.2. COMET: Detailed Formulas

y(t)yi (t).

1

A consistent estimate of C can be formed from the data and used in (18) to obtain a large-sample ML estimator of u and a. Note that the cost function in (18) can be rewritten in an alternative form. Consider the normalized cost function (where we have omitted the dependence on u, a for notational convenience)

where the weighting matrix, C21, is the inverse of the (asymptotic) covariance matrix of the residuals, g˜ 5 gˆ 2 g. In order to use EXIP an expression for the matrix C is required. Below, a straightforward derivation is given.

1

(20)

C 0 E5(gˆ 2 g)(gˆ 2 g)*6 5 JCJ*.

(gˆ 2 g(u, a)) TC21 (gˆ 2 g(u, a))

rˆi 5

RjiR,

N

where # denotes the Kronecker matrix product. The covariance matrix of gˆ 2 g is obtained easily by using the relationship in (17):

The Markov-like estimates of u and a are obtained by fitting the ‘‘data,’’ gˆ , to the model (17) in a weighted least squares sense

0 g˜ T (u, a)C21g˜ (u, a),

1

which results in

Since rˆ is constructed from a Hermitian matrix, it is clear that the mapping between rˆ and gˆ is one-toone. Given rˆ, gˆ is uniquely defined and vice versa. Thus the transformation matrix J is invertible. The model of the covariance matrix in (2) can now be expressed in terms of g as g(u, a) 5 J vec (R(u, a)) 5 Jr(u, a).

i

s51

5 rir*j 1

where J is an 3 matrix whose entries are 0, 1, 1 or 6j/2. The real-valued entries of g ˆ are obtained ⁄2 ˆ. by forming linear combinations of the elements in R ˆ are real and for the The diagonal elements of R off-diagonal elements we have m2

1

ˆ 21r˜(u, a), 5 r˜*(u, a)W

(19)

t51

(26)

ˆ. ˆ 5R ˆT#R where W The cost function in (25) or (26) is sometimes referred to as the generalized least squares criterion [25,45]. This cost function is known to be a large-

Since the signals are independent from snapshot to snapshot and circularly Gaussian distributed, the following relation is easily established

190

sample approximation to the maximum likelihood criterion [16,19]. In the previous literature, the asymptotic distribution of the estimates minimizing (26) has been derived by means of a rather involved analysis and then compared to that of the maximum likelihood estimates. The derivation above establishes the large-sample equivalence of the minimizer of (26) and the maximum likelihood estimate in a straightforward manner based on first principles. Also, the algebraic equivalence between the cost function explicitly using real-valued quantities (18) and the cost function in (23) is another useful byproduct of the above derivation. Since r(u, a) is a linear function of a (for given u)—see (5)—and since F(u) must have full column rank (by the identifiability assumptions guaranteeing the unique determination of (u, a) from given R), the minimization of (26) with respect to a is immediate,

Substituting (28) into (27) leads to the concentrated function ˆ 21/2 P'ˆ 21/2 ˆ 21/2rˆ, W uˆ 5 arg min rˆ*W W F(u) u

(30)

where '

PWˆ 21/2F(u) 5 I 2 PWˆ 21/2F(u) ˆ 21/2F(u)[F*(u)W ˆ 21F(u)] 21F*(u)W ˆ 21/2. 5I 2 W

(31)

Note that (30) can be further simplified by observing that

ˆ 21 (rˆ 2 r(u, a)) min (rˆ 2 r(u, a))*W

ˆ 2T/2 # R ˆ 21/2 ) vec (R ˆ) ˆ 21/2rˆ 5 (R W

(32)

ˆR ˆ 21/2 ) 5 vec (I), ˆ 21/2R 5 vec (R

(33)

which results in the following cost function whose minimizer yields the COMET estimate of the signal parameter vector,

a

ˆ 21/2rˆ 2 W ˆ 21/2F(u)a\2 (27) 5 min \W a

ˆ 21rˆ, ˆ 21F(u)]21F*(u)W ⇒a ˆ 5 [F*(u)W

uˆ 5 arg max vec* (I) PWˆ 21/2F(u) vec (I) u

(28)

ˆ 21 )F(u) 5 arg max vec* (R

ˆ 21/2 is where \·\ denotes the Euclidean vector norm, W 21 ˆ a Hermitian square root factor of W , and u is yet to be determined. Note that in general the parameter vector a is real-valued. Optimality of the parameter estimate obtained from (26) can only be ensured if the estimate a ˆ in (28) is also real-valued. This fact is not obvious at first since the quantities in (28) are complex-valued, and in effect was not recognized in previous works which ended up with estimation procedures more complicated than necessary. Note, however, that since J is invertible, (28) can be rewritten as ˆ J*) 21Jrˆ, ˆ J*)21JF] 21(JF)*(JW aˆ 5 [JF)*(JW

(34)

u

ˆ 21 ), ˆ 21F(u)] 21F*(u) vec (R 3 [F*(u)W

(35)

where in the last equation we have used the fact that ˆ 21/2vec (I) 5 (R ˆ 2T/2 # R ˆ 21/2 ) vec (I) W ˆ 21 ). 5 vec (R

(36)

ˆ must be computed but not the Thus, the inverse of R square-root factor. Once uˆ is obtained as above, large-sample ML estimates of m and s can be obtained with (28) (estimates of m and s may be of interest in some applications).

(29)

where Jrˆ 5 gˆ is real and so is JF (since Jr 5 (JF)a must be real for any real-valued a). It remains to ˆ )J*. To this end, let R ˆ be ˆ J* 5 J(R ˆT#R examine JW fixed and let 5z(t)6 be a sequence of independent random variables that are circularly Gaussian disˆ. tributed with mean zero and covariance matrix R N Also, let rˆ 5 vec (St51 z(t)z*(t)/N). The previous calculations, leading to (14) and (23), imply that Jrˆ is real-valued, and also that the covariance matrix of ˆ )J*. Hence the latter matrix ˆT#R this vector is J(R must be real.

4. STATISTICAL ANALYSIS AND CRB Derivation of the limiting distribution of (uˆ , aˆ ) is straightforward. Since u is usually the parameter vector of most interest, we focus on the analysis of uˆ in what follows. The compact form of the COMET cost function in (35) is used here to derive a compact closed-form expression for the covariance of uˆ .

191

The consistency of uˆ is easily established. As N = `, uˆ converges to the set of global minima of the asymptotic function corresponding to (35), which can be written as

' 21/2 PW F(u0 )a0 21/2F(u) W

5

00

00

where a(u) represents a deterministic array response vector. In this case (3) becomes Rs (u, m) 5 A(u)Rx (m)A*(u).

The derivative of vec (R) in (39) is then easily evaluated. Consider the derivative with respect to the pth element of u,

2

,

(37)

r

­ vec (R) ­up

where (u0, a0 ) denote the true values of (u, a). To prove the consistency of uˆ we need to show that the set of global minimizers of (37) contains a single point: u0. The set in question is seen to be defined by the solutions to the equation F(u)a 5 F(u0 )a0 5 r

CRB

5N

1

­ vec (R) * ­u

2

21/2

W

(38)

­ vec (R) ­u

­A

1­u R A*2 1 vec 1AR x

­A* x

p

5 (AcRcx # I)

­up

2

­ vec (A) ­up

1 (I # ARx )

­ vec (A*) ­up

­ vec (A) ­up

1 (I # ARx )Lmn

­ vec (Ac ) ­up

,

(41)

where (·) c denotes the complex conjugate, and the permutation matrix Lmn is such that for any m 3 n matrix A, vec (AT ) 5 L mn vec (A).

(42)

The matrix expression corresponding to (41) can be written as

' PW 21/2F(u)

3 W21/2

5 vec

5 (AcRcx # I)

for some vector a. Under the assumption of parameter identifiability discussed in Section 2, the only solution to (38) is u0 (and a0 ), which establishes the consistency of the estimates. By the properties of (consistent) ML estimators and the EXIP, the asymptotic distribution of uˆ is Gaussian with covariance matrix that equals the CRB. In Appendix A the following matrix expression for the CRB on the signal parameters is derived 21

(40)

,

­ vec (R)

(39)

­u

5 (AcRcx # I)D 1 (I # ARx )LmnDc,

(43)

where the (mnu ) 3 nu matrix D is given by D 5 ­ vec (A)/­u. An explicit expression for F(u) will depend on the parameterization of Rx (m); however, in general we have

which provides a useful means of assessing the accuracy of the obtained estimates. Remark 1. In the COMET approach to estimating (u, m, and s) there is no guarantee that the ˆ (derived from uˆ , m ˆ s and Q ˆ , and s ˆ ) are estimated R positive semi-definite, for a given sample length, because no constraint has been imposed on m and s. Now, since the COMET estimate is consistent (as ˆ will be positive definite for a ˆ s and Q shown before), R ‘‘sufficiently large’’ N. Hence, COMET is also a large-sample realization of the ML method which constrains m and s to lie in the set 5m, s0Rs and Q . 06.

F(u)a 5 vec (A(u)Rx (m)A*(u)) 1 Ss

(44)

5 (Ac (u) # A(u))vec (Rx (m)) 1 Ss

(45)

5 [(Ac (u) #A(u))K S]

m

3s4 ,

(46)

where K is a fixed matrix that depends on the parameterization of Rx (m). Note that a priori information such as uncorrelatedness between sources is easily incorporated (see Example 2.2) in this parameterization.

We conclude this section by considering a wellknown special case of our signal model. EXAMPLE 4.1 (Point source model in array processing). As in Example 2.1, let A(u) 5 [a(u1 ) · · · a(unu )]

192

Since COMET is asymptotically equivalent to ML, one can approximate (49) by

A compact matrix expression for the CRB (39) is now easily obtained. When R 5 ARxA* 1 s2I with a general Hermitian signal covariance Rx, a compact expression for the CRB is available [13]. Showing that (39) reduces to this expression under the previous special parameterization is an interesting exercise. j

l(uˆ , aˆ ) where uˆ , a ˆ are the COMET estimates.

(50)

In order to show that the above approximation is suitable, consider the following simple calculation. ˆ Let bˆ T 5 [uˆ T aˆ T ] be the COMET estimate and b be

5. DETECTING THE NUMBER OF SIGNALS

the ML estimate. Then In order to apply the estimation techniques described previously, the number of parameters, also called the ‘‘model order,’’ must first be determined. This is often called the detection problem. The number of free parameters of the model under discussion is nu 1 nµ 1 ns.

ˆˆ )(b ˆˆ )1 ­ ˜l(b ˆ 2 bˆˆˆ ) ˜l (b ˆ ) 5 ˜l (b b ˆˆ )*­ ˜l (bˆˆ )(b ˆˆ ) 1 · · · ˆ 2b ˆ 2b 1 21(b bb ˆ 5 l˜(bˆ ) 1 O(1/N 2 ),

(47) ˆ where we have made use of the fact that ­b˜l(bˆ ) 5 0 and bˆ 2 bˆˆ 5 O(1/N ). The notation O(·) is an in probability version of the conventional notation [46]. ˆ ) provides an approximaIt follows from (51) that ˜l (b ˆˆ ). Apparently, even an tion of order O(1/N 2 ) for ˜l (b

Often, we wish to determine the number of signal parameters, nu. When nµ and ns are a priori specified (or are a function of nu ) we say that we have a signal detection problem. The parameter nu is then related to the number of signals present in the observation process. However, one can also pose the problem of determining an appropriate structure for Rs and Q (i.e., estimating nµ and ns, provided the structures tested for these matrices are nested). Below, we show how the COMET parameter estimates can be used in a standard detection scheme for model order selection. Two detection methods based directly on the minimum of the COMET cost function will also be derived.

approximation of order O(1/N 3/2) would be sufficient ˆ to preserve the asymptotic distribution of ˜l (bˆˆ ). A test is formulated for the hypothesis that nu 1 nµ 1 ns parameters are adequate to describe the observed data against the alternative that R has an arbitrary structure. Thus, the hypotheses are that y(t) is N(0, R), i.e., Gaussian distributed with zeromean and covariance matrix R, with

5.1. Likelihood Ratio Test

H0: R 5 Rs (u, m) 1 Q(s),

Since the ML problem under discussion is a standard one, the detection question can be solved by generalized likelihood ratio test (GLRT) rules. The negative log-likelihood function of the observed data sequence is given by (within an additive constant; see (11))

3

4

ˆ6 . log 0R(u, a)0 1 Tr 5R21(u, a)R

H1: R is arbitrary. Under H1, the ML estimate of R is simply given by ˆ [25]. Under H0, the the sample covariance matrix, R COMET estimates are used in (48). The resulting GLRT (the difference between the minima of the negative log-likelihood functions (48) under H1 and H0 ) is given by

(48)

5

l(u, a) , N

(51)

˜l(u, a)

ˆ 02 Nm. GLRT 5 l(uˆ , aˆ ) 2 N log 0R

Note that the structured covariance matrix estimate under H0 requires nu 1 nµ 1 ns parameters and that the unstructured covariance matrix estimate under H1 requires m 2 parameters. Standard likelihood

For GLRT the quantity of interest is min l(u, a). u,a

(52)

(49)

193

parameterized models. The test statistic

ratio test theory then states that 2 3 GLRT , Asx2 (m 2 2 nu 2 nµ 2 ns ),

COMET-AIC (nu 1 nµ 1 ns )

(53)

ˆ 21/2 P'ˆ 21/2 W ˆ 21/2rˆ1 2(nu 1 nµ 1 ns) (55) 5 min r*W W F(u) u

i.e., the test statistic, 2GLRT, is asymptotically chi-squared distributed with m 2 2 nu 2 nµ 2 ns degrees of freedom under H0. The GLRT detection scheme sequentially tests different null hypotheses against H1 until one is accepted. In each case, the COMET estimates have to be computed and the test statistic is compared against a threshold. The threshold is chosen according to the tail area of the chi-square distribution (with the appropriate number of degrees of freedom) to get a certain probability of correct detection (typically a value slightly larger than 0.9). If the threshold is exceeded, a new hypothesis is chosen and the procedure repeated.

is evaluated for increasing model orders. This involves computing the COMET estimates for each case. The model order that minimizes COMET-AIC is taken to be the estimate of the number of parameters. Unlike the chi-squared tests of GLRT and COMET, COMET-AIC based on (55) requires parameter estimation for an over-fitted model which is a disadvantage. However, a threshold does not have to be chosen which can be a distinct advantage; see Section 6.4.

5.2. COMET Based Detection

6. ARRAY SIGNAL PROCESSING APPLICATIONS

An alternative to GLRT is to derive the (asymptotic) distribution function for the minimum value of the COMET criterion (35) and base the detection test on this distribution. In Appendix B we show that when the correct number of parameters is assumed, i.e., the parameters nu, nµ, and ns match the number of parameters in R,

The COMET estimate is in general obtained by minimizing a multidimensional non-linear cost function. Before presenting the array signal processing applications, some implementational aspects of COMET are presented in the following.

ˆ 21/2 P'ˆ 21/2 ˆ 21/2rˆ W min rˆ*W W F(u)

6.1. Calculating the COMET Estimate

u

, Asx2 (m 2 2 nu 2 nµ 2 ns ).

Since the COMET criterion is a continuous function of the parameters, a Newton-type method is a good choice for performing the minimization in (35), at least locally [48,49]. In this method the parameters are updated according to

(54)

This result can be used to estimate (nu, nµ, n s ) or only nu (signal detection) when nµ and ns are a priori specified. The minimum value of the COMET cost function is evaluated for increasing model orders. For each model order, the minimum is tested against a threshold that is chosen according to the tail area of the chi-square distribution with the appropriate number of degrees of freedom. If the threshold is exceeded, the model order is increased and the procedure is repeated.

uk11 5 uk 2 µk [V9(uk )] 21V8(uk ),

(56)

where uk is the estimate at the kth iteration, µk is the step length, V 9 is the second derivative of the criterion function, and V 8 is the gradient. A frequent problem with this procedure is that the second derivative matrix (the so-called Hessian) is not always positive definite, especially as we move away from the optimum. An alternative is to use an approximation to the Hessian. It is the natural to consider the asymptotic instead of the exact Hessian, which results in the so-called scoring method. This has several advantages: the asymptotic Hessian is guaranteed to be positive (semi-) definite, it is computationally less costly, and it only involves first-order derivatives of the covariance matrix. The scoring

5.2.1. COMET-AIC Detection A modification of detection schemes which are based on explicit knowledge of the limiting distributions and which does not involve the subjective choice of a threshold is Akaike’s information criteria (AIC) [47]. It can be formulated for COMET cost function in (54) and involves an additional term that increases with the model order to penalize over-

194

Note that only the second term above depends on u1. Thus the problem is reduced to

recursion can be written as uk11 5 uk 2 µkH21(uk )g(uk ).

(57) uˆ 1 5 arg max vec* (I)

The approximation of the Hessian matrix, H, and the gradient vector, g, required by the algorithm are available as by-products of the statistical analysis in Section 4; see (B.4) and (B.5). There are several ways of choosing the step length, µk, in (57), for example, a line search may be performed at every iteration to find the µk that provides the largest decrease in the cost functions. See [48,49] for more details and alternatives. It should be noted that the asymptotic Hessian is proportional to the asymptotic covariance matrix of the estimation errors in (39). Thus, the scoring algorithm provides an estimate of the accuracy of the parameter estimates at no additional cost. Furthermore, when computing the gradient vector, aˆ is readily available; see (28). The scoring technique above converges rapidly to the local optimum in the attraction basin of which it is initialized. One way of initializing a multidimensional search is by relaxing the optimization problem, i.e., the multidimensional problem is broken into a series of one-dimensional search problems. First, u is chosen to be a scalar parameter corresponding to one signal.1 After this parameter is determined, it is fixed and the optimization is performed over one additional u parameter. This is repeated until the correct number of parameters is reached, after which the rapidly converging scoring technique is preferably employed. For the case when the unknown parameter enters in only one column of F(u), the relaxation technique above can be implemented efficiently [50]. Let F 5 [F f1 ] where the column vector f1 depends on the parameter u1 to be optimized over and F contains the fixed parameters. The projection matrix in the COMET cost function (35) may then be decomposed using the projection matrix update formula (see, e.g., [50]):

u1

3

1

ˆ 21/2P'ˆ ˆ 21/2f f*1W W 1 W21/2F

ˆ 21/2P'ˆ 21/2 W ˆ 21/2f1 f*1W W F

vec (I)

(59)

'

5 arg max u1

ˆ 21/2 P ˆ 21/2 vec (I)0 2 0f*1W W F ˆ 21/2 P'ˆ 21/2 W ˆ 21/2f f*1W 1 W F

.

(60)

The matrices above need to be computed only once; thereafter only vector–matrix products are required to evaluate the cost function. In specific applications, by taking the structure of F into account, efficient implementation of the one-dimensional maximization in (60) can often be obtained. Finally, it is our experience that the optimization of the COMET criterion in general requires some degree of regularization of the sample covariˆ should be replaced by R ˆ 1 lI, ance matrix; i.e., R where l is the regularization parameter, especially ˆ exhibits a large eigenvalue spread. The local when R properties of the criterion function (around the optimum) which determine the accuracy of the parameter estimates are mainly determined by the domiˆ . These are only to a small nant eigenelements of R degree affected by the regularization. However, the cost function is better conditioned when applying regularization. By this we mean that the search algorithm is less likely to converge to a local minimum, and its convergence speed is increased. A further discussion on ill-conditioning in generalized least squares can be found in [51].

6.2. Direction Estimation in Underwater Acoustics An array of hydro-acoustic sensors is often employed for source localization in underwater applications, which is a challenging problem primarily due to the difficulty of accurately modeling the propagation medium [52]. The background noise is often spatially correlated along the array due to, for example, surface noise and engine noise. Also, the ambient noise field may give rise to spatially correlated noise whose covariance structure will depend on the array configuration [53]. Additionally, any modeling errors present may be viewed as a con-

PWˆ21/2F(u) 5 PWˆ21/2F ' ˆ 21/2f f* W ˆ 21/2P'ˆ 21/2 PWˆ21/2FW 1 1 W F

' ˆ 21/2f1f*1W ˆ 21/2P'ˆ 21/2 PWˆ21/2FW W F

. (58)

1 Sometimes one signal is parameterized by more than one parameter; see for example the application in Section 6.3

195

tribution to the noise and, in general, they will be spatially correlated as well. If accurate information about the noise covariance is available, the measured data can be pre-whitened leading to the wellstudied spatially white noise problem. However, this is rarely possible due to the non-stationarity of noise. The sensitivity of high-resolution parametric estimation methods to spatially correlated noise has been reported for example in [54,55]. Several approaches have been proposed to address the source localization problem in the presence of an unknown noise covariance [8,56–63]. It is obviously not possible to perform source localization based on second-order statistics in completely unknown noise fields. Therefore, any estimation method requires some prior knowledge on the noise covariance structure. Herein, we will examine a class of noise covariance models which are well suited to underwater applications and the COMET estimator. We reiterate that the COMET technique can be used with any linear parameterization of the spatial noise covariance matrix. The performance of the covariance matching method will be evaluated for different model orders on experimental data collected from a full-scale hydro-acoustic array deployed in the Baltic Sea. As will be shown, a significant reduction in the estimation bias and total root-mean-squared estimation error can be achieved by using an appropriate noise model.

Assume that the spatial power density of the noise is given by p(u) where u is the spatial parameter associated with the angle of incidence. Usually, the received signals are bandpass filtered around the center frequency of interest. Therefore, the background noise can be viewed as a narrowband signal when modeling its propagation across the array. The spatial covariance matrix of the noise is then given by the integral Q5

e

2p

0

a(u)p(u)a*(u) du.

(61)

The sources that we are attempting to localize have an impulsive spatial power density; however, the noise will be assumed to have a smooth spatial power density function. Let the noise spatial power density be parameterized in terms of its Fourier coefficients `

o (u

p(u) 5 u01

k

cos (ku) 1 vk sin (ku)).

(62)

k51

If p(u) is sufficiently smooth, the coefficients of the higher order terms decrease as k increases. Hence by choosing K in (63) sufficiently large the following parameterized model of the spatial noise covariance is obtained K

Q5

osQ, k

k

(63)

k50

6.2.1. Noise Covariance Modeling

where K is even, s2k 5 uk, s2k11 5 vk, and (for k $ 0):

In most array signal processing methods, the spatial noise covariance matrix Q must be known up to a scalar. However, in many cases this is an unrealistic assumption and the noise covariance and signal parameters must be estimated simultaneously. Herein, we will consider a linearly parameterized spatial noise covariance matrix (see (6)) with the unknown parameters s. This parameterization allows for incorporation of prior knowledge on the background noise and interference into the model. For example, for a uniform linear array, it is reasonable to assume that the noise covariance is Toeplitz and the covariances between well-separated sensors are zero. This structure may easily be expressed as in (6). Also, modeling different noise powers in different sensors is easily incorporated in this model. In the experiment presented the set of known matrices, Qi, is obtained by using knowledge of the sensor array response, a(u), and assuming that the spatial power of the noise field is smooth; see also [20,23,26, 34,64].

Q2k 5

e

Q2k11 5

2p

0

e

a(u) cos (ku)a*(u) du

2p

0

a(u) sin (ku)a*(u) du.

(64) (65)

This modeling approach is attractive in that it captures the two extremes in terms of spatial power distribution. If power is received from a very small angular region it is considered to be a signal of interest and is modeled as a point source. The direction to the source will then be estimated. On the other hand, a spatially diffuse source (power received from a large angular region) is modeled as noise.

6.2.2. An Underwater Experiment We next present some results of underwater source localization based on the noise covariance matrix modeling above and the COMET method. The experimental data were collected with seven hydrophones 196

source measured relative to array broadside. It follows from (64)–(65) and (66) that the elements of Qk are given by

uniformly spaced along a horizontal line approximately 0.9 m apart and 30 m under the sea surface. The data were collected in the Baltic Sea under calm summer conditions and in an area with a horizontally stratified sea bottom. The experiment was conducted and the data were made available by the Swedish National Defense Research Establishment. A transmitter was placed at 18 known locations about 1300 meters from the array. Most signal power was concentrated at 285 Hz. The data were collected by coherent receivers and bandpass filtered around 285 Hz, and complex (I/Q) samples were formed by Hilbert transformation. At each location, 10 batches of 50 samples each were collected. We assume that the sensors are omni-directional (in the horizontal plane) and uniformly spaced. If additionally the propagation medium is assumed to be isotropic the array response is given by

a(u) 5

3

1 e

j2pDl sin u

···

4

,

il Q 2k 5

e

il Q 2k11 5

2p

0

e

e j2pDl(i2l) sin u cos (ku) du

2p

0

e j2pD l(i2l) sin u sin (ku) du.

(67) (68)

In Fig. 3, the estimated locations for the 18 different experiments are displayed. The COMET estimator is used assuming spatially white noise with equal power in all sensors (hence K 5 0 in (63)). The sample statistics for each location are based on the estimates obtained from the 10 batches available. The straight line shows the direction to the source as obtained from optical measurements that are quite accurate. The zero degree represents array broadside. Note that the standard deviations of the estimates are often quite small, whereas the biases can be significant. This is most likely due to modeling errors. Any estimator assuming a white noise model will show comparable performance to that displayed in Fig. 3. Figure 4 shows the same case as above but now seven real noise covariance matrix parameters are

(66)

e j2pDl(m21) sin u

where Dl is the element spacing in wavelengths, m is the number of sensors, and u is the angle to a far field

FIG. 3. Mean of the estimated u 6 one standard deviation (in degrees) versus true direction of arrival, when assuming an identity noise covariance matrix.

197

FIG. 4. Mean of the estimated u 6 one standard deviation (in degrees) versus true direction of arrival, when using seven noise covariance parameters.

6.3. Direction and Angular Spread Estimation in Wireless Communications

used to model the background noise as in (67)–(68). The overall bias has decreased compared to the first case, whereas the standard deviation has increased slightly in a few cases. This was expected since more parameters are estimated compared to the white noise case. In Fig. 5 the overall bias and root-meansquared error is displayed as a function of the number of noise parameters. As the number of noise parameters increases, the performance improves up to a point. When using too many parameters, problems with identifiability may arise, degrading the accuracy of the estimates. Furthermore, the minimization problem can become ill-conditioned and searching for the global minimum is difficult. Determining the appropriate noise model or, in the present case, the appropriate model order is an intricate problem. Our experience is that the detection methods presented in Section 5 result in model order estimates that are too large for practical purposes. Model perturbations are most likely the dominant source of errors here, whereas the detection methods are tailored to errors caused by additive noise. In [20], the asymptotic distribution of the estimates is analyzed as a function of model perturbations, which may provide some guidelines in choosing an appropriate noise model.

Antenna arrays at the base stations of mobile communication systems may be used to increase system capacity and improve coverage [65–67]. When designing a communication system using antenna arrays, it is essential to have an accurate model of the channel taking the spatial dimension into account. The point source model is often inappropriate owing to multipath propagation. Non-line-of-sight propagation and local scattering at the mobiles is commonplace, leading to a stochastic channel. The goal of this section is to describe a channel model and estimate its parameters from experimental data using the COMET technique.

6.3.1. Channel Modeling In [19,68], a model of the Rayleigh fading due to local scattering is developed taking the spatial dimension into account. The propagation between the mobile and the array is modeled as a superposition of a large number of rays originating from local scatterers in the vicinity of the mobile. The array output vector, assuming a single mobile, can then be ex-

198

FIG. 5. The average bias and root-mean-squared error of the direction estimate, over the set of 18 experiments, versus the number of noise covariance parameters.

channel. These variations are referred to as Rayleigh fading and have been extensively studied for the scalar channel. Below we develop a vector channel model for Rayleigh fading. We will assume independent scatterers with random phase; thus the expected value of the channel is zero:

pressed as y(t) 5

o g (t)s(t 2 t (t))a(h 1 f (t)) 1 e(t), n

n

n

(69)

n

where s(t) is the complex envelope of the transmitted signal, a(h) is the array response vector, gn (t) are the complex gains, tn (t) are the relative time delays, fn (t) is the angular deviation of the nth multipath component from the nominal direction h to the mobile, and e(t) is additive spatially white noise. We will assume that the relative time delays for different propagation paths are small compared to the inverse of the bandwidth of the communication signal. This is a reasonable assumption if the scattering originates in the vicinity of the mobile. In this case, the time delay can be modeled as a phase shift that can be incorporated in gn (t), and we have

E5v(t)6 5

(71)

n

Assuming a large number of scatterers the secondorder moment of the channel transfer vector can be modeled as E5v(t)v*(t)6

5o o 5o

5E

y(t) 5 s(t)

o

n

n

n

gn (t)a(h 1 fn (t)) 1 e(t)5 v(t)s(t) 1 e(t).

o E5g (t)6a(h 1 f (t)) 5 0.

(72)

0gn (t)0 2a(h 1 fn (t))a*(h 1 fn (t))

(73)

E50gf (t)0 26a(h 1 f)a*(h 1 f) df.

(74)

5E

(70)

6 6

gn(t)g*m (t)a(h 1 fn (t))a*(h 1 fm(t))

m

n

n

<

The channel transfer vector, v(t), is sometimes referred to as the spatial signature of the source. Even small movements (on the order of a wavelength, about 0.35 m in the present case) of the mobile will change the scattering, causing large variations of the

e

p

2p

We will assume that the random process 0gf (t)0 2 is stationary during the data collection process that is described in what follows. This is reasonable if the

199

mobile moves within the local scattering region, i.e., the distance and angle h to the mobile are approximately constant. If the angular distribution of scatterers is assumed to be a Gaussian function then the expected value of the angular power distribution, E50 gf (t)0 26, is a Gaussian function of f, i.e., E50 gf (t)0 2 6 5

1

Î2psh

2

e2f

2

/2sh

,

the uplink, the covariance matrix of the downlink channel can be estimated. This allows the construction of transmit weights that either minimize or maximize the expected received signal power at a specific mobile. Therefore, given estimates of h and sh it may be possible to form downlink spatial beamformers that direct or suppress energy in certain directions.

(75)

6.3.2. Channel Estimation The COMET technique will now be used to estimate the parameters in the array covariance model (77). The primary purpose is to estimate sh to gain insight into the propagation conditions that may be expected in mobile communications; see also [19]. Although the accuracy of the parameter estimate obtained with COMET is excellent, computationally less complex estimators may be formulated for realtime angular spread estimation. In [70] standard rooting techniques such as root-MUSIC, ESPRIT, and MODE are modified to provide estimates of the angular spread. Note from (77) that the covariance matrix depends linearly on the signal power, p, and noise variance, s2e . Thus,

leading to E5v(t)v*(t)6 <

1

Î2psh

e

p

2p

2

e 2f

2

/2sh

a(h 1 f)a*(h 1 f) df.

(76)

The parameter, sh, is referred to as the angular spread. Now, assume an m element uniform linear array with element spacing Dl in wavelengths, and a constant modulus signal s(t) with power p. Then the signal received at the array has zero mean and the covariance matrix [19] E5 y(t)y*(t)6 5 R(h, sh, p, s2e ) < pa(h)a*(h) ( B(h, sh) 1 s2e I

(77)

a(h) 5 [1, e j2pDl sin h, . . . , e j2pDl(m21) sin h ] T

(78)

2 2 sh

5B(h, sh )6 kl 5 e 22[pDl(k2l)]

cos2 h

,

F(h, sh ) 5 [vec (a(h)a*(h) ( B(h, sh )) vec (I)] a5

(79)

p

3s 4 . 2 e

The COMET estimates of u 5 [h from (35).

where s2e is the variance of the additive noise and ( denotes elementwise multiplication. Note that the covariance matrix of the channel is parameterized by the nominal direction to the mobile, h, and the angular spread, sh. The angular spread is a critical parameter when attempting to exploit the spatial dimension in a wireless communication system. When using base station antenna arrays to increase system capacity, spatially selective reception and transmission must be achieved. For the uplink (mobile to base), there are numerous ways to estimate the channel and also achieve spatial noise and interference rejection [69]. However, since current cellular systems in general use different carrier frequencies for the uplink and the downlink (base to mobile), the downlink channel will not be the same as the uplink channel (the channel is wavelength dependent; see (78)–(79)). Consequently, it is more difficult to achieve accurate spatial discrimination on the downlink. However, if the nominal direction and angular spread of a mobile can be determined for

(80) (81)

sh ] T are obtained

6.3.3. Experimental Results In this section, the COMET estimator is applied to data collected from an experimental array used in a number of field tests. The data have been collected and made available by Ericsson Radio Systems [71]. A uniform linear array with 10 dipole antennas spaced 0.4 wavelengths apart at 870 MHz was mounted on the roof of a 30-m-high building in Kista, a suburb of Stockholm. The broadside of the array was directed toward the mobile, which was transmitting a sine wave. Data from three different experiments were collected and processed. Two of the mobile sites (referenced as Ring and Sja¨ll) were in a suburban environment and one (Skog) was in a forest. In all cases the mobile was positioned about 1 km from the array and there was no direct line of sight to the base station. In each experiment, five trials were conducted 20 m apart along a straight line. For each trial, 16 data sets were collected. A data set consists of 170 samples of the array output

200

FIG. 6. The results from the five trials at the mobile site ‘‘Sja¨ll.’’ For each trial, 17 estimates of h and sh are computed. Shown is the mean of the nominal direction estimates plus/minus the mean of the angular spread estimates.

during the collection/observation of which the mobile was stationary. The mobile was then moved 0.1 m and the next data set collected. For each trial, 17 COMET estimates of the nominal direction, h, and the angular spread, sh, were obtained and sample statistics of these estimates were formed. Each estimate of h and sh was obtained from a sample covariance matrix formed by taking 10 vector samples from each data set in a trial (10 3 16 data sets 5 160 array output samples). Each data set may be viewed as measuring one realization of the channel. Thus, the sample covariance matrix formed corresponds to 16 different channel realizations. Note that this was done with the purpose of emulating a case in which the mobile moves during the data collection. (A data set of 160 samples, made up as previously described, corresponds to a distance of 1.5 m traveled by the mobile.) The results obtained by applying the COMET method to the experimental data are shown in Figs. 6, 7, and 8. The figures show the mean of nominal direction estimates plus/minus the mean angular spread estimates for each trial. Table 1 displays the sample statistics based on the 17 estimates from each trial. It is possible to see a trend in the direction

estimates corresponding to the 80-m movement (approximately 5° at the distance of 1 km) of the mobile for Ring and Sja¨ll, where the route was approximately parallel to the array. There is no such trend in the results for Skog, where the mobile moved toward the array. Note that the angular spread estimates are similar for all five trials within an experiment. This is to be expected since the environment causing the local scattering at the mobile does not change much during the 80-m movement. The angular spreading is largest (approximately 3.5°) for Skog, probably due to the diffuse scattering caused by the forest. The results from Ring show angular spreading of the same order as that for Skog, whereas the angular spread is about 2° for Sja¨ll. Among the three experiments, the conditions for Sja¨ll were ‘‘closest’’ to a line-of-sight environment resulting in smaller and more stable angular spread estimates. The standard deviations of the nominal direction and angular spread estimates are very small, indicating a small amount of noise as well as stationary measurement conditions. To ease the interpretation of the results, we note that an angular spread of 3° at 1 km distance corresponds to a scattering radius of approximately 50 m (or 180 wavelengths at 870 MHz).

201

FIG. 7. The results from the five trials at the mobile site ‘‘Skog.’’ For each trial, 17 estimates of h and sh are computed. Shown is the mean of the nominal direction estimates plus/minus the mean of the angular spread estimates.

6.4. Estimation and Detection in Radar Signal Processing

ploiting the structure of the received signals. The return signals from different pulses are assumed coherent (a phase reference is required over multiple pulses as will become clear shortly). Herein, we do not restrict the return signals to be coherent from pulse to pulse. We apply the detection schemes of Section 5 to a situation with severe jamming and show that the COMET estimator allows accurate detection. In general, the channel outputs are passed through a filter matched to the pulse shape used; each sample of the filter outputs corresponds to a time delay, i.e., a certain range or range bin. Let x(t) denote the complex baseband antenna array output for a certain range bin where a target is assumed to be present. Our purpose here is to investigate methods for detecting whether a signal (target) is really present in this range bin. Several return pulses are sampled and x(t) denotes the array output snapshot corresponding to the tth pulse, t 5 1, . . . , N. The antenna output is modeled as [72,76]

In a pulsed radar system, the reflected energy from transmitted pulses is used to detect, localize, and possibly identify targets. The pulses are focused spatially by using an antenna array and the return signal is sampled (after appropriate filtering) at different times corresponding to different ranges (range binning). The return signal is often obtained from a beamformer focused in the direction of interest, the output of which is often referred to as the main channel. These radar systems are sensitive to jamming signals that may be present in the sidelobes of the array. To handle this problem, one or more reference signals can be obtained from adjacent sensors or from a beamformer focused away from the direction of interest. The reference signal(s) are obtained from the auxiliary channel(s). If the jamming signal is strong relative to the return signal, the auxiliary channel may provide good estimate of the interference which may be used to cancel the jamming signal in the main channel [72,73]. This technique may experience problems when the return signal from the target is too strong or the jammer is too close to the target. In [74,75], a model-based estimation and detection scheme is introduced, ex-

x(t) 5 a(h)s(t) 1 e(t),

(82)

where a(h) is the array response vector, h is the direction to the source, s(t) is the source signal after 202

FIG. 8. The results from the five trials at the mobile site ‘‘Ring.’’ For each trial, 17 estimates of h and sh are computed. Shown is the mean of the nominal direction estimates plus/minus the mean of the angular spread estimates.

where h is the direction to the target, relative to array broadside. In general, both the azimuth and elevation of a target must be determined. However, for simplicity, we assume that the array and the (far field) target belong to the same elevation plan, and hence consider the azimuth as the only unknown. In radar systems, the number of antenna elements is usually quite large. In general only a small number of channels can be digitized, allowing more advanced signal processing. The main channel provides a spatially focused beam which is a linear combination of the antenna outputs. The auxiliary channels are also obtained as linear combinations of the antenna outputs and often correspond to just individual antennas. Thus, the digitized signals can be modeled as (ignoring quantization effects)

matched filtering, and e(t) is the additive noise. Different assumptions can be made on the source signal. In [74,75], it was modeled as a complex sinusoid with an unknown Doppler frequency. Here, the source signal is assumed to have phase and amplitude that vary randomly from pulse to pulse: s(t) 5 bt e jf t.

(83)

Thus, phase coherence between pulses is not required. In the following examples, we consider the case of a target signal, an uncorrelated jammer, and receiver noise.

6.4.1. A Simulation Example We present an example comparing the detection performances of the methods in Section 5. A uniform linear array of 30 omni-directional antennas, spaced half wavelength apart, is assumed. Thus the array response is given (cf. (66)) by

3 4

y(t) 5

e jp sin h , · · · e j29p sin h

h

(85)

where the first column of the m 3 l matrix Th contains the beamforming weights for the main channel, y1 (t), focused in direction h and the subsequent columns are the weights for the l 2 1 auxiliary channels, z(t). The main channel is formed by using a Chebyshev beamformer steered toward array

1

a(h) 5

y1 (t)

3 z(t) 4 5 T*x(t),

(84)

203

broadside h 5 0° with a 3-dB beamwidth of 6°. The first four elements in the array are chosen as the auxiliary channels; thus, l 5 5. A target signal from direction h1 and with power p1 is present as well as a jamming signal from direction h2 and with power p2. The covariance model assumed for y(t) has the structure R 5 E5y(t)y*(t)6 5 p1T*ha(h1 ) a*(h1 )Th 1 p2T*ha(h2 )a*(h2 )Th 1 s2Q,

(86)

where Q 5 T*hTh is an l 3 l known matrix, and s2 is the unknown noise variance. The COMET method will be used to detect the presence of a target and estimate its parameters, in the presence of a jamming signal. By applying COMET to the covariance model above, a cost function depending only on h1 and h2 is obtained. We will investigate the ability to detect a target at the array broadside (h1 5 0°) when a jammer is present in the first sidelobe, which occurs at 7.5°; thus, h2 5 7.5°. The first sidelobe provides 237 dB of power suppression compared to the main lobe. The jammer power is 10 dB above the noise power and the power of the target signal will be measured relative to the noise. For each detection trial, N 5 256 vector samples are collected and processed. The COMET estimator is applied to the 256 samples of the five channel data under the assumption that no signal is present. Only the noise vari-

FIG. 9. Probability of correct detection. The solid line corresponds to the COMET detection method, the dashed line to GLRT, and the dotted line to COMET-AIC. A jammer is present in the first sidelobe of the main channel with a power of 10 dB above the noise power.

ance, s2, is estimated in this case. If this hypothesis is rejected, one signal is estimated and so on until a hypothesis is accepted. The threshold is chosen, according to the distribution of the limiting cost function derived in Section 5, such that the probability of false detection is 0.005 when no signal is present. Also, the COMET-AIC criterion is evaluated for increasing number of signals until an increase in cost is obtained. The minimum provides an estimate of the number of signals as described in Section 5. Detection based on the GLRT criterion using the COMET parameter estimates is also performed using the same threshold as above. In Fig. 9 the probability of correct detection is displayed versus the signal-to-noise ratio of the target signal. Figure 10 displays the root-meansquared error of the direction estimate to the target source versus the signal-to-noise ratio. Only the cases where the correct number of sources where detected are included. For all these cases, the estimated direction to the target is close to the true direction. Thus, when the correct number of sources is detected, the direction estimate is quite accurate. In this example, the detection performances of the methods are comparable. The main advantage of COMET as compared to maximum likelihood is the simplified estimation procedure. The detection procedures are coupled with the estimation of the signal parameters providing robustness against any jamming signal and any possible pointing errors. Of course, as with any model-based technique, the

TABLE 1 Means and Standard Deviations of the Nominal Direction and Angular Spread Estimates for the Five Trials in Each Experiment Experiment

Trial

Nominal direction (°)

STD

Angular spread (°)

STD

Sja¨ll

1 2 3 4 5

20.35 0.91 1.3 2.2 4.3

0.10 0.031 0.061 0.036 0.032

2.3 1.8 2.2 2.3 2.3

0.098 0.046 0.049 0.041 0.03

Skog

1 2 3 4 5

0.52 0.81 21.7 1.1 20.28

0.078 0.12 0.081 0.071 0.095

4.0 3.3 3.8 2.9 3.5

0.095 0.11 0.094 0.068 0.11

Ring

1 2 3 4 5

0.049 0.075 0.023 0.063 0.034

4.2 3.7 2.7 4.0 3.0

0.042 0.041 0.022 0.067 0.035

212 28.4 26.5 25.7 25.0

Note. Each direction and angular spread estimate is based on 160 vector samples. The sample statistics above are based on 17 parameter estimates in each trial.

204

In this example, 16 array output snapshots (N 5 16) are available for detecting the target signal. This is repeated 404 times on independent data sets. In Fig. 11 the probability of correct detection with COMET-AIC is displayed versus the target signal-tonoise ratio. In Fig. 12 the corresponding bias and RMS values of the direction estimates for the target signal are displayed. Note that accurate signal detection and estimation is achieved when the SNR is greater than 5 dB with only 16 array snapshots and in the presence of a 25 dB stronger jammer. The chi-squared-based COMET and GLRT detection methods perform very poorly and overestimate the number of signals when the threshold is chosen according to the limiting distribution function of the minimized criteria as in the previous example. Therefore, these results are not included in the figure. This is most likely due to modeling errors and also the small number of data samples. The performance of COMET and GLRT can of course be improved by adjusting the threshold appropriately; this is, however, avoided in the COMET-AIC procedure.

FIG. 10. Bias (x) plus/minus the standard deviation of the direction estimate to the target signal. Only the cases where SNR $ 220 dB are displayed.

robustness against perturbations in the array response model must be considered in a practical system as observed in the next example.

7. CONCLUSIONS

6.4.2. A Radar Experiment Here we attempt to detect and estimate a target signal in data collected by an experimental radar system in the presence of a strong jammer. These data were collected under the ARPA/NAVY Mountaintop Program and the reader is referred to [77] for more details concerning the experiment. A uniform linear antenna array which is horizontally oriented with respect to the ground is used. The array is composed of 14 horizontally polarized elements with an element spacing of 0.333 m. All elements are equipped with coherent receivers and thus no prior beamforming is performed. A remotely sited jammer at approximately 42° relative to the array broadside and with transmit frequency 435 MHz provided the only actual signal which is present in the data. Also, scattering of the jamming signal caused by the environment may be present. A target signal was artificially introduced in the data at array broadside (0°), for a specific range bin. It is a constant modulus signal with random phase and its power was varied. Furthermore, temporally and spatially white Gaussian noise was added to the data scaled such that the jammer power is 30 dB above the noise power. Thus, a signal described by (82) was added to the collected data. The jammer power is estimated as the average power at an element of the original data set. In the estimation process, the array is assumed to be an ideal ULA and the noise is assumed to be spatially white.

Signal detection and parameter estimation from sensor array measurements are problems that arise in several engineering applications. Examples include target localization with phased array radars, detection and classification of sources with underwater arrays, and signal separation in wireless communications using antenna arrays. When digital multi-

FIG. 11. Probability of correct detection for COMET-AIC. A jammer is present at 42° with a power of 30 dB above the noise power.

205

in wireless communications. The angular spread of sources in typical mobile communication environments was estimated. Finally, a radar application was considered. Using a COMET based scheme, accurate signal detection and estimation could be achieved.

APPENDIX A: COMPACT EXPRESSION FOR THE CRB Below, we provide a simple derivation of Cov (uˆ ) 5 CRB. The COMET estimate is given by the minimizing argument of (30):

FIG. 12. Bias (x) plus/minus the standard deviation of the direction estimate to the target signal. Only the cases where SNR $ 1 dB are displayed.

ˆ 21/2rˆ. ˆ 21/2P'ˆ 21/2 W f(u) 5 rˆ*W W F(u)

(A.1)

Expanding the derivative of the cost function around the uˆ gives

channel data are available, advanced estimation and detection schemes can be implemented improving performance significantly compared to single channel or analog processing methods. We have presented a class of covariance matching estimation techniques that are well suited for several array processing problems with spatially correlated noise and structured signal covariance matrix. The technique can easily accommodate knowledge about uncorrelated signals and spatial correlation structure of the noise. Typically, the well-known subspace-based methods are not applicable to these cases. In situations where the noise and signal covariance matrices are parameterized linearly, the COMET approach can provide a more attractive solution than the maximum likelihood estimator. We derived the COMET cost function from first principles and at the same time established that COMET is a large-sample approximation of the maximum likelihood method. Identifiability conditions from second-order statistics were given under the most general conditions. A compact expression for the estimation error covariance matrix was also provided. The detection problem was addressed and two tests based on COMET were presented that are naturally associated with the estimation procedure. Three array processing applications have been presented where the COMET approach showed good promise. An underwater source localization problem was examined. By appropriately modeling the spatial covariance of the background noise, a significant reduction in bias was demonstrated. The COMET approach was then used to estimate a channel model

0 5 ­u f(uˆ ) 5

­f(uˆ ) ­u

. ­u f(u)0u5u0 1 [lim ­uu f(u)] N=`

0

(uˆ 2 u0)

u5u0

⇒ uˆ 2 u0. 2H21g,

(A.2)

(A.3)

where H 5 lim ­uu f(u)0 u5u0 N=`

and g 5 ­u f(u)0 u5u0.

(A.4)

Straightforward differentiation of (A.1) gives '

ˆ 21/25P ˆ 21/2 W ˆ 21/2­ F(F*W ˆ 21F)21 ­k f (u) 5 2rˆ*W k W F ˆ 21/2 1 (· · ·)*6W ˆ 21/2rˆ 3 F*W ˆ 21/2­ Faˆ 6 ˆ 21/2 P'ˆ 21/2 W . 22 Re 5rˆ*W k W F ˆ 21/2P'ˆ 21/2 W ˆ 21/2­ Faˆ , 5 22rˆ*W k W F

(A.5)

where ­k 5 ­/­uk, (· · ·)* denotes a term equal to the conjugate transpose of the previous term, and a ˆ 5 a ˆ (u) is given by (28). Note that the derivative above can be shown to be real-valued just as aˆ is shown to be real-valued in Section 3. The limiting second 206

APPENDIX B: ASYMPTOTIC DISTRIBUTION OF COMET COST FUNCTION

derivative is given by lim ­kpf(u) 5 2r*W21F(F*W21F) 21

N=`

Let 3 ­pF*W

21/2 ' PW21/2FW21/2­kFa ' PW21/2FW21/2­pFa.

21/2

5 2a*­kF*W

dk 5 ­kF(u)a 5 ­kr,

(B.1)

D 5 [d1 · · · dnu ].

(B.2)

(A.6) and let

Next, note that the expression above for ­kf (u), when evaluated at the true signal parameters, can be significantly simplified since

According to the statistical analysis previously done ˆ 21/2 5 (rˆ 2 r)*W ˆ 21/2P'ˆ 21/2 W ˆ 21/2 ˆ 21/2P'ˆ 21/2 W rˆ* W W F W F ' PW21/2FW21/2

21/2

5 (rˆ 2 r)*W

uˆ 2 u0 . 2H21g,

1 o(1/ÎN)

5 O(1/ÎN).

(B.3)

where (see (A.6), (A.8)) (A.7) '

H 5 2D*W21/2PW21/2FW21/2D

(B.4)

By making use of (A.7), one can write (asymptotically) '

­kf(u) . 22(rˆ 2

g 5 2D*W21/2PW21/2FW21/2 (rˆ 2 r).

' r)*W21/2PW21/2FW21/2­kFa '

5 22a*­kF*W21/2PW21/2FW21/2(rˆ 2 r).

Using (B.3) we obtain

(A.8)

f(u0 ) . f(uˆ ) 1 ­u f(uˆ )(u0 2 uˆ )

Observing that

1

lim NE5­k f(u)­p f(u)6

N=`

'

5 4a*­kF*W21/2PW21/2FW21/2­pFa

'

2

3

4

(uˆ 2 u0 )* lim ­uu f(u0 ) (uˆ 2 u0 ) 1 2

N=`

(uˆ 2 u0 )*H(uˆ 2 u0 ),

(A.10)

ˆ 21/2P'ˆ 21/2 W ˆ 21/2rˆ 2 1g*H21g f(uˆ ) . rˆ*W 2 W F(u)

Note that ­pFa evaluated at a0 may be written as

. (rˆ 2 r)*W21/2PW21/2F(u)W21/2(rˆ 2 r)

­up

a5

­r ­up

5

­ vec (R) ­up

(B.6)

where we also have used the fact that f 8(uˆ ) 5 0. The following asymptotic approximation for f (uˆ ) is obtained by combining (B.6) and (B.3),

[CRB21 ] kp 5 Na*­kF*W21/2PW21/2F(u)W21/2­pFa.

1

5 f(uˆ ) 1

(A.9)

and using (A.3) it is readily concluded that

­F

(B.5)

'

'

.

2 2(rˆ 2 r)*W21/2PW21/2F(u)

(A.11)

'

3 W21/2DH21D*W21/2PW21/2F(u)W21/2(rˆ 2 r)

Hence, the expression in (A.10) can be rewritten in matrix form as in (39).

˚ W21/2 (rˆ 2 r), 5 (rˆ 2 r)*W21/2P 207

(B.7)

REFERENCES

where ˚ 5 P' 21/2 P W F(u) '

'

2 2PW21/2F(u)W21/2DH21D*W21/2PW21/2F(u) .

1. Schmidt, R. O. Multiple emitter location and signal parameter estimation. IEEE Trans. Antennas and Propagation 34 (1986), 276–280. 2. Roy, R., and Kailath, T. ESPRIT—Estimation of signal parameters via rotational invariance techniques. IEEE Trans. Acoust. Speech Signal Process. ASSP-37 (1989), 984–995. 3. Stoica, P., and Sharman, K. A novel eigenanalysis method for direction estimation. Proc. IEE F (1990), 19–26. 4. Stoica, P., and Sharman, K. Maximum likelihood methods for direction-of-arrival estimation. IEEE Trans. Acoust. Speech Signal Process. ASSP-38 (1990), 1132–1143. 5. Viberg, M., Ottersten, B., and Kailath, T. Detection and estimation in sensor arrays using weighted subspace fitting. IEEE Trans. Signal Process. SP-39 (1991), 2436–2449. 6. Stoica, P., and Nehorai, A. MUSIC, maximum likelihood and Crame´r-Rao bound. IEEE Trans. Acoust. Speech Signal Process. ASSP-37 (1989), 720–741. 7. Stoica, P., and Nehorai, A. MUSIC, maximum likelihood and Crame´r–Rao bound: Further results and comparisons. IEEE Trans. Acoust. Speech Signal Process. ASSP-38 (1990), 2140– 2150. 8. Wu, Q., and Wong, K. M. UN-MUSIC and UN-CLE: An application of generalized canonical correlation analysis to the estimation of the directions of arrival. IEEE Trans. Signal Process. 42 (1994). 9. Stoica, P., and Ha¨ndel, P., and Nehorai, A. Improved sequential MUSIC. IEEE Trans. Aerospace Electron. Systems AES-31 (1995), 1230–1239. 10. Haardt, M., and Nossek, J. A. Unitary ESPRIT: How to obtain increased estimation accuracy with a reduced computational burden. IEEE Trans. Signal Process. SP-43 (1995), 1232– 1242. 11. Bo¨hme, J. F. Estimation of spectral parameters of correlated signals in wavefields. Signal Processing 11 (1986), 329–337. 12. Stoica, P., and Nehorai, A. Performance study of conditional and unconditional direction-of-arrival estimation. IEEE Trans. Acoust. Speech Signal Process. ASSP-38 (1990), 1783–1795. 13. Ottersten, B., Viberg, M., Stoica, P., and Nehorai, A. Exact and large sample ML techniques for parameter estimation and detection in array processing. In Radar Array Processing (Haykin, Litva, and Shepherd, Eds.). Springer-Verlag, Berlin, 1993, pp. 99–151. 14. Stoica, P., and Moses, R. Introduction to Spectral Analysis. Prentice Hall, Upper Saddle River, NJ, 1997. 15. Forster, P., and Vezzosi, G. Approximate ML after spheroidal beamforming in the narrowband and wideband cases. InProc. EUSIPCO, Grenoble, France, 1988, pp. 307–310. 16. Kraus, D., and Bo¨hme, J. F. Asymptotic and empirical results on approximate maximum likelihood and least squares methods for sensor array processing. In Proc. ICASSP 90, Albuquerque, NM, 1990, pp. 2795–2798. 17. Kraus, D. Approximative Maximum-Likelihood-Scha¨tzung and verwandte Verfaren zur Ortung und Signalscha¨tzung mit Sensorgruppen. Ph.D. thesis, Aachen, Germany, 1993. 18. Rao, C. R., Zhang, L., and Zhao, L. C. Multiple target angle tracking using sensor array outputs. IEEE Trans. Aerospace Electron. Systems AES-29 (1993), 268–271. 19. Trump, T., and Ottersten, B. Estimation of nominal direction of arrival and angular spread using an array of sensors. Signal Process. 50 (1996), 57–69.

(B.8)

˚ is idempotent. First note Next, we establish that P that ' ˚ 2 5 P' 21/2 1 4PW21/2F(u)W21/2DH21D* W21/2 P W F(u) '

'

3 PW21/2F(u)W21/2DH21D*W21/2PW21/2F(u) '

'

2 4PW21/2F(u)W21/2DH21D*W21/2PW21/2F(u) '

'

5 PW21/2F(u) 2 2PW21/2F(u)W21/2D '

3 (2H21 2 2H21D*W21/2PW21/2F(u)W21/2DH21) '

3 D* W21/2PW21/2F(u) .

(B.9)

Using (B.4) we have '

2H21 2 2H21D*W21/2PW21/2F(u)W21/2DH21 5 H21 (2I 2 HH21 )5 H21.

(B.10)

˚ . Note also that the rank of P ˚ is given ˚25P Thus, P by ˚ ) 5 Tr 5P ˚ 6 5 Tr 5P' 21/2 6 rank (P W F(u) '

'

2 Tr 5PW21/2F(u)W21/2DH21D* W21/2PW21/2F(u) 6 5 m 2 2 nµ 2 ns 2 nu,

(B.11)

' 2 since PW 21/2F(u) has rank m 2 nµ 2 ns, D has rank nu, and H has full rank under the identifiability conditions. In (B.7), the vector W21/2(rˆ 2 r) is AsN (0, I) and ˚ is idempotent with rank equal to m 2 2 nµ 2 since P ns 2 nu, the minimum of the cost function is asymptotically chi-square distributed with m 2 2 nu 2 nµ 2 ns degrees of freedom. To summarize, we have proved that when the correct number of parameters is assumed, i.e., the parameters nu, nµ, and ns match the number of parameters in R,

f(uˆ ) , Asx2 (m 2 2 nu 2 nµ 2 ns ).

(B.12)

208

20. Go¨ransson, B. Direction estimation in unknown noise fields by spatial noise suppression. In IEEE/IEE Workshop on Signal Processing in Multipath Environments, Glasgow, Scotland, Apr. 1995. 21. Wax, M., Sheinvald, J., and Weiss, A. J. Detection and localization in colored noise via generalized least squares. IEEE Trans. Acoust. Speech Signal Process. ASSP-44 (1996), 1734–1743. 22. Chonavel, T., and Loubaton, P. Estimation of the autocorrelation coefficients of band-limited signals. Application to the source localization problem. In Proc. EUSIPCO, 1992, pp. 1825–1828. 23. Friedlander, B., and Weiss, A. Direction finding using noise covariance modeling. IEEE Trans. Signal Process. 43 (1995), 1557–1567. 24. Giannakis, G. B., and Halford, S. D. Asymptotically optimal blind fractionally spaced channel estimation and performance analysis. IEEE Trans. Signal Process. 45 (1997), 1815–1830. 25. Anderson, T. W. An Introduction to Multivariate Statistical Analysis, 2nd ed. Wiley, New York, 1984. 26. Vanpoucke, F., and Paulraj, A. A harmonic noise model for direction finding in colored ambient noise. IEEE Signal Process. Lett. SPL-2 (1995), 135–137. 27. Go¨ransson, B. Direction finding in the presence of spatially correlated noise fields. In Proc. EUSIPCO, EURASIP, Sept. 1994, pp. 1321–1324. 28. Valaee, S., Champagne, B., and Kabal, P. Parametric localization of distributed sources. IEEE Trans. Signal Process. SP-43 (1995), 2144–2153. 29. Meng, Y., Stoica, P., and Wong, K. M. Estimation of the directions of arrival of spatially dispersed signals in array processing. IEE Proc. Radar Sonar Navig. 143 (1996), 1–9. 30. Borgiotti, G. V., and Kaplan, L. J. Superresolution of uncorrelated interference sources by using adaptive array techniques. IEEE Trans. Antennas and Propagation 27 (1979), 842–845. 31. Stoica, P., and So¨derstro¨m, T. On reparameterization of loss functions used in estimation and the invariance principle. Signal Process. 17 (1989), 383–387. 32. Schmidt, R. O. A Signal Subspace Approach to Multiple Emitter Location and Spectral Estimation. Ph.D. thesis, Stanford University, Stanford, CA, 1981. 33. Stoica, P., and So¨derstro¨m, T. On array signal processing in spatially correlated noise fields. IEEE Trans. Circuits Systems II 39 (1992), 879–882. 34. Vanpoucke, F. Algorithms and Architectures for Adaptive Array Signal Processing. Ph.D. thesis, Katholieke Universiteit Leuven, 1995. 35. Wax, M., and Ziskind, I. On unique localization of multiple sources by passive sensor arrays. IEEE Trans. Acoust. Speech Signal Process. ASSP-37 (1989), 996–1000. 36. Gershman, A. B., Ermolaev, V. T., and Flaksman, A. G. Angular adaptive resolution of uncorrelated sources. Radiophys. Quantum Electron. 31 (1988), 673–677. 37. Swingler, G. N. An approximate expression for the Crame´r–Rao bound on direction-of-arrival estimation with small ring arrays. J. Acoust. Soc. Amer. 94 (1993), 1400–1403. 38. Jansson, M., Go¨ransson, B., and Ottersten, B. Analysis of a subspace-based spatial frequency estimator. In Proceedings of ICASSP, 1997, pp. 4001–4004. 39. Wax, M. On unique localization of constrained-signals sources. IEEE Trans. Acoust. Speech Signal Process. ASSP-40 (1992), 1542–1547.

40. Linebarger, D. A., Sudborough, I. H., and Tollis, I. G. Difference bases and sparse sensor arrays. IEEE Trans. Inform. Theory IT-39 (1993), 716–721. 41. Bo¨hme, J. F. Separated estimation of wave parameters and spectral parameters by maximum likelihood. In Proc. ICASSP 86, Tokyo, Japan, 1986, pp. 2818–2822. 42. Jaffer, A. G. Maximum likelihood direction finding of stochastic sources: A separable solution. In Proc. ICASSP 88, New York, 1988, Vol. 5, pp. 2893–2896. 43. Stoica, P., and Nehorai, A. On the concentrated stochastic likelihood function in array signal processing. Circuits Syst. Signal Process. 14 (1995), 669–674. 44. Graham, A. Kronecker Products and Matrix Calculus with Applications. Ellis Horwood, Chichester, 1981. 45. Browne, G. W. Generalized least squares estimates in the analysis of covariance structures. South African Statist. J. 8 (1974), 1–24. 46. Mardia, K. V., Kent, J. T., and Bibby, J. M. Multivariate Analysis. Academic Press, London, 1979. 47. Akaike, H. A new look at the statistical model identification. IEEE Trans. Automat. Control. AC-19 (1974), 716–723. 48. Gill, P. E., Murray, W., and Wright, M. H. Practical Optimization. Academic Press, London, 1981. 49. Dennis, J. E., and Schnabel, R. B. Numerical Methods for Unconstrained Optimization and Nonlinear Equations. Prentice Hall, Englewood Cliffs, NJ, 1983. 50. Ziskind, I., and Wax, M. Maximum likelihood localization of multiple sources by alternating projection. IEEE Trans. Acoust. Speech Signal Process ASSP-36 (1988), 1553–1560. 51. So¨derkvist, I. On algorithms for generalized least squares problems with ill-conditioned covariance matrices. Comput. Statist. 11 (1996), 303–313. 52. Gershman, A. B., Turchin, V. I., and Zverev, V. A. Experimental results of localization of moving underwater signal by adaptive beamforming. IEEE Trans. Signal Process. 43 (1995). 53. Thalham, R. J. Noise correlation functions for an isotropic noise field. J. Acoust. Soc. Amer. 69 (1981), 213–215. 54. Li, F., and Vaccaro, R. J. Performance degradation of DOA estimators due to unknown noise fields. IEEE Trans. Signal Process. SP-40 (1992), 686–689. 55. Viberg, M. Sensitivity of parametric direction finding to colored noise fields and undermodeling. Signal Process. 34 (1993), 207–222. 56. Paulraj, A., and Kailath, T. Eigenstructure methods for direction-of-arrival estimation in the presence of unknown noise fields. IEEE Trans. Acoust. Speech Signal Process. ASSP-34 (1986), 13–20. 57. Le Cadre, J. P. Parametric methods for spatial signal processing in the presence of unknown colored noise fields. IEEE Trans. Acoust. Speech Signal Process. ASSP-37 (1989), 965– 983. 58. Stoica, P., Viberg, M., and Ottersten, B. Instrumental variable approach to array processing in spatially correlated noise fields. IEEE Trans. Signal Process. SP-42 (1994), 121–133. 59. Viberg, M., Stoica, P., and Ottersten, B. Array processing in correlated noise fields based on instrumental variables and subspace fitting. IEEE Trans. Signal Process. 43 (1995), 1187–1199. 60. Kay, S. M., and Nagesha, V. Maximum likelihood estimation of signals in autoregressive noise. IEEE Trans. Signal Process. 42 (1994), 88–101.

209

61. Nagesha, V., and Kay, S. M. Maximum likelihood estimation for array processing in colored noise. IEEE Trans. Signal Process. SP-44 (1996), 169–180. 62. Wong, K. M., Reilly, J. P., Wu, Q., and Qiao, S. Estimation of the directions of arrival of signals in unknown correlated noise, parts I and II. IEEE Trans. Signal Process. SP-40 (1992), 2007–2028. 63. Stoica, P., Viberg, M., Wong, K. M., and Wu, Q. Maximumlikelihood bearing estimation with partly calibrated arrays in spatially correlated noise fields. IEEE Trans. Signal Process. SP-44 (1996), 888–899. 64. Bo¨hme, J. F., and Kraus, D. On least squares methods for direction of arrival estimation in the presence of unknown noise fields. In Proc. ICASSP 88, New York, 1988, pp. 2833– 2836. 65. Andersson, S., Millnert, M., Viberg, M., and Wahlberg, B. An adaptive array for mobile communication systems. IEEE Trans. Vehicular Technol. 40 (1991), 230–236. 66. Zetterberg, P., and Ottersten, B. The spectrum efficiency of a base-station antenna array system for spatially selective transmission. IEEE Trans. Vehicular Technol. 44 (1995), 651–660. 67. Naguib, A. F., Paulraj, A., and Kailath, T. Capacity improvement with base-station antenna arrays in cellular CDMA. IEEE Trans. Vehicular Technology. 43 (1994), 691–698. 68. Zetterberg, P. Mobile Cellular Communications with Base Station Antenna Arrays: Spectrum Efficiency, Algorithms and Propagation Models. Ph.D. thesis, Kungliga Tekniska Ho¨gskolan, Stockholm, Sweden, 1997. 69. Ottersten, B. Array processing for wireless communication. In 8th IEEE Signal Processing Workshop on Statistical Signal and Array Processing, Corfu, Greece, 1996, pp. 466–473. 70. Bengtsson, M., and Ottersten, B. Rooting techniques for estimation of angular spread with an antenna array. In Proceedings of IEEE Vehicular Technology Conference 97, 1997, pp. 1158–1162. 71. Forsse´n, U., Karlsson, J., Johannisson, B., Almgren, M., Lotse, F., and Kronestedt, F. Adaptive antenna arrays for GSM900/DCS1800. In Proceedings of IEEE Vehicular Technology Conference 1994, pp. 605–609. 72. Farina, A. Antenna-Based Signal Processing Techniques for Radar Systems. Artech House, Norwood, MA, 1992. 73. Applebaum, S. P. Adaptive arrays. IEEE Trans. Antennas and Propagation AP-24 (1976), 650–662. ¨ . A comparison of 74. Viberg, M., Ottersten, B., and Erikmats, O model-based detection and adaptive sidelobe cancelling for radar array processing. In Proc. Nordic Antenna Symposium, Eskilstuna, Sweden, 1994. 75. Swindlehurst, A. L., and Stoica, P. Maximum likelihood methods in radar array signal processing. Proc. IEEE, 86 (1998), 421–441. 76. Haykin, S. Radar array processing for angle of arrival estimation. In Array Signal Processing (S. Haykin, Ed.). Prentice Hall, Englewood Cliffs, NJ, 1985. 77. Titi, W., and Marshall, D. F. The ARPA/NAVY mountaintop program: Adaptive signal processing for airborne early warning radar. In Proc. ICASSP 96, Atlanta, Georgia, 1996.

¨ RN OTTERSTEN was born in Stockholm, Sweden, in 1961. BJO He received his M.S. in electrical engineering and applied physics from Linko¨ping University, Linko¨ping, Sweden, in 1986. In 1989 he received his Ph.D. in electrical engineering from Stanford University, Stanford, CA. Dr. Ottersten has held research positions at the Department of Electrical Engineering, Linko¨ping University; the Royal Institute of Technology, Stockholm; and the Information Systems Laboratory, Stanford University. During 1996/1997 Dr. Ottersten was Director of Research at ArrayComm, Inc., in California. In 1993, he received the Signal Processing Society Paper Award. In 1991 he was appointed professor of signal processing at the Royal Institute of Technology (KTH), Stockholm, and he is currently head of the Department of Signals, Sensors & Systems at KTH. His research interests include wireless communications, stochastic signal processing, sensor array processing, and time series analysis. PETRE STOICA received his M.Sc. and D.Sc., both in automatic control, from the Bucharest Polytechnic Institute, Romania, in 1972 and 1979, respectively. In 1993, he was awarded an honorary doctorate by Uppsala University, Sweden. Since 1972, he has been with the Department of Automatic Control and Computers at the Polytechnic Institute of Bucharest, Romania, where he holds the position of professor of system identification and signal processing. He spent 1992, 1993, and the first half of 1994 with the Systems and Control Group at Uppsala University in Sweden, as a guest professor. In the second half of 1994 he held a Chalmers 150th anniversary visiting professorship with the Applied Electronics Department at Chalmers University of Technology, Gothenburg, Sweden. Currently he is affiliated with the Systems and Control Group of Uppsala University. He is also a corresponding member of the Romanian Academy. RICHARD ROY received his B.S. EE and B.S. physics degrees from MIT in 1972; his M.S. physics degree from Stanford in 1973, and his M.S. EE and Ph.D. EE degrees from Stanford in 1975 and 1987, respectively. His Ph.D. research was in the field of digital signal processing, specifically estimation and parameter identification. Dr. Roy has been involved in the intelligence industry while with ESL/TRW, the oil-well services business, and the aerospace industry over the past 20 years. From 1975 through 1993, he was also a member of the Information Systems Laboratory at Stanford University as a graduate student and a research associate, and was also a consultant to industry, working in the areas of parameter estimation and system identification with applications to control and communications systems during that time. In April of 1992 he founded ArrayComm, Inc., a Californiabased corporation involved in the development of intelligent antenna products based on its proprietary SDMA technology for wireless telecommunication applications. The company currently has strategic relationships with several of the world’s largest telecommunications firms and is widely recognized as one of the leaders in its field. Dr. Roy currently serves as chief technical officer of ArrayComm. He is widely published internationally and holds several patents in the area of intelligent antenna technology and wireless telecommunications.

210