S-100 44 STOCKHOLM. Maximum Likelihood Estimation of Nominal Direction of Arrival and ... in mobile communication applications to increase the capacity of ...
Maximum Likelihood Estimation of Nominal Direction of Arrival and Angular Spread Using an Array of Sensors T~onu Trump and Bjorn Ottersten January 1995 Also in Proceedings of COST 229 Adaptive Systems, Intelligent Approaches, Massively Parallel Computing and Emergent Techniques in Signal Processing and Communications, October, 1994, Vigo, Spain
IR{S3{SB{9501
ROYAL INSTITUTE OF TECHNOLOGY Department of Signals, Sensors & Systems Signal Processing S-100 44 STOCKHOLM
KUNGL TEKNISKA HÖGSKOLAN Institutionen för Signaler, Sensorer & System Signalbehandling 100 44 STOCKHOLM
MAXIMUM LIKELIHOOD ESTIMATION OF NOMINAL DIRECTION OF ARRIVAL AND ANGULAR SPREAD USING AN ARRAY OF SENSORS T~onu Trump and Bjorn Ottersten
Signal Processing, Royal Institute of Technology (KTH), S{100 44 Stockholm, Sweden.
ABSTRACT
The problem of estimating nominal direction of arrival and angular spread of a source surrounded by a large number of local scatterers using a linear array of sensors is addressed. The maximum likelihood (ML) estimator is considered together with an estimator based on least squares t of sample covariance to theoretical covariance matrix derived from a parametric model of the situation. Analysis of the least squares based algorithm is carried out. It is shown that although the least squares estimate is not ecient, it is less computationally demanding than the ML approach and provides good enough estimates to initialize the ML estimator. Finally the described method is applied to experimental data.
1. INTRODUCTION The problem of signal parameter estimation in sensor array processing has received much attention for several years. Many algorithms have been proposed for estimation of directions of arrival (DOA) of signals impinging on an array of sensors [5]. The formulation of classical DOA estimation problem, however, requires the number of signal sources to be less than the number of sensors in the array in order to reach identi ability. An interesting alternative approach has been proposed in [4], where knowledge about transmitted signals is exploited in addition to the knowledge about the sensor array and as a result the number of signals can exceed the number of sensors. Unfortunately in communication applications the transmitted signals are not generally known a priori. Moreover, in some radio propagation scenarios, especially in urban areas when there is no line of sight between the source and the receiving antenna array, the array receives many T~onu Trump is on leave from Dept. of Radio and Communication Engineering, Tallinn Technical University, Tallinn, Estonia This work was supported by NUTEK.
rays from the vicinity of source and the assumption on number of rays becomes restrictive. In this paper we have made an attempt to model the scattering situation by a reasonable spatial distribution function and estimate the parameters of this distribution instead of marginal directions of arrival. Selection of the distribution function does not appear to be critical as long as we consider relatively small spreads around the nominal DOA. The particular distribution function chosen herein is the Gaussian distribution as it was used in [2] and [9]. Another possible alternative would be to use uniform distribution [6, 8]. It has recently been proposed to use antenna arrays in mobile communication applications to increase the capacity of cellular system [1, 7]. The capacity increase can be obtained either by hosting multiple mobiles per channel per cell or by decreasing the frequency reuse distance. A recent capacity study [9] indicates that for transmitting base-station antenna arrays the knowledge of both, nominal direction to the mobile and the angular spread characterizing the scattering situation around the mobile are of paramount importance.
2. DATA MODEL Assume a large number of rays impinging on a uniform linear array (ULA), all originating from a vicinity of a source and caused by independent re ections around the source. The noise-free signal received by k{th sensor is X xk (t) = s(t) gn ej n ej2(k?1) sin(+ n) ; (1) n
where s(t) is the transmitted signal, gn ej n is the complex ray gain, is the nominal direction to source, n is a small deviation from the nominal direction and is the array element separation in wavelengths. Assuming that the distribution of rays around the nominal direction is normal N (0; ), it can be shown
that k; l{th element of the correlation matrix of the array outputs is given by
E [xk xl ] S [J0 (2(k ? l)) + 2
1 X
m=1
+ 2j
(2)
J2m (2(k ? l))e?2m2 2 cos(2m)]
1 X
m=1
J2m?1 (2(k ? l))e?0:5(2m?1)2 2
sin((2m ? 1))] + 2 kl ; where S is the signal power, Jm denotes the order m Bessel function of the rst kind and kl is the Kronecker delta. If the spread is small (sin ) we can simplify the above expression to E [xk xl ] Sej2(k?l) sin e?2[(k?l) ]2 cos + 2 kl : (3) The correlation matrix can thus be written as Rx S A()A () B(; ) + 2I = S R + 2 I; (4) where denotes the elementwise product, A is the array response in nominal DOA A = [1; ej2sin ; : : : ; ej2(m?1)sin ]T ; (5)
m is the number of sensors and the kl{th element of the m m matrix B given by 2 2 2 (B)kl = e?2[(k?l)] cos : (6) Note that if there is more than one independent source, the correlation matrix includes sum over the 2 ). corresponding Si R(i ; ;i
3. PARAMETER ESTIMATION The maximum likelihood (ML) estimate of the parameter vector = [S 2 2 ]T can be obtained by minimizing the negative log{likelihood function ^ = arg min L(S; 2 ; 2 ; ); (7)
where
L(S; 2 ; 2 ; ) = log jRxj + TrfR?x 1 R^ N g; (8) R^ N is the sample covariance R^ N = N1 PNi=1 x(ti )x (ti ) and N is the number of snapshots.
It can be easily veri ed that for an unambiguous array case the problem is identi able i.e. Rx (1 ) = Rx(2 ) if and only if21 = 2, whenever the number of sensors m 2 and < 1.
This is a four dimensional nonlinear optimization problem which can be solved by a Newton search ^k+1 = ^k ? k H?1 L0 ; (9)
where L0 and H are the gradient and the Hessian respectively and k is a steplength. The i{th element of the gradient is given by
L0i = TrfR?x 1 @@Rix (I ? R?x 1R^ N )g
(10)
and the ij {th element of the Hessian is 2 Hij = Trf[R?x 1 @@ iR@xj ? R?x 1 @@Rix R?x 1 @@Rjx ] (I ? R?x 1 R^ N ) + R?x 1 @ Rx R?x 1 @ Rx ]R?x 1 R^ N g: (11)
@i
@j
The derivatives of the correlation matrix R are provided in the Appendix. Finding derivatives of Rx , given respective derivatives of R is elementary. The search method will converge to the global optimum if initialized close enough to the true parameters. Simulations indicate that for most of practical purposes the \close enough" initial estimate can be obtained from an estimator based on a least squares (LS) t of sample correlation to the model based theoretical correlation matrix ; = arg min Trf(R^ N ? S^R(; ) ? ^ 2 I)2 g ;
= arg min V (; ); (12) ; ^ (13) ^2 = m1 TrfR^ N g ? S; ^ ^ S^ = TrfRTrNfRRg2?g ?TrmfRN g : (14) This estimator needs a two{dimensional search compared to four{dimensional search required by the ML approach. Moreover, by our experience the initialization of the LS estimator is less critical than that of ML. In this paper the ESPRIT algorithm is used to obtain initial estimate of . The spread variance 2 is initialized to zero. Because the LS algorithm is computationally less demanding than the ML approach, it can be of interest by itself. We therefore present the analysis of the LS algorithm in the next section.
? 4TrfRx @@Rx gE [TrfR^ N @@Rx g]
4. ANALYSIS OF THE LS ALGORITHM Consider the criterion function in the form l =k R^ N ? Rx() k2F :
(15)
The sample covariance R^ N tends to true covariance Rx() as N tends to in nity which implies consistent estimates. Since ^N maximizes the criterion we have by mean value theorem that 0 = lN0 (^ ) = lN0 (0 ) + l00 ( )~;
(16)
where ~ = ^ ? 0 is the estimation error and is a vector laying between ^ and 0 . Due to the consistency of the estimate ! 0 if N ! 1 and consequently ~ ?H?1(0 )l0 (0 ); (17) where H is the limiting Hessian. An element of the gradient
@ l = 2Trf(R ? R^ ) @ Rx g x N @ @i N i
(18)
can be recognized as a sum of i.i.d. zero mean variables and therefore we have by the central limit theorem that
p
where
NlN0 (0 ) AsN(0; Q);
(19)
Q = Nlim N E [lN0 (lN0 )T ]: !1
(20)
Hence, the asymptotic distribution of the estimation error is given by
p
N ~ AsN(0; C);
where
(21)
C = H?1QH?1:
(22) The ij {th element of the limiting Hessian can be found as @ Rx @ Rx + (R ? R^ ) @ 2 Rx g Hij = Nlim 2Tr f x N @ @ !1 @ @ j
i
i
j i
Rx @ Rx g: (23) = 2Trf @@ j @i Next, nd the ij {th element of Q. Consider E [ @@ lN @@ lN ] = 4E [Trf(Rx ? R^ N ) @@Rx g i j i @ R Trf(Rx ? R^ N ) @ x g] (24) j Rx gTrfRx @ Rx g = 4TrfRx @@ @j i
j
? 4E [TrfR^ N @@Rx g]TrfRx @@Rx g i j @ R @ R + 4E [TrfR^ N @ x gTrfR^ N @ x g] i j
Since E [R^ N ] = Rx , the rst two terms cancel each other. For the fourth term we have Q4 = E [TrfR^ N @ Rx gTrfR^ N @ Rx g]
@i @j X Rx xm(t) = N12 E [ xl (t) @@ i lm tlm X xo ( ) @@Rjx xp ( )] op op X 1 = N2 E [xl (t)xm (t)xo ( )xp ( )] @R @ Rt lmop x x @i op @j lm :
Application of the Gaussian moment factoring theorem gives now X Q4 = N12 (E [xl (t)xm (t)]E [xo ( )xp ( )] t lmop
+E [xl (t)xp ( )]E [xo ( )xm (t)])
@R @R x x @i op @j lm Rx gTrfRx @ Rx g = TrfRx @@ @j i R x Rx @ R x g + N1 TrfRx @@ @j i
Substituting into (24) we nally obtain E [ @ l @ l ] = 4 TrfR @ Rx R @ Rx g:
@i N @j N
N
x
@i x @j
(25)
Using (25) and (23) in (22) gives the asymptotic covariance of the LS estimate. One can see that this di ers from the Cramer{Rao lower bound, given by the inverse of the Fisher information matrix, the ij {th element of which equals to FIMij = N TrfR?x 1 @ Rx R?x 1 @ Rx g: (26)
@i
@i
5. EXPERIMENTAL RESULTS AND CONCLUSIONS In the simulation study we compared the performance of the ML and LS algorithms to the theoretical results
1.4
B B 1.2 B BB 1.0 B r
JJ b
JHH LS aa? a JJ aPPP 0.6 J PPPX @@ X ML PPXX? 0.4 ``` `hhhhhh 0.2 0.8
each trial the mobile was stationary and 177 samples of relative amplitude and phase were collected. The typical results obtained by applying the ML method to experimental data are shown in Figure 2. The estimation was performed on correlation matrices, obtained by averaging over each 1.5 meter segment. It is possible to see a trend in the DOA estimates corresponding to the 80m movement (approximately 5 at the distance of 1km) of the mobile. The estimated an-
b
b
b
b
b
r
b
b
r
r
r
r
r
r
b
r
b
-5
r
-6
10 20 30 40 50 60 70 80 90 100
Figure 1: RMS value of DOA estimate (in degrees) versus number of snapshots N . presented above. The performance of the ML algorithm was found to t the asymptotic results well for as small number of snapshots as 50 to 100. The performance of the LS algorithm appears to be worse than that of the ML algorithm as predicted by the theoretical analysis. As an example the e ect of number of snapshots, N , on the RMS value of DOA estimates is plotted in Figure 1. The uniform linear array, consisting of 10 elements is used. The received signal consists of 30 components with random amplitudes and phases. The DOA of each component signal is selected randomly under Gaussian distribution with mean zero and standard deviation, corresponding to 5 angular spread. The RMS values are based on 200 experiments. Finally the ML algorithm was applied to experimental data obtained from eld tests of radio propagation in suburban areas, performed by Ericsson Radio Systems AB [3]. The experimental setup was as follows. A 10 element ULA with 0.4 wave{length sensor spacing at 870 MHz was set 2 meters above the roof of an 30 meters high oce building. The broadside of the array was directed toward the mobile source, which was transmitting a sine wave. Data from 3 different experiments was used, during two of which the mobile was located in a suburban and during one it was in a rural area, inside a forest. In all experiments the mobile was placed about 1 km away from the array. In each experiment the mobile was moved 80 meters along a straight line and there were no line of site between the mobile and the receiving antenna array most of the time. The signal was collected with 10 cm intervals along 5 equidistant 1.5 m segments of the line. During
b
-7
b
-8 -9 -10
b
b b
1
2
3 a)
4
5
4 b
3 b
b
2 b
1
1
2
b
3 b)
4
5
Figure 2: Estimated nominal DOA, , (a) and angular spread, , (b) versus trial number using experimental data. gular spread was around 3 degrees which corresponds to scattering area with radius equal to about 50m around the mobile.
6. APPENDIX Here we present the formulae for derivatives of correlation matrix R: Denoting p = 2(k ? l), it can be easily found from (2) that the kl-th elements of the derivative matrices are given by
@R 1 X ?2m2 2 sin(2m)] = ? 2 2 mJ ( p ) e 2 m @ kl m=1 1 X + j 2 (2m ? 1)J2m?1 (p) m=1
e?0:5(2m?1)2 2 cos((2m ? 1))]
2(R)kl [(0:5p)2 2 cos() sin()
+ j 0:5p cos()];
! 1 X 2 2 @R = ? 2 2m2 J2m (p)e?2m cos(2m)] @ 2 kl m=1 1 2 X ? j 2 (2m 2? 1) J2m?1 (p) m=1
e?0:5(2m?1)2 2 sin((2m ? 1))]
?2(R)kl [(0:5p)2 cos2 ()]: The second derivatives are
@ R2 1 X 2 2 = ? 2 4m2 J2m (p)e?2m cos(2m)] @2 kl m=1 1 X ? j 2 (2m ? 1)2 J2m?1 (p) m=1
e?0:5(2m?1)2 2 sin((2m ? 1))]
(R)kl [p4 2 cos2 () sin2 () ! @ R2 @ ( 2 )2 kl = +
! @ R2 @@ 2 kl = ?
+p2 2 sin2 () + j 2p3 2 cos2 () sin() +j 0:5p sin()]; 1 X 2 2 2 4m4 J2m (p)e?2m cos(2m)] m=1 1 4 X j 2 (2m 4? 1) J2m?1 (p) m=1 e?0:5(2m?1)2 2 sin((2m ? 1))] (R)kl [2(0:5p)4 cos4 ()]; 1 X 2 2 2 4m3 J2m (p)e?2m sin(2m)] m=1 1 3 X j 2 (2m 2? 1) J2m?1 (p) m=1 ? 0 e :5(2m?1)2 2 cos((2m ? 1))]
(R)kl [?0:5p4 2 cos3 () sin() +(2(k ? l))2 sin() cos() ?j 0:5p3 cos3 ()]:
7. REFERENCES [1] S. Anderson, M. Millnert, M. Viberg and B. Wahlberg , \An Adaptive Array for Mobile Communication Systems," IEEE Trans. on Vehicular Technology, vol.40, pp.230{236, no. 1, Feb. 1991. [2] F. Adachi, M.T. Feeney, A.G. Williamson and J.D. Parsons, \Crosscorrelation Between the Envelopes of 900MHz Signals Received at a Mobile Radio Base Station Site," IEE Proc. Vol.133, Pt.F. Oct. 1986. [3] U. Forssen, J. Karlsson, B. Johannisson, M. Almgren, F. Lotse and F. Kronestedt, \Adaptive Antenna Arrays for GSM900/DCS1800," Proc. 43rd IEEE Vec. Tech. Conf., pp.605{609, Stockholm, 1994. [4] J. Li, B. Halder, P. Stoica, M. Viberg and T. Kailath, \Decoupled Maximum Likelihood Angle Estimation for Signals with Known Waveforms," submitted to IEEE Trans. on Signal Processing. [5] B. Ottersten and M. Viberg, Sensor Array Signal Processing, In preparation. [6] J. Salz and J.H. Winters, \E ect of Fading Correlation on Adaptive Arrays in Digital Wireless Communications," Proc. ICC'93, pp. 1768{1774, Geneva, May 1993. [7] S.C.Swales, M.A.Beach, D.J.Edwards and J.P. McGeehan, \The Performance Enhancement of Multibeam Adaptive Base-Station Antennas for Cellular Land Mobile Radio Systems," IEEE Trans. on Vehicular Technology, vol. 39, pp.56{ 67, no. 1, Feb. 1990. [8] Y. Yamada, K. Kagoshima and K. Tsunekawa, \Diversity Antennas for Base and Mobile Stations in Land Mobile Communication Systems," IEICE Trans. Vol.E 74, Oct. 1991. [9] P Zetterberg and B. Ottersten, \The Spectrum Eciency of a Basestation Antenna Array System for Spatially Selective Transmission," Proc. 43rd IEEE Vec. Tech. Conf., pp. 1517{1521, Stockholm, 1994.
