September 2017 doc.: IEEE 802.

11-17/1387R0

IEEE P802.11
Wireless LANs

High-Accuracy Indoor Geolocation using Collaborative Time of

Arrival (CToA) - Whitepaper
Date: 2017-09-07

Author(s):
Name Affiliation Address Phone email
94 Em Hamoshavot Rd.
Leor Banin Intel Corporation leor.banin@intel.com
Petah Tikva, Israel, 49527
94 Em Hamoshavot Rd.
Ofer Bar-Shalom Intel Corporation ofer.bar-shalom@intel.com
Petah Tikva, Israel, 49527
94 Em Hamoshavot Rd.
Nir Dvorecki Intel Corporation nir.dvorecki@intel.com
Petah Tikva, Israel, 49527
94 Em Hamoshavot Rd.
Yuval Amizur Intel Corporation yuval.amizur@intel.com
Petah Tikva, Israel, 49527

Abstract
Collaborative time of arrival (CToA) is the next generation, indoor geolocation method, which is
designed for enabling scalability of the existing IEEE802.11/Wi-Fi-based, geolocation systems. The
technique leverages on the IEEE802.11 fine timing measurements (FTM) capabilities, enabled in state-
of-the-art Wi-Fi chipsets, and supports two concurrent operation modes; the CToA “client-mode”
enables “GPS-like” operation indoors, and allows an unlimited number of clients to privately estimate
their position and navigate indoors, without exposing their presence to the network. The CToA “network-
mode” is designed for large-scale asset-tracking applications, and enables a centric positioning server to
pinpoint objects equipped with wireless, Wi-Fi-based, low-power electronic tags (e-Tags).

The CToA method is broadcast-based and operates over an un-managed network, built out of cheap,
unsynchronized units called “CToA broadcasting stations” (bSTA). The bSTAs, which are stationed at
known locations, periodically broadcast a unique beacon transmission and publish its time of departure
(ToD). Neighbor bSTA units and clients that receive the beacon broadcast, measure and log its time of
arrival (ToA). Every bSTA publishes its most recent timing measurement log as part of its next beacon
broadcast. CToA clients combine their own ToA measurements with those published by the bSTAs, in
order to estimate and track their location.
CToA e-Tag clients act similar to bSTAs, and simply wake-up sporadically to broadcast a CToA beacon.
The ToA of that broadcast is measured and logged by the receiving bSTAs similarly to beacons broadcast
by other bSTAs. The timing measurement report is then delivered to a centric positioning server that can
estimate and track the location of numerous CToA-based e-Tags, simultaneously.

The paper outlines the principles of the CToA method and the mathematical background of the position
estimation algorithms. In addition, performance examples as well as theoretical analysis of the expected
positioning accuracy are provided.

Submission page 1 Banin et al., Intel Corporation

 1

High-Accuracy Indoor Geolocation using

Collaborative Time of Arrival
Leor Banin, Ofer Bar-Shalom, Nir Dvorecki, and Yuval Amizur

Abstract—Collaborative time of arrival (CToA) is the next WLAN-based location technology relies, to a great extent, on
generation, indoor geolocation protocol, which is designed for the WLAN received signal strength indicator (RSSI) infras-
enabling scalability of the existing IEEE802.11/Wi-Fi-based, ge- tructure. The RSSI is a measure of the RF energy received by
olocation systems. The protocol leverages on the IEEE802.11 fine
timing measurements (FTM) capabilities, enabled in state-of- the station. WLAN stations estimate the RSSI of the beacons
the-art Wi-Fi chipsets, and supports two concurrent operation broadcast by access points (AP), and use this metric to sort
modes; the CToA “client-mode” enables “GPS-like” operation between the APs based on their signal quality and proximity.
indoors, and allows an unlimited number of clients to privately The RSSI metric is measured in units of [dBm], and in general
estimate their position and navigate indoors, without exposing is inversely proportional to the logarithm of the squared
their presence to the network. The CToA “network-mode” is
designed for large-scale asset-tracking applications, and enables distance between the transmitter and the receiver [16]. RSSI-
a centric positioning server to pinpoint objects equipped with based mobile device positioning exists in two main flavors:
wireless, Wi-Fi-based, low-power electronic tags (e-Tags). path-loss models, and “fingerprinting”. Path-loss models relate
The CToA protocol is a broadcast-based protocol that operates the received signal power to the propagation distance. A set
over an un-managed network, built out of cheap, unsynchronized of RSSI measurements obtained from different WLAN APs
units called “CToA broadcasting stations” (bSTA). The bSTAs,
which are stationed at known locations, periodically broadcast in the vicinity of the client station, enable it to estimate its
a unique beacon transmission and publish its time of departure position via trilateration methods [18], [19]. While this method
(ToD). Neighbor bSTA units and clients that receive the beacon is relatively simple to implement, it is prone to yield fairly
broadcast, measure and log its time of arrival (ToA). Every bSTA inaccurate positioning results due to the large variations in the
publishes its most recent timing measurement log as part of its RSSI measurements [6]. The alternative approach is to corre-
next beacon broadcast. CToA clients combine their own ToA
measurements with those published by the bSTAs, in order to late the RSSI measurements against a pre-calibrated database
estimate and track their location. CToA e-Tag clients act similar of RSSI “fingerprints”, measured over a pre-defined grid and
to bSTAs, and simply wake-up sporadically to broadcast a CToA stored in a server. The fingerprint approach provides better
beacon. The ToA of that broadcast is measured and logged by the accuracy compared to the path-loss based RSSI. However, as
receiving bSTAs similarly to beacons broadcast by other bSTAs. the method’s accuracy is sensitive to even minor changes in
The timing measurement report is then delivered to a centric
positioning server that can estimate and track the location of the propagation channel (e.g., a placement of a new sales-stand
numerous CToA-based e-Tags, simultaneously. in a shopping mall), frequent re-calibrations and updates of
The paper outlines the principles of the CToA protocol and the the fingerprints database are required. The high-maintenance
mathematical background of the position estimation algorithms. incurred by this type of positioning systems obviously limits
In addition, performance examples as well as theoretical analysis their scalability.
of the expected positioning accuracy are provided.
Index Terms—geolocation, Indoor navigation, fine timing mea- Facing the limited positioning accuracy enabled by
surement, FTM, time delay estimation, Maximum likelihood
RSSI/path-loss based location technologies and the limited
estimation, WLAN, Wi-Fi, IEEE 802.11
scalability of fingerprint-based systems, industry vendors be-
gan seeking alternative WLAN-based positioning technolo-
I. I NTRODUCTION gies, which will enable to achieve higher positioning accuracy.
HE challenge of accurate indoor location and navigation
T has been attracting an increasing amount of attention
since the mid 1990’s. Cultivated by the cellular revolution and
Taking advantage of the high bandwidth supported by the
WLAN systems (ranging between 20-160 MHz), the approach
pursued was geolocation based on time-delay estimation [11].
the U.S. federal communication committee (FCC) enhanced Though the early releases of the IEEE802.11TM standard
911 services (E911) [1], indoor location has ignited a rapid included means for time delay estimation, the timing resolution
development of mobile location technologies. The ubiquity of enabled by these mechanisms was in the microseconds level
IEEE802.11TM wireless local area network (WLAN) technol- - too coarse for any practical indoor positioning applications.
ogy in mobile devices, which to date, has already reached an High-accuracy positioning in a dense multipath environment
attach-rate of 100% in the smart-device segment [2], facilitated imposed several hardware design changes in the existing
the development of WLAN-based indoor location systems. WLAN chipsets, in order to increase the timing resolution
Due to the lack of standard infrastructure for high-resolution from the microseconds level to the nanosecond level (or
timing measurement capabilities in its early releases, existing even sub-nanosecond level). The solution that was endorsed
by the IEEE802.11TM group, was a novel time-delay based
L. Banin, O. Bar-Shalom, N. Dvorecki and Y. Amizur are with Intel’s
Location Core Division, 94 Em Hamoshavot Rd., Petah Tikva 49527, Israel. ranging protocol called “fine-timing measurement” (FTM).
Corresponding author’s e-mail: oferbarshalom@gmail.com. The FTM protocol enables a WLAN station to measure its
 2

distance w.r.t. another station1 . The range measurement is combines these parameters along with its own estimated ToA,
based on high-resolution, time delay estimation, which also t2 , and measured ToD, t3 values, to obtain a range estimate
accounts for the latency imposed by the chipset hardware. The w.r.t. the responding station3 . The FTM protocol has first
hardware-imposed latency (e.g. the receive/transmit filters’ appeared in the 2016 release of the IEEE802.11TM standard
group-delay and other hardware latencies), is measured and [23] (formerly known as IEEE802.11REVmc). Followed the
pre-calibrated by the chipset in order to reach the required standard release, the Wi-Fi alliance - the organization that
timing resolution. Obtaining an accurate time delay estimate promotes Wi-Fi technology - has announced in February 2017
in a dense-multipath environment is typically implemented on a Wi-Fi location certification program to certify WLAN
using some super-resolution method , which are applied to the devices complying with the FTM protocol.
estimated channel response, [9], [10]. FTM is a point-to-point Being a P2P, single-user protocol, the FTM protocol is lim-
(P2P), single-user protocol, which includes an exchange of ited in scenarios where an extremely large number of users are
multiple message frames between an initiating WLAN station requesting positioning services simultaneously. Provided no
(STA) and a responding STA. The initiating STA attempts to FTM message transactions are lost on the way due to temporal
measure its range w.r.t. the responding station (e.g., WLAN AP channel interruptions, the initiating station should be able to
or a dedicated FTM responder). The FTM message sequence obtain a range estimate w.r.t. the responding station within
chart is illustrated in Fig. 1. The time of flight (ToF) between about 30ms. Hence, obtaining its position, which involves
ranges estimation towards 3 additional stations, should ideally
take about 100-120ms. This implies that each FTM responder
may be able to serve about 30 client stations per second.
Clearly, with more and more navigating stations attempting
/ŶŝƚŝĂƚŝŶŐ ZĞƐƉŽĚŝŶŐ
to execute FTM sessions, the collision likelihood increases,
^d ^d
which effectively decreases the number of stations that can be
serviced. Consider for example, large stadiums hosting rock
concerts or major sports events. In such occasions it is easy
to imagine tens of thousands of users navigating throughout
the stadium area using location-based services. Servicing all
these users might require to deploy a network of thousands of
”ϭϬŵƐ;ƌĞĐŽŵŵĞŶĚĞĚͿ FTM responders around the stadium. The protocol described
ƚϭ in the sequel, dubbed “collaborative time of arrival” (CToA),
ƚϮ is aimed to provide a more cost-effective solution for such use
ƚϯ
cases.
ƚϰ aϮϬŵƐ
Paper Organization: The remainder of the paper is orga-
nized as follows. Section II gives an overview of the CToA
protocol and its challenges. The mathematical model for the
positioning problem considered, is formulated in section III.
This section is divided into three parts; the first part, ad-
dressed in section III-A, outlines the measurement models
and maximum-likelihood position estimators for client-mode
CToA in the absence of clock-skewness. These measurement
models are then used in section III-B for obtaining measures
Fig. 1. FTM Protocol Message Flow Example for the expected positioning accuracy. The third part, which
is outlined in section III-C, introduces the effect of the clock-
the two stations is calculated using (1), drifts on the measurement models. This section also details
(t4 − t1 ) − (t3 − t2 ) the Kalman filter algorithm that is executed by the client
ToF = (1) device and used for estimating and tracking all the time-
2
where t1 denotes the time of departure (ToD) measured by varying parameters in the system. System-level simulation
responding station, and t4 denotes the time of arrival (ToA), performance are described in section IV. In the appendix
which is estimated by the responding station. The values of t1 we derive the Cramér-Rao lower bounds for the positioning
and t4 are reported back to the initiating station2 after the com- problem and the associated concentration ellipses, both of
pletion of the FTM measurement phase. The initiating station which are used in the main text to illustrate the theoretical
system performance discussed in section III-B.
1 Notice that FTM only enables to measure the range between two stations.
Notation: We use lower-case letters to denote scalars,
Obtaining a position estimate based on multiple range measurements is out
of the standard scope. However, the standard does define mechanisms for the
lower-case, boldface letters to denote vectors, and upper-
responding stations to provide their location information (such as, absolute or case, boldface letters to denote matrices. We further
relative position coordinates, floor level etc.), in an information element (IE),
called location configuration information (LCI). The LCI of the responding 3 The exchange of the FTM measurement message and its acknowledgement
stations may be used by the initiating station to estimate its absolute or relative (ACK) frame, which has to be sent out after exactly a short inter-frame spacing
position. (SIFS) of 16µs, is assumed to finish within a short period, during which the
2 The values of t and t are reported at picosecond granularity. clocks of the two stations do not drift appreciably.
1 4
 3

use the following nomenclature throughout the paper. • Network-Mode CToA - designed to enable a network
{·}T transpose administrator to simultaneously track the position of a
diag{x} N × N diagonal matrix, large number of clients. This mode is useful for large
whose diagonal is the vector x scale asset tracking (e.g., using eTags4 ), fleet manage-
IN N × N identity matrix ment, law-enforcement, etc. CToA clients in operating
1N N × 1 vector of 1’s in the network-mode do not listen for CToA beacons,
0N N × 1 vector of 0’s but only transmit CToA beacons (at rather low rate), in
⊗ Kronecker product order to enable the network administrator to track their
∥x∥ norm √∑ x,
√ of the vector position. The sporadic, short transmissions executed by
i.e. xT x = 2
i xi such devices enables them to operate for long periods
argmin∥y(x)∥ search for the value of x that minimizes using small, coin-cell batteries.
x
the norm of y(x) As bSTAs activity in the two modes is identical, the CToA
n ∼ N (µ, σ) Gaussian-distributed noise with network can support both modes simultaneously.
mean, µ, and standard-deviation, σ. Since the protocol is based on broadcasting of CToA
E{.} Expectation operator beacons, it is uni-directional; the CToA beacons are unac-
c The speed of light, knowledged if not received. Yet, the fact that beacons can
c = 2.99792458 · 108 m/s. be received by multiple neighbor bSTAs in the vicinity of the
broadcasting bSTA, gives a level of redundancy and immunity
II. S CALING - UP THE FTM P ROTOCOL against frame losses. Each CToA beacon is associated with a
A. CToA Overview unique packet identification index (PID), which is assigned to
CToA is a geolocation protocol designed for scaling up the it by the broadcasting unit. The PID is typically implemented
number of clients that could be serviced simultaneously. This as a running counter, and is independently maintained by every
can be achieved through the use of broadcast approach rather bSTA. Each CToA beacon has a ToD time-stamp (measured by
than a P2P or a point to multi-point ranging approach. The pro- the broadcasting unit), and multiple ToA measurements - all of
tocol operates over an un-managed network of unsynchronized which are associated with the same PID. The PID enables the
and independent units called “CToA broadcasting stations” CToA clients (operating in “client-mode”), or the positioning
(bSTA), which together form a high-precision, geolocation server (in “network-mode”) to associate between the ToD and
network. The bSTAs, which are deployed at known locations, its corresponding ToA measurements, collected either by the
are implemented using either standard WLAN APs that have client itself or different bSTAs.
the ability to measure accurate ToA or network-detached, In addition to the ToD, some of the CToA beacons broadcast
FTM-responders. by every bSTA include also a data log of the timing mea-
According to the CToA protocol, the bSTA units serve several surements collected by the bSTA during the past n-seconds.
purposes. Every bSTA: This data log is called “CToA location measurement report”
(CLMR). The timing measurements included in the CLMR
• Periodically broadcasts a CToA “beacon” and measures
reports are used by the cSTA, which combines them with
the ToD of that beacon.
its own ToA measurements in order to estimate its position.
• Listens for CToA beacons broadcast by its neighbor
Although the CLMR logs are maintained by each of the
bSTAs, and measures their ToA.
bSTAs independently, the protocol also enables the CLMR
• Maintains a log of its current ToD and ToA measure-
logs broadcast by one bSTA to be aggregated by its neighbors,
ments, and publishes its most recent measurements log
thereby providing an immunity mechanism against “hidden-
as part of its next CToA beacon broadcast.
nodes” in the wireless network.
• Periodically announces its location as part of its CToA
Figures 2-3 illustrate an example of CToA beacon broadcast
beacon broadcasts.
and its reception. The example assumes a CToA network,
The CToA protocol supports two modes: a client-mode and which consists of 3 bSTA units and a single cSTA. These units
a network-mode. While both modes rely on the same protocol are assigned with (simplified) medium-access control (MAC)
principles, they are targeted towards different usage models: addresses: 10:01, 10:02 and 10:03, while the cSTA has a MAC
• Client-Mode CToA - may be visioned as the indoor address of 55:55. As depicted in Fig. 2, at some point indicated
counterpart of the global navigation satellite systems by ToD time-stamp of “199678” (measured in picoseconds and
(GNSS). It is designed for enabling an unlimited number referenced to the time-base of bSTA#1), bSTA#1 broadcasts
of clients to estimate their location and navigate, simul- a CToA beacon associated with PID “1551”. The ToD and
taneously, while maintaining their privacy. The CToA the PID are logged in the CLMR log maintained by bSTA#1.
client stations (cSTAs) only listen to the bSTA broadcasts. This log also includes the MAC address of the broadcasting
Once a cSTA receives a broadcast, it measures its ToA bSTA. The beacon propagates, and as illustrated in Fig. 3, it is
and combines it with the ToD/ToA measurements log received by the neighbor bSTAs (#2 & #3), and by the cSTA
published by the bSTA in the CToA beacons, in order
4 An “e-Tag” refers to an electronic tag - a small, wireless-enabled device
to determine its position. Since the cSTAs do not trans-
that is attached to a larger object, and enables to remotely monitor various
mit, their presence is not exposed and their privacy is measures related to the object including its location. In the context of this
maintained. paper the e-Tag is assumed to be Wi-Fi enabled.
 4

- all of which update their CLMR logs accordingly: bSTA#2 broadcasts), will be fed into a Kalman filter that will produce
measured the ToA of the beacon to be “329673” (referenced an updated estimate of the cSTA’s current position.
to its own time-base) and updates that value in its CLMR log,
along with it MAC address as the receiving unit. Similarly do B. The CToA Protocol
bSTA#3 and the cSTA, which estimate the ToAs of that same
The CToA protocol follows the principles of the channel
beacon to be “341006”, and “133564”, respectively.
sounding mechanism introduced for IEEE802.11ac standard
(a.k.a very high-throughput/VHT) [23]. The channel sounding
dŽĞĂĐŽŶƌŽĂĚĐĂƐƚ protocol, which was originally proposed for determining the
optimal beamforming weights at the transmitter side [17],
relies on a transmission of a null-data packet (NDP), which
consists of only a known sequence of OFDM symbols, but
with no data payload. The transmission of the NDP is pre-
ceded by a packet called NDP announcement (NDPA), which
W/ dŽ dyD dŽ ZyD informs the receivers of an NDP frame that is about to be
ϭϱϱϭ ϭϵϵϲϳϴ ϭϬ͗Ϭϭ
transmitted after a standard, short inter-frame spacing/SIFS
of 16µs from the end of the NDPA packet. The NDPA
contains information for the receiver for estimating the channel
ď^dηϭ response.
DϭϬ͗Ϭϭ CToA relies on a similar protocol structure; the CToA
beacons broadcast by the bSTAs consist of an NDPA frame
announcing the upcoming transmission of the NDP frame,
ď^dηϮ which is used by the receivers for measuring its ToA. As
ď^dηϯ
DϭϬ͗ϬϮ mentioned above, the NDPA also includes data that enables
DϭϬ͗Ϭϯ
the receiving clients to estimate their location. The protocol
messaging sequence for the client-mode is illustrated in Fig. 4.
dŽůŝĞŶƚ;Đ^dͿ
Dϱϱ͗ϱϱ As shown, cSTAs in this mode only receive, but do not trans-
mit any data, maintains their privacy. The protocol messaging
Fig. 2. CToA Beacon Broadcast Example

dŽĞĂĐŽŶZĞĐĞƉƚŝŽŶ Đ^d;ĐůŝĞŶƚͿ ď^dŶ ď^dũ ď^dŵ

^/&^
dŽŶŝ

dŽŝ
dŽũŝ

^/&^
dŽũů
ď^dηϭ W/ dŽ dyD dŽ ZyD
W/ dŽ dyD dŽ ZyD
DϭϬ͗Ϭϭ ϭϱϱϭ ϭϵϵϲϳϴ ϭϬ͗Ϭϭ ϯϮϵϲϳϯ ϭϬ͗ϬϮ dŽů dŽŶů dŽŵů
ϭϱϱϭ ϭϵϵϲϳϴ ϭϬ͗Ϭϭ ϯϰϭϬϬϲ ϭϬ͗Ϭϯ

W/ dŽ dyD dŽ ZyD

ď^dηϮ Fig. 4. The CToA Protocol - Client-Centric Operation
ϭϱϱϭ ϭϵϵϲϳϴ ϭϬ͗Ϭϭ ϭϯϯϱϲϰ ϱϱ͗ϱϱ
ď^dηϯ DϭϬ͗ϬϮ
DϭϬ͗Ϭϯ sequence for the network-mode CToA is illustrated in Fig. 5.
dŽůŝĞŶƚ;Đ^dͿ
Dϱϱ͗ϱϱ C. The Un-managed Operation of the CToA Network
Fig. 3. CToA Beacon Reception Example Besides using unsynchronized broadcasting units, the CToA
network is also un-managed in the sense no that coordination
The CLMR logged by each of the bSTAs will be broadcast or scheduling protocol between the bSTAs (or the receiving
as part of their next CToA beacon broadcast event, and the clients) is required for its operation. Ideally, if the bSTAs are
CLMR logs will be updated accordingly. As will be explained implemented as dedicated units, the CToA network could be
in sectionIII-C, the entries of the CLMR log table will undergo allocated with a unique operation channel (with bandwidth of
a time-stamp matching step, in which the ToA measurements 20MHz, 40MHz or 80MHz, depending on the spectrum band
will be matched with their corresponding ToD value. This in which the system operates). Any bSTA can first scan the
information, combined with the position information of the spectrum to detect prior CToA broadcast activity, and once
bSTAs, (which is also provided as part of their CToA beacon detected - the bSTA can contend on accessing that channel
 5

WŽƐŝƚŝŽŶŝŶŐ^ĞƌǀĞƌ each of them is set to operate at a different channel (aka, its

“native channel”). To enable an un-managed CToA network to
coexist with a live Wi-Fi data network, the scheme illustrated
in Fig. 6 may be used. According to this scheme, the AP/bSTA
periodically (e.g., once every few minutes) scans the spectrum
to detect CToA broadcast activity. Following the scan results,
the AP/bSTA announces to its associated STAs, on a short
Đ^d;ĞͲdĂŐͿ ď^dŶ ď^dũ ď^dŵ “un-availability window” (such window should typically last 1
^/&^ millisecond or less). During that window, the bSTA/AP hops
between the native channels of its neighbor bSTA/APs. On
^/&^ dŽũŝ each of these channels it broadcasts a short CToA beacon,
dŽdžŬ
which only includes its ToD but lacks the CLMR. When the
dŽŶdžŬ
dŽũdžŬ dŽŵdžŬ
AP returns to its native channel, it broadcasts a longer CToA
beacon that includes both its ToD and its recent CLMR.
This process is depicted in Fig. 6. This example illustrates
^/&^ the frequency hopping process of two AP-based bSTAs (#1
dŽũů and #2), in a CToA network that consists of a total of 4 AP-
dŽŶů dŽŵů based bSTAs. As shown, while bSTA#1 goes off of its native
channel, it hops between the channels occupied by AP/bSTAs
#2,#3 and #4, and on each channel it broadcasts a short bea-
Fig. 5. The CToA Protocol - Network-Centric Operation con. The broadcasting process, including the standard WLAN
channel arbitration, lasts about 200µs while the beacon itself,
>ĞŐĞŶĚ͗ (which consists of NDPA-SIFS-NDP) consumes about 100µs.
WŶŽŶͲĂǀĂŝůĂďŝůŝƚLJǁŝŶĚŽǁŽŶ͞ŶĂƚŝǀĞ͟ĐŚĂŶŶĞů
When bSTA#1 returns to its native channel, it broadcasts a
͞EĂƚŝǀĞ͟WĐŚĂŶŶĞůʹdŽZyнWZy
longer beacon that includes also the CLMR log. A similar
dŽĞĂĐŽŶďƌŽĂĚĐĂƐƚďLJď^dηϭ
process is executed by bSTA#2.
dŽĞĂĐŽŶďƌŽĂĚĐĂƐƚďLJď^dηϮ By scanning the medium, the CToA clients can detect the
ŚĂŶŶĞůĂƌďŝƚƌĂƚŝŽŶ CToA activity and infer the broadcast periodicity of the bSTAs.
ĨƌĞƋƵĞŶĐLJ
Once the cSTA figures out this information, it simply needs
to hop between the native channels used by the bSTA/APs,
and collect the CLMRs broadcast on each of these channels.
Under the assumption of 5Hz beacon broadcasting rate, the
ď^dηϮ

ď^dηϮ

unavailability of the AP for its associated STAs is about

0.5% or less. Accordingly, within 200 milliseconds the client
should be able to gather sufficient information for estimating
ď^dηϯ

ď^dηϯ

its position.
ď^dηϰ

ď^dηϰ

D. The CToA Time-Tracking Challenge

The general approach described above, of synchronizing a
network using time-stamped broadcast transmissions, is well-
known [20]. However, whilst most network-based applications
ď^dηϭ

ď^dηϭ

(e.g., audio or video distribution) would be satisfied with

microsecond or even millisecond-level synchronization, ac-
ƚŝŵĞ curate geolocation applications require sub nanosecond-level
accuracy. Attaining such a high timing accuracy in a network
ϴϬϮ͘ϭϭĂnjEW
;dŽͲŽŶůLJͿ
^/&^ EW built out of unsynchronized, and independent broadcasting
ƚŝŵĞ units, which rely on low accuracy oscillators as their clock
source, is extremely challenging.
Fig. 6. bSTA Frequency Management with AP-based bSTAs To understand the challenge, let us compare the problem of
client geolocation in a CToA network to a client geolocation in
a GNSS network. GNSS networks implement a similar time-
for broadcasting its beacons, and listen to broadcasts of its stamped broadcasts approach for enabling an unlimited num-
neighbor during the remaining time. ber of client receivers to navigate simultaneously, worldwide.
A more challenging case is when the bSTA functionality Yet, there are two fundamental differences between the time-
is implemented as part of a standard Wi-Fi AP. In such case, tracking of receivers in a GNSS network compared to CToA
the AP is obligated to provide data transaction services to its network.
associated STAs. To ensure high-data throughput and capacity, 1) The GNSS network is fully synchronized, whereas the
the Wi-Fi network typically consists of a grid of APs, where CToA network is not.
 6

2) The broadcasts in a GNSS network are received simul- under the idealistic assumption that the bSTA clocks do not
taneously, while in a CToA network the broadcasts are drift over time, such that their offsets w.r.t. to the client’s
staggered in time. clock are time-invariant. Under this assumption we derive the
Let’s delve into these two differences and explore them in a position estimators for two cases:
st
bit more details. In GNSS networks the satellite vehicles (SV) • ”1 Fix” - corresponds to the scenario, at which the client
are synchronized using on-board atomic clocks, which have a first attempts to sync and estimate its position.
frequency stability of approx. 10−14 . This frequency stability • “Tracking” - corresponds to the scenario where the client
translates into a clock-drift of roughly 1 nanosecond per day already has an estimate of its position and the bSTA
[21], (which is equivalent to a ranging error of about 35cm - timing-related parameters, and continues tracking them
an error, which is further corrected by the GNSS system). using a Kalman filter.
Since GNSS networks are fully-synchronized, in terms of These measurement models are used for developing approxi-
timing parameters, the GNSS client receiver needs to estimate mate performance bounds that can predict the expected posi-
only the offset and the drift between its internal clock, and tioning accuracy. Then, in section III-C, we define the Kalman
the GNSS network clock. The GNSS client receiver’s clock filter that enables the client to simultaneously estimate and
is typically generated using a crystal oscillator/XO with a track its own location coordinates, as well as the clocks
frequency stability in the order of 10−6 , which is commonly parameters the bSTAs, (both of which it receives directly, as
expressed in units of parts-per-million/ppm). Tracking these well as of bSTAs received indirectly via other bSTA).
parameters (along with additional system states such as po-
sition and velocity), is done using a Kalman filter algorithm
A. CToA MLE Solution in the Absence of Clock Drift
[22].
In the CToA network, since the bSTAs are unsynchro- We shall now derive the maximum-likelihood estimates
nized, each bSTA contributes a clock-offset and drift that (MLE) of the client position under the assumption that the
need to be estimated and tracked. Furthermore, the different clocks of the client and the bSTA do not drift over time, so
medium-access control (MAC) methods used by GNSS and that the clock offsets are fixed. For simplicity, the derivation
CToA impose an additional challenge. In GNSS networks, assumes a horizontal position only, (which is typically of most
the multiplexing at the code space (CDMA) or the frequency interest in indoor-positioning scenarios). The extension to 3-
space (FDMA), ensures that broadcast transmissions from dimensional positioning is straightforward.
all SVs are received simultaneously at the client. On the We define a “measurement” as the time-of-flight (ToF)
other hand, CToA relies on the “listen-before-talk” MAC of a broadcast transmission between two endpoints, A and
of the IEEE802.11TM , which effectively results in timing B. The transmitting endpoint, A, measures the broadcast’s
measurements being staggered in time. Given that a typical time of departure (ToD), while the receiving endpoint, B,
Wi-Fi XO has an accuracy of ±25ppm, consecutive timing measures its time of arrival (ToA). Both timing measurements
measurements taken of the same broadcasting source may are referenced to a 3rd party clock, and thus have offsets
accumulate significant time-drift [14]. This effectively means marked by νA and νB , respectively.
that while one bSTA clock offset is being measured, the other z , ToFAB = (ToAB + νB ) − (ToDA + νA ) (2)
bSTAs clock offsets keep on drifting apart. As an example,
consider two bSTA broadcast timing measurements taken by By denoting the coordinates vectors of the endpoints as qA and
a static receiver, while the beacons are being broadcast at rate qB , respectively, and ignoring any non-line of sight (NLoS)
of 5Hz. The drift of the second ToD time-stamp accumulates to timing biases, the ToF between the two endpoints may be
up to 5µs (w.r.t. its nominal value). This effectively translates expressed as,
into a ranging error of: ≈ 5µs · 3 · 108 m/s = 1500m! Clearly, 1
the clock skewness poses a major challenge for the receiver, ToFAB ∼ = ∥qB − qA ∥ ≡ ToAB − ToDA (3)
c
which requires the application of filtering techniques for
Combining (2) with (3) we obtain the definition of a noiseless
tracking these changes over time. As will be described in the
ToF measurement as,
sequel, the CToA client uses a Kalman filter for tracking the
various timing-related parameters as well as it own location. 1
z̃ , ∥qB − qA ∥ + νB − νA (4)
Assuming the rate at which the clock offsets vary in time is c
slow enough, a beacon broadcasting rate of 3-5Hz by each If the 3rd party also acts as the receiving endpoint, then νB = 0
bSTA is sufficient for the cSTAs to accurately track the clock and the noiseless ToF measurement is defined as,
behavior of the bSTAs.
1
z̄ , ∥qB − qA ∥ − νA (5)
c
III. CT OA P ROBLEM F ORMULATION & P OSITIONING 1) MLE Solution for CToA Client’s First Fix: Assume that a
A LGORITHMS single CToA client station (cSTA), located at: p = [x, y]T , at-
In the following section, the mathematical background of tempts to estimate its position using time-delay measurements
the CToA client position estimation is established. To facilitate it gathers from M CToA bSTAs, whose locations are known to
the explanation we split the derivation into two parts; first, the cSTA, where the mth bSTA is located at qm = [xm , ym ]T .
in section III-A, we address the position estimation problem The cSTA collects two types of time delay measurements:
 7

• bSTAi → bSTAj measurements, where i, j ∈ Let ei denote an M × 1 vector of zeros, whose ith entry is 1.
1 . . . M , ∀i ̸= j . These time-delays are measured by the Using this notation,
bSTAs and published in their CToA beacon broadcast.
The cSTA collects L measurements of this type, where νi = eTi ν (8)
the lth measurement is denoted by z̃l . Each measurement νj − νi = (eTj − eTi )ν (9)
is subjected to additive measurement error denoted by
ñl ∼ N (0, σ̃). where, ν , [ν1 , . . . , νM ]T . Now, the timing measurements
• bSTAi → cSTA, where i ∈ 1 . . . M . These time-delays
may be recast as,
are measured by the cSTA itself, and the cSTA collects L 1
z̄ℓ ∥p − qi ∥ − eTi ν + n̄ℓ
= (10)
measurements of this type, where the ℓth measurement is c
denoted by z̄ℓ . Each measurement is subjected to additive z̃l = ToFij + (eTj − eTi )ν + ñl ,
measurement error denoted by n̄ℓ ∼ N (0, σ̄). Typically, 1
σ̄ > σ̃. ≈ ∥qj − qi ∥ + (eTj − eTi )ν + ñl (11)
c
Let νi denote the (unknown) offset between the cSTA & We may further define the following vectors and matrices,
bSTAi clocks. A single bSTAi → cSTA, ToA measurement
of the ℓth broadcast made by the ith bSTA, may be modeled z̄ , [z̄1 , . . . , z̄L ]T ,
as, z̃ , [z̃ , . . . , z̃ ]T ,
[1 ] L
1 z̄
z̄ℓ = ∥p − qi ∥ − νi + n̄ℓ , ℓ = 1, . . . , L (6) z ,
c z̃
Similarly, the measurement of the lth broadcast that is received d¯ℓ (p) , ∥p − qi ∥
by the jth bSTA, may be modeled as, d˜l , ∥qi − qj ∥
d̄(p) , [d¯1 (p), . . . , d¯L (p)]T
z̃l = ToFij − νi + νj + ñl
, [d˜1 , . . . , d˜L ]T
d̃
= ToAj − ToDi − νi + νj + ñl (7) [ ]
d̄
=
1
∥qj − qi ∥ + βij − νi + νj + ñl , l = 1, . . . , L d(p) , (12)
c d̃

Notice that the ToF (scaled by the propagation speed, c), Ē , [−ei,1 , . . . , −ei,L ]T
may represent a biased version of the range between the Ẽ , [(ej,1 − ei,1 ), . . . , (ej,L − ei,L )]T
[ ]
two bSTAs. This may happen due to some obstruction in the Ē
propagation medium, resulting in a non-line of sight (NLoS) E ,
Ẽ
link between the two endpoints. The scalar βij > 0 represents
n̄ , [n̄1 , . . . , n̄L ]T
the NLoS ranging bias between the ith and jth bSTAs.
ñ , [ñ1 , . . . , ñL ]T
However, since the position of the bSTAs is assumed to be [ ]
known a priori, βij can be easily estimated and eliminated. For n̄
n ,
the sake of derivation clarity, it shall be further assumed that ñ
βij = 0, ∀i, j. An example for timing measurements collected Using the definitions of (12), we may recast (10)-(11) as,
under the assumptions that the bSTA clocks are stable and do
not drift over time is illustrated in Fig. 7. z = c−1 d(p) + Eν + n (13)
Assuming that the measurement noise is Gaussian-distributed
FORFNRIIVHWVWLPHGLIIHUHQFH
with the mean and covariance as follows,
E{n} = 0
E67$% [ 2 ]
σ̄ IL 0
E{nnT } = ,W (14)
0 σ̃ 2 IL
PHDVXUHPHQWQRLVH
Then the maximum likelihood estimate (MLE) of the cSTA
E67$$
position vector, p̂, may be obtained as,
p̂ = argmin(z − c−1 d − Eν)T W−1 (z − c−1 d − Eν) (15)
p,ν

The estimate of the clock offsets vector may be found using

weighted least-squares (WLS) criteria,
F67$
ν̂ = (ET W−1 E)−1 ET W−1 (z − c−1 d) (16)
Define,
Fig. 7. Timing Measurements with Stable Clock B , [W−1 − W−1 E(ET W−1 E)−1 ET W−1 ] (17)
 8

Then, by substituting (16) back in (15) we get, concentration ellipse at that position (namely, 2σ), as the
−1 −1 measure of accuracy. The bound is calculated using,
p̂ = argmin(z − c d) B(z − c
T
d) (18) √
p Bound @ 95% = κ · λmax (23)
The nonlinear minimization problem in (18) can be solved via
where κ is calculated using (77) and λmax is evaluated for the
2-dimensional grid search (or 3-dimensional search, in case of
measurement models (13).
3-positioning), over all the possible locations.
The bound was evaluated over a pre-defined grid covering
2) MLE Solution for a CToA Client in “Tracking” Mode:
the office floor with a resolution of 0.5m×0.5m. The theo-
Once the CToA client receiver has converged to the true values
retical error was computed under the assumption is that at
of the bSTA clock offsets and continuously tracks them, these
each grid point the CToA client receives timing measurements
clock offsets may be considered as “known” (up to some
from only the 4 nearest bSTAs. Fig. 8 illustrates the expected
estimation error). In such case it would be reasonable to
accuracy for a “1st -Fix” scenario, while Fig. 9 corresponds
assume that z ≃ z̄, d ≃ d̄. Define,
to the “tracking” scenario. The measurement noise standard
ζ , z̄ − Ēν̂ (19) deviations assumed in the “1st -Fix” scenario were σ̄ =1.5m,
and σ̃ =0.6m. For the “tracking” scenario the assumed residual
The resulting measurement model in this case may be recast error was σr =0.3m.
as,
ζ = c−1 d̄(p) + ň (20)
The additive noise vector, ň, is assumed to be Gaussian
distributed with the following properties:
E{ň} = 0
√
E{ňň }
T
, σ̄ 2 + σr2 · IL (21)
where σr corresponds to the standard deviation residual esti-
mation error of the clock offsets. The client position MLE in
this case is obtained by minimizing the following cost function
2
p̂ = argmin ζ − c−1 d̄(p) (22)
p

where again, the nonlinear minimization problem in (22) can

be solved via grid search over the position coordinates space.

B. Approximate CToA Performance Analysis Fig. 8. Accuracy of CToA at “1st -Fix” Mode in a Typical Office Environment
In order to obtain theoretical performance bounds, the
skew-less measurement models derived in section III-A were
used for calculating the respective Cramér-Rao lower bounds
(CRLB). The derived bounds are affected mainly on the
geometrical properties of the network deployment, as well as
the additive noise levels. Yet, these bounds ignore propagation
models which may take into account the type of materials
through which the signals propagate.
To illustrate the expected positioning accuracy, the bounds
are depicted as “heat-maps”, where the colors are mapped to
the location according to the size of the minimal theoretical po-
sitioning error predicted by the CRLB. Under the assumption
that the additive measurement noise is Gaussian-distributed,
for each location on a given grid of locations one can calculate
the concentration ellipse that defines the smallest area at which
the CToA receiver is contained with a given probability. As
explained in Appendix C, the concentration ellipse is tightly
related to the position estimation error covariance matrix
predicted by the CRLB, which is derived in Appendix A. Fig. 9. Accuracy of CToA at Tracking Mode in a Typical Office Environment
Fig. 8- 9 describe the geometric-dependent accuracy of CToA
in a typical office network deployment, where the magenta As can be seen, once the CToA client has its EKF con-
rings denote the position of the bSTAs. The color mapping verged, the positioning error over the entire area drops signifi-
corresponds to the size of the major axis of the 95% percentile cantly compared to the “1st -Fix” scenario. It can be further be
 9

shown that positioning accuracy of the “1st -Fix” resembles to clock skew model, the instantaneous value of the clock-offset
the accuracy of a hyperbolic, time-difference of arrival (TDoA) of the nth bSTA is calculated using,
based, positioning system. When the client converges its time
tracking, the performance is similar to the performance that
can be achieved using ToA-based positioning system (namely,
an FTM-based system in which the client estimates ranges νn (ti ) = νn (ti−1 ) + ν̇n × ∆t (24)
(or round-trip time/RTT) to individual bSTAs). The following
proposition proves that the asymptotic accuracy of the client-
based CToA system is equivalent to “tracking” mode and
hence to the achievable RTT accuracy. where νn (ti−1 ) corresponds to the previous estimated value of
Proposition 1 (Asymptotic CToA Performance): The asymp- the clock offset, (whose νn (t0 ), is its initial value), ν̇n denotes
totic positioning accuracy for a CToA client in “1st -Fix” mode, the clock-skew (or the change rate of the clock offset), and
approaches the accuracy attained in “tracking” mode, given ∆t = ti − ti−1 .
N → ∞ replicas of the bSTA→bSTA timing measurements. In the following section we outline the algorithm, which
Proof: See Appendix B. enables the client station to estimate and track its location. The
Kalman Filter is the optimal estimate for linear system models
with additive independent white noise in both the transition
C. Coping with Clock Drifts Using Kalman Filtering and the measurement system models. Yet, in many systems,
including navigation systems, the measurement model is not
linearly dependent in the parameters of interest. In such cases,
FORFNRIIVHWVWLPHGLIIHUHQFH the extended Kalman filter (EKF), which is the nonlinear
FORFNVNHZIUHTXHQF\GLIIHUHQFH
PHDVXUHPHQWQRLVH
version of the Kalman filter, is widely used [22], [24]. In
the EKF, the state transition and observation models are not
required to be linear functions of the states, but instead, may
E67$
be only differentiable functions. In the client-mode CToA, the
EKF is executed by the client and is used by the client to
estimate and track its own position coordinates, as well as
E67$
the timing parameters of the stray bSTA units, from which
it receives the measurement broadcasts. In the network-mode
CToA, the EKF is executed at a centric positioning server,
connected to the CToA network, and is used for tracking the
F67$ position of multiple clients simultaneously, as well as tracking
the timing parameters of all the network bSTAs. Consequently,
the system states tracked by the EKF in each of the modes
are different; the network-mode EKF needs to track position
Fig. 10. Timing Measurements with Skewed Clocks
per client, as well as timing parameters of all network bSTA.
In the following section we focus on the client-mode CToA
The analysis outlined in the former section ignored any EKF.
clock skews, which result from the XO’s frequency deviations.
Such deviations may be caused due to multiple effects such as: 1) CToA EKF System Model: The EKF is described by
ambient temperature changes, phase noise, thermal noise, ag- two equations: a system model equation and an observation
ing, and so on. To illustrate the timing measurements collected (measurement) model with additive noise. The system model
under clock skews, consider the example depicted in Fig. 10. is defined by the following recursive equation,
In this example there is a single client station (cSTA) and two
bSTAs marked #A and #B. Each of the bSTAs has an initial
timing offset w.r.t. the cSTA clock, which is denoted by νA (t0 )
and νB (t0 ), respectively. In addition, the bSTAs XOs’ output xk = Fk xk−1 + wk , k≥0 (25)
frequencies are skewed, such that the clock offsets drift over
time at rates, which are denoted by ν̇A and ν̇B , respectively.
To facilitate the explanation, the assumption in this illustration
is that the skew rate is time-invariant. Yet, in practice it may where the index k denotes the discrete time-step. The vector
change over time (e.g., due to ambient effects). Each bSTA xk denotes an N × 1 states vector, which describes the
measures the broadcast events timing (ToD or ToA) relative parameters being estimated and tracked by the filter. The
to its native time base. Once the measurements are conveyed states vector for the client-mode CToA consists of the client’s
to the cSTA for enabling it to compute its location, then the position coordinates and per-bSTA clock parameters (clock
cSTA needs to resolve the instantaneous clock offset of the offset and clock offset change rate (or skew)). This vector is
bSTA, which is associated with that measurements and is a defined as follows. The size of the EKF state vector is thus:
function of the offset drift rate. Under the assumed 1st -order N = 3 + 2M , where M denotes the number of bSTAs being
 10

received by the cSTA (both directly and indirectly)5 . where,

T
 
pk , [xk , yk , zk ] σ̃x2 0 0
νk , [ν1,k , . . . , νM,k ]
T Qp,k , 0 σ̃y2 0  (33)
T 0 0 σ̃z2
ν̇ k , [ν̇1,k , . . . , ν̇M,k ]
[ ]T In general, the determination of the noise variance values of
xk , pTk , ν Tk , ν̇ Tk (26)
Qk , is challenging, and is often resorted to some heuristic
The state-vector xk is associated with a covariance matrix, methods. Commonly it is assumed that most of the clock
deviation is dictated by the clock skew, rather than clock mea-
Pk = E{(xk − x̄k )(xk − x̄k )T } (27) surement noise. The values of {σ̃x2 , σ̃y2 , σ̃z2 } are determined
where x̄k , E{xk }. When the filter is initialized the state- according to the motion assumptions of the cSTA device (e.g.,
covariance matrix is assumed to be, pedestrian, vehicle/drone etc.).
  2) CToA EKF Measurement Model: The measurement
P̃p,0 0 0 model is defined as,
P0 =  0 2
σν,0 IM 0  (28)
0 0 2
σν̇,0 IM zk = h(xk ) + vk (34)
2
where σν,0 , σν̇,0 denote the initial values for the standard where zk is a J × 1 vector of measurements, in which each
deviations of the clock offsets and drifts, and P̃p,0 denotes entry corresponds to a ToF measurement that follows the def-
the initial value of states covariance matrix, given by inition in (2). The vector h(x) , [h1 (x), h2 (x), . . . , hJ (x)]T ,
 2  denoting the nonlinear measurement model vector transfer
σx,0 0 0 function, and vk denotes the additive measurement noise that
P̃p,0 ,  0 2
σy,0 0  (29) has the following statistical properties:
2
0 0 σz,0
E{vk } = 0
The initial values of the standard deviations constructing
the initial states covariance matrix are commonly determined E{vk vkT } 2
= Rk = σ m I
empirically. E{vk vjT } = 0, ∀k ̸= j (35)
The dynamic system-model linear transfer function is denoted E{vk xTk } = 0, ∀k
by Fk , an N × N block-diagonal matrix defined as follows,
  E{wk vjT } = 0, ∀k, j
I3 0 0
Fk ,  0 IM ∆tIM  (30) As discussed in sec. III-A, there are two types of transfer
0 0 IM functions, which depend on the type of the measurement
(bSTAi → cSTA or bSTAi → bSTAj ). From (10)-(11) it
where ∆t corresponds to the elapsed time between two con- is easy to see that the corresponding measurement transfer
secutive discrete time steps. functions are given by,
The vector wk denotes a random N × 1 model noise vector,
which described the uncertainties in the system model and has 1
hℓ (xk ) = ∥pk − qi ∥ − eTi ν k (36)
the following statistical properties: c
1
hl (xk ) = ∥qj − qi ∥ + (ej − ei )T ν k (37)
E{wk } = 0 c
E{wk wkT } = Qk Since the measurement transfer function, h(·) is nonlinear, it
E{wk wjT } = 0, ∀k ̸= j (31) cannot be applied to estimate the measurements covariance
E{wk xTk } = 0, ∀k matrix directly. Instead we linearize h(·) by replacing it
with its first order Taylor series expansion, calculated around
In the CToA EKF system model, the process noise, wk is x̂k|k−1 :
assumed to be distributed as, wk ∼ N (0, Qk ), where system
model noise covariance matrix, Qk , is a block-diagonal matrix h(xk ) ∼
= h(x̂k|k−1 ) + Hk · (xk − x̂k|k−1 ) (38)
given by, where the notation, x̂n|m represents the estimate of x at time
 
Qp,k 0 0 n given observations up to and including time m ≤ n. The
Qk = ∆t ·  0 σ̃ν2 IM 0  (32) matrix Hk denotes the Jacobian of the measurement model
0 0 2
σ̃ν̇ IM function vector h(·), which is a J × N matrix defined as,
 ∂h1 ∂h1 ∂h1 
5 In a GPS system, since the entire GPS network is synchronized and uses ∂x1 ∂x2 · · · ∂x N
 .. 
highly-stable atomic clocks, the receiver has to track only its clock offset Hk ,  ... ..
. . 
and drift w.r.t. the GPS system clock. Furthermore, the SVs orbital motion
generates a substantial Doppler offset on the GPS carrier frequency, which ∂hJ ∂hJ
· · · ∂hJ
enables to estimate the GPS receiver 3-dimensional speed. Thus, the EKF 
∂x1 ∂x2 ∂xN
state vector in the GPS receiver typically includes a total of 8 states: 3 states 
∂hi 
for the receiver position, 3 states for the 3-dimensional receiver speed and 2 [Hk ]ij ≡  (39)
states for the clock model. ∂xj 
x=x̂k|k−1
 11

The Jacobian is obtained by calculating the partial derivatives trajectory was generated using a light detection and ranging
of (36)-(37). Equations (40)-(41) define the corresponding (LIDAR)-based ground truth tool [13]. The LIDAR system
lines of the matrix Hk . used integrates a 270o laser scanner, which uses a dedicated
[ ] map and laser measurements to estimate its position. The
(pk − qn )T
[Hk ]ℓ = , −eTi , 0TM (40) LIDAR output is a series of position reports generated at a
c∥pk − qn ∥
[ ] rate of 20Hz with an accuracy of 10-30 cm. The map is
[Hk ]l = 0T3 , (ej − ei )T , 0TM (41) obtained in advance by performing a survey of the venue using
The EKF is a recursive estimator, in which only the es- the LIDAR, during which a structure map is created using
timated state from the previous time step and the current a simultaneous localization and mapping (SLAM) algorithm.
measurement are required for the computation of the estimate This is a one-time procedure, and the generated map is then
for the current state. The state of the filter is represented by used in subsequent sessions for localization of the device.
two variables: the vector x̂k|k , which denotes the a posteriori Fig. 11 depicts the reference client trajectory (denoted by p
state estimate at time k given observations up to and including and marked by red dots), and the position estimates (denoted
at time k, and Pk|k , the a posteriori error covariance matrix. by p̂ and marked by line-connected blue dots). The location
The CToA EKF algorithm is summarized in Algorithm 1. of the bSTAs is marked by the red rings. The simulation
generated independent clock sources for each of the bSTAs.
Initialize These clocks were generated with an accuracy of ±10ppm
1. Use (16) and (18) for obtaining and with Gaussian-distributed clock noise with zero mean
an estimate for p̂0 and ν̂ 0 . and standard deviation of ±10−9 (1 ppb). The bSTAs were
2. Set x̂0 = [p̂T0 , ν̂ 0T , 0TM ]T set to broadcast CToA beacons at 2Hz. Each bSTA exchange
Predict measurements its 4 nearest neighbor bSTAs. On every point
EKF time as well as EKF states, are predicted according along the trajectory, the CToA client used measurement from
to the ToA of the received NDP packet. the 4 nearest bSTAs. The cumulative distribution function
Predicted state estimate: (CDF) of the positioning errors is depicted in Fig. 12. As can
be seen, the estimation accuracy is equal or better to 1.5m for
x̂k|k−1 = Fk x̂k−1|k−1 (42) 95% of the position estimates along the trajectory.
Predicted covariance estimate:
Pk|k−1 = Fk−1 Pk−1|k−1 FTk−1 + Qk−1 (43)
Update
The measurements included in the LMR conveyed by the
packet are updated according to the
new EKF predicted time.
Innovation (measurement residual):
ỹk = zk − h(x̂k|k−1 ) (44)
Innovation covariance:
Sk = Hk Pk|k−1 HTk + Rk (45)
Near-Optimal Kalman gain:
Kk = Pk|k−1 HTk S−1
k (46)
Updated state estimate:
Fig. 11. CToA bSTA Deployment a Typical Office Environment
x̂k|k = x̂k|k−1 + Kk ỹk (47)
Updated estimate covariance:
Pk|k = (I − Kk Hk ) Pk|k−1 (48) V. C ONCLUSIONS
Collaborative time of arrival (CToA), which is the next
Algorithm 1: CToA Client-Mode EKF Algorithm
generation, indoor geolocation protocol was presented. The
protocol is designed for enabling scalability of the existing
IEEE802.11/Wi-Fi-based, geolocation systems. The protocol
IV. CT OA P ERFORMANCE IN AN I NDOOR N ETWORK uses beacon broadcast-based fine-time delay measurements,
In the following section we outline some simulation perfor- collected by both clients and un-managed bSTA units, in
mance of an indoor CToA network. To analyze the accuracy of order to concurrently enable an unlimited number of clients
the client position estimates, the EKF position estimates were to privately navigate indoors, and enable a positioning server
compared against a “ground truth” trajectory. The reference to track a plethora of clients of a different type (e-Tags).
 12

100
Client Positioning Error CDF where d̄˙ x , d̄˙ y denote the vectors containing the partial deriva-
tives w.r.t. the client’s position coordinates, which are given
90
by,
80
∂ d¯i (xi − x)
70
= −√ (52)
60
∂x (xi − x)2 + (yi − y)2
¯
CDF [%]

∂ di (yi − y)
50 = −√ (53)
∂y (xi − x)2 + (yi − y)2
40

30 Using (51), the FIM elements can be found as,

20
2 ˙T ˙
Jxx = 2c−2 ḋTx W−1 ḋx = d̄ d̄
10 c2 σ̄ 2 x x
2
0
0 1 2 3 4 5 6 7 8 Jxy = Jyx = 2c−2 ḋTx W−1 ḋy = 2 2 d̄˙ Tx d̄˙ y
c σ̄
Positioning Error [m]
2
Jyy = 2c−2 ḋTy W−1 ḋy = 2 2 d̄˙ Ty d̄˙ y (54)
Fig. 12. CToA Client Positioning Error CDF in a Typical Office Environment c σ̄
Jxνi = 2c−1 ḋTx W−1 Eei
Jyνi = 2c−1 ḋTy W−1 Eei
Due to the infrequent nature of the beacon broadcasts,
and the skewness of the clocks used by the network units, Jνi νj = 2eTi ET W−1 Eej
the estimation of the position, as well as the clock-related
parameters is implemented using a Kalman filter, which is Define,
executed by every CToA client independently (or by the [ ]
Jxx Jxy
positioning server, for tracking the position of the CToA e- Jpp , (55)
Jyx Jyy
Tags). System simulations indicate that the network is capable [ ]
Jxν
of reaching a positioning accuracy of about 1.5m in a typical Jpν , (56)
office environment in 95% of the cases. These results also Jyν
match the theoretical performance predicated by the CRLB. The FIM is given by,
[ ]
A PPENDIX Jpp Jpν
J= (57)
JTpν Jνν
A. CToA Cramér Rao Lower Bound
We shall now derive the Cramér-Rao lower bound (CRLB) The CRLB is obtained by inverting the complete FIM.
for the CToA method in the absence of clock drifts. The [ ( ) ]
T −1
CRLB provides a lower bound on the covariance matrix of −1 Jpp − Jpν J−1 J
νν pν ()
J = (58)
any unbiased estimator. () ()
1) CRLB CToA Client in “1st Fix” Mode: Since the ob-
servations vector, m, is distributed m ∼ N (µ, W), the ijth The bound on the position coordinates is given by the top-left
entry of the Fisher information matrix (FIM) may be obtained block of J−1 ,
as [29], [ ]
T −1
{ } { T } C1ppFix = Jpp − Jpν J−1
st
νν Jpν (59)
∂W −1 ∂W ∂µ ∂µ
Jij = tr W−1 W +2 W−1 (49)
∂ψi ∂ψj ∂ψi ∂ψj 2) Approximate CRLB for CToA client in “Tracking” Mode:
where ψi is the ith element of the unknown parameters vector, When the EKF is converged and the bSTA clock offsets are
ψ , [p , ν ] . Noticing that W is free of any unknown
T T T known (up to some residual error), and are being continuously
parameters, then tracked, then µ ≃ c−1 d, and ψ , p.
{ T } Consequently,
∂µ −1 ∂µ
Jij = 2 W (50) { T }
∂ψi ∂ψj 2 ∂µ ∂µ
Jij ≃ 2 2 (60)
For, µ , c−1 d + Eν, the partial derivatives w.r.t. the client’s c (σ̄ + σr2 ) ∂ψi ∂ψj
position coordinates are given by, The partial derivatives are obtained using (51), and the CRLB
∂µ on the position coordinates estimation error is obtained by,
= c−1 ḋx ≡ c−1 [d̄˙ Tx 0TL ]T
∂x [ ]
1 Jyy −Jxy
∂µ
= c−1 ḋy ≡ c−1 [d̄˙ Ty 0TL ]T
Tracking
Cpp ≃
(51) J J − Jxy 2 −Jxy Jxx
∂y [ xx 2 yy ]
∂µ σxx σxy
= Eei = 2 (61)
∂νi σxy σyy
 13

B. Proof of Proposition 1 where,

Assume that every CToA broadcast includes N replicas of { }
R = p̂ : (p̂ − p)T Σ−1 (p̂ − p) ≤ κ (70)
the timing measurements collected by the bSTA. If the clock
offsets were time-invariant then one could define, As shown in [27, eq. (39)], for the 2-D case,
[ ] [ 2 ]
Ē σ̄ IL 0 Pin = 1 − exp(−κ/2)
Ě , , W̌ , (62) (71)
1N ⊗ Ẽ 0 σ̃ 2 IN ×L
Further we have,
Next, from (54) we have, [ ]
{ } 2
σxx σxy
Jνi νj = 2eTi ET W−1 Eej (63) Σ = E (p̂ − p)(p̂ − p)T =
σxy 2
σyy
(72)
Then, under (62), Jνi νj becomes,
The eigenvalues of Σ can be found by solving: det{Σ−λI} =
Jˇνi νj = 2eTi ĚT W̌−1 Ěej 0.
[ ] [ ]
[ −2 T ] Ē 1 2 √
−2
T
= 2ei σ̄ Ē σ̃ 1N ⊗ Ẽ
T T Eej λ1 = σ + σyy + (σxx − σyy ) + 4σxy
2 2 2 2 2 (73)
1N ⊗ Ẽ 2 xx
( ) [ ]
= 2eTi σ̄ −2 ĒT Ē + N · σ̃ −2 ẼT Ẽ ej 1 2 √
λ2 = σ − σyy + (σxx
2 2 − σ 2 )2 + 4σ 2 (74)
2N σ̃ −2 eTi ẼẼT ej . 2 xx yy xy
≈ (64)
N →∞

Recall that from (59) we have,

[
C1ppFix = Jpp − Jpν J−1
st ]
T −1 [ J
νν Jpν (65)
[
Hence, under N → ∞, J̌pν J̌−1 νν J̌pν → 0.
T

-
Thus, given enough bSTA→bSTA measurements (equivalent
1st Fix
to an EKF in “tracking” mode), Cpp ≈ CTracking
pp (up to
additive noise level scaling). J
This concludes the proof.

C. On the Relation between the Concentration Ellipse and the  NO

CRLB
 NO
In [27], Torrieri derived the theory for bounding Gaussian-
distributed estimation errors. In the geolocation problem con-
sidered, the parameters of interests include the 2-D coordinates Fig. 13. Concentration ellipse and coordinate axes
vectors, which are defined as: p ∈ R2 , p = [x, y]T . The
following appendix is based on the derivation in [27], and As depicted in Fig. 13, assuming that the principle axes of
outlines the relation between the concentration ellipse, which the concentration ellipse lie on the axes ξ1 , ξ2 , which are
provides a measure of accuracy for an unbiased estimator of counterclockwise rotated w.r.t. to axis system γ1 , γ2 by an
a 2-D Gaussian vector, and the CRLB. angle ϑ, then:
Let p̂ denote an unbiased estimate of p, then given that the [ ] [ ]
additive noise can be modeled as Gaussian, the probability ξ1 γ1
= AT (75)
density function (pdf) of the estimation error is given by, ξ2 γ2
[ ]
exp − 12 (p̂ − p)T Σ−1 (p̂ − p) where,
f (p̂|p) = n√ (66) [ ]
(2π) 2 det(Σ) cos ϑ − sin ϑ
A= (76)
where, sin ϑ cos ϑ
{ } is an orthogonal matrix (with eigenvectors as columns). Since
Σ , E (p̂ − p)(p̂ − p)T (67)
Σ is symmetric positive-definite matrix, and thus, so is Σ−1 ,
The loci of constant density is defined as, the matrix AT Σ−1 A is diagonal, provided that |ϑ| ≤ π4 .
(p̂ − p)T Σ−1 (p̂ − p) , κ (68) It can be shown that for a consecration ellipse defined by
γ T Σ−1 γ, √the principle√axes of the concentration ellipse are
The scalar κ is a constant that determines the size of the n- given by 2 κλ1 and 2 κλ2 . In order to obtain a bound on
dimensional region enclosed by the surface, which in the 2-D the maximal positioning error with a given probability, Pin ,
case is an ellipse. The probability that p̂ is contained inside at a specific position, one needs to evaluate √the size of half
that region is given by, of the ellipse’s major axis that is given by κλmax , where
∫ ∫ ∫ λmax is obtained by (73) and,
Pin = · · · f (p̂|p)dp̂ (69)
κ = −2 ln(1 − Pin ) (77)
R
 14

To summarize, the concentration ellipse is defined by 3 pa- [22] P. D. Groves, Principles of GNSS, Inertial, and Multisensor Integrated
rameters: its major axis, its minor axis and its rotation angle Navigation Systems, Artech House Boston-London, 2008.
[23] IEEE Std 802.11TM -2016 (Revision of IEEE Std 802.11-2012) - Part
- all of which can be obtained from the theoretical covariance 11: Wireless LAN Medium Access Control (MAC) and Physical Layer
matrix defined by the CRLB. This concludes the appendix. (PHY) Specifications, IEEE 802.11 Working Group, December 7th,
2016.
[24] S. J. Julier and J. K. Uhlmann, “Unscented filtering and nonlinear
R EFERENCES estimation,” Proc. of the IEEE, vol. 92, no. 3, pp. 401-422, Mar. 2004.
[25] G. A. Terejanu, “Extended Kalman Filter Tutorial,” available online:
[1] “Wireless E911 Location Accuracy Requirements,” Federal Communi- https://www.cse.sc.edu/ terejanu/files/tutorialEKF.pdf
cations Commission, PS Docket No. 07-114, Feb. 3, 2015. [26] D. Taniuchi, X. Liu, D. Nakai and T. Maekawa, “Spring Model
[2] “Mobile Device Feature Attach Rate and Penetration,” ABI Research, Based Collaborative Indoor Position Estimation With Neighbor Mobile
August 14, 2014 Devices,” in IEEE Jour. of Selected Topics in Signal Processing, vol.
[3] J. Zheng, C. Wu, H. Chu and P. Ji, “Localization Algorithm Based on 9, no. 2, pp. 268–277, Mar. 2015.
RSSI and Distance Geometry Constrain for Wireless Sensor Network,” [27] D. J. Torrieri, “Statistical Theory of Passive Location Systems,” IEEE
2010 International Conference on Electrical and Control Engineering, Trans. on Aerospace Systems, vol. AES-20, no. 2, pp. 183-198, Mar.
Wuhan, 2010, pp. 2836-2839. 1984.
[4] J. Torres-Sospedra et al., “Comprehensive analysis of distance and sim- [28] Q. M. Chaudhari, E. Serpedin and K. Qaraqe, “On Maximum Likeli-
ilarity measures for Wi-Fi fingerprinting indoor positioning systems,” hood Estimation of Clock Offset and Skew in Networks With Expo-
Expert Systems with Applications, vol. 42, no. 23, pp. 9263-9278, Dec. nential Delays,” IEEE Trans. on Signal Processing, vol. 56, no. 4, pp.
2015. 1685-1697, April 2008.
[5] A. T. Parameswaran, M. I. Husain, S. Upadhyaya, “Is RSSI a Reliable [29] H. L. Van Trees, Optimum Array Processing. Part IV of Detection,
Parameter in Sensor Localization Algorithms: An Experimental Study,” Estimation, and Modulation Theory, John Wiley & Sons, New York,
2009 2002.
[6] E. Elnahrawy, X. Li and R. P. Martin, “The limits of localization using
signal strength: a comparative study,” 2004 First Annual IEEE Commu-
nications Society Conference on Sensor and Ad Hoc Communications
and Networks, IEEE SECON 2004., 2004, pp. 406-414.
[7] K. Wu, J. Xiao, Y. Yi, D. Chen, X. Luo and L. M. Ni, “CSI-Based
Indoor Localization,” IEEE Trans. on Parallel and Distributed Systems,
vol. 24, no. 7, pp. 1300-1309, July 2013.
[8] Y. Shu et al., “Gradient-Based Fingerprinting for Indoor Localization
and Tracking,” IEEE Trans. on Industrial Electronics, vol. 63, no. 4,
pp. 2424-2433, April 2016.
[9] X. Li and K. Pahlavan, “Super-resolution TOA estimation with diversity
for indoor geolocation,” IEEE Trans. on Wireless Communications, vol.
3, no. 1, pp. 224-234, Jan. 2004.
[10] F. X. Ge, D. Shen, Y. Peng and V. O. K. Li, “Super-Resolution Time
Delay Estimation in Multipath Environments,” IEEE Trans. on Circuits
and Systems I: Regular Papers, vol. 54, no. 9, pp. 1977-1986, Sept.
2007.
[11] L. Banin, U. Schatzberg, Y. Amizur, “Next Generation Indoor Position-
ing System Based on WiFi Time of Flight,” 26th International Tech-
nical Meeting of the Satellite Division of The Institute of Navigation,
Nashville TN, Sept. 16-20, 2013.
[12] U. Schatzberg, L. Banin and Y. Amizur, “Enhanced WiFi ToF indoor
positioning system with MEMS-based INS and pedometric informa-
tion,” 2014 IEEE/ION Position, Location and Navigation Symposium -
PLANS 2014, Monterey, CA, 2014, pp. 185–192.
[13] L. Banin, U. Schatzberg and Y. Amizur, “WiFi FTM and Map Infor-
mation Fusion for Accurate Positioning,” 2016 International Confer-
ence on Indoor Positioning and Indoor Navigation (IPIN), Alcalá de
Henares, Spain, 4-7 October 2016.
[14] H. Kim, X. Ma and B. R. Hamilton, “Tracking Low-Precision Clocks
With Time-Varying Drifts Using Kalman Filtering,” IEEE/ACM Trans.
on Networking, vol. 20, no. 1, pp. 257–270, Feb. 2012.
[15] X. Cao et al., “Joint Estimation of Clock Skew and Offset in Pairwise
Broadcast Synchronization Mechanism,” IEEE Trans. on Comm., vol.
61, no. 6, pp. 2508–2521, June 2013.
[16] T. S. Rappaport, Wireless Communications: Principles and Practice,
Prentice Hall PTR, 1996.
[17] E. Perahia and R. Stacey, Next Generation Wireless LANs: 802.11n
and 802.11ac, and Wi-Fi Direct, Cambridge University Press, 2nd Ed.,
2013.
[18] A. J. Weiss, “On the Accuracy of a Cellular Location System Based
on RSS Measurements,” IEEE Trans. on Vehicular Technology, vol. 52,
no. 6, pp. 1508–1518, Nov. 2003.
[19] S. Mazuelas et al., “Robust Indoor Positioning Provided by Real-Time
RSSI Values in Unmodified WLAN Networks,” IEEE Jour. of Selected
Topics in Signal Processing, vol. 3, no. 5, pp. 821–831, Oct. 2009.
[20] J. Elson, L. Girod and D. Estrin, “Fine-Grained Network Time Syn-
chronization using Reference Broadcasts,” ACM SIGOPS Operating
Systems Review - OSDI ’02, vol. 36, no. SI, pp. 147-163, Winter
2002.
[21] M. A. Lombardi, T. P. Heavner and S. R. Jefferts, “NIST Primary
Frequency Standards and the Realization of the SI Second,” The
Journal of Measurement Science, vol. 2, no. 4, pp. 74-89, Dec. 2007.