On Optimal Dimension Reduction for Sensor Array Signal Processing

(1)

On Optimal Dimension Reduction for Sensor Array Signal Processing

SOREN ^ANDERSON

Department of Electrical Engineering Linkoping University

S-581 83 Linkoping, Sweden

Revised version

Abstract

The computational complexity for direction-of-arrival estimation using sensor arrays increase very rapidly with the number of sensors in the array. One way to lower the amount of computations is to employ some kind of reduction of the data dimension. This is usually accomplished by employing linear transformations for mapping full dimension data into a lower dimensional space. Dierent approaches for selecting these transformations have been proposed. In this paper, a transformation matrix is derived that makes it possible to theoretically attain the full-dimension Cramer-Rao bound also in the reduced space.

A bound on the dimension of the reduced data set is given, above which it is always possible to obtain the same accuracy for the estimates of the source localizations, using the lower-dimension data, as that achievable by using the full dimension data. Furthermore, a method is devised for designing the transformation matrix. Numerical examples, using this design method, are presented, where the achievable performance of the (optimal) Weighted Subspace Fitting method with full dimension data is compared to the performance obtained with reduced dimension data. The problem of estimating parameters of sinusoidal signals from noisy data is also addressed by a direct application of the results derived herein.

1

(2)

\On Optimal Dimension Reduction for

Sensor Array Signal Processing", 26 pages.

Keywords: array processing, dimension reduction, Cram er-Rao bound, cisoids-in-noise

List of Figures

1 Evaluation of the performance measure ( k ), (28 ), over the interval

^;

3 : 5 3 : 5 ] for the transformation matrix obtained from design DOAs

^;

3

^;

1 : 5 0 1 : 5 3 ]

^T

. Scenario: ULA, SNR=3 dB, m = 100 N = 20. 17 2 CRB's and RMS error of DOA estimate vs. DOA of emitter using the

transformation matrix obtained from design-DOAs

^;

3

^;

1 : 5 0 1 : 5 3 ]

^T

. Scenario: ULA, SNR=3 dB, m = 100 N = 20. : : : : : : : : : : : : : : 18 3 RMS error of DOA estimate vs. emitter DOA of left source using the

transformation matrix obtained from design-DOAs

^;

3

^;

1 : 5 0 1 : 5 3 ]

^T

. Scenario: ULA, Symmetrically placed emitters, m = 100 N = 20. : : : : 19 4 RMS error of

^;

0 : 2 DOA estimate vs. out-of-sector DOA using the the

transformation matrix obtained from design-DOAs

^;

1 : 5 0 1 : 5 ]

^T

. Sce- nario: ULA, SNR=3 dB, One out-of-sector emitter, m = 100 N = 10. : 20 5 RMS error of cisoid frequency estimates vs. SNR. Scenario: 100 obser-

vations, i.e., m = 100 N = 1. : : : : : : : : : : : : : : : : : : : : : : : : 21 6 RMS error of cisoid frequency estimates vs. number of data points. Sce-

nario: SNR = 3 dB, N = 1. : : : : : : : : : : : : : : : : : : : : : : : : : 22

2

(3)

1 Introduction

Sensor array signal processing deals with methods for processing measurements of an array of sensors, located in a waveeld at dierent points in space. Typical applications with radiating sources are radio telescopy, where radio sources are measured by means of antenna groups and passive (listening only) sonar, where hydrophone arrays are either towed behind a vessel or dropped into the ocean to measure ship or submarine noise.

The problem in these applications is to detect and estimate incoming signals in order to determine, e.g., the locations of the emitting sources.

The interest in performance improvement for Direction-Of-Arrival (DOA) estimation algorithms can result in antenna arrays composed of a large number of sensor elements.

Since the computational requirements are directly aected by the dimension of the collected data, the burden increases rapidly with the number of sensors. This increase is often proportional to the cube of the dimension, m , of the Element SPace (ESP), i.e., of the dimension of the data. This is the case, e.g., for high-resolution methods based on eigendecomposition of the array covariance matrix, which requires

^O

( m

³

) operations,

1]. Therefore, in order to secure the computability (in terms of time) of the algorithms it is useful to introduce some sort of mapping that reduces the dimension of the data set, from m to r ( r << m ), before applying signal processing algorithms. Another advantage is that if the transformation matrix is implemented using analog technology, there is a need for only r analog-to-digital converters, in contrast to the m ones needed for ESP computations. The space into which the set of full dimensional data is mapped is referred to as the Reduced Dimension BeamSpace (RDBS).

One approach to reduce the amount of computations is to employ a linear (matrix) transformation for mapping the full dimension ESP data into the lower dimensional RDBS, and then apply a signal processing algorithm to this new set of data. The design of the matrix transformation is guided by some subjective criterion as, for instance:

selection of spatial sector, 2]-3], maximization of average signal-to-noise ratio within a specied sector, 4, 5], or minimization of output interference power (assuming a priori knowledge of the interference scenario), 6]. In 7], the design is made by analytic deter-

3

(4)

mination of a transformation matrix through minimization of the RDBS signal-to-noise ratio (SNR) at which two closely spaced emitters can be resolved. The derivation is carried out using the framework of the MUSIC algorithm, 8].

Apart from a low computational complexity, other aspects of performance may be improved due to the reduction of data dimension. These include (for spatial sector selection methods):

- Sources outside the selected sector are ltered out and not detected, thereby the dimension of the parameter space is reduced.

- Lower sensitivity to assumptions on spatial whiteness of the additive noise as pointed out in 2].

- Finally, as indicated by simulations in 4], and shown in 7, 9], the resolution threshold for the MUSIC algorithm can be lower in RDBS than in ESP.

The main purpose of this paper is to analyze under which circumstances it is possible to retain the property of asymptotically achieving the full-dimension CRB, while reducing the dimension of the data set used for the DOA estimation. Somewhat sur- prisingly, it is possible to obtain, in a very simple way, analytical expressions for the requirements to be met by an optimal transformation matrix. These expressions rely on knowledge of the functional form of the array (calibrated array of arbitrary geometry) and the true locations of the sources present. This last condition is, however, not limiting the possibilities of using the results derived herein to outline a design method.

If a spatial sector is specied, in which sources of interest are assumed to be located, it is possible to, in a close-to-optimal way, estimate these DOAs. Numerical examples are presented to illustrate the design method. An alternative approach would be to rst employ beamforming for obtaining crude estimates of the emitter locations as the peak locations in the beamformer output, followed by the design method proposed herein.

As far as the author is aware, there has not been any discussion in the literature of how the variance of the DOA estimates improves, when the dimension of RDBS is increased. This means that, in this respect, there is really no guidelines available for selection of the RDBS dimension. Herein, it is shown that if the RDBS dimension is no less than twice the number of signals, the ESP CRB can always be attained.

4

(5)

2 Data Model

Consider an array of arbitrary geometry consisting of m sensor elements. Assume that d < m wavefronts impinge on the array, due to d far eld sources emitting narrowband signals. The output of the sensor array can then be modeled as a weighted superposition of the d wavefronts, corrupted by additive sensor noise, uncorrelated with the emitter signals. That is, the output of the i

^th

sensor is represented by

x

i

( t ) =

^X^d

j⁼¹

a

i

(

j

) s

j

( t ) + n

i

( t ) : (1) Here, a

i

(

j

) is a complex scalar representing the propagation delay of the j

^th

emitter signal and the gain and phase adjustments by the i

^th

sensor. In matrix notation, by stacking the x

i

:s into a column vector, we obtain

x

( t ) =

^a

(

¹

) :::

^a

(

d

) ] s ⁽ ^t ^{) +} n ⁽ ^t ^{) =}

^A

⁽

⁰

⁾ s ⁽ ^t ^{) +} n ⁽ ^t ⁾ ⁽²⁾

where

⁰

is a d -dimensional parameter vector corresponding to the true signal parameters. The array response vector

^a

(

j

) = a

¹

(

j

) ::: a

m

(

j

) ]

^T

contains the sensor responses to a unit wavefront from the direction

j

.

The additive noise, n ( t ), is modeled as a zero-mean stationary, temporally and spatially white Gaussian random process, i.e.,

E n ⁽ ^t ⁾ n

⁽ ^s ^{)] =}

²^Im

ts

: (3) The symbol (

)

denotes the Hermitean transpose and

ts

is the Kronecker delta. The noise is uncorrelated with the signal waveforms, which are modeled as zero-mean stationary, Gaussian random processes with covariance

E s ⁽ ^t ⁾ s

⁽ ^s ^{)] =}

^ts^S

^: ⁽⁴⁾

The covariance matrix of the array output is then given by

R

xx

= E x ⁽ ^t ⁾ x

⁽ ^t ^{)] =}

^A

⁽

⁰

⁾

^SA

⁽

⁰

^{) +}

²^I^m

^: ⁽⁵⁾

The ESP data consists of N samples of the array outputs,

^x

(1) :::

^x

( N ). If we introduce the m

r matrix

^T

, d < r

m , and the mapping

^x ^7! ^z

=

^T^x

, from

5

(6)

ESP to RDBS, a new set of observations result, consisting of the r -dimensional vectors

fz

( t )

^g^N_t⁼¹

. Then, in RDBS, the representation of the array observations will be

z

( t ) =

^T^A

(

⁰

) s ( t ) +

^T

n ( t ) : (6) The covariance matrix of the RDBS sensor noise process equals

^T

²^Im^T

=

²^T^T

. It will be required that

T

=

^Ir

(7)

implying that the RDBS sensor noise is spatially white whenever the ESP noise is so.

Furthermore, it is assumed that the array manifold has no ambiguites, i.e., any matrix

^a

(

¹

) :::

^a

(

m

)] has full rank, m , for distinct

i

:s. If d

⁰

denotes the rank of the signal covariance matrix,

^S

, unique estimates of the DOAs are ensured by also assuming that 2 d < r

^;

d

⁰

, 10]. We also regard the parametrization of the array response vector,

^a

( ), as being known.

6

(7)

3 RDBS Cramer-Rao Bound

The previously proposed transformation matrix design methods have mainly been focused on obtaining as high as possible RDBS SNR or as low as possible resolution threshold in RDBS. Apart from the results (improved resolution threshold in RDBS for the MUSIC algorithm) in 7], there have not been any real guidelines for how much the ESP dimension actually can be reduced without loss in performance measured by, for example, estimation accuracy. Furthermore, there have not been any discussion if there is a possibility of achieving the ESP CRB, denoted CRB

^ESP

, which is the best one could hope for, by using the RDBS data set.

Since there exist estimation methods, e.g. the WSF method, 11]-12], that indeed attains the CRB

^ESP

, the same method can be used in RDBS to reach the CRB

^RDBS

(

^T

).

Therefore, our aim in this section is to show that there exists a transformation matrix that makes CRB

^RDBS

(

^T

) = CRB

^ESP

. Of course, the expressions for the CRB's are to be evaluated at

⁰

, the true DOAs. However, the result is still useful for deriving practical design methods, as demonstrated in the next section.

The emitter waveforms are modeled as zero-mean, stationary Gaussian random processes in Section 2, which means we should consider the expression for the stochastic

¹

Cramer-Rao bound. This bound is under the previously stated assumptions as follows,

13].

Theorem 1

E

^

^;⁰

^

^;⁰^T

CRB

^ESP

(8) where

CRB

^ESP

=

²

2 N

Re

^D^P^?^A^D^SA^R^;1^xx^AS^T^;1

(9) Here, the notation

^A ^B

denotes the Schur product of the matrices

^A

and

^B

, i.e., (

^A^B

)

ij

=

^Aij^Bij

, and where

P

?

A

=

^I ^;^A

(

^A^A

)

^;1^A

(10)

1OrunconditionalCRB.

7

(8)

is the orthogonal projection onto the null space of

^A

, and

D

=

^d

(

¹

) :::

^d

(

d

) ] (11)

d

(

j

) = @

@

j^a

(

j

) j = 1 ::: d : (12) The expression for CRB

^RDBS

(

^T

) is obtained by substitution of

^A

,

^D

and

^R^xx

by

^T^A

,

T

D

and

^R^zz

=

^T^R^xx^T

=

^T^ASA^T

+

²^Ir

, respectively. The design method to be presented later is based on the following theorem, which is the main result of this paper.

First, dene

^<^fM^g

to be the range space of the matrix

^M

.

Theorem 2 Assume that r

2 d and that

^T

is an m

r matrix satisfying

^T^T

=

^Ir

. Then, if

<f

^A

(

⁰

)

^D

(

⁰

)]

^g^<^fT^g

(13) the relation CRB

^RDBS

(

^T

⁰

) = CRB

^ESP

(

⁰

) holds true. The symbol \

" above should be read \is contained in".

Proof

Introduce rst the notation

F

=

^D^T ^Ir^;^T^A

(

^A^T^T^A

)

^;1^A^T^T^D

(14) and

G

=

^SA^T ^T^ASA^T

+

²^Ir

;1

T

AS

: (15)

Here, the dependence on

⁰

for

^A

and

^D

has been suppressed for notational convenience.

Then the expression for the CRB in RDBS reads CRB

^RDBS

=

²

2 N

h

Re

^F ^G^Tⁱ^;1

(16) Starting with the expression for

^F

we have

F

=

^D^T ^Ir^;^T^A

(

^A^T^T^A

)

^;1^A^T^T^D

=

^D^T ^T^T ^;^T^A

(

^A^T^T^A

)

^;1^A^T^T^D

=

^D^T^T^Im^;^A

(

^A^T^T^A

)

^;1^A^T^T^D

(17)

8

(9)

where the requirement that

^T^T

=

^Ir

is utilized. Now, making use of the assumption (13), there exists a

^T

such that the relations

TT

A

=

^A

TT

D

=

^D

(18)

are satised. Clearly, this means that

^F

=

^D^P^?^A^D

, which is the rst part of the CRB expression in ESP. For the second part, we notice rst that

T

ASA

T

+

²^Ir

;1

=

^;2^Ir^;

^;2^T^AS

^;2^A^T^T^AS

+

^Id

;1

A

T

^;2

(19) by making use of the matrix inversion lemma. Inserting (19) into (15) gives

G

=

^S^A^T^T^AS

^;2^;

^;2^S^A^T^T^AS

^;2^A^T^T^AS

+

^Id

;1

A

TT

AS

^;2

=

^S^A^AS

^;2^;

^;2^S^A^AS

^;2^A^AS

+

^Id

;1

A

AS

^;2

=

^S^A^ASA

+

²^Im

;1

AS

=

^S^A^R^;1^xx^AS

(20)

where the relation

^T^T^A

=

^A

was used to obtain the second equality and the matrix inversion lemma applied to get the third one. This completes the proof.

²

A matrix

^T

of size m

2 d that makes the relations (13) and (18) satised can be obtained by using the Singular-Value-Decomposition (SVD), 1]. Let

^A^D

] =

^U^V

(21)

be the SVD of the matrix

^A^D

], where

^U

is m

2 d ,

= diag

¹

:::

²d

]

(

¹

²

:::

²d

0), and

^V

is 2 d

2 d . Since

^U

spans

²

the range space of

^A ^D

] and

^U^U

=

^I²d

, one may take

^T

=

^U

. A computationally slightly more ecient way to calculate

^T

is to use Gram-Schmidt orthogonalization on the matrix

^A ^D

]. An advantage with the SVD approach is, however, that it suggests a natural way to reduce the dimension further simply by selecting only the r ( d < r < 2 d ), rst columns of

^A^D

]. This will be motivated in Section 4.

2The matrixÛÛ is in fact a projection matrix, projecting onto the range space of Â^D].

9

(10)

A further note is that we can compute a matrix

^T

=

^U

, as in (21) above, using a number r > 2 d of \design DOAs"

and still be able to optimally estimate the d true source locations. That is, provided that the true DOAs

⁰

are included in

, since we then have

^T

(

)

^T

(

)

^A

(

⁰

)

^D

(

⁰

)] =

^A

(

⁰

)

^D

(

⁰

)]. Moreover, if the array response vectors,

^a

( ), are smooth functions of it is likely that also -values close to the design DOAs can be \near optimally" estimated. The design method presented in Section 4 is based on this observation.

A corollary to the theorem is that also the ESP deterministic

³

Cramer-Rao bound,

14], is preserved in RDBS. The underlying assumption is in this case that the signals are assumed to be non-random, i.e. to be the same in all realizations. This is natural for the problem of estimating the frequencies of cisoids-in-noise, for instance. The deterministic CRB can be attained only if the number m is large. This is due to the fact that the number of parameters (more signal waveform parameters to estimate) grows as more snapshots, N , is taken. The stochastic CRB, however, is achieved by the WSF method if only the number of snapshots is large. The expression for the deterministic bound is

E

^

^;⁰

^

^;⁰^T

CRB

^ESP^DET

=

²

2 N

Re

^D^P^?^A^D^S

^

^T^;1

(22) where

^

S

= 1 N

N

X

t⁼¹^s

( t )

^s

( t ) (23)

and, as readily can be seen, (13) will guarantee that also in this case we will have CRB

^ESP^DET

= CRB

^RDBS^DET

(

^T

).

Since the problem of estimating the frequencies of complex sinusoids in noise can be formulated as a large- m (ULA and N = 1) problem, there is a possibility of performing the estimation in RDBS in an ecient manner by using the deterministic ML method in RDBS ecient here alludes to estimation accuracy as well as to computational requirements. In 15], the consistency and asymptotic (for m

^! ¹

) eciency of the ML estimates for sine wave parameters was established for the single-experiment case ( N = 1). A numerical example of this application is given in Section 6.

3Sometimes referred to asconditionalCRB.

10

(11)

4 Transformation Matrix Design

In this section, a design method is outlined for RDBS estimation of DOAs in a prespecied spatial sector. The performance of the proposed method is evaluated in Section 6.

Design Steps.

1. Specify rst an interval,

a

b

], within which sources of interest are assumed to be located.

2. Compute the beamwidth

⁴

,

BW

, of the array, and construct the \design vector"

of DOAs as

=

a

+

BW

a

+ 2

BW

:::

b^;

BW

b

]

^T

. Let the number of design DOAs be denoted r= 2, assuming that the sector is chosen in such a way that r is even. Then r = 2(

b ^;

a

) =

^BW

.

3. Compute the SVD of

^A

(

)

^D

(

)] as

^U^V

=

^A

(

)

^D

(

)].

4. Take as the transformation matrix (of dimension m

r )

^T

=

^U

.

5. If possible, reduce the dimension of the transformation matrix by evaluating the performance measure (see (28) and the discussion below) over the spatial interval of interest.

Remarks: Note that the specication of the spatial sector in step 1 above also can include several separated sectors. Moreover, if sectors are specied beforehand, the corresponding

^T

's can be precomputed and stored (or implemented in analog hardware).

But, if the sector specication is made by rst nding the peak locations of the spatial spectrum of a beamformer output, the sector choices will be data dependent. Hence, in this case, the transformation matrix must be computed \on-line". Furthermore, the spacing between the design DOAs in step 2 are chosen somewhat arbitrarily. For obtaining enough sector coverage at the end points of the interesting interval it seems necessary

4The beamwidth of a ULA with half-wavelength element spacing is approximately 2=mradians.

11

(12)

to have a slightly smaller spacing there. To examine if a particular spacing is acceptable one can, preferably, evaluate the performance measure (28).

It is seen from Theorem 2 that if the matrix

^T^T

^A

(

⁰

)

^D

(

⁰

)] approximates the matrix

^A

(

⁰

)

^D

(

⁰

)] well enough, we should be very close to asymptotically attain the ESP CRB also in RDBS. The design procedure guarantees that, due to the

BW

- distance between the \design DOAs", the approximation is very good over the whole interval. However, if the spatial sector is large compared to the beamwidth, the RDBS dimension may be unacceptably high. To further reduce the dimension, it is proposed to use only the k \largest" left singular vectors of

^A

(

⁰

)

^D

(

⁰

)]. One motivation for this is that

^Uk

=

^U

(: 1 : k ), the notation meaning the k rst columns of

^U

, solves the following problem:

Uk

= arg min

_rank

(T)=k^k^T^T

^A

(

)

^D

(

)]

^;

^A

(

)

^D

(

)]

^k²

(24) This can be seen as follows. It is well-known that

Bk

= arg min

_rank

(B)=k^k^B^;

^A

(

)

^D

(

)]

^k²

(25) where

^Bk

=

^P_ki⁼¹

i

u

i

v

_i

is the expression for the best rank- k approximation of the matrix

^A

(

)

^D

(

)]. Here, u

i

and v

i

are the left and right singular vectors of

^A

(

)

^D

(

)], corresponding to the i

^th

largest singular value, respectively. The optimal criterion value is

k⁺¹

. If we take

^T

=

^Uk

and use the fact that

UV

=

^X²^d

i⁼¹

i

u

i

v

_i

=

^A

(

)

^D

(

)] (26) we obtain

TT

^A

(

)

^D

(

)] =

^X^k

i⁼¹

i

u

i

v

_i

=

^Bk

: (27)

Now, if the ( k +1)

^th

singular value is small enough compared with the k larger ones, it should be possible to delete r

^;

k columns in the optimal m

r matrix

^T

, and thus reduce the RDBS dimension, hopefully without too much loss in the RDBS optimal performance.

As a measure of closeness between the ESP and RDBS performance over a spatial interval, we take the following:

( k ) =

R^b

^a

CRB

^ESP

( ) d

R^b

^a

CRB

^RDBS

(

^Tk

) d

!

1=²

k = d + 1 ::: r : (28)

12

(13)

Here, the index k of

^Tk

indicates the RDBS dimension. Then, 0

( k )

1, with the upper bound beeing desirable to attain. By depicting the values of ( k ) obtained for dierent k 's, it is easy to see for which dimensions the performance measure is suciently close to 1. To make the measure easy to compute, it can be approximated by taking sums over closely spaced 's instead of using integrals. A further note is that this criterion only considers a single emitter traversing the interval of interest. This is, however, no limitation here since we are only interested in making sure that the transformation matrix

\covers" the whole interval, which will be guaranteed if the index is close enough to 1.

13

(14)

5 Source localization using WSF

The WSF method, 11]-12], is a member of the same general class of subspace t- ting algorithms as the deterministic ML method. In this class, WSF is optimal in the sense that it asymptotically achieves the lowest estimation error variance for large N . Assuming that the emitter waveforms are generated according to a Gaussian distribution, WSF then gives the lowest asymptotic estimation error variance of any unbiased estimator, i.e., it asymptotically attains the stochastic CRB, (9). In the following some of the details of the method are reviewed.

The covariance matrix of the array output is given by equation (5). Let d

⁰

be the rank of

^S

. Then, clearly, we have d

⁰

d . Let

R

xx

=

^X^m

i⁼¹

i^ei^ei

=

^Ess^Es

+

^Enn^En

(29) be the eigendecomposition of

^R^xx

, where

¹

²

:::

m

. The matrix

^Es

contains the d

⁰

eigenvectors of

^R^xx

corresponding to the largest eigenvalues, i.e.,

^Es

=

^e¹

:::

^ed⁰

]. The range space of

^Es

is called the signal subspace. It's orthogonal com- plement is the noise subspace, which is spanned by the columns of

^En

=

^ed⁰⁺¹

:::

^em

].

Furthermore, the smallest eigenvalue of

^R^xx

has multiplicity m

^;

d

⁰

and is equal to the noise variance. Thus,

n

=

²^I

. It is also well-known that the signal subspace is a subset of

^<^fA

(

⁰

)

^g

, i.e.,

^<^fEs^g ^<^fA

(

⁰

)

^g

with equality if and only if the signals are non-coherent, i.e., when d

⁰

= d .

The eigendecomposition of the sample covariance matrix, ^

^R^xx

, is dened in a similar fashion as (29),

^

R

xx

= 1 N

N

X

t⁼¹^x

( t )

^x

( t ) = ^

^Es

^

_s^E

^

_s

+ ^

^En

^

_n^E

^

_n

: (30) The WSF criterion is stated as follows. Find the best least squares t of the two subspaces

^

Es

and

^A

(

), i.e.,

^

= arg min

Q

k^E

^

_s^W¹⁼²^;^AQ^k²_F

= argmax

Tr

^h^P^A^E

^

_s^W^E

^

_sⁱ

(31) Here,

^W

is a d

⁰

d

⁰

Hermitean positive denite weighting matrix for which the optimal choice is

^W

= (

s^;

²^I

)

²^;1_s

.

^P^A

is a projection matrix, projecting onto the range space of

^A

, i.e.,

^P^A

=

^A

(

^A^A

)

^;1^A

.

14

(15)

Techniques for numerically calculating the estimates and detecting the number of emitters present can be found in in 16]. The location estimation of sources in RDBS is easily accomplished using this estimation algorithm, simply by substitution of the array response vector: from

^a

( ) to

^T^a

( ).

15

(16)

6 Numerical Examples

In the examples to follow, the employed WSF-algorithm is supplied with a detection step, based on the asymptotic

²

-distribution of the normalized WSF cost function, which estimates the number of emitters present in the actual scenario, see 16, 17]. The only a priori knowledge used for the detection and DOA estimation is the functional form of the array response vector and the spatial sector.

In the rst three examples, the noise power will be

²

= 1 and the emitter waveforms are Gaussian distributed. Furthermore, a uniform linear array with half-wavelength element spacing is employed. The estimated root-mean-square (RMS) error curves depicted in the gures are based on 200 independent trials for each DOA considered and N = 20 independent snapshots are used in each trial. The last example consists of an application of the results derived herein to the problem of estimating the frequencies of complex sinusoids from noisy data.

Example 6.1

This rst example demonstrates the design procedure by means of a simple scenario. The sector of interest is assumed to be

^;

3 : 5 3 : 5 ] and the scenario constitutes an m = 100 element ULA as well as an emitter signal for which we have SNR=3 dB. Following the design steps of Section 4, we compute the array beamwidth to approximately 1.14 degrees. Take

BW

= 1 : 5 , which means we can choose the vector of design DOAs as

=

^;

3

^;

1 : 5 0 1 : 5 3 ]

^T

. This choice of

BW

is slightly larger than the beamwidth, but will turn out to be small enough anyway. Then we obtain, through the SVD of

^A

(

)

^D

(

)], a transformation matrix of dimension 100

10. To see if it is possible to make further dimension reduction and still attain the ESP CRB over the interval of interest, we take a look at the (approximation of the) performance measure ( k ), see (28), computed over the interval

^;

3 : 5 3 : 5 ]. In Figure 1 it can readily be seen that it should be sucient to have a dimension ranging from about seven to nine for the RDBS. This means that we can delete at least the last column in

^T

, corresponding to the smallest singular value. We select here the RDBS dimension to be nine. (The eect

16

(17)

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

3 4 5 6 7 8 9 10

RDBS dimension

Performance Index

Figure 1: Evaluation of the performance measure ( k ), (28), over the interval

^;

3 : 5 3 : 5 ] for the transformation matrix obtained from design DOAs

^;

3

^;

1 : 5 0 1 : 5 3 ]

^T

. Scenario: ULA, SNR=3 dB, m = 100 N = 20.

of removing one column of the transformation matrix is mainly that the side-lobes of the beampattern are reduced.) The agreement between the CRB's is demonstrated by Figure 2, where the CRB's are depicted over the interval of interest for a RDBS dimension of nine. The dierence between the ESP and RDBS CRB's in Figure 2 is very small and, hence, we can regard a dimension of nine as sucient. The corresponding curves for RDBS dimension 10 is not shown since they, for all practical purposes, are identical to those of Figure 2. However, at the end points of the interval of interest, the RDBS CRB for dimension 10 is a little closer to the ESP CRB than what is the dimension-nine RDBS CRB. Further reduction of the RDBS dimension makes this dierence more noticeable.

Finally, Figure 2 also depicts the RMS-error in RDBS and the CRB's versus direction- of-arrival for one single source for an RDBS dimension of nine. Clearly, for DOAs in the interval

^;

3 : 5 3 : 5 ], the best possible performance (ESP CRB) can be considered

reached.

²

17

(18)

0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02

-5 -4 -3 -2 -1 0 1 2 3 4 5

CRB’s and RMS error (degrees)

Direction-of-arrival

Dimension of RDBS=9 CRB (ESP):

CRB (RDBS):

RMS WSF (RDBS):

Figure 2: CRB's and RMS error of DOA estimate vs. DOA of emitter using the transformation matrix obtained from design-DOAs

^;

3

^;

1 : 5 0 1 : 5 3 ]

^T

. Scenario: ULA, SNR=3 dB, m = 100 N = 20.

Example 6.2

The scenario for this example is the same as for Example 6.1, but here we consider two sources, symmetrically located with respect to array broadside. Figure 3 depicts the ESP and RDBS CRB's and the RDBS RMS error for the estimate of the direction to the left source as the angle between the emitters decreases. That is, we are interested in comparing the experimental estimation accuracy (RMS) for two sources in the interval- of-interest with the optimal ESP CRB. Note that the \Direction-of-Arrival axis" in Figure 3 is half the separation between the emitters. We can see that, for this case, the ESP and RDBS CRB-curves are almost identical in the interval of interest. Moreover, the lower bound on the estimation error accuracy is nearly attained, as is demonstrated

by the depicted RMS-error curve.

²

Example 6.3

18

(19)

0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 0.055

-5 -4.5 -4 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0

RMS error (deg)

Direction-of-arrival

CRB (RDBS):

RMS WSF (RDBS):

Figure 3: RMS error of DOA estimate vs. emitter DOA of left source using the transformation matrix obtained from design-DOAs

^;

3

^;

1 : 5 0 1 : 5 3 ]

^T

. Scenario: ULA, Symmetrically placed emitters, m = 100 N = 20.

This example examines what happens when a source is located close to \the edge" of, or outside, the sector-of-interest. The reason is that it is essential to know what this means for the estimation accuracy of a source inside the sector. The interval of interest is assumed to be

^;

2 : 5 2 : 5 ] and the transformation matrix is designed by using the design DOAs

^;

1 : 5 0 1 : 5 ], giving an RDBS dimension of six. An emitter within the interval of interest is placed at

^;

0 : 2 and the theoretical as well as the experimental estimation accuracy of this DOA are examined, in Figure 4, as a function of another emitter whose DOA varies from 2 to 6 . The DOA of this second emitter is referred to as the out-of- sector DOA. As depicted in Figure 4, the experimental RMS-curve for the

^;

0 : 2 emitter is very close to the optimal CRB as long as the second emitter is located within the sector-of-interest. This behavior is expected from the theory. However, when the second emitter slides out of the sector-of-interest, the estimation accuracy of the

^;

0 : 2 emitter rapidly degrades. The performance is especially bad for out-of-sector DOAs at which the RDBS CRB and the ESP CRB for the

^;

0 : 2 emitter dier. The conclusion is then

19

(20)

0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 0.022 0.024

2 2.5 3 3.5 4 4.5 5 5.5 6

RMS error (deg)

Out-of-Sector DOA

CRB (RDBS):

RMS WSF (RDBS):

Figure 4: RMS error of

^;

0 : 2 DOA estimate vs. out-of-sector DOA using the the transformation matrix obtained from design-DOAs

^;

1 : 5 0 1 : 5 ]

^T

. Scenario: ULA, SNR=3 dB, One out-of-sector emitter, m = 100 N = 10.

that a study is needed that centers on the issues of reducing the sensitivity to out-of- sector emitters while still trying to retain \as much as possible" of the (ESP) optimal

estimation accuracy.

²

Example 6.4

This last example is a direct application of the results derived herein to the problem of estimating the frequencies of complex sinusoidal signals (cisoids) from noisy data.

Assume that the collected data points have been generated as y ( t ) =

^X^d

k⁼¹

k

e

ⁱ⁽^!^k^t⁺^'^k⁾

+ n ( t ) t = 0 ::: m

^;

1 (32) where the noise process,

^f

n ( t )

^g

, is a temporally white complex Gaussian, 18], process with variance

²

= 1 and the '

k

:s are constants with values in the interval 0 2 ) radians. The amplitudes,

k

, are real-valued. The !

k

:s are all assumed to be distinct

20

(21)

0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016

-10 -5 0 5 10 15 20

RMS error (rad)

SNR Dimension of RDBS=6

CRLB (ESP) For both emitters:

RMS WSF (RDBS) Emitter @0.06 rad:

Figure 5: RMS error of cisoid frequency estimates vs. SNR.

Scenario: 100 observations, i.e., m = 100 N = 1.

from one another and satisfy 0 < !

k

< . The latter requirement is due to the sampling theorem. We then can write

x

=

2

6

4

y (0) ...

y ( m

^;

1)

3

7

5

=

^A

(

^!

)

^s

+

ⁿ

(33)

where

!

= !

¹

::: !

d

]

^T

(34)

s

=

^h

¹

e

^i'¹

:::

d

e

^i'^dⁱ^T

(35)

A

(

^!

) = a ( !

¹

) ::: a ( !

d

)] (36) a ( ! ) =

^h

1 e

^i!

::: e

^i!⁽^m^;1)ⁱ^T

(37)

n

= n (0) ::: n ( m

^;

1)]

^T

: (38) From (37), the array response vector can be seen to have the same form as that of a ULA of dimension m , and from (38),

ⁿ

satises (3). Since the emitted waveforms,

^s

, is

21

(22)

0 0.005 0.01 0.015 0.02 0.025

40 60 80 100 120 140 160

RMS error (rad)

Dimension of RDBS Dimension of RDBS=6

CRLB (ESP) For both emitters:

Figure 6: RMS error of cisoid frequency estimates vs. number of data points.

Scenario: SNR = 3 dB, N = 1.

a sequence of numbers, the deterministic CRB, (22), is applicable. The ML-estimator is as follows, 14].

^

!

= argmax

_!

Tr

^h^P^A

(

^!

) ^

^Rxx

i

(39)

Here, ^

^Rxx

is the sample covariance matrix, (30). The criterion (39) is only a slight modication of the WSF-criterion, which means that the method employed earlier can be used also for calculating the deterministic ML-estimates. The number of cisoids will in this example be assumed known. When this is not the case, the number of signals can be estimated from the set of data, see e.g. 19]. In Figure 5, the RMS error curves for the frequency estimates of two cisoids, at 0.06 and 0.15 rad., respectively, are depicted versus SNR here dened as

_k²

=

²

. The RMS error curves are based on 200 independent trials. The number of data points is m = 100 and the RDBS transformation matrix is computed to cover the frequency interval ranging from 0.03 to 0.33 radians, which makes an RDBS dimension of 6 necessary. Note that this means that the ML-estimation of the cisoid frequencies is based on only six (complex-valued) RDBS observations. However, as expected from the theory, the CRB is nearly attained for the frequencies at hand,

22

(23)

as demonstrated in Figure 5. Figure 6 depicts the RMS error for the same frequencies as above, but here versus the number of ESP data points, m . The SNR is in this case

xed to 3 dB. The RDBS dimension is 6, the same for all m . Note that this means that the width of the beampattern of the transformation matrix decreases as the number of data points grows when the \design frequencies" are constant for all m . This is because the beamwidth of the ULA decreases with increasing m . The frequency sector varies, approximately, from 0 0 : 36] for m = 40 to 0 : 06 0 : 29] for m = 150. As can be seen, the estimation accuracy tends to the CRB as the ESP dimension grows.

²

23

(24)

7 Conclusions

A high dimension of the RDBS is not necessary to obtain accurate estimates of the source localizations. As demonstrated here, it is possible to make a signicant reduction of the ESP-dimension without loosing any accuracy at all. There also exists a bound for the RDBS dimension above which it is always possible to retain the property of optimally estimating the directions-of-arrival. This bound was shown to be 2 d , i.e., twice the true number of sources present. Whether or not this is the lowest possible bound is, however, still an open question. Moreover, a practical design strategy for a transformation matrix was devised that achieve near optimal performance, in the sense of ESP CRB. The optimal transformation matrix is easily obtained from a singular value decomposition, assuming knowledge only of the functional form of the array response.

This optimal transformation is thus independent of emitter signal correlations as well as of the powers of the sources. Furthermore, any ecient estimation method, i.e. any method that attains the stochastic (or deterministic) Cramer-Rao bound, can utilize this transformation matrix for obtaining optimal dimension reduction. The problem of estimating the frequencies of complex sinusoids was also addressed. This problem can be formulated as one with high ESP dimension, equal to the number of collected data points, which makes the results derived herein applicable.

The possibility of this large reduction in data dimension may seem to be a bit contra- dictive at a rst glance, since all the information present in ESP cannot be present also in RDBS. However, the criterion considered here is mainly to retain the ESP estimation accuracy in a spatial sector. This ensures that the informative part of the reduced data set is as large as possible, i.e., the information contained in the selected spatial sector in ESP is preserved. In addition to this, the sector criterion also de-emphasizes the requirement to know the true directions-of-arrival in order to make the RDBS inherit the optimal estimation-accuracy properties from ESP.

24

(25)

References

1] G.H. Golub and C.F. Van Loan, Matrix Computations, 2 nd edition, Johns Hopkins University Press, Baltimore, MD., 1989.

2] P. Forster and G. Vezzosi, \Application of Spheroidal Sequences to Array Process- ing", ICASSP 87, pages 2268{2271, April 1987.

3] K. Buckley and X.L. Xu, \Spatial-Spectrum Estimation in a Location Sector", IEEE Trans. ASSP, pages 1842{1852, Nov. 1988.

4] B.D. Van Veen and B. Williams, \Structured Covariance Matrices and Dimen- sionality Reduction in Array Processing", Proc. 4:th Workshop Spectr. Est. and Modeling, pages 168{171, Aug. 1988.

5] B.D. Van Veen and B. Williams, \Dimensionality Reduction in High Resolution Direction of Arrival Estimation", Asilomar Conf., Pacic Grove, California, Oct.

1988.

6] B.D. Van Veen and R.A. Roberts, \Partially Adaptive Beamformer Design Via Output Power Minimization", IEEE Trans. ASSP, ASSP-35 (11):1524{1532, Nov.

1987.

7] H.B. Lee and M.S. Wengrovitz, \Resolution Threshold of Beamspace MUSIC For Two Closely Spaced Emitters.", IEEE Trans. ASSP, ASSP-38 (9):1545{1559, Sept.

1990.

8] R. O. Schmidt, A Signal Subspace Approach to Multiple Emitter Location and Spectral Estimation, PhD thesis, Stanford Univ., Nov. 1981.

9] X.L Xu and K.M. Buckley, \Statistical Performance Comparison of MUSIC in Element-Space and Beam-Space", Proc. ICASSP 89 Glasgow, pages 2124{2127, May 1989.

10] M. Wax and I. Ziskind, \On Unique Localization of Multiple Sources by Passive Sensor Arrays", IEEE Trans. on ASSP, ASSP-37 , July 1989.

25

On Optimal Dimension Reduction for Sensor Array Signal Processing