Joint Estimation of 2-D Arrival Angles, Propagation Delays, and Doppler Frequencies in Wireless Communications Martin Haardt,1 Christopher Brunner,1, 2 and Josef A. Nossek 2 1. Siemens AG, Mobile Networks, OEN MN P 36 Hofmannstr. 51, D-81359 Munich, Germany Phone: +49 (89) 722-29480, Fax: +49 (89) 722-44958 E-mail:
[email protected]
Abstract – In wireless communications, the estimation of the dominant directions of arrival at the base station facilitates the use of advanced mobile positioning algorithms, the implementation of beamforming schemes for smart antennas, and the development of directional channel models. In this paper, we present an efficient high-resolution technique to estimate the multipath propagation structure of the wireless channel in terms of the two-dimensional arrival angles (azimuth and elevation), the corresponding propagation delays, and the Doppler shifts of the dominant wavefronts. Unlike the short-term fading, these four parameters are stationary over long time intervals and can, therefore, be estimated with a high accuracy. In the proposed estimation procedure, the chip waveform of the sounding sequence and the correlation matrix of the (colored) additive noise are taken into account. To achieve these goals, 4-D Unitary ESPRIT is applied to the measurements of a uniform rectangular array in the space-frequency domain. Automatic pairing of the 4-D estimates is obtained by computing the simultaneous Schur decomposition of four real-valued, non-symmetric matrices. This efficient closedform procedure neither requires any search nor any other heuristic pairing strategy and can estimate the parameters of significantly more wavefronts that there are antenna elements. In a direct sequence CDMA system, an extension of this channel sounding technique can be used to estimate the channel parameters of each user.
1. Introduction Spatio-temporal radio channel measurements that include the dominant directions of arrival (DOAs), the Doppler frequencies, and the corresponding propagation delays are desired in a variety of applications including sonar, radar, seismic exploration, and mobile communication. In mobile communication, they facilitate the design of efficient smart antenna concepts which exploit the inherent spatial diversity of the radio channel to achieve a significant increase of the spectral efficiency. Two-dimensional (2-D) antenna array geometries at the base station of a mobile communication system allow the estimation of azimuth as well as elevation angles of the dominant wavefronts that impinge on the array. Taking the elevation into account, facilitates the reduction of the near-far problem in cellular radio systems. High-resolution channel sounding algorithms that reveal the multipath propagation structure in terms of 1-D arrival angles and propagation delays (2-D channel sounding) have been presented in [6, 8]. An extension to the joint estimation of 2-D arrival angles and propagation delays (3-D channel sounding) has been provided in [2]. In these references, no appreciable Doppler shift was assumed. Here, we extend these results to the 4-D case by jointly estimating the 2-D arrival angles and propagation delays along with the corresponding Doppler shifts of the dominant multipath components. Thereby, the time-variance of the channel parameters is reduced, and the result-
2. Institute for Network Theory and Circuit Design Munich University of Technology, D-80290 Munich, Germany Phone: +49 (89) 289-28511, Fax: +49 (89) 289-28504 E-Mail:
[email protected] ing estimation accuracy is greatly enhanced if there is an appreciable Doppler shift. In a direct sequence CDMA system, extensions of these multidimensional channel sounding techniques can be used to estimate the channel parameters of each user [2, 9]. 2. 4-D Data Model To measure the propagation characteristics of the radio channel, assume that the transmitter emits a sounding signal s(t) consisting of a repeatedly transmitted (spreading) sequence w(t), where
s(t) =
1 X
n=;1
w(t ; nT )
and
w(t) =
N X c
m=1
dm p(t ; mTc):
(1)
The spreading sequence w(t) is composed of Nc (complex) chips dm , 1 m Nc , weighted by the chip waveform p(t) 2 R. If the chips are of duration Tc , the spreading sequence w(t) is of length T = Nc Tc . Moreover, we define another sequence bm of length Nc such that the sequences dm and bm satisfy the following cross-correlation properties
8N < c rdb (m) = d(n+m)modN bn : n=0 0 Nc ;1
X
c
if m = l Nc with l ZZ otherwise,
2
where the overbar denotes complex conjugation. Several choices for the sequences dm and bm are, for example, discussed in [5].1 The (mobile) sounding unit is assumed to be in the far field of the antenna array so that the impinging wavefronts are approximately planar. For simplicity, we assume that a uniform rectangular array (URA) of omnidirectional sensors is located at the receiver. Extensions to other 2-D centro-symmetric array geometries with a dual invariance structure and arbitrary (but identical) sensor characteristics are straightforward [1]. Assume that the URA contains M1 sensors in x-direction and M2 sensors in y -direction such that the noise-corrupted array measurements at time t are contained in the matrix 2D (t), which is of size M1 M2 . Due to multipath propagation, d narrowband planar wavefronts of wavelength are impinging on the URA. Thus, there are d distinct multipaths, each parameterized in terms of the corresponding azimuth i , elevation i , propagation delay i , Doppler frequency fDi , and complex amplitude i (fading). In the baseband, the noise-corrupted measurements
X
1 The sequences dm and bm might, for example, be realized as shifted m-sequences such that dm bm 8m [5]. In this case, the correlators of Figure 2 are matched filters which provide the optimum signal to noise ratio (SNR) at their outputs.
=
3. Efficient High-Resolution Parameter Estimation
z
3.1. Generation of the Temporal Invariance Structure
f
i
Array
y
ui
i
180 < i 0 90 x(t) = vecfX2D (t)g may be modeled as Figure 1: Definitions of azimuth (; ( i ).
x(t)
= = (4)
where i A2D
d X i=1
A2D
(1) (2) )
a(i
i ej
i
(4)
i
180 ) and elevation
t s(t ; i ) + n(t)
1
(1) (2) ) 2
(1) (2) ) = a((2) ) a((1) )
i
(r) ) =
a(i
i
1 e
j
(r)
i
i
:::
where
e (Mr ;1) j
(r)
i
iT
(3)
r = 1 2. Furthermore, the diagonal matrices = diagfi gdi=1 and (t) = diagfe i t gdi=1 contain the complex amplitudes i and Doppler shifts e i t of the d multipath components, respectively. (4)
j
(4)
j
Finally, the vector
s(t) = s(t ; 1 )
s(t ; 2 ) : : :
s(t ; d ) T
(4)
contains d delayed versions of the sounding signal s(t). Here, i denotes the unknown propagation delay of the ith propagation path. Recall that the proposed 4-D channel sounder jointly estimates the 4-D parameters of the dominant wavefronts, i.e., their 2-D spatial (1) (2) frequencies i and i as well as their propagation delays i and Doppler frequencies fDi . The phase changes due to the Doppler frequencies over the duration T = Nc Tc of the spreading sequence w(t) are so small that they
t
can be neglected in practice [7]. Therefore, the Doppler shifts e i during the duration of the nth spreading sequence in (2) may be apj
proximated as e
j
(4)
i
(nT +T=2) :
Z =!
(4)
0
1 X
rpp ( ; nT ):
n=;1
(5)
rsc ( )e; n! d: 0
j
;=!
0
;
f g
h
= Nc Tc = !20
Then it is well known that the Fourier series of rsc (
contains the 2-D directional information. Here, vec A denotes a vector-valued function that maps an m n matrix A into an mndimensional column vector by stacking the columns of the matrix. (1) The spatial frequencies in x-direction i = (2=)x ui and the (2) spatial frequencies in y -direction i = (2=)y vi are scaled versions of the corresponding direction cosines ui = cos i sin i and vi = sin i sin i of the ith source relative to the x- and y -axes, 1 i d as illustrated in Figure 1. Moreover, x and y are the distances between the sensors in x- and y -direction. Note that the 2-D array steering vectors can be written as Kronecker products of 1-D Vandermonde vectors, a(i
Rsc (n) = 2!0
(2)
(2) a((1) d d )
a(2
1
according to
Its Fourier series will be denoted as
= 2f i , and the 2-D array steering matrix D
c
rsc ( ) = E s(t + )c(t) = Nc
(t) s(t) + n(t)
= a((1) (2) )
N X
bm p(t ; mTc ) m=1 is a periodic repetition of rpp ( ) with period T c(t) =
vi x
g
Let rpp ( ) = E p(t + )p(t) be the autocorrelation function of the chip waveform p(t) such that the cross-correlation function of the periodic sounding signal s(t) defined in (1) with the correlator sequence
Unit Ball
(3)
; i ) is given ; ;
by e j n!0 i Rsc (n) = ej n i Rsc (n) where i = 2i =T and the propagation delay i < T . If rsc ( ) and rsc ( i ) are sampled above Nyquist rate, this property still holds for the DFTs of the resulting sequences. Assume that the array measurements x(t) are sampled at rate fs = MTcc , where the number of samples per chip Mc is an integer. If one samples above Nyquist rate, the correlator (5) can also be realized in the discrete time domain. Using Nc Mc samples of the correlator output, we perform an (Nc Mc )-point DFT for each antenna as illustrated in Figure 2. Note that only the spectral lines that correspond to the M3 largest values of Rsc (n) are computed. Thereby, the influence of the additive noise and the computational complexity can be reduced significantly. Let each row of C(M1 M2 ) (Nc Mc ) contain Nc Mc samples of the output cor (n) of the corresponding antenna after the correlator. The transformation from the space-time domain to the space-frequency domain is achieved by post-multiplying cor (n) by the (Nc Mc ) M3 DFT matrix W such that 3D (n) = cor (n)W: It is well known that each column of W is of the form (3)
X
2
X X
X
= 1 e; `!0 e; 2`!0 : : : e; (Mc Nc ;1)`!0 T : Here, the columns of W compute M3 (Nc Mc ) spectral lines cenw`
j
j
w
j
tered at DC according to W=
Nc Mc ;(M3 ;1)=2
w0
;
w(M3 1)=2
where we have invoked the wrap-around property of the DFT matrix. The construction of the space-frequency data matrix 3D (n) is also illustrated in Figure 2. Due to the described Fourier properties, the data model in (2), and the approximation at the end of Section 2, 3D (n) can be expressed as
X
X X3D(n) = d X i=1
i e
(6) j
(4)
i
(nT +T=2) a((1) (2) ) a((3) )T + N
i
N
i
i
3D
(n)
where 3D (n) contains the additive noise in the space-frequency (3) domain, the Vandermonde vector a(i ) CM3 is defined in (3), and = diag Rsc (n) (nM=3Nc1)M=c2 (M3 1)=2 :
2
f
g ; ; ;
A/D
baseband I & Q demod.
Mc samples / chip
3.2. Joint 4-D Frequency Estimation
Nc Mc samples
correlator
DFT
Mc
baseband I & Q demod.
A/D samples / chip
x11x12 x1M3 Nc Mc samples
correlator
DFT
x21x22 Mc samples / chip
baseband I & Q demod.
A/D
correlator
x2M3 Nc Mc samples
Mc samples / chip A/D
uniform rectangular array of size
M1 M2
X3D (n) =
x11 x21 x31 x41
x12 x22 x32 x42
correlator
x1M3 x2M3 x3M3 x4M3
x31x32 x3M3 Nc Mc samples DFT
x41x42
x4M3
()
Figure 2: Generation of the 3-D invariance structure of X3D n using a uniform rectangular array of size M1 M2 . Here, M3 Nc Mc spectral lines are calculated in the frequency domain. The 4-D invariance structure is realized by combining matrices X3D n from M4 different periods of the spreading sequence.
=2 2 ()
X
Apart from this diagonal matrix , 3D (n) may be interpreted as a single snapshot of a (virtual) uniform cube array of the size M1 M2 M3 corresponding to the nth repetition of the spreading sequence in (1). Assuming that the 4-D channel parameters do not change significantly during several periods of the spreading sequence w(t), we, therefore, define the matrix
X4D =
vec X
3D
(n1 );1
X
(n2 );1 vec X3D (nM4 );1
vec
3D
which is of size (M1 M2 M3 ) M4 .2 Using equation (6) and ele-
mentary properties of the Kronecker product, the matrix expressed as
X4D =
X4D can be
d X (2) (3) (4) T i0 a((1) i i i ) a(i ) + N i=1
4D
(7)
(4) ) 2 CM4 is defined in (3),
where the Vandermonde vector a(i
(1) (2) (3) ) = a((3) ) a((2) ) a((1) )
a(i
i0 = i e
j
(4)
i
i
T=2 , and
N = vec N 4D
i
i
3D
(n1 );1
i
N
i
3D
In the sequel, we use N realizations of the four-dimensional vector x4D = vec 4D at N different points in time. Moreover, 4-D smoothing [3] (using Msubr elements in the rth direction, 1 r 4, and a total of L subarrays) is applied to decorrelate the wavefronts and to reduce the computational complexity.
fX g
2 Note that the estimation accuracy of the Doppler shifts is increased if n1 n2 : : : nM4 do not denote subsequent repetitions of the spreading sequence, e.g., they may be chosen such that ni+1 = ni + N .
0
0 fN g
f0 0 g
f
;
g
3.3. Estimation of the Amplitudes If the array geometry is known and the array is calibrated, the 2-D steering matrix A2D is completely specified by the estimated spatial (1) (2) frequencies i and i . Moreover, the diagonal matrices and (t) in (2) are completely determined by the estimated temporal (3) (4) frequencies i and i , respectively. Fast fading influences the amplitudes over short intervals, whereas DOAs, Doppler frequencies, and propagation delays remain constant over longer periods. We, therefore, estimate the amplitudes over Q sequences, where Q might be significantly smaller than M4 . To this end, let us collect QNc Mc samples of x(t) as defined in (2). Then the resulting matrix XQ C(M1 M2 ) (QNc Mc ) may be expressed as
2
= x(t1 ) x(t2 ) x(tQ ) = A2D S + NQ where S = (t1 )s(t1 ) (t2 )s(t2 ) (tQ )s(tQ ) XQ
(n2 );1 vec N3D (nM4 );1
vec
fX g
f0 0 g
DFT
baseband I & Q demod.
Note that the arrival angles, propagation delays, and Doppler frequencies of each multipath component in (2) are relatively stationary, whereas the complex amplitude i of each path is nonstationary and subject to (Rayleigh) fading. Therefore, it is the goal of the 4-D channel sounder to provide high-resolution estimates of the d (1) (2) (3) (4) four-dimensional frequency quadruplets i , i , i , and i along with their correct pairing. A very simple and efficient solution is provided by 4-D Unitary ESPRIT. Let x4D and n4D denote smoothed versions of x4D = vec 4D and n4D = vec 4D , respectively. Then estimates of the smoothed signal covariance maH and of the corresponding covariance matrix trix Rxx = E x4D x4D 2 Rnn = E n n H can be obof the noise (plus interference) N 4D 4D tained as explained in [2, 9]. Note that we do not have to estimate 2 . Both covariance matrices are used as an input the scalar factor N for 4-D Unitary ESPRIT as presented in [1, 3]. In this application, 4-D Unitary ESPRIT is modified slightly to take into account a noise 2 Rnn that is not necessarily equal to the scaled correlation matrix N identity matrix, i.e., colored noise. This modification has been described in [2]. Automatic pairing of the 4-D estimates is achieved by computing the simultaneous Schur decomposition (SSD) of four realvalued, non-symmetric matrices. This closed-form procedure neither requires any search nor any other heuristic pairing strategy [1, 3]. By jointly estimating the two-dimensional arrival angles, the Doppler frequencies, and the corresponding propagation delays, the proposed 4-D channel sounder is able to resolve more dominant wavefronts than there are antenna elements. The maximum number of wavefronts that can be estimated via 4-D Unitary ESPRIT is given by dmax = min m1 m2 m3 m4 2LN where mr = M (Msubr 1)=Msubr 1 r 4, [3]. Moreover, it avoids difficulties in cases where two or more wavefronts have equal arrival angles, propagation delays, or Doppler frequencies.
(8)
and
NQ contains the corresponding noise samples. Hence, the linear minimum variance unbiased estimate of the amplitudes = 1 2 : : : d ]T is given by
=
;BH R;1 B;1 BH R;1 vecfXQg NN NN
where RNN = E vec NQ vec NQ H B = ST A2D , and the Khatri-Rao product of two matrices denotes the columnwise Kronecker product.
f f g f g g
power versus azimuth and delay
0
0
−5
−5
−10
−10
−10
−15 −20 −25
normalized power [dB]
0
−15 −20 −25
−15 −20 −25
−30
−30
−35
−35
−35
−40 1
−40 5
−40 200
0.5
−30
1
−1
azimuth
delay in seconds
u−axis
−5
x 10
−5
−200
0 −5
−10
−10
−10
−25
normalized power [dB]
0 −5
normalized power [dB]
0
−20
−15 −20 −25
−20 −25
−30
−35
−35
−35
−40 1
−40 5
−40 200
−30
5
1 0.5
0
0
0
−0.5
100
5
delay in seconds
u−axis
azimuth
−5
x 10
−5 delay in seconds
3 2
−100
−6
−1
4
0
0
−6
x 10
−0.5 −1
v−axis
x 10 delay in seconds
−15
−30
0.5
1
estimated power versus Doppler shift and delay
−5
−15
−6
0
Doppler shift in Hz
delay in seconds
estimated power versus azimuth and delay
estimated power in the u−v plane
3 2
−100
−6
−1
v−axis
5
0
−6
x 10
−0.5
4
0
0
0
−0.5
100
5
0.5
0
power [dB]
power versus Doppler shift and delay
−5
normalized power [dB]
normalized power [dB]
power in the u−v plane
−200
1
x 10
0
Doppler shift in Hz
−6
delay in seconds
Figure 4: Power of the impinging wavefronts as a function of the 2-D arrival angles (direction cosines), the propagation delays, and the corresponding Doppler shifts in the u-v plane (on the left side), in the azimuth-delay (- ) plane (in the middle), and in the Doppler-delay (Df- ) plane (on the right side). The 3-D plots on the top show the power of the impinging wavefronts as predicted via ray tracing, whereas the bottom row visualizes the 4-D estimates obtained by applying the proposed 4-D channel sounder (simulation parameters: M1 M2 , x : , M c , Nc , M3 , M4 , y Msub1 Msub2 , Msub3 , Msub4 , Tc : s, fc GHz, N ,Q , SNR dB).
=
=4
= 10
=4
= 0 244
=2
4. Simulations In the simulations, a sounding signal characterized by a square-root raised-cosine spectrum ( = 0:22) is used. First we compare the performance of the 2-D version of the presented channel sounding technique with the performance of JADE as described in [8]. The RMS errors of the azimuths (left side) and the corresponding propagation delays (right side) averaged over the two wavefronts are shown in Figure 3. Note that the 2-D version of the presented chanRMSE of the azimuths in degrees
1
10
Nc Nc Nc Nc
= 15 = 15 = 63 = 63
Mc Mc Mc Mc
RMSE of the delays in Tc
1
10
=2 =4 =2 =4
Nc Nc Nc Nc
= 15 = 15 = 63 = 63
Mc Mc Mc Mc
=2 =4 =2 =4
0
10 0
10
−1
RMSE
RMSE
10 −1
10
= 5 5 = = 0 45 = 20 = 1 = 20
= 10
= 63
= 51
=4
nel sounder provides more accurate results than JADE. Moreover, we present simulation results that use realistic 2-D arrival angles derived from a ray tracing program for macrocell environments [4]. They are based on a three-dimensional topographical model of downtown Munich as explained in [2]. The top row of Figure 3 depicts the power of the impinging wavefronts as a function of the 2-D arrival angles (direction cosines), the propagation delays, and the corresponding Doppler shifts in the u-v plane (on the left side), in the azimuthdelay (- ) plane (in the middle), and in the Doppler-delay (fD - ) plane (on the right side) as predicted via the ray tracing tool. The bottom row of the same figure shows the estimation results obtained with the proposed channel sounder based on 4-D Unitary ESPRIT. Clearly, the dominant 2-D arrival angles, the propagation delays, the corresponding Doppler shifts, and the amplitudes of the dominant wavefronts are estimated correctly. References
−2
10
[1] M. Haardt, Efficient One-, Two-, and Multidimensional High-Resolution Array Signal Processing, Shaker Verlag, Aachen, Germany, 1996, ISBN 3-8265-2220-6. −2
[2] M. Haardt, C. Brunner, and J. A. Nossek, “Efficient high-resolution 3-D channel sounding”, in Proc. 48th IEEE Vehicular Technology Conf. (VTC ’98), Ottawa, Canada, May 1998.
10
−3
10
o x
JADE 2-D Channel Sounder
o x
−3
10
10
[3] M. Haardt and J. A. Nossek, “Simultaneous Schur decomposition of several non-symmetric matrices to achieve automatic pairing in multidimensional harmonic retrieval problems”, IEEE Trans. Signal Processing, vol. 46, pp. 161–169, Jan. 1998.
JADE 2-D Channel Sounder
−4
20
30
40
SNR [dB]
10
10
20
30
40
SNR [dB]
Figure 3: RMS error (RMSE) of the azimuths (left side) and the corresponding propagation delays (right side) averaged over both wavefronts (1 1 : Tc , and 2 2 : Tc ) as a function of the SNR (averaged over 500 trials) for different values of Nc and Mc using a uniform linear array of 5 elements ( x : , N ). The circles (o) denote the results obtained via JADE (m1 m2 L Nc , see [8]), whereas the crosses (x) denote joint 2-D estimation of the spatial and temporal frequencies via the 2-D version of the presented channel sounder (M3 , Msub1 , Msub3 ).
= 10
= 13
= 15 = 6 1
= 0 5 = 10 =4 =1 =
= 10
=2
=5
[4] T. K¨urner, D. C. Cichon, and W. Wiesbeck, “Concepts and results for 3D digital terrain-based wave propagation models: An overview”, IEEE Journal on Selected Areas in Communications, vol. 11, pp. 1002–1012, Sept. 1993. [5] U. Martin, Ausbreitung in Mobilfunkkan¨alen: Beitr¨age zum Entwurf von Meßger¨aten und zur Echosch¨atzung, Ph. D. dissertation, University of Erlangen, Erlangen, Germany, Oct. 1994, in German. [6] U. Martin, “A directional radio channel model for densely build-up urban area”, in Proc. 2nd European Personal Mobile Communications Conference, pp. 237–244, Bonn, Germany, Sept. 1997. [7] K. I. Pedersen, B. H. Fleury, and P. E. Mogensen, “High resolution of electromagnetic waves in time-varying radio channels”, in Proc. 8th IEEE Int. Symp. on Personal, Indoor and Mobile Radio Commun. (PIMRC), Helsinki, Finland, Sept. 1997. [8] A. J. van der Veen, M. C. Vanderveen and A. Paulraj, “Joint angle and delay estimation using shift-invariance techniques”, IEEE Trans. Signal Processing, vol. 46, pp. 405–418, Feb. 1998. [9] M. D. Zoltowski, J. Ramos, C. Chatterjee, and V. Roychowdhury, “Blind adaptive 2D rake receiver for DS-CDMA based on space-frequency MVDR processing”, IEEE Trans. Signal Processing, June 1996, submitted for publication.