Wide-band maximum likelihood direction finding by using tree-structured EM algorithm

Wide-Band Maximum Likelihood Direction Finding

by

Using

Tree-Structured EM Algorithm

Nail

Cadalli

and Orhan Arikan

Department of Electrical and Electronics Engineering, Bilkent University,

TR.-06533

Ankara,

TURKEY.

Phone:

90-3 12-26643

07,

Fax:

90-3 12-2664

126,

e-mail: cadalliaee

.

bilkent

.

edu.

tr

,

oarikan@ee

.

bilkent

.

edu.

tr

Abstract

A t r e e sti-uctiired E x p e c t a t i o n M a x i m i z a t i o n ( E M )

~ l ~ g o r L t h l l 7 ~ ‘ i s p r o p o s e d a n d applied t o t h e wide-band a n -

gle of ur,rzvul e s t i n a a t i o n .

It

naay be s e e n a s a g e n - r.ruliznt~ion O I L EM u s i n g the i d e a s of C a s c a d e

E M

al-

q o r i t l ~ m an d S p a c e A l t e r n a t i n g G e n e r a l i z e d EM algo- i~ithna. Also, for p a s s i v e d a t a a c q u i s i t i o n , r o b u s t a n d e f i c z e n t a l t e r n a t i v e s for t h e e s t i m a t i o n of t h e s o u r c e .sr,qnals a r e i n iiestigated.

1.

Introduction

Tn many da.ta. acquisition systems, reception d a t a acquired by a.n array of sensors are processed t o ob- t,aiii informa.tion a b o ut t h e source locations. When t,he sources are located relatively far away from the sensors. only t,he direction of arriva.ls of th e acquired source signals c,an be reliably obtained. Although t he Maxiriiuiii likelihood (ML) estimation provides more accurate estimates for t h e direction of arrivals, due t o t,lie higher computational cost of obtaining the ML esti- mates, i t has not found much use in practice. However,

by exploiting t h e superposition property of th e d a t a x q u i s i t i on system, th e complexity of the ML estima- t,iou ca.n be grea.tly reduced by using th e Expectation R~Iasimization ( EM ) algorithm [11 2. 51. In EM for- rnalism, the observation, i n c o m p l e t e d a t a is obtained

via a. many-to-one mapping from t he c o m p l e t e d a t a

spa.ce t,lia.t includes signals which we would obtain as

the sensor o ut p ut s if we were able to observe the effect

of‘ each source separately. T h e EM algorithm iterates between estimating t h e log-likelihood of t he complete data. using t h e incomplete d a t a an d th e current param- eter est8iriia,tes (E-step) a n d maximizing t he estimated log-likelihood function to obtain t he updated parame- ter est>iina.tes (R/I-step). Under mild regularity condi-

tions, th e iterations of t he EM algorithm converges t o

a stationary point of t he observed log-likelihood func- tion, where a t each iteration th e likelihood of tlie es- timated parameters is increased [ll]. In this study, a

tree structured hierarchy is used for the description of relation between t he complete d a t a space a n d the ob- servations. Within this hierarchy it is possible t o com- bine in one algorithm t he ideas of th e Cascade EM a nd Space Alternating Generalized EM algorithms [ 3 , 71.

For the estimation of unknown signals arriving from different directions t o a passive array, alternative regu- larized estimation schemes t o t he common least squares solution are investigated. For this purpose two differ- ent methods are used. T h e first one is an adaptive Tikhonov type regularized least-squares

(RLS)

estima- tion method, which is coinputationally intensive a.nd th e second one is a n averaged least-squares estimation (LSSET) method over a set of angles in a neighbor- hood of the nominal angles. I t has been demonstrated t h a t when RLS or LSSET methods are used in tlie es- timation of the received signals, t he EM algorithm has better convergence behavior.

2.

Signal Model

For t he case of M sources with direction of arrivals

81, 1

5

1

5

M , t he measured signal at t h e i’th sensor of an array with P sensors is

M

Yi(t) = x a i ( t , 01)

*

S i ( t - Ti(O1))

+

% ( t )

1=1

1

5

i

5

P

,

t=O,T,2T

,...,

(N-1)T (1)

where s , ( t ) is t h e wide-band signal of t h e I’th source, u i ( t ) is t he 0 mean spatially and temporally white Gaussian noise at th e

i’th

sensor, ri(0) is t h e time delay of the source signal from th e direction B as i t propa- gates to the i’th sensor relative to t8he phase center of

554

thr c?.rray. n i ( t , 8) is the trime domain function for the

p i l i of' the i't81i sensor which is dependent on frequency

i>nd the directmion of a.rriva.1, 8. T h e frequency domain ~~rl)resrnt~ation of (1) is,

where F is the DFT size which is chosen sufficiently la.rge and Y i ( t ) , & ( k , 8 ) , S l ( k ) ve U i ( k ) are the ti:a.ns- Corms of y j ( t ) , a i ( $ ,

e ) ,

s j ( t ) and u ; ( t ) respectively. Let the following definitions be made (' is the transpose 0pera.tor) rl; ~! b ( R : , O ) = [ & ( k , Q ) e - - T " ' 2-k r P ( @ ) & ( k ,

q e - j F T . ] "

= \ b I ( k , O ) .

.

.bp(k,O)]' B(k, 0 ) = [ b ( k : , 8,) . b ( k ,

en,)]

S ( k ) = [ S l ( t ) . '

.

S & f ( k ) ] ' Y ( k ) = [ Y 1 ( k ) " . k p ( k ) ] '

[Tsing these ckfinit,ions (2) becomes

Y(X-)

=

B ( t , O ) S ( k )

+

U ( k ) 0

5

k

< F

_{( 3 )}

This final coiiipa.ct forin of the measurement relation,

whicb is t,lre same as the signal model of the Cramer- R.a.o Lower Bound formula in [SI, will be used in our deriva.tions.

3.

Wide-Band

EM

Algorithcn

Since, the measurement noise if; modeled as nor-

mally tlistribnt,ed additive noise, the ]probability clen- sity of' the observations are Gaussia>n. Hence, the log-

liltelihood fiinction of tthe observations has the following f'aiiiiliar form ( t is the conjugate trainspose operator)

,

F - 1

L ( 0 ,

S ; Y ) =

-

[ Y ( k ) - B(k, @ ) S ( k ) ] +

k = l l

[ Y ( k ) - B(k, @)S(k)l (4) hi order t,o fincl the ML estimate, likelihood function

of' t8he obaerva,tions should be maximiszed with respect t.o 0 and S ( k ) . However, the direct maximization of t,liis function is not only computationally demanding but, a.lso clue t,o the local inaxima structure of the: I.ike- lihood fuiict,ioii it4 ia not guaranteed to converge t o the glohal maxima.. T h e Especta,t,ion Maximization

(EM)

nietbocl of olit,a.ining the ML estimate overcomes this

difficulty by a.n iterative search in much lower dimen- sional parameter spaces [l]. T h e EM method requires the identiihcation of so ca.lled complete da ta space. In our application the commonly used complete da.ta is

Xl(k) = [ X 1 l ( k ) . . . X p l ( k ) ] ' which is the signal that would be observed a.t the sensors if we were able to see the effect of

I'th

source only. Then the many-to-one ma.pping for all sources from the complete d at a space t o the incomplete d ata space can be written as

A4

Yl(k) = C X l ( k ) O < k < F (5)

f=1

The mean. of the complete d at a Xl(k) is b(H, & ) S j ( k ) and it is normally distributed. T h e log-likelihood func- tion of the complete data is

F - 1 IV

C,(@,

S ;

:X)

= -

C

IlXi(k)

-

b(k,

Q5(k)l12 (6)

k = O I = 1

Here, the observed signal is decomposed t o M con- stituents. Therefore t o estimate 01 and S i ( k ) , only X l ( k ) is used besides the observation. At the n'th it- eration of the EM algorithm expectation step condi- tionally estimates the likelihood of the complete data

L,(O,

S

1

On, Sn,). Maximization step then finds the irmximizer of the estimated likelihood a.nd assigns to

Byt1.

To

find b ( k , B l ) S [ ( k ) it is sufficient t o know

Xl(k),

therefore

in

expectation step

X r ( k )

is estimated.

It

can be shown t h a t ,

([GI,

p. 164),

X ; ( k )

=:

E{Xi(k)l6;2, S p ( k ) , Y ( k ) }

1

=

b ( k ,

ep)sp(k)

+

,,[y(kj - ~ ( k ,

on)sn(k))]

O < k < F (7)

In maximization step complete d at a likelihood which is formed by using Xy(k) is maximized with respect t o

0,

and S,(:k). The

01

update is found as

e;+'

=

, ar gm ax 0

(8) where there is two maximization problems inside one another. If S l ( k ) is unkno~711 they must be simnltane- ously solved. For a given 0 value, the solution of the inner maximization is

Si(11)

= [b(k,

e)bt(k, O)]-'bt(k, Q)Xl(k) (9)

-

b t ( k , W i ( k )

-

Ilbik,

Q)Il2

Iiivrtiiig thiq ~ s p r e s s i o n into ( 8 ) and solving for the oiitci iiicixiniizatioii

@if+'

i> found. For t h a t maximiza-

1 1 0 1 1 liiitv~r search may he used Finally, a t the n ' t h

itei,itioii of tlic EM algorithm the update formulas are

follow% F - l o f ' + i = ai'g inax (r b=O (10) b t ( k . @)Xi" (k)Xltn(k)b(k, 8) p ( k ,

@)l12

If

,SI(

k ) is 1;iiown. as in active array applications, (8) is

lucecl t,o oiie maximization problem and there

wiiiaiiis no iierd for (1 1). If (10) and (11) are run

logetlir,r, i.c. .. in the ca.se of unknown source signals,

o;"''

sliould I W close t o t,rue direction values for

s;+'

t,o c o n ~ w g e t,o t,rur signa.1 waveforms.

i\lt,er ( L O ) , @"+' is available. If i t is inserted into ( 3 ) , S " + ' ( k ) call he solved for by using a number of al- t,clrna.t i\cs. For iiistaiice the least squares

(LS)

solution

is ;IS follows,

S ( k ) = c w y niin IIY(k) - B(k, 0 ) S ( k ) 1 I 2

S(kl

= [Bt(kl @jB(k, 0j]-'Bt(k,0)Y(k)(12)

Rr.girlariza.t,ion may be applied on the LS solution which is called reguhrized least squares,

RLS,

S(X,j = [Bt(k* O ) B ( k , 0 )

+

pI]-'Bt(ll-, O)Y(b) (13) I t is iiiiportant. t.o chose p in the regularization and i t c m I>(, clioseii optimally [4, 91. Another alternative in

sourc(> signal c->st,ima.tsion may be the following which \vi11 hi1 referred to a.s LSSET solution, where

K is

a set

ol angles iii a. neighborhood of 0 ,

ISM algorithm starts with n = 0 a t which time 0' is ;ivailalilP obtaincd by using a rough estimation. To find

Xj'(k) in ( 7 ) ,

,S'p

is needed and it is estimated by one

of (,he niet,hotls iiient,ioiied above. EM shows niono-

i oiiic increase oC t,he likelihood and it!s convergence is-

slit's I i ; \ ~ e 11c:eii investiga.ted [ l , 111.

structured as a binary tree as shown in Figure 1.

Yi,j

(k)

is the intermediate incomplete d a t a between the observation Y(k) and the complete d a t a X t ( k ) ' s .

In this setting EM a.lgorithm is run for two sources a t a time using the intermediate data. a t the joint node of two leaves. This provides a n update for t h e corre- sponding DOA and source signals. T h e value of the intermediate d a t a is found by using, in (3), t h e origi- nal observation Y(k) and the current source signal es- timates other than the ones which are t o he updated by the current run. For instance, t o ruii

EM

algorithm for

X,(k)

and

X,(k)

we form the required incomplete d a t a as

where Yl,l(k} is found by using Y(k) and the cur- rent estimates for the last there source signals in ( 3 ) .

Yz,a(k) may be found similarly and EM algorithm is

run for t h a t branch too. This may be repeated a num- ber of times and then by using the updates obtained for the first 4 source signals and

DOA's,

branch of

Y l , ? ( k ) may be processed. Tile idea of putting inter-

Y1.I

i L index Y i , j

level

Figure 1. An example for the tree structure. mediate d a t a mappings between Y ( k ) and X l ( k ) ' s can be associated with t h a t of t h e Cascade EM, CEM, al- gorithm but here there is more t h a n one intermediate d a t a space. Due t o the limited space, t h e generalization of CEM t o multiple levels is not presented here. T h e tree structure may also be associated with Space Al- ternating Generalized EM algorithm in t h e sense t h a t not all of the parameters are updated at a time. Also EM is run on a more noisy data reducing the informa- tion content of intermediate observations a n d this is reported t o speed the convergence [3].

5 .

Simulations and Conclusions

4 .

Tree-Structured

EM

In tlric. section wcve will use a different mapping from

t h e complete d a t a to the incomplete d a t a which is

Observation signals are obtained by simulation of a

linear array of sensors. T h e number of signals are as- sumed t o be known since there are studies in detection

[IO] The source signals are taken as coherent pulse modulated chirp signals with band width comparable t,o the center frequency. Noise is assumed to be in- tlelwndeiit ideiit(ica1 Gaussian distributed. First the

EM (10-3)

E

Tree-EM

DOA estimation BIIW

I " " " '-7 5.3 5.5 5.5 5.5 6.3 6.1 4.6 2.2

oo121

0 01

1

E I 0.008 - E 0 0 0 6 -

a

0 0.004 - 0 0 0 2 - 01 " ' I

.

" ' 0 2 4 6 0 10 12 14 16 18 20 lleralmn number

T h e DOA error norms for iterations of original EM and tree-structured EM dgorithms are shown in the next table. The original EM algorithm, could not converge t o true DOA values. Furthermore, it diverges from the initial angle values. But within t,lie saint. number of t80tal it e d i o n s the tree-structured EM converges with much lower DOA error to 0 = [ 3 5 . 3 -50.0 -20.0 50.71'.

By this study, an improvement on ER4 algorithm is realized not only by using robust signal estimation schemes but also by changing the d a t a mapping of the original algorithm.

References

Figure 2. Averaged traces of DOA 'error norms.

, -A

- 4 -2 0 2 4 6

SNR

-901

-6

Figure 3. DOA error variances.

signal estimation alternatives are insert*ed in the

EM[

al- gorithm and their relative performances are compared. EM algorit,hm is run for two sources impingingfrorn

35'

a.nd -50' a t a n SNR level of OdB. The averaged traces of error norm of DOA estimation, which describes the c,onvergence behaviours, can be seen in Figure 2 where

EM, LS, RLS, LSSET refers t o ( l l ) , ( l a ) , (13), (1.4:) re- spec,tively. The DOA error variance together with the

CRLB for ea.ch a.lternative is plotted in Figure 3. For

LSSET,

K

consists of 5 a.ngles in a 1' neighborhood of

hhe c,urreiit DOA. This figures out to be computation- ally less complex than

RLS.

To compare the tree-structured Eh4 algorithm. with the original EM algorithm four sources from directions

0 = [35'

-

50' - 20" 50'1' a.re usedl at SNR=lOdB.

Initial DOA's a.re given as 0 0 = [33' --48' - 18' 48'1'.

A. P. Dempster,

N.

M. Laird, a.nd D. B. Rubin. Max-

imum likelihood from incomplete data via the EM al-

gorithm. J . Roy. Stat. Soc., B-39:1-37, 1977.

M. Feder and E. Wkinstein. Pa,rameter estimation of

superimposed signals using the EM algorithm. I E E E

Trans:. Acoust., Speech, Signal Processing, 36(4):477-

489, April 1988.

J. A . Fessler and A. 0 . Hero. Space a.lternat- ing generalized expectation-maxiinizatioii algorithm.

IEEE' Trans. S i g n d Processing, 42( 10):2664-2677,

Oct. 1994.

A. Hoer1 and W. Kennard. Rigde regression: bia.sed

estimation for non-orthogoid problems. Technonset-

M. I. Miller a.nd D. R. Fuhrmann. hlaximam- likelihood narrow-band direction finding and the EM algorithm. I E E E Truns. Acoust., Speech, Siynul Pro-

cessing, 38(9):1560-1577, Sep. 1990.

A. Papoulis. Probability, Random Variables and

Stochastic Processes. McCraw-Hill, 3rd edition, 1991.

M. Segal and E. Weinstein. The cascade em algorithm.

Proc. I E E E , 76(10):1388-1390, Oct. 1988.

P. Stmoica and A. Nehorai. MUSIC, Maximum Likeli- hood and Cramer-Rao Bound. I E E E Trans. Acoust., Speech, Signal Processing, 37(5):720-741, May 1989. H. D. Vinod. Survey of ridge regression anti re-

lated techniques for improvement,s over ordinay least, squares. Review of Econon2ic.s and S'trrtistics, 60:121-

131, 1978.

H. PJang and M . Kaveh. Coherent, signal-suhspa.ce processing for the detect,ioii a.nd estimat,ion of angles

of ar:rival of multiple wide-hand sources. I E E E Trans. Acoust., Speech, Signal Proces.sin,g, ASSP-33(4):823-

831, Aug. 1985.

C. F. J. Wu. On the convergence properties of the EM algorithm. The Arinols of Statistics, 11(1):95-103, 1983.

~ Z C Y , 12:55-68, 1970.