Direct Location Estimation using Indoor Wireless Communication Systems Peter Voltz Polytechnic Institute of NYU
Overview • Localization by TOA Estimation • Maximum Likelihood Line of Sight TOA Estimation – Analysis – Performance examples and comparison
• Direct Location Estimation – Analysis – Performance examples and comparison
Multipath Issue – Suppose we transmit this pulse s (t )
t (ns)
And the channel impulse response is
Then the received signal looks like
• Time of Arrival (TOA) Estimation in multipath • OFDM Signal Structure • Maximum Likelihood estimation of the LOS delay only
• Frequency Domain (OFDM subcarriers) Technique • Uses only Receiver FFT outputs • Comparison with approximate Maximum likelihood estimation of individual paths (MODE)
OFDM Parameters • • • •
64 Subcarriers 52 Occupied Subcarriers 4 µs Symbol Interval 0.8 µs Guard Interval
The MODE Method for Channel Time Delay Estimation P. Stoica and K.C.Sharman, “Maximum Likelihood Methods for Direction-of-Arrival Estimation,” IEEE Trans. On Acoustics, Speech and Signal Processing, Vol. 38, No.7, July 1990. R. Wu, J. Li, Z-S Liu, “Super Resolution Time Delay Estimation via MODE-WRELAX,” IEEE Trans. On Aerospace and Electronic Systems, Vol. 35, No. 1, January 1999.
The MODE Method for Channel Time Delay Estimation Our OFDM symbol duration is T + TG , where TG is the guard interval, and T is the receiver integration time over which the subcarriers are orthogonal. A single symbol of the transmitted OFDM signal with N subcarriers and data vector d = [d0 , d1,..., d N −1 ] can be expressed as
s (t ) =
N −1
1 ∑ T k =0
2π kt dk e T j
;
−TG < t < T .
We assume here that the symbol begins at t = −TG seconds so that a receiver would begin its T - second integration period at t = 0 if there were no time delay between transmitter and receiver.
Assume, on the other hand, that the signal is received after passing through a multipath channel with impulse response h(t ) =
L −1
∑ ai δ (t − τ i )
i =0
and frequency response H (ω ) =
L −1
∑
ai e − j ωτ i
i =0
in which we assume that 0 ≤ τ 0 ≤ τ 1 ≤ ⋅ ⋅ ⋅ ≤ τ L −1 ≤ TG . Then, on the interval 0 < t < T , the received signal, including complex receiver noise with power density N 0 , is
L −1
N −1
1 r (t ) = ai ∑ ∑ T i =0 k =0
2π k (t −τ i ) dk e T + n(t ) j
or, 2π L −1 − j kτi 1 r (t ) = d k ∑ ai e T ∑ T k =0 i =0 N −1
j 2π k t T e + n(t ) .
After the standard receiver sampling and FFT processing, the k th element of the output vector is
2π L −1 − j kτ i y k = d k ∑ ai e T i =0
+ n k ,
where nk is complex gaussian noise with variance σ 2 = N 0 . If d k = 1 then (4) is effectively
yk =
L −1
∑ ai e
i =0
or
−j
2π kτ i T
+ nk
2π yk = H T
k + nk .
Suppose we measure every mth subcarrier. Then L −1
ymk = ∑ ai e i =0
−j
2π m kτi T
+ nmk
Which may be written in vector form as y = Aa + n With
1 e− jθ 0 2 e− jθ 0 A = e− jθ 0 3 ⋅ ⋅ − j θ 0 M −1 e
( (
(
) )
)
1
⋅ ⋅
e − j θ1
⋅ ⋅
− j θ1 2
⋅ ⋅
(e ) − jθ 3 (e ) 1
⋅ ⋅
(e )
⋅ ⋅
⋅ ⋅
− j θ1 M −1
e − j θ L −1 2 e − j θ L −1 2π − j θ L −1 3 ; θ i = mτ i e T ⋅ ⋅ − j θ L −1 M −1 e 1
( (
(
The probability density function of y is
) )
)
(
1
p ( y ) = N 2 N exp − y − A a π σ
2
/σ 2
)
And the maximum likelihood estimate of the channel amplitudes and time delays (angles) is found by minimizing the quantity Q = y − Aa
2
.
For any choice of the angles (which determine A) the solution for a is well known to be
( ) *
aˆ = A A
−1 *
A y
and the resulting error is
Q = y − A aˆ *
2
( )
= y − A A* A
( ) *
Q = y I − A A A
−1 *
A y
( )
−1 * 2
A y = I − A A* A
−1 *
2
A y
.
Next we must minimize this over the angles. This is a nonlinear minimization problem. It can be parameterized differently in terms of the following polynomial L
(
b( z ) = b0 z L + b1 z L −1 + . . . + bL = b0 ∏ z − e − j θ l
)
l =1
Since b(z ) has all its zeros on the unit circle its coefficients must satisfy the constraints
bl = bL∗ −l ,
l = 0,1, ... , L
Define the M × ( M − L) matrix B b0 ⋅ ⋅ B = bL 0
⋅
⋅
0 ⋅ b0 ⋅ ⋅ ⋅ bL
Since the zeros of b(z ) are e − j θ l it follows that
B* A = 0 .
Equation (15) shows that the space spanned by the columns of B is the orthogonal complement of the space spanned by the columns of
( ) *
−1 *
A. Now the matrix I − A A A A is the projection matrix onto the orthogonal complement of the space spanned by the columns of A. Due to (15), another way to express this projection matrix is
( ) *
B B B *
−1 *
B . Therefore we can re-write Q in equation (11) as
( ) *
Q= y B B B
−1 *
B y .
Minimizing Q is still a nonlinear optimization problem in the bl . It can be approached in an iterative way as follows, however.
First, let β be a vector containing the real and imaginary parts of the first half of b (the independent parameters). Since B∗ y is a linear function of the elements of β , we can re-write it as B* y = S β
For a suitable matrix S . Then (16) can be re-written * *
( ) *
Q=β S B B
−1
Sβ
Now we iterate as follows:
.
1.Solve for βˆ which minimizes Q = β * S * W S β , where W = I in the first iteration. The solution is the real eigenvector corresponding to the smallest eigenvalue of Q = S * W S . 2.Let Bˆ be the matrix B populated suing the values of βˆ in the vector b . 3.Solve for βˆ which minimizes Q = β * S * W S β , where now
( ) ˆ*
Q = B Bˆ
−1
.
4.Iterate steps 2 and 3 several times.
Effect of Transmit and Receiver Filters The above development depends on the channel impulse response being simply a finite number of impulses. In reality there are transmit and receive filters which causes the intervening channel from OFDM generation to reception to lose this property. But the approach still works provided the cascade of transmit and receive filters has unity gain and linear phase over the band containing the pilot subcarriers. To see this note first that regardless of the intervening channel, if the guard interval is longer than the overall channel delay spread, then the OFDM outputs will be (assume unity pilot symbols): 2π yk = H eff T
k + nk
where H eff (ω ) = H (ω ) H T (ω ) H R (ω )
And since, on the band of interest, H T (ω ) H R (ω ) = e − j ω t0
we have H eff (ω ) = H (ω ) e
− j ω t0
L −1 = ∑ ai e − j ωτ i i =0
− jω t 0 e
H eff (ω ) =
L −1
∑
ai e − j ω (τ i + t0 )
i =0
So that yk =
L −1
∑ ai e
−j
2π k (τ i + t 0 ) T
+ nk
i =0
Which is the same as (5) but with an added delay in each impulse.
Maximum Likelihood Method for LOS Delay Only P. Voltz and D. Hernandez, “Maximum Likelihood Time of Arrival Estimation for Real-Time Physical Location Tracking of 802.11a/g Mobile Stations in Indoor Environments,” Proceedings of PLANS 2004 (Position Location and Navigation Symposium), April 26-29 2004, Monterey, CA.
OFDM Signal: N −1
1 s (t ) = ∑ T k =0
2π kt dk e T
− TG < t < T
;
j
Multipath Channel: h(t ) =
L −1
∑
;
ai δ (t − τ i )
0 ≤ τ 0 ≤ τ1 ≤ ⋅ ⋅ ⋅ ≤ τ L −1 ≤ TG
i =0
Received signal: r (t ) =
L −1
N −1
1 ai ∑ ∑ T i =0 k =0
2π k (t −τ i ) T dk e + n(t ) j
2π L −1 − j k τi 1 T r (t ) = d k ∑ ai e ∑ T k =0 i =0 N −1
j 2π k t T + n(t ) e
After the standard receiver sampling and FFT processing, the k th element of the output vector is 2π L −1 − j kτ i T yk = d k ∑ ai e i =0
where nk
+ nk
is complex gaussian noise with variance σ 2 = N 0
2π L −1 − j k (τ 0 +τ i ) yk = d k ∑ ai e T + nk i =0 2π 2π −j kτ 0 L −1 −j kτ i T = dk e T ∑ ai e i =0
yk = H k d k e
where L −1
−j
Hk
H k = ∑ ai e i =0
2π kτ0 T + nk
−j
is given by 2π kτ i T
;
τi = τ 0 +τi + nk
k = 0 → N −1
With the zero-delay frequency response defined as L −1
H 0 (ω ) = ∑ ai e − jω τ i i =0
We see that 2π Hk = H0 k T
is the channel zero-delay frequency response at the frequency of the k th
subcarrier,
ωk =
2π k T
Impulse Response h(t )
t (ns)
Zero delay Impulse Response h0 (t )
t (ns)
CHANNELS WITH KNOWN STATISTICS Define h = [ H 0 , H1, .. , H N −1 ]T
Assume that h is a zero mean, complex gaussian vector with
{ }
known covariance matrix K h = E hh∗ . Then y = GD h + n 2π 2π − j 2π τ − j 2 τ − j ( N −1)τ 0 0 0 G = diag 1, e T , e T , ... , e T
D = diag{d 0 , d1, d 2 , ... , d N −1}
,
n=
a complex gaussian vector independent of h with
{ }
K n = E nn∗ = σ 2 I
Note we assume that both h and n have zero mean. Now, for a given value of τ 0 , y is a zero mean complex gaussian vector, and has the probability density function given below p( y τ 0 ) =
1
π N Det ( K y )
exp(− y ∗ K −y 1 y ) ,
where K y = GDK h D∗G ∗ + σ 2 I
The maximum likelihood estimate of τ 0 is the value that maximizes p( y τ 0 ) . We can show that Det( K y ) does not depend on τ 0 . In fact for any matrices A and B of appropriate size, Det ( I + AB ) = Det ( I + BA )
Det ( K y ) = Det (GDK h D∗G ∗ + σ 2 I ) = σ 2 N Det ( I + GDK h D∗G ∗ / σ 2 ) = σ 2 N Det ( I + G *GDK h D∗ / σ 2 )
and G *G = I
, so that
Det ( K y ) = σ 2 N Det ( I + DK h D∗ / σ 2 )
; not dependent on τ 0
Next, for the exponent in p( y τ 0 ) , we must say a little more about the matrix K h. It can be very poorly conditioned if the delay spread of the channel is small. In this case, if K h has rank p < N , it can be factored as K h = RR*
where R is an N × p matrix of rank p . Then K y can be written
(
K y = σ 2 I + GDRR* D∗G ∗ = σ 2 I + GDRR* D∗G ∗ / σ 2
)
and, using the well known fact that
( I + AB )−1 = I − A( I + BA)−1 B when either inverse exists, we can write K −y1 as (use (use A = GDR / σ and B = R* D∗G * / σ ) −1 * ∗ * * ∗ * 1 GDR R D G GDR R D G K −y 1 = I− I+ σ σ σ σ σ2
Now, since G *G = I this reduces to
(
)
−1 * ∗ * 1 −1 2 * ∗ Ky = I − GDR σ I + R D DR R D G 2
σ
.
.
Using this in p( y τ 0 ) , the maximum likelihood solution for time delay τ 0 is the maximizer of *
(
2
* ∗
Q (τ 0 ) = y GDR σ I + R D DR
)
−1 * ∗ *
R D G y
.
Q (τ 0 ) = y*GFG* y
(
2
* ∗
F = DR σ I + R D DR
)
−1 * ∗
R D
.
Note that the matrix F can be pre-computed if σ 2 , K h and D are known. Note also that a one-dimensional search must be performed to find the maximum.
Consider a single path channel h(t ) = a0 δ (t − τ 0 ) H k = a0 e
−j
2π kτ 0 T
1 1 h = a0 ⋅ = a0 u ⋅ 1
= a0
y = G h + n = a0 1,
2π 2π − j τ 0 − j 2τ 0 e T ,e T ,
⋅ ⋅ ⋅,
{ } { 2 } uu∗ = σ h2 uu∗
K h = E h h* = E a0
R =σh u
(
2
*
F = Rσ I +R R
)
−1 *
R =
1
σ 2 / σ h2 + N
u u∗
T 2π − j ( N −1)τ 0 e T
+
Q (τ ) = y*G (τ ) FG * (τ ) y =
y*G (τ )uu∗G * (τ ) y
σ 2 / σ h2 + N
Ignoring the constant in the denominator ∗ *
Q (τ ) = u G (τ ) y
2
Now 2π − j 2π τ 0 − j 2π 2τ 0 − j ( N −1)τ 0 u∗G * (τ ) = 1, e T , e T , ⋅ ⋅ ⋅,e T
u∗G * (τ ) y = a0
N −1 j 2π k (τ −τ 0 ) e T
∑
k =0
+
N −1
∑ nk
2 2π j kτ e T
k =0
2π sin N (τ − τ 0 ) T + η (τ ) u∗G * (τ ) y = a0 2π sin (τ − τ 0 ) T
2
COMPARISON OF ML AND MODE • For the multipath channel model the MODE technique amounts to an approximate maximum likelihood estimator for all the amplitudes and time delays. • Apparently a high dimensional search problem, but MODE is an approximate ML computation. Reduces computational complexity significantly, yet gives very good performance. • In simulating the MODE technique the number of desired paths to be estimated for best average performance was 5. • The first step after applying MODE was to reject any delays that fell outside a chosen time window and then choose the smallest delay within the acceptable window as the TOA estimate. In the simulations presented here, the time window was within 100ns from the expected signal arrival time.
• Fig. 1 - Each randomly chosen channel has 5 impulses. • Time delays between successive impulses are independent poisson with mean value 50ns. • Amplitudes are complex gaussian random variables with exponentially decaying RMS values. RMS amplitude of the last impulse is 20% of the first. • For ML technique, 50 channels randomly chosen to estimate correlation matrix • For both techniques 500 channels randomly chosen to determine error statistics for TOA estimate . • The ordinate in Fig. 1 contains values of TOA error, and the abscissa is the percent of random channels for which the TOA is less than or equal to the ordinate. The Signal-to-noise ratio (SNR) is 20dB.
Results 20 18 16
TOA error (ns)
14 12 ML
10
MODE
8 6 4 2 0 0
20
40
60
80
100
Percentile
Fig. 1. 5 impulses and mean separation 50ns
Fig. 2 is similar to Fig. 1 with the exception that the channel now contains 50 random impulses with a mean time separation of 5ns instead of 50ns.
40 35
TOA error (ns)
30 25 ML
20
MODE
15 10 5 0 0
20
40
60
80
100
Percentile
Fig. 2. 50 impulses and mean separation 5ns
Trend: ML technique usually showed better performance than MODE, and the difference in performance was greater When the channel was more complex, i.e. more delay paths.
In some cases MODE did better when the number of paths was small, but the difference in performance in those cases was slight.
COMPARISON USING RAY-TRACE DATA
a)
b)
Fig. 3. a) Overhead Representation of Building Model for Ray-Tracing Environment (Blue represents walls, shelves, and other RF-influencing materials;Green crosses represent RF transmitter positions; Red dots represent simulated mobile positions b) Three-Dimensional View.
Building layout modeled after typical supermarket with approximate Dimensions of seventy-five by one hundred feet, with shelves of varying heights (three to eight feet) and outer walls six meters in height. RF-influencing materials incorporated into the model categorized by a relative dielectric constant. This parameter was set to a value of ten for shelving and six for outer walls. Five access points 5.5 meters in height. Multiple mobile terminal positions 1.5 meters in height. Ray-traced impulse response profiles used in testing the algorithm. For each transmitter, impulse responses were generated by ray-tracing for roughly 500 receiver locations.
20 18 16
TOA error (ns)
14 12 ML
10
MODE
8 6 4 2 0 0
20
40
60
80
100
Percentile
Fig. 4. TOA error with Ray Trace Data
Fig. 4 - The ML correlation matrix estimate was obtained by averaging 100 randomly selected channels, and 500 channels were again used to compile performance statistics. The SNR is again 20dB here. Performance in this case is nearly identical. The reason for this is that the ray-tracing data for this particular scenario showed that most of the channel impulse responses had a small number of significant paths. As indicated earlier, the two techniques are comparable in this kind of environment.
Direct Location Estimation recall
y = GD h + n
ith BS : ( xi , yi , h) MS : ( x, y,0)
2π − j ( N −1)τ 0 − j 2Tπ τ 0 − j 2Tπ 2τ 0 G = diag1, e ,e , .. . , e T
= diag{1, e
−j
2π τ 0 ( x, y ) T
,e
−j
2π 2τ 0 ( x , y ) T
,..., e
−j
2π ( N −1)τ 0 ( x , y ) T
τ 0 ( x, y ) = d / C = ( x − xi ) 2 + ( y − yi ) 2 + h 2 / C
}
• Likelihood Function M
(1)
(2)
p ( y , y ,L , y
(M )
1 ( i )* ( i ) −1 ( i ) ( x, y )) = ∏ N exp( − y ( K y ) y ) (i ) i =1 π Det (K y ) M
1
=
π
M
NM
(i ) Det ( K ∏ y )
exp(−∑ y ( i )* (K (yi ) ) −1 y ( i ) ) i =1
i =1
( xˆ , yˆ ) = arg max log{ p(y (1) , L, y ( M ) | ( x, y ))} ( x, y ) M
= arg max{− ∑ y ( i ) (K (yi ) ) −1 y (i ) } ( x, y )
i =1
*
• Cost Function (2-dim.) M
( xˆ, yˆ ) = arg max{∑ ( x, y )
i =1
1
(i )*
(i ) (i ) (i ) * (i ) y G ( x , y ) F G ( x , y ) y } 2
σi
where
F (i ) = DR( i ) (σi2I + R ( i )*D*DR( i ) )−1 R (i )* D* K h(i ) = R (i ) (R (i ) )*
TOA-LS Method • Lease Square Method ( xi , yi , h) given: AP coordinate τi TOA estimation goal: estimate MS location (x, y)
di = ( x − xi ) + ( y − yi ) + h = τ i ⋅ C 2
2
2
2
2
2
2
xi yi x y − xi x − y i y + + + + + h 2 = (τ i C ) 2 2 2 2 2 R2 2
• In matrix form
x1 x 2 M xM
y1 y2 M yM
− 0.5 ( x12 + y12 ) / 2 − (τ 1C ) 2 − h 2 x 2 2 2 2 − 0.5 ( x2 + y2 ) / 2 − (τ 2C ) − h y = M M R 2 2 2 2 − 0.5 ( xM + yM ) / 2 − (τ M C ) − h
Aθ = b
• Lease Square Solution
ˆθ = (ATA)−1ATb
TOA-WCLS Solution • Weighted Constraint LS (ref [1]) – Constrained Optimization Problem
θˆ = arg min ( Aθ − b) T W(Aθ − b) θ
subject to
q θ + θ Pθ = 0 T
where
T
0 q = 0 − 1
R = x2 + y2 1 0 0 P = 0 1 0 0 0 0
• Lagrangian L(λ , θ) = ( Aθ − b) T W(Aθ − b) + λ (q T θ + θ T Pθ) – Fix λ , solve for minimizer θˆ ˆθ = ( A T WA + λP ) −1 ( A T Wb − λ q) 2
∂L(λ , θ) =0 ∂θ
– Solve for
λ
from
∂L(λ , θ) =0 ∂θ θ =θˆ
Results • Five BS geometry is considered in a square. • BS coordinates [5,10]m, [5,50]m, [80,20]m, [10,75]m, and [90,90]m and MS coordinates [x,y] = [20,20]m. • The two-dimensional maximization for the direct method is performed using a built-in MATLAB function fminsearch in 100m by 100m square. • Channel impulse responses for each of the five BS’s are generated randomly, with the time delays representing a poisson process with average time between points equal to 10 ns and an exponential power delay profile with the power decaying by the fraction .03 over a 200ns delay spread.
Results • Average SNR =
1 M
M
∑ SNR i =1
σ is chosen such that 2 i
• Weighting Matrix W wii =
SNRi M
∑ SNR
i
i =1
i
SNRi ⋅ d i2 = k
Results (cont.) Mean S quare Error 90 TOA-LS TOA-WCLS Direct Es t.
80 70 60 50 E S M
40 30 20 10 0
5
10
15
20 25 Average S NR (dB)
30
35
40
Results (cont.) Mean S quare Error 20 TOA-LS TOA-WCLS Direct Es t.
18 16 14 ) B d( E S M
12 10 8 6 4 2
5
10
15
20 25 Average S NR (dB)
30
35
40
Results (cont.) 90% P e rcentile Error 14 TOA-LS TOA-WCLS Direct Es t.
12
) r et e m ( r or r E n oi t a m it s E
10
8
6
4
2
5
10
15
20 25 Average S NR (dB)
30
35
40
Results (cont.) • Direct Estimation Cost Function x 10
5
12 10 8 6 4 100 2 0 20
50 40 60 80 100
0
Results (cont.)
Direct Es timation (Average S NR=20dB) 100 90 80 70
(x,y)=(20,20)
60 50 Y 40
r=1.40348m (50% err) center=(20.0945,19.7638 )
30 20 10 0
0
10
20
30
40
50 X
60
70
80
90
100
Results (cont.) TOA-WCLS (Average S NR=20dB) 100 90 80 70 60 50 Y 40 30
r=1.8149m (50% error) Center=(19.9565,19.5865 )
20 10 0
0
10
20
30
40
50 X
60
70
80
90
100
Results (cont.) TOA-LS (Average S NR=20dB)
100 90 80 70 60 50 Y 40 30
r=2.1878m (50% error) Center=(19.7311,19.8105 )
20 10 0 0
10
20
30
40
50 X
60
70
80
90
100
