Method and system for estimating position

Abstract

A method of estimating a position of an object is disclosed. The method comprises receiving a plurality of signals transmitted between each of a respective plurality of reference stations of a known position and the object. The method further comprises, for each received signal: processing the signal using a matrix featuring a characteristic power delay profile and a characteristic signal waveform to provide a processed signal, and calculating a direct estimate of the position of the object based, at least in part, on the processed signals.

Description

RELATED APPLICATIONS

This application is a National Phase of PCT Patent Application No. PCT/IL2013/050811 having International filing date of Oct. 3, 2013, which claims the benefit of priority under 35 USC §119(e) of U.S. Provisional Patent Application No. 61/709,361 filed on Oct. 4, 2012. The contents of the above applications are all incorporated by reference as if fully set forth herein in their entirety.

FIELD AND BACKGROUND OF THE INVENTION

The present invention, in some embodiments thereof, relates to a positioning technique and, more particularly, but not exclusively, to a method and system for estimating a position in a multipath environment.

Locating a target in a wireless system involves the collection of location information from radio signals propagating between the target and a number of base stations at known position.

The ability to determine the location or position of a wireless device has become increasingly desirable as an ever increasing number of people carry wireless devices, such as mobile phones, pagers, wireless email/Internet devices, and communications radios on a daily basis. In most cases, because individuals carry the wireless devices on or about their person, it is a reasonable assumption that the position of a particular wireless device is also that of its owner. As a result, locating an individual may be accomplished by locating the wireless device they carry.

With the proliferation of private and public Wi-Fi networks in recent years, several positioning systems based on Wi-Fi networks have been introduced. In a wide-area Wi-Fi positioning system, the location and characteristics of Wi-Fi access points are used to locate Wi-Fi enabled mobile devices.

Time of arrival (TOA) based positioning techniques rely on measuring the propagation times between the target and the base stations. These measurements are based on estimation of the first arrival time of the transmitted signal. In case of multipath channel, where the transmitted signal is reflected by objects, walls, cars, buildings, people etc., the received signal contains several overlapping replicas of the transmitted signal and with the addition of noise accurate estimation of the first arrival time is a considerable challenge.

TOA estimation methods based on passing the received signal through a matched filter (MF) whose output peak epoch is the TOA estimate are disclosed in W. Chung and D. Ha, “An accurate ultra wideband (UWB) ranging for precision asset location,” in Proc. IEEE Conf. Ultrawideband Syst. Technol. (UWBST), Reston, Va., November 2003, pp. 389-393; and K. Yu and I. Oppermann, “Performance of UWB position estimation based on time-of-arrival measurements,” in Proc. IEEE Conf. Ultrawideband Syst. Technol. (UWBST), Kyoto, Japan, May 2004, pp. 400-404.

An alternative solution where a search for a finite number of MF output peak epochs is performed and then the TOA is estimated as the earliest peak epoch is disclosed in V. Lottici, A. D'Andrea, and U. Mengali, “Channel estimation for ultrawideband communications,” IEEE J. Select. Areas Commun., vol. 20, no. 9, pp. 1638-1645, December 2002. Other solutions are based on comparing the received signal energy to a threshold and then performing a forward or backwards search with a heuristic TOA estimation criterion [L. Stoica, A. Rabbachin, and I. Oppermann, “A low-complexity noncoherent IR-UWB transceiver architecture with TOA estimation,” IEEE Trans. Microwave Theory Tech., vol. 54, pp. 1637-1646, June 2006; I. Guvenc, Z. Sahinoglu, and P. V. Orlik, “TOA Estimation for IR-UWB Systems With Different Transceiver Types,” IEEE Trans. Microwave Theory and Tech., vol. 54, no. 4, April 2006; D. Dardari, C.-C. Chong, and M. Z. Win, “Analysis of threshold-based TOA estimator in UWB channels,” in Proc. Eur. Signal Process. Conf. (EUSIPCO), Florence, Italy, September 2006; A. Chehri, P. Fortier and P. M. Tardif, “Time-of-Arrival Estimation For IR-UWB Systems Based on Two Step Energy Detection” Biennial Symposium on Commun., pp. 369-373, June 2008.

Other approaches include the “frequency-domain super-resolution TOA estimation” which includes Multiple Signal Classification (MUSIC) [X. Li, “Super-Resolution TOA Estimation With Diversity for Indoor Geolocation,” IEEE Trans. Wireless Commun., vol. 3, no. 1, January 2004], Estimation of Signal Parameters via Rotational Invariance Techniques (ESPRIT) [R. Roy, A. Paulraj and T. Kailath, “Estimation of Signal Parameters via Rotational Invariance Techniques-ESPRIT”, IEEE Trans. Acoust., Speech, Signal Process., vol. 37, no. 7, pp. 984-995, July 1989.] and other algorithms [R. J. J. Vidal, M. Najar, “High resolution time-of-arrival detection for wireless positioning systems”, In IEEE Vehicular Technology Conference, Vancouver, Canada, pp. 2283-2287, September 2002; M. Navarro and M. Najar, “Joint Estimation of TOA and DOA in IR-UWB”, IEEE Workshop on Signal Processing Advances in Wireless Commun., Helsinki, Finland, pp. 17-20, June 2007]. Win and Scholtz [M. Z. Win and R. A. Scholtz, “Characterization of ultra-wide bandwidth wireless indoor communications channels: A communication theoretic view,” IEEE J. Sel. Areas Commun., vol. 20, pp. 1613-1627, December 2002] introduced the generalized maximum likelihood (GML) solution. The method relies on the assumption that the number of multipath components is finite and known and all the multipath coefficients and their arrival times are jointly estimated.

An iterative approximation to the GML was introduced by [J. Y. Lee and R. A. Scholtz, “Ranging in a dense multipath environment using an UWB radio link,” IEEE J. Select. Areas Commun, vol. 20, no. 9, pp. 1677-1683, December 2002.]. The multipath coefficients are estimated one by one and removed from the received signal sequentially until the strongest remaining coefficient is below a preset threshold. The TOA is determined according to the earliest arrival path.

European Publication No. EP1596217 discloses a method of detecting the time of arrival (TOA) of the first-arrival-path signal in a transmission system where wide band signals are transmitted. A parameter related to the energy of a plurality of samples of a received signal resulting from a multipath propagation of a transmitted signal is evaluated. For each TOA candidate sample, the observation period is divided into a corresponding pair of complementary windows, where the first window in each pair comprises the samples preceding the candidate sample and the other one comprises the remaining samples in the period. The signal energy in both windows of each window pair is evaluated, and the sample giving rise to windows whose energies maximize a function related to the energy ratio of the two windows in each pair is chosen as the estimated TOA sample.

Additional background art includes I. Guvenc and C. C. Chong, “A Survey on TOA Based Wireless Location and NLOS Mitigation Techniques”, IEEE Comm. Surveys and Tutorials, vol. 11, no. 3, pp. 107-124, 2009; K. W. K. Lui, H. C. So and W. K. Ma, “Maximum A Posteriori Approach to Time-of-Arrival-Based Localization in Non-Line-of-Sight Environment”, IEEE Trans. on Vehicular Technology, vol. 53, no. 3, pp. 1517-1523, Match 2010; P. C. Chen, “A non-line-of-sight error mitigation algorithm in location estimation”, in Proc. IEEE Int. Conf. Wireless Commun. Netw. (WCNC), vol. 1, pp. 316-320, September 1999; A. Al-Jazzar and J. J. Caffery, “ML and Bayesian TOA location estimators for NLOS environments”, in Proc. IEEE Veh. Technol. Conf. (VTS), vol. 2, pp. 1178-1181, September 2002; Y. T. Chan, H. Y. Chin Hang and Pak-chung Ching, “Exact and Approximate Maximum Likelihood Localization Algorithms”, IEEE Trans. on Vehicular Technology, vol. 55, no. 1, pp. 10-16, January 2006; S. Gezici and H. V. Poor “Position Estimation via Ultra-Wide-Band Signals”, Proceeding of the IEEE, vol. 97, no. 2, pp. 386-403, February 2009; A. J. Weiss, “Direct Position Determination of Narrowband Radio Frequency Transmitters”, IEEE Signal Processing Letters, vol. 11, no. 5, May 2004; P. Closas, C. Fernandez-Prades and J. A. Fernandez-Rubio, “Maximum Likelihood Estimation of Position in GNSS”, IEEE Signal Processing Letters, vol. 14, no. 5, May 2007; P. Closas, C. Fernandez-Prades and J. A. Fernandez-Rubio, “Cramer-Rao Bound Analysis of Positioning Approaches in GNSS Receivers”, IEEE Transactions On Signal Processing, vol. 57, no. 10, October 2009; P. Closas, C. Fernandez-Prades, D. Bernal and J. A. Fernandez-Rubio, “Bayesian Direct Position Estimation”, IEEE Aerospace Conference, 2010; C. Yen and P. J. Voltz, “Indoor positioning based on statistical multipath channel modeling,” EURASIP Journal on Wireless Communications and Networking 2011; K. Papakonstantinou and D. Slock, “Direct Location Estimation for MIMO Systems in Multipath Environments”, IEEE Global Communications Conference (GLOBECOM). November 2008; and M. Eric and D. Vucic, “Direct position estimation of UWB transmitters in multipath conditions”, IEEE International Conference on Ultra-Wideband (ICUWB 2008), vol. 1. 2008.

SUMMARY OF THE INVENTION

According to an aspect of some embodiments of the present invention there is provided a method of estimating a position of an object. The method comprises: receiving a plurality of signals transmitted between each of a respective plurality of reference stations of known position and the object. The method further comprises, for each received signal: processing the signal using a matrix featuring a characteristic power delay profile and a characteristic signal waveform to provide a processed signal. The method further comprises calculating a direct estimate of the position of the object based, at least in part, on the processed signals.

According to some embodiments of the invention the matrix is based on an autocorrelation matrix of the received signal and is a function of the position of the object.

According to some embodiments of the invention the method further comprises partitioning the signal into at least a first segment containing noise and a second segment containing at least a multipath signal corresponding to a plurality of paths between the object and a respective reference station, wherein the calculation of the direct estimate is based, at least in part, on the segments.

According to some embodiments of the invention the calculating the estimate comprises combining, for each signal, contributions from the first segment and the second segment in a statistically independent manner.

According to an aspect of some embodiments of the present invention there is provided a system for estimating a position of an object. The system comprises a plurality of reference stations of know positions; and a data processor, configured for executing the method as described above and further exemplified below.

According to an aspect of some embodiments of the present invention there is provided a computer software product. The product comprises a computer-readable medium in which program instructions are stored, which instructions, when read by a data processor, cause the data processor to receive a plurality of signals transmitted between each of a respective plurality of reference stations of know position and an object and execute the method as described above and further exemplified below.

According to an aspect of some embodiments of the present invention there is provided a method of estimating a time-of-arrival (TOA) of a signal transmitted between an object and a reference station. The method comprises: receiving a multipath signal having a plurality of components respectively corresponding to a plurality of paths between the object and the reference station; processing the signal using a matrix featuring a characteristic power delay profile and a characteristic signal waveform to provide a processed signal; and calculating an estimate of the TOA of the signal based, at least in part, on the processed signal.

According to some embodiments of the invention the matrix is based on a noise variance.

According to some embodiments of the invention the matrix is based on an autocorrelation matrix of the signal and is a function of the time of arrival of the received signal.

According to some embodiments of the invention the processing comprising modeling the signal using a function characterized by the matrix, the function being selected from the group consisting of a Gaussian function, a Lorentzian function, a hyperbolic secant function, a logistic distribution, a Fourier transform and a wavelet transform.

According to some embodiments of the invention at least one of the signals is a diversity signal having diversity measurements, and the method further comprises processing the diversity measurements, wherein the estimate is based on the processed diversity measurements.

According to some embodiments of the invention the diversity comprises at least one of: time diversity, frequency diversity, interference diversity, antenna polarity diversity and space diversity.

According to some embodiments of the invention the space diversity is formed by a system selected from the group consisting of Multiple-Input Multiple-Output (MIMO) system, a Single-Input Multiple-Output (SIMO) system and a Multiple-Input Single-Output (MISO) system.

According to some embodiments of the invention the diversity is formed by frequency hopping.

According to some embodiments of the invention the diversity measurements comprises data received from a plurality of statistically independent channels of a respective reference station.

According to some embodiments of the invention the diversity measurements comprises data received from a plurality of correlated channels of a respective reference station.

According to some embodiments of the invention the method further comprises partitioning the signal into at least a first segment containing noise and a second segment containing at least a multipath signal corresponding to a plurality of paths between the object and a respective reference station, wherein the calculation of the TOA estimate is based, at least in part, on the segments.

According to some embodiments of the invention the partitioning comprises employing an energy search procedure.

According to some embodiments of the invention the calculating the estimate comprises combining, for each signal, contributions from the first segment and the second segment in a statistically independent manner.

According to some embodiments of the invention a signal transmitted between the object and each of the at least two reference stations is a two-way signal, and wherein the matrix features a difference between a signal transmission time of the object and a signal transmission time of a respective reference station.

According to some embodiments of the invention at least two of the reference stations are synchronized, and the object is unsynchronized with any reference station.

According to some embodiments of the invention a signal transmitted between the object and each of the at least two reference stations is a one-way signal, and wherein the plurality of reference stations comprises at least three reference stations.

According to some embodiments of the invention the estimate is in three-dimensions, and wherein the plurality of reference stations comprises at least four reference stations.

According to some embodiments of the invention at least two of the reference stations are unsynchronized, wherein the object is unsynchronized with any reference station, and wherein a signal transmitted between the object and each of the at least two reference stations is a two-way signal.

According to some embodiments of the invention at least two of the reference stations are synchronized, wherein the object is unsynchronized with any reference station, and wherein a signal transmitted between the object and each of the at least two reference stations is a two-way signal.

According to some embodiments of the invention the reference stations are located at indoor positions, and the method is utilized for estimating a distance between a reference station and an object.

According to some embodiments of the invention the object is a moving object.

According to some embodiments of the invention the reference stations are located at outdoor positions, and the method is utilized for estimating an outdoor distance of the object.

According to some embodiments of the invention the object is an inventory item.

According to some embodiments of the invention the object is a human or an animal.

According to some embodiments of the invention the object is a vehicle.

According to some embodiments of the invention the object is a mobile phone.

According to some embodiments of the invention the reference stations are Access Points, transmitting and receiving over a wireless network according to the LAN protocol, and wherein the object is a client of the wireless network.

According to some embodiments of the invention the reference stations system are satellites are transmitting according to a Global Navigation Satellite System (GNSS) protocol, and the object comprises a receiver for GNSS signals, where the GNSS is selected from the group consisting of a Global Positioning System (GPS), a Global Orbiting Navigation Satellite System (GLONASS) and a Galileo navigation system.

According to some embodiments of the invention the signals are in the time domain.

According to some embodiments of the invention the signals are in the frequency domain.

According to some embodiments of the invention the calculation of the direct estimate and/or TOA estimate comprises using the matrix for constructing a polynomial objective function, and searching for frequency at which the polynomial objective function has an extremum.

According to some embodiments of the invention the matrix is the same for all received signals.

According to some embodiments of the invention the matrix is different for at least two of the received signals.

According to some embodiments of the invention the characteristic power delay profile is calculated based on an estimate of a signal-to-noise ratio.

According to some embodiments of the invention the method further comprises updating the characteristic power delay profile using data collected from signal transmissions between the reference stations and other objects.

According to some embodiments of the invention the method further comprises, obtaining at least one additional power delay profile, wherein the processing the signal is executed separately for each power delay profile.

According to some embodiments of the invention the method comprises scoring each power delay profile and selecting a power delay profile based on the score, wherein the calculation of the direct estimate and/or the TOA estimate is based on the selected power delay profile.

According to some embodiments of the invention the method further comprises jointly estimating a plurality of pairs of candidate power delay profiles and respective candidate direct estimations of the position of the object, scoring each pair of the plurality of pairs, and selecting the direct estimation of the position of the object based on the score.

According to some embodiments of the invention the method further comprises jointly estimating a plurality of pairs of candidate power delay profiles and respective candidate TOAs, scoring each pair of the plurality of pairs, and selecting the TOA estimate based on the score.

According to some embodiments of the invention the method further comprises, obtaining at least one additional signal waveform, wherein the processing the signal is executed separately for each signal waveform.

According to some embodiments of the invention the method comprises scoring each signal waveform and selecting a signal waveform based on the score, wherein the calculation of the direct estimate and/or the TOA estimate is based on the selected signal waveform.

According to some embodiments of the invention the method further comprises jointly estimating a plurality of pairs of candidate signal waveforms and respective candidate direct estimations of the position of the object, scoring each pair of the plurality of pairs, and selecting the direct estimation of the position of the object based on the score.

According to some embodiments of the invention the method further comprises jointly estimating a plurality of pairs of candidate signal waveforms and respective candidate TOAs, scoring each pair of the plurality of pairs, and selecting the TOA estimate based on the score.

According to some embodiments of the invention the method further comprises: jointly estimating a plurality of sets of: candidate power delay profiles, candidate signal waveforms, and respective candidate direct estimations of the position of the object, scoring each set of the plurality of sets, and selecting the direct estimation of the position of the object based on the score.

According to some embodiments of the invention the method further comprises: jointly estimating a plurality of sets of: candidate power delay profiles, candidate signal waveforms, and respective candidate TOAs, scoring each set of the plurality of sets, and selecting the TOA estimate based on the score.

According to some embodiments of the invention the method further comprises obtaining a model describing a relation between the power delay profile and a distance from the object to a respective reference station, wherein the matrix is expressed in terms of the model.

According to an aspect of some embodiments of the present invention there is provided a system for estimating a position of an object, comprising: a plurality of reference stations of know positions; and a data processor, configured for executing the method as described above and further exemplified below to provide a TOA of signals transmitted between the object and each of the reference stations and for calculating the position based on the TOA.

According to an aspect of some embodiments of the present invention there is provided a computer software product, comprising a computer-readable medium in which program instructions are stored, which instructions, when read by a data processor, cause the data processor to receive a plurality of signals transmitted between each of a respective plurality of reference stations of know position and an object and execute the method as described above and further exemplified below.

According to an aspect of some embodiments of the present invention there is provided a method of estimating a position of an object. The method comprises comprising: receiving a plurality of signals transmitted between each of a respective plurality of reference stations of known position and the object; for each signal, estimating a time-of-arrival (TOA) of the signal using the method as described above and further exemplified below, thereby providing a plurality of TOA estimates; and calculating an indirect estimate of the position of the object based, at least in part, on the plurality of TOA estimates.

According to some embodiments of the invention the calculation of the indirect estimate comprises utilizing at least one procedure selected from the group consisting of a weighted least-squares procedure, a non-weighted least-squares procedure and a Kalman filter procedure.

According to an aspect of some embodiments of the present invention there is provided a method of estimating a distance between an object and a reference station, comprising: receiving a two-way multipath signal having a plurality of components respectively corresponding to a plurality of paths between the object and the reference station; processing the signal using a matrix featuring a difference between a signal transmission time of the object and a signal transmission time of the reference station; and calculating an estimate of the propagation delay of the signal based, at least in part, on the processed signal.

Unless otherwise defined, all technical and/or scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which the invention pertains. Although methods and materials similar or equivalent to those described herein can be used in the practice or testing of embodiments of the invention, exemplary methods and/or materials are described below. In case of conflict, the patent specification, including definitions, will control. In addition, the materials, methods, and examples are illustrative only and are not intended to be necessarily limiting.

Implementation of the method and/or system of embodiments of the invention can involve performing or completing selected tasks manually, automatically, or a combination thereof. Moreover, according to actual instrumentation and equipment of embodiments of the method and/or system of the invention, several selected tasks could be implemented by hardware, by software or by firmware or by a combination thereof using an operating system.

For example, hardware for performing selected tasks according to embodiments of the invention could be implemented as a chip or a circuit. As software, selected tasks according to embodiments of the invention could be implemented as a plurality of software instructions being executed by a computer using any suitable operating system. In an exemplary embodiment of the invention, one or more tasks according to exemplary embodiments of method and/or system as described herein are performed by a data processor, such as a computing platform for executing a plurality of instructions. Optionally, the data processor includes a volatile memory for storing instructions and/or data and/or a non-volatile storage, for example, a magnetic hard-disk and/or removable media, for storing instructions and/or data. Optionally, a network connection is provided as well. A display and/or a user input device such as a keyboard or mouse are optionally provided as well.

BRIEF DESCRIPTION OF THE DRAWINGS

Some embodiments of the invention are herein described, by way of example only, with reference to the accompanying drawings. With specific reference now to the drawings in detail, it is stressed that the particulars shown are by way of example and for purposes of illustrative discussion of embodiments of the invention. In this regard, the description taken with the drawings makes apparent to those skilled in the art how embodiments of the invention may be practiced.

In the drawings:

FIG. 1 is a schematic illustration of a baseband system model, used according to some embodiments of the present invention;

FIG. 2 is a schematic illustration of a section of a signal, according to some embodiments of the present invention;

FIG. 3 is a is a schematic illustration of a partitioned observation vector, according to some embodiments of the present invention;

FIG. 4 shows a signal used in computer simulations executed according to some embodiments of the present invention;

FIG. 5 shows root mean square error (RMSE) as obtained by computer simulations executed according to some embodiments of the present invention, using the signal shown in FIG. 4; Ray signal to noise ratio (RSNR) of 30 dB;

FIG. 6 shows RMSE as obtained by computer simulations executed according to some embodiments of the present invention, using the signal shown in FIG. 4; RSNR of 45 dB;

FIG. 7 shows a signal used in additional computer simulations executed according to some embodiments of the present invention;

FIG. 8 shows performances of an estimator according to some embodiments of the present invention, evaluated using the signal shown in FIG. 7 and a channel modeled as white complex Gaussian process with exponential decaying variance;

FIG. 9 shows performances of an estimator according to some embodiments of the present invention, evaluated using the signal shown in FIG. 7 for IEEE 802.15.4a in residential indoor environment channel model;

FIG. 10 shows performances of an estimator according to some embodiments of the present invention, evaluated using the signal shown in FIG. 7 (single channel) for IEEE 802.15.4a in indoor office environment channel model;

FIG. 11 shows performances of an estimator according to some embodiments of the present invention, evaluated using the signal shown in FIG. 7 (five channels) for IEEE 802.15.4a in indoor office environment;

FIG. 12 shows a signal used in further computer simulations executed according to some embodiments of the present invention;

FIG. 13 shows the root mean square error (RMSE) as a function of the signal to noise ratio (SNR), as obtained by computer simulations performed according to some embodiments of the present invention for four reference stations and a target located at (25,25,25)m;

FIG. 14 shows the root mean square error (RMSE) as a function of the signal to noise ratio (SNR), as obtained by computer simulations performed according to some embodiments of the present invention for six reference stations and a target located at (25,25,25)m;

FIG. 15 shows the root mean square error (RMSE) as a function of the signal to noise ratio (SNR), as obtained by computer simulations performed according to some embodiments of the present invention for four reference stations and a target located at (25,25,47.5)m;

FIG. 16 shows the root mean square error (RMSE) as a function of the signal to noise ratio (SNR), as obtained by computer simulations performed according to some embodiments of the present invention for six reference stations and a target located at (25,25,47.5)m;

FIG. 17 is a flowchart diagram of a method of estimating a position of an object according to some embodiments of the present invention;

FIG. 18 schematic illustration of a system for estimating a position of an object, according to some exemplary embodiments of the present invention;

FIG. 19 shows functions that can be used for approximating a power delay profile, according to some embodiments of the present invention;

FIG. 20 shows additional functions that can be used for approximating a power delay profile, according to other embodiments of the present invention;

FIG. 21 shows additional functions that can be used for approximating a power delay profile, according to other embodiments of the present invention;

FIG. 22 shows a signal waveform composed by Sinc function multiplied by a Hamming window, which can be used according to some embodiments of the present invention; and

FIG. 23 shows a Dirichlet signal waveform, which can be used according to some embodiments of the present invention.

DESCRIPTION OF SPECIFIC EMBODIMENTS OF THE INVENTION

The present invention, in some embodiments thereof, relates to a positioning technique and, more particularly, but not exclusively, to a method and system for estimating a position in a multipath environment.

Before explaining at least one embodiment of the invention in detail, it is to be understood that the invention is not necessarily limited in its application to the details of construction and the arrangement of the components and/or methods set forth in the following description and/or illustrated in the drawings and/or the Examples. The invention is capable of other embodiments or of being practiced or carried out in various ways.

Position estimation methods can be partitioned into two classes: direct and indirect position estimation.

In the direct position estimation (DPE) the target position is estimated in one step, directly from the received signals (the raw measurements) from all base stations together. The received signals are transferred to the server for joint processing.

The indirect position estimation (IPE) is done in two stages. In a first stage, a few position related parameters are estimated at each base station independently. In a second stage, the position is estimated based on the parameters obtained in the first stage. Typically, the position estimation in the second stage is the position that minimizes the sum of squared errors between the measured distances (obtained in the first step) and the distances corresponding to the estimated position. If a figure of merit is provided with each distance measurement then the position can be estimated by weighted least squares [Gezici et al., supra].

It was found by the present inventors that indirect positioning requires passing less data to the location estimation processor.

As used herein, “multipath” refers to the reception of a line of sight signal as well as a non-line of sight signal transmitted from a signal source.

Multipath signals may result, for example, from the reflection of a signal from a nearby reflector, such as the ground, objects, buildings, face or the surface of a body of water. Multipath signals may also result when signals are significantly refracted. In general a non-reflected or refracted (straight line) signal is referred to as “line-of-sight signal,” whereas a reflected or refracted signal is referred to as a “diverted signal.” The term “multipath signal” refers to a signal which includes a line-of-sight signal and at least one diverted signal.

Since the diverted signal does not travel along the line-of-sight, it always arrives at the receiver later than the line-of-sight signal.

Referring now to the drawings, FIG. 17 is a flowchart diagram of a method of estimating a position of an object according to some embodiments of the present invention. It is to be understood that, unless otherwise defined, the operations described hereinbelow can be executed either contemporaneously or sequentially in many combinations or orders of execution. Specifically, the ordering of the flowchart diagrams is not to be considered as limiting. For example, two or more operations, appearing in the following description or in the flowchart diagrams in a particular order, can be executed in a different order (e.g., a reverse order) or substantially contemporaneously. Additionally, several operations described below are optional and may not be executed.

The object can be a cellular phone, an inventory item, a human being, an animal, a vehicle, a robot and the like. The method is suitable for both indoors and outdoors applications, and can be used for estimating the position of a stationary or moving objects.

The method begins at 10 and continues to 11 at which the method receives a plurality of signals transmitted between each of a respective plurality of reference stations of known position and the object. Optionally and preferably the method receives at least one signal from each reference station. Thus, for example, when the method employs K reference stations, the method receives at least K signals. The reference stations can be positioned indoors or outdoors.

Many types of reference stations are contemplated. For example, in some embodiments of the present invention one or more (e.g., all) of the reference stations is an Access Point (AP) 180B communicating with the object over a wireless link. Any AP known in the art can be employed. Representative examples include, without limitation, an IEEE 802.11 Wireless Local Access Network such as, but not limited to, IEEE 802.11a-g, IEEE 802.15.4, and derivatives thereof; cellular-based data systems including, without limitation, the “TIA/EIA-95-B Mobile Station-Base Station Compatibility Standard for Dual-Mode Wideband Spread Spectrum Cellular System” (the IS-95 standard), the standard offered by a consortium named “3rd Generation Partnership Project” (3GPP) and embodied in a set of documents including Document Nos. 3G TS 25.211, 3G TS 25.212, 3G TS 25.213, and 3G TS 25.214 (the W-CDMA standard), the standard offered by a consortium named “3rd Generation Partnership Project 2” (3GPP2) and embodied in “TR-45.5 Physical Layer Standard for cdma2000 Spread Spectrum Systems” (the IS-2000 standard), and the high data rate (HDR) system that conforms to the TIA/EIA/IS-856 standard (the IS-856 standard).

In some embodiments of the present invention one or more (e.g., all) of the reference stations is a satellite transmission system, which can transmit according to a Global Navigation Satellite System (GNSS) protocol, such as, but not limited to, a Global Positioning System (GPS) protocol, a Global Orbiting Navigation Satellite System (GLONASS) protocol and a Galileo navigation system protocol.

In any event, the object includes a transmitter and/or receiver configured for communicating with each of the reference stations.

Each received signal can be a multipath signal. The received signal can be in the time domain or the frequency domain, as desired.

When the signal is the time domain, the signal is optionally a complex sampled signal represented by a vector, wherein each component of the vector can be a signal intensity value and phase associated with a different time point within the signal. Optionally, but not necessarily, the components are equally spaced in terms of the corresponding time points.

When the signal is the frequency domain, the signal is optionally a complex sampled signal represented by a vector, wherein each component of the vector can be a signal intensity value and phase associated with a different frequency within the signal. Optionally, but not necessarily, the components are equally spaced in terms of the corresponding frequencies. Alternatively, the components can be spaced according to a non-linear scale, such as, but not limited to, a logarithmic scale. A representative example of a signal in the frequency domain is, without limitation, a Fast Fourier Transform (FFT) or Discrete Fourier Transform (DFT) vector.

While the embodiments below are described with a particular emphasis to signals in the time domain, it is to be understood that more detailed reference to signals in the time domain is not to be interpreted as limiting the scope of the invention in any way, and one of ordinary skills in the art, provided with the details described herein would know how to adjust the procedures of the present invention for the case in which the signal is represented in the frequency domain.

Throughout this specification, underlined or bold symbols represent vector or matrix quantities.

The method optionally and preferably continues to 12 at which each received signal is processed using a matrix R featuring a characteristic power delay profile σ_h²(t) and a characteristic signal waveform and optionally a noise variance.

The power delay profile is typically a channel average power as a function of the time delay and it can be expressed by σ_h²(t)=E{|h(t)|²} where E{x} is the expectation of the quantity x, and h(t) is the multipath channel signal. The power delay profile expresses the statistical characteristics of the multipath channel and is typically dependent on the specific multipath environment, and may also vary within a specific environment according to the object positions. For example, for a given reference station positions the power delay profile may differ for different target positions.

The power delay profile can be provided in more than one way. In some embodiments of the present invention it is provided from external source, for example, a user input, in some embodiments of the present invention it is a predetermined profile, and in some embodiments of the present invention it is estimated by the method.

The power delay profile can be initialized per a specific location environment and kept fixed. In some embodiments of the present invention the power delay profile can be adaptively varied to better match the multipath statistical characteristics of the specific environment and of the specific object position.

Optionally and preferably the power delay is estimated independently per each link between reference station and object. Estimating the power delay profile is preferably performed utilized a priori information such as SNR information or information pertaining to the expected range between the object and the reference station. Alternatively, the power delay profile can be estimated jointly and incorporeally with the position estimation.

In some embodiments of the present invention two or more power delay profile candidates are employed, and a separate matrix is calculated for each delay profile candidate. The signal processing as further detailed hereinbelow can be executed for each matrix. In some embodiments, the power delay profile candidates are scored during the processing. In these embodiments, the power delay profile score is used for selecting one power delay profile from the candidates (e.g., the power delay profile with the highest score) and the position is estimated using the selected power delay profile.

Also contemplated are embodiments in which a model describing the relation between the power delay profile and the distance from the object to a respective reference station is obtained. In these embodiments, the matrix R is expressed in terms of this model. A representative example of this embodiment is provided in Example 8 of the examples section that follows.

Two or more power delay profile candidates can also be employed according to some embodiments of the present invention in a joint estimation of candidate positions and candidate power delay profiles. In these embodiments, a score is optionally and preferably given to each pair of candidate power delay profile and candidate position. The estimated position of the object can then be selected based on the score. Typically, but not necessarily, the candidate position that is associated with the higher score is selected as the estimated position of the object.

A more detailed description of an embodiment which employs a plurality of power delay profile candidates is described hereinunder.

The present inventors have found several functions that can be used to approximate the power delay profile. The advantageous of these functions is that they depend on one or two unknown parameters which are relatively easy to estimate.

Representative examples of functions which can be used as approximations to the power delay profile will now be provided, with reference to FIGS. 19-21. The power delay profile is initiated at the TOA, τ.

In some embodiments of the present invention the power delay profile is approximated by
σ_h²(t)≈c(t)p(t),
where c(t)=0 for t<τ, c(t)=1 at t=τ, c(t)=λ for t>τ and p(t)=0 for t<τ, p(t)=e^{−(t−τ)/Tp} for t≧τ.

Thus, for a given TOA the approximated channel power delay profile is dependent on the two parameters λ, T_p.

The λ parameter is the multipath density (for example, the probability of a multipath arrival within an infinitesimal interval that follows the TOA and does not include the first arrival). In some embodiments of the present invention λ is approximated as constant for t>τ. The value of λ increases as the multipath environment has more dense reflectors and possibly reduces as the SNR increases.

The T_P parameter relates to the SNR. As the SNR increases T_P reduces. Hence the power delay profile of the present embodiments is determined by the SNR which is specific per object position and by a measure of the environment reflectors density which is general parameter that is typically common to most of the target positions. The SNR and measure of reflectors density can be estimated as known in the art. FIG. 19 shows an example for c(t) 100, p(t) 101, σ_h²(t) 102 and the TOA 103 for a typical indoor multipath environment.

In some embodiments of the present invention the power delay profile is approximated as
σ_h²(t)≈c(t)j(t)
where c(t) is as defined above and j(t)=1 for τ≦t<τ+T_M and j(t)=0 otherwise, where T_M is the average channel delay spread that includes the energetic part of the multipath. For example T_M can be the multipath duration from the TOA until the multipath power decays by 8 dB. T_M can be estimated per specific received signal and/or according to the SNR. FIG. 20 shows an example of c(t) 110, j(t) 111, σ_h²(t) 112 and the TOA 113.

In some embodiments of the present invention the power delay profile is approximated as
σ_h²(t)≈z(t)
where z(t)=1 τ≦t<τ+T_M and z(t)=0 otherwise, where T_M was described above. This approximation is particularly useful for environment crowded with reflectors resulting in a dense multipath arrival times (for example, an office or other industrial indoor environment). FIG. 21 shows an example σ_h²(t) 120 for highly dense arrival times with TOA marked by 121.

The characteristic signal waveform typically describes the shape of the result of convolution between the transmitted signal by the reference station and/or object and the receiver filters. In case the convolution result is repetitive then the characteristic signal waveform is optionally and preferably represented based on one repetition of the convolution. The characteristic signal waveform can be provided from external source, for example, a user input, or it can be predetermined. The present inventors found that it is not mandatory to use the exact shape of the signal waveform to construct the matrix R.

Thus, in various exemplary embodiments of the invention an approximation of the signal waveform is employed. The present embodiments also contemplate use of data received from other objects in the same environment for the purpose of updating the power delay profile and/or signal waveform. In these embodiments, the method adaptively learns over time the parameters that determine the power delay profile and/or waveform, particularly, but not exclusively, the parameters that depend on the environment and that are approximately the same for all objects in the environment. In various exemplary embodiments of the invention the method calculates the aforementioned λ parameter based on data pertaining to the transmission between the reference stations and other objects in the same environment.

In some embodiments of the present invention two or more signal waveform candidates are employed, and a separate matrix is calculated for each signal waveform candidate. The signal processing as further detailed hereinbelow can be executed for each matrix. In some embodiments, the signal waveform candidates are scored during the processing. In these embodiments, the signal waveform score is used for selecting one signal waveform from the candidates (e.g., the signal waveform with the highest score) and the position is estimated using the selected signal waveform. The signal waveform scoring is optionally and preferably performed independently from the position estimation.

Two or more signal waveform candidates can also be employed according to some embodiments of the present invention in a joint estimation of candidate positions and candidate signal waveforms. In these embodiments, a score is optionally and preferably given to each pair of candidate signal waveform and candidate position. The estimated position of the object can then be selected based on the score of the pairs. Typically, but not necessarily, the candidate position that is associated with the higher score is selected as the estimated position of the object.

A more detailed description of an embodiment which employs a plurality of signal waveform candidates is described hereinunder.

The present embodiments also contemplate use of a plurality of candidate signal waveforms as well as a plurality of candidate power delay profiles. This can be done in more than one way.

In some embodiments of the present invention each candidate signal waveform is assigned with a signal waveform score and each candidate power delay profile can be assigned with a candidate power delay profile score. In these embodiments, the power delay profile score and signal waveform scores can be used for selecting one power delay profile and one signal waveform from the candidates (e.g., the power delay profile and signal waveform with the highest scores) and the position is optionally and preferably estimated using the selected power delay profile and signal waveform.

In some embodiments of the present invention the method a joint estimation of candidate position, candidate signal waveform and candidate power delay profile is executed. In these embodiments, a score is optionally and preferably given to each set of candidate signal waveform, candidate position and candidate power delay profile. The estimated position of the object can then be selected based on the score of the sets. Typically, but not necessarily, the candidate position that is associated with the higher score among the sets is selected as the estimated position of the object.

Generally, the estimation accuracy improves as the approximation is more accurate. For example, in some embodiments of the present invention the characteristic signal waveform describes a Sinc function multiplied in the time domain by a window centered around the Sinc peak. Examples of such window functions are rectangular window, Hanning window, Hamming window, Blackman window, Kaiser window or Chebchev window. FIG. 22 presents an example of a signal waveform 130 composed by Sinc function multiplied by a Hamming window and FIG. 23 presents an example of Dirichlet signal waveform 140.

The matrix R can be based on an autocorrelation matrix of the signal and is preferably a function of the position of the object. A representative and non-limiting example of an autocorrelation matrix suitable for the present embodiments is provided in the Examples section that follows (see Appendix B of Example 1).

When the signal is in the frequency domain, the autocorrelation matrix is different than when the signal is in the time domain. The relation between the autocorrelation matrices can be formulated as follows. Let y be the received vector in the time domain, and let z be its DFT vector. Thus, z=F y where F is a DFT matrix corresponding to the vector y. The autocorrelation of z is given by E{zz^H}=F E{yy^H}F^H, where E{yy^H} is the autocorrelation of the received signal in the time domain which as explained is obtained by the power delay profile and the transmitted signal waveform. The matrix F can be used, according to some embodiments of the present invention for constructing an objective function, which can be expressed as a sum of powers, for example, a polynomial function. This objective function can be used in subsequent operations for estimating the position of the object. Additional examples of calculations in the frequency domain are provided in Example 6 of the Examples section that follows.

The method proceeds to 13 at which a direct estimate of the position of the object is calculated based, at least in part, on the processed signals.

The term “direct position estimate” can be better understood by first considering the term “indirect position estimate.”

As used herein, an “indirect position estimate” is a procedure in which the position of an object is estimated based on one or more previously estimated time-of-arrival or range values, wherein the position of an object is estimated only after the one or more time-of-arrival or range values is/are estimated and explicitly expressed.

As used herein, “direct position estimate” refers to any procedure suitable for estimating the position of an object from signals transmitted between the object and a plurality of reference station, with the provision that the procedure used for estimating the position is other than an indirect position estimate.

In some embodiments of the present invention the direct position estimate is a procedure in which the position is estimated from the received signal directly without applying preprocessing in a manner that position information is lost.

The estimation typically employs a maximum likelihood (ML) procedure. Broadly speaking, the ML procedure involves the calculation of argmin_p{S(p)} or argmax_p{S(p)}, where S is a function of the (unknown) position p. In various exemplary embodiments of the invention S comprises contributions from each signal in a statistically independent manner. Denoting the sampled signal for the kth reference station by the vector y^k, a representative and non-limiting example for the function S can be written as:

$S=\sum _{k=0}^{K-1}\ln \left(R^k\left({\tau }^k\right)\right)+{\left(y^k\right)}^H{\left(R^k\left({\tau }^k\right)\right)}^{-1}{y}^k$
where the superscript H denotes the conjugate transpose operation, τ^k is given by:

${\tau }^k=\frac{1}{c}\sqrt{{\left(p_{\mathrm{ref}}^k-p\right)}^T\left(p_{\mathrm{ref}}^k-p\right)},$
c is the propagation speed of the signal (approximately the speed of light), p^k_ref is the (known) position of the kth reference station, and the superscript T denotes the transpose operation.

The matrix R^k optionally and preferably is calculated using a noise variance σ_v². The noise variance can be obtained in more than one way. In some embodiments of the present invention σ_v² is calculated long before the signal arrives. In some embodiments of the present invention σ_v² is calculated based on a segment of the signal as further detailed hereinbelow.

In some embodiments of the present invention σ_v² is obtained in the following manner.

Let y be the received vector in the time domain which is given by
y=Cz+n
where n is an additive noise vector with length N, z is a vector of L multipath samples, and C is a convolution matrix of the signal waveform. Specifically, each column in C is a cyclic shifted representation of the signal waveform.

Using singular value decomposition C=UΩV^T. Let w be the projection of y to the subspace that spans the multipath z:
w=U^Ty+n=U^TUΩΩ^Tz+U^Tn=ΩV^Tz+U^Tn

The dimension of w is N but only the L first elements contain the signal and the remaining N-L elements contain the noise. Let {tilde over (y)}={w₀, w₁, . . . , w_L-1} be the first L elements of w. Thus, {tilde over (y)} contains all the signal energy and L noise samples. The noise variance can therefore be estimated as

${\sigma }_v^2=\frac{{y}^2-{\tilde{y}}^2}{N-L}$
where ∥x∥ represents the norm or vector x. The normalization is by N-L since ∥{tilde over (y)}∥² contains also L noise samples.

The calculation of the correlation matrix optionally and preferably includes estimating signal attenuation. This can be done by calculating an overall factor α for multiplying the power delay profile. Let E_s be the known signal energy for α=1, then α is given by

$\alpha =\frac{{\tilde{y}}^2-{\sigma }_v·L}{E_s}$

In some embodiments of the present invention the processing operation 12 comprises partitioning each signal into at least a first segment containing noise and a second segment containing at least a multipath signal corresponding to a plurality of paths between the object and a respective reference station. Thus, for example, denoting the sampled signal for a given reference station by the vector y=[y₀, y₁, . . . , y_L-1], where L is the length of the vector, the method can partition y into [y₀^T, y₁^T]^T, where y₀ is a vector representing the first segment, y₁ is a vector representing the second segment. The length (number of components) of y₀ is denoted N_τ, and the length of y₁ is therefore L−N_τ. The signal can be partitioned, for example, by employing a coarse energy search procedure. Typically, the energy corresponding to the second segment is higher than the energy corresponding to the first segment.

In some embodiments of the present invention the second segment y₁ is further divided into a transient sub-segment y_1,0, and a generally steady state sub-segment y_1,1. The length of y_1,0 can equal the number d of samples within the time interval T_D which is the time that includes the energetic part of the signal waveform, for example about 98% of the signal waveform is included in duration T_D. The length of y_1,1 can be L−N_τ-d.

The method proceeds by defining a characteristic multidimensional width for each of the segments. Typically, the widths are defined by modeling each segment using a localized function. Representative examples of localized functions suitable for the present embodiments, including, without limitation, a Gaussian function, a Lorentzian function, a hyperbolic secant function, a logistic distribution, a Fourier transform and a wavelet transform. In various exemplary embodiments of the invention a Gaussian function is employed both for y₀ and for y₁.

It was found by the present inventors that the characteristic multidimensional width of the first segment y₀ can be based on the noise variance σ_v² which can be calculated, for example, by averaging the energy over the entire length of y₀ or portion thereof. Alternatively, σ_v² can be predetermined. For different reference stations, different noise width σ_v² can be defined. Alternatively, the same noise width can be employed for two or more (e.g., all) reference stations.

The characteristic multidimensional width of y₁ can be defined using an autocorrelation matrix R₁ of y₁ as a function of the (unknown) position p of the object. In various exemplary embodiments of the invention different autocorrelation matrices are defined for different reference stations. Yet, use of the same autocorrelation matrix for two or more, e.g., all the reference stations, not excluded from the scope of the present invention. The multidimensional width of y₁ can be defined as the determinant of the autocorrelation matrix R₁.

In embodiments in which the segmentation operation is employed, the calculation of direct estimate of the position is also based on the segments y₀ and y₁. For example, in these embodiments the function S can be written as:

$S=\sum _{k=0}^{K-1}{N}_{\tau }^k\ln \left({\sigma }_v^2\right)+\frac{1}{{\sigma }_v^2}{y_0^k}^2+\ln \left(R_1^k\left({\tau }^k\right)\right)+{\left(y_1^k\right)}^H{\left(R_1^k\left({\tau }^k\right)\right)}^{-1}{y}_1^k.$
where R₁^k(τ^k)=E{y₁^k(y₁^k)^H|τ^k} is the autocorrelation matrix of y₁^k and is a function of the power delay profile, the signal waveform and possibly the noise variance, σ_v², and possibly the signal attenuation factor α.

When the sub-segments y_1,0 and y_1,1 are defined, the expression for S can be further simplified as

$S=\sum _{k=0}^{K-1}\gamma +\frac{1}{{\sigma }_v^2}{y_0^k}^2+{\left(y_{1,0}^k\right)}^H{\left(R_{1,0}^k\left({\tau }^k\right)\right)}^{-1}{y}_{1,0}^k+{{\tilde{y}}_{1,1}^k}^2$
where R_1,0^k(τ^k)=E{y_1,0^k(y_1,0^k)^H|τ^k}, {tilde over (y)}_1,1^k is the result of a linear convolution between y_1,1^k and 1/√{square root over (λ_i)} where λ_i are the eigenvalues of R_1,1^k(τ^k)=E{y_1,1^k(y_1,1^k)^H|τ^k=0} and

${\gamma }^k=N_{\tau }^K\left(\ln \left({\sigma }_v^2\right)-\frac{1}{L_k}\sum _{i=0}^{L-1}\ln \left({\lambda }_i\right)\right)$
where L_k is the length of the received vectors y^k.

According to some embodiments of the present invention at least one of the signals is a diversity signal having diversity measurements, and the method further comprises processing the diversity measurements. In these embodiments, the estimate is based on the processed diversity measurements.

As used herein, a “diversity measurement” refers to an output of a single channel, and “diversity measurements” refer a plurality of outputs of a respective plurality of channels, wherein no two channels are identical thereamongst. Two or more of the channels may or may not be correlated.

It is recognized that each channel provides output in the form of a signal. Thus, although the term “diversity signal” is in the singular form, it is used for representing a set of signal components, wherein each signal component corresponds to a single channel output.

According to some embodiments of the present invention the diversity comprises at least one of: time diversity, frequency diversity, interference diversity, antenna polarization diversity and space diversity.

For example, in some embodiments, the signals are transmitted by a Multiple-Input Multiple-Output (MIMO) system or a Single Input Multiple Output (SIMO) system or a Multiple Input Single Output (MISO) system which generate space diversity using multiple antenna elements. A MIMO, MISO or SIMO system employs multiple transmitters and multiple receivers that yield multiple channels. The correlation between the channel realizations depends on the space distance between the antennas at each side (transmitter and receiver). As the antennas are more separated in space the correlation is weaker. When the antennas are separated by more than half the carrier frequency wavelength the channels can be assumed statistically independent.

In some embodiments of the present invention the diversity is formed by frequency hopping. Frequency hopping is a known technique utilized to achieve frequency diversity and/or interference diversity. The multipath varies with changes in the carrier frequency. Thus, transmitting in several carrier frequencies (frequency hopping) yields multiple multipath realizations. The larger the difference between the carrier frequencies the weaker is the correlation between the received realizations.

The diversity measurements can comprise data received from a plurality of statistically independent channels of a respective reference station. A representative example of use of such diversity measurements is described in Example 1 of the Examples section that follows. The diversity measurements can alternatively comprise data received from a plurality of correlated channels of a respective reference station. A representative example of use of such diversity measurements is described in Example 3 of the Examples section that follows.

In some embodiments of the present invention the position estimation is performed using a set of power delay profile candidates. For example, each individual power delay profile candidate can correspond to a different multipath environment (e.g. indoor, outdoor etc). This is particularly useful when the specific power delay profile environment in which the object and reference stations placed is unknown. In these embodiments scenario a joint direct position estimation and power delay profile can be employed.

For example, let Σ={σ₀²(t), σ₁²(t), . . . , σ_J-1²(t)} be a set of J power delay profile candidates. The estimation typically employs a maximum likelihood (ML) procedure. Broadly speaking, the ML procedure involves the calculation of arg min_p,σi2(t)S(p,σ_i²(t)) where S is a function of the (unknown) position p and the power delay profile σ_i²(t). Denoting the sampled signal for the kth reference station by the vector y^k, a representative and non-limiting example for the position estimation, {circumflex over (p)}, in this case is:

$\hat{p}=\arg {\min }_p\sum _{k=0}^{K-1}{\min }_{{\sigma }_i^2\left(t\right)\in \Sigma }\ln \left(R^k\left({\tau }^k,{\sigma }_i^2\right)\right)+{\left(y^k\right)}^H{\left(R^k\left({\tau }^k,{\sigma }_i^2\right)\right)}^{-1}{y}^k$
where in the inner minimization σ_i²(t) can be chosen from the set Σ, the superscript H denotes the conjugate transpose operation, τ^k is given by:

${\tau }^k=\frac{1}{c}\sqrt{{\left(p_{\mathrm{ref}}^k-p\right)}^T\left(p_{\mathrm{ref}}^k-p\right)},$
c is the speed of the signal propagation (approximately the speed of light), p^k_ref is the (known) position of the kth reference station and the superscript T denotes the transpose operation.

In some embodiments of the present invention there exists at least one received signal for which the power delay profile used to process the signal is chosen from a set of power delay profile candidates and then the processed signal is used for the position estimation. Denote by {circumflex over (σ)}_k²(t) the chosen power delay candidate for the kth received vector denoted by y^k. A representative example of choosing {circumflex over (σ)}_k²(t) is given by

${\hat{\sigma }}_k^2\left(t\right)=\arg {\max }_{{\sigma }_i^2\left(t\right)\in \Sigma }{\int }_{\tau }\exp \left(-\ln \left(R^k\left({\tau }^k,{\sigma }_i^2\right)\right)-{\left(y^k\right)}^H{\left(R^k\left({\tau }^k,{\sigma }_i^2\right)\right)}^{-1}{y}^k\right)d\tau .$
Then the position estimation can be estimated as:

$\hat{p}=\arg {\min }_p\stackrel{K-1}{\sum _{k=0}}\ln \left(R^k\left({\tau }^k,{\hat{\sigma }}_k^2\right)\right)+{\left(y^k\right)}^H{\left(R^k\left({\tau }^k,{\hat{\sigma }}_k^2\right)\right)}^{-1}{y}^k$

In some embodiments of the present invention the position estimation is performed using a set of signal waveform candidates. In these embodiments the signal waveform can also be incorporated into the direct position estimation.

For example, let χ={x₀(t), x₁(t), . . . , x_I−1(t)} be a set of I signal wave form candidates x_i(t)_i. A representative and non-limiting example for position estimation in this embodiment is given by:

$\hat{p}=\arg {\min }_p\sum _{k=0}^{K-1}{\min }_{x_i\in x,{\sigma }_i^2\left(t\right)\in \Sigma }\ln \left(R^k\left({\tau }^k,{\sigma }_i^2,x_i\right)\right)+{\left(y^k\right)}^H{\left(R^k\left({\tau }^k,{\sigma }_i^2,x_i\right)\right)}^{-1}{y}^k$
where the inner minimization is on the set of signal waveform candidates and power delay profile candidates.

It is to be understood that it is not intended to limit the scope of the present invention to the calculation of position estimation using the chosen power delay candidate, since the power delay can be provided as output, without calculating position estimation, or it can be used for calculating other physical quantities such as, but not limited to, TOA. In some embodiments of the present invention both the position estimation and the TOA are calculated using the chosen power delay.

The present inventors contemplate a scenario in which two or more of the reference stations (e.g., all) are synchronized, and the object is synchronized with the synchronized reference stations. In some of these embodiments the signal transmitted between the object and each of the synchronized reference stations is a one-way signal. A representative example of position estimation in this scenario is provided in Example 2 of the Examples section that follows. In other embodiments the signal transmitted between the object and each of the synchronized reference stations is a two-way signal.

As use herein, “one-way signal” refers to a signal transmitted by a first entity and received by a second entity, wherein the second entity does not transmit a response signal back to the first entity.

As use herein, “two-way signal” refers to a pair of signals, one transmitted by a first entity and received by a second entity, and another transmitted by the second entity (typically, but not necessarily, in response to the first signal) and received by the first entity.

In the above description each of the first and second entities can be the object or a reference station. When the first entity is the object the second entity is a reference station, and when the first entity is a reference station the second entity is the object.

The present inventors also contemplate a scenario in which two or mote the reference stations (e.g., all) are synchronized, and the object is unsynchronized with any reference station. In these embodiments the signal transmitted between the object and each of the synchronized reference stations is a one-way signal, and the number of the reference stations is at least d+1, where d is the dimensionality of the estimated position. Specifically, when it is desired to have a two-dimensional position estimate, the number of the reference stations is at least 3, and when it is desired to have a three-dimensional position estimate, the number of the reference stations is at least 4.

The present inventors also contemplate a scenario in which two or more of the reference stations are unsynchronized both thereamongst and with the object. In these embodiments, the signal transmitted between the object and each of the unsynchronized reference stations is a two-way signal.

Representative examples of position estimators in various synchronization scenarios are provided in Examples 4 and 5 of the Examples section that follows.

Reference is now made to FIG. 18 which is a schematic illustration of a system 50 for estimating a position of an object 52, according to some exemplary embodiments of the present invention. System 50 comprises a plurality of reference stations 54 of known positions. Reference stations 54 can be positioned indoors or outdoors. In the representative illustration of FIG. 18, an indoor environment 56 with a plurality of obstacles 58 distributed therein, so that the signals (not shown) transmitted between object 52 and stations 54 are generally multipath signals.

System 50 also comprises a data processor 60 configured for executing the method described herein. Data processor 60 can be a general purpose computer or dedicated circuitry. Data processor 60 communicates with stations 54 and/or object 52 and receives data pertaining to the multipath signals. Data processor 60 can be located near or in one of the reference stations 54 or within object 52 or anywhere in environment 56 or at a remote location, as desired.

The signals can be transmitted from reference stations 54 and received by object 52 or from object 52 and received by reference stations 54. In some embodiments of the present invention one-way signals are employed, wherein object 52 transmits signals and at least one reference station 54 receives signals, or alternatively at least one reference station 54 transmits signals and the object 52 receives signals. In some embodiments of the present invention two-way signals are employed, wherein object 52 both receives and transmits signals, and at least one of reference stations 54 both receives and transmits signals.

According to an aspect of some embodiments of the present invention there is provided a method of estimating a TOA of a signal transmitted between an object and a reference station.

The method comprises receiving a multipath signal having a plurality of components respectively corresponding to a plurality of paths between the object and the reference station. The method further comprises processing the signal using a matrix featuring a characteristic power delay profile and a characteristic signal waveform and optionally a noise variance to provide a processed signal, as further detailed hereinabove.

In some embodiments of the present invention two or more power delay profile candidates are employed, as further detailed hereinabove. In some embodiments of the present invention the power delay profiles are scored during the processing. In these embodiments, the power delay profile score is used for selecting one power delay profile from the candidates (e.g., the power delay profile with the highest score) and the TOA is estimated using the selected power delay profile.

Also contemplated are embodiments in which a model describing the relation between the power delay profile and the distance from the object to a respective reference station is obtained, as further detailed hereinabove. A representative and non limiting example for these embodiments is provided in the Examples section that follows (see Example 8).

Two or more power delay profile candidates can also be employed according to some embodiments of the present invention in a joint estimation of candidate TOA/distance and candidate power delay profiles, as further detailed hereinabove.

In some embodiments of the present invention two or more signal waveform candidates are employed, and a separate matrix is calculated for each signal waveform candidate, as further detailed hereinabove. Thus, signal processing can be executed for each matrix. In some embodiments, the signal waveform candidates are scored during the processing. In these embodiments, the signal waveform score is used for selecting one signal waveform from the candidates (e.g., the signal waveform with the highest score) and the TOA is estimated using the selected signal waveform. The signal waveform scoring is optionally and preferably performed independently from the TOA estimation.

Two or more signal waveform candidates can also be employed according to some embodiments of the present invention in a joint estimation of candidate TOA/distance and candidate signal waveforms, as further detailed hereinabove.

The present embodiments also contemplate use of a plurality of candidate signal waveforms as well as a plurality of candidate power delay profiles, as further detailed hereinabove. For example, a joint estimation of candidate TOA, candidate signal waveform and candidate power delay profile can be executed. In these embodiments, a score is optionally and preferably given to each set of candidate signal waveform, candidate TOA and candidate power delay profile. The estimated TOA of the object can then be selected based on the score of the sets. Typically, but not necessarily, the candidate TOA that is associated with the higher score among the sets is selected as the estimated TOA.

The method, according to some embodiments of the present invention, can comprise receiving more than one received signal. In these embodiments, the TOA estimation is performed by utilizing the plurality of received signals. Optionally and preferably one or more of the received signals is a diversity signal as further detailed hereinabove.

The method further comprises calculating an estimate of the TOA of the signal based, at least in part, on the processed signal.

For example, the estimation of a TOA, τ, can employ an ML procedure wherein argmin_τ{S_t(τ)}, where S_t is a function of the (unknown) TOA. A representative example for the function S_t can be written as:

$S_t\left(\tau \right)=\sum _{k=0}^{K-1}\ln \left(R^k\left(\tau \right)\right)+{\left(y^k\right)}^H{\left(R^k\left(\tau \right)\right)}^{-1}{y}^k.$
where y⁰, y¹, . . . , y^K−1 are received vectors at either the object or the reference station and R^k(τ) is the matrix featuring a characteristic power delay profile, a characteristic signal waveform and optionally a noise variance. In various exemplary embodiments of the invention the vectors y⁰, y¹, . . . , y^K−1 are statistically independent.

In some embodiments of the present invention the method comprises partitioning the signal into at least a first segment containing noise and a second segment containing at least a multipath signal corresponding to a plurality of paths between the object and a the reference station, as further detailed hereinabove. In some embodiments of the present invention the method further comprises defining a characteristic multidimensional width for each segment, as further detailed hereinabove. In these embodiments the estimation of a TOA, τ, is based on the segments. A representative example for the function S_t in these embodiments can be written as:

$S_t\left(\tau \right)=\sum _{k=0}^{K-1}{N}_{\tau }\ln \left({\sigma }_v^2\right)+\frac{1}{{\sigma }_v^2}{y_0^k}^2+\ln \left(R_1\left(\tau \right)\right)+{\left(y_1^k\right)}^H{\left(R_1\left(\tau \right)\right)}^{-1}{y}_1^k.$

This expression for S_t(τ) is preferred from the standpoint of accuracy. When a reduced complexity is desired, S_t(τ) is optionally and preferably expressed in terms of the sub-segments y_1,0 and y_1,1. This can be done, for example, by partitioning the matrix R₁ to form a block matrix having a first on-diagonal block R_1,0 that corresponds to y_1,0 and a second on-diagonal block R_1,1 that corresponds to y_1,1. The sub-segments y_1,0 and y_1,1 and the corresponding on-diagonal blocks R_1,0 and R_1,1 can then be used to define an expression for S_t(τ). For example,

$S_t\left(\tau \right)=\sum _{k=0}^{K-1}\gamma +\frac{1}{{\sigma }_v^2}{y_0^k}^2+{\left(y_{1,0}^k\right)}^H{R}_{1,0}^{-1}\left(\tau -N_{\tau }\right)y_{1,0}^k+{{\tilde{y}}_{1,1}^k}^2$
where, ∥{tilde over (y)}_1,1∥² is given by the approximation y_1,1^HR_1,1⁻¹(N_τ)y_1,1≈∥{tilde over (y)}_1,1∥², γ is defined as

$\gamma =N_{\tau }\left(\ln \left({\sigma }_v^2\right)-\frac{1}{L}{\sum }_{i=0}^{L-1}\ln \left({\tilde{\lambda }}_i\right)\right),$
and {tilde over (λ)}_i are the eigenvalues of R_1,1(N_τ=0).

In some embodiments of the present invention the TOA estimation is performed using a set of power delay profile candidates. For example each individual power delay profile corresponds to a different multipath environment (e.g. indoor, outdoor etc). This is particularly useful when the specific power delay profile environment in which the object and reference station are placed is unknown. In these embodiments a joint TOA and power delay profile estimation can be employed.

For example, let Σ={σ₀²(t), σ₁²(t), . . . , σ_J-1²(t)} be a set of J power delay profile candidates. The estimation typically employs a maximum likelihood (ML) procedure. Broadly speaking, the ML procedure involves the calculation of arg min_τ,σi2(t)S_t(τ,σ₁²(t)) where S is a function of the (unknown) position p and the power delay profile σ_i²(t). Denoting the sampled signal for the kth reference station by the vector y^k, a representative and non-limiting example for the TOA estimation, {circumflex over (τ)}, in this case is given by

$\hat{\tau }=\arg {\min }_{\tau }\stackrel{K-1}{\sum _{k=0}}{\min }_{{\sigma }_i^2\left(t\right)\in \Sigma }\ln \left(R\left(\tau ,{\sigma }_i^2\right)\right)+{\left(y^k\right)}^H{\left(R\left(\tau ,{\sigma }_i^2\right)\right)}^{-1}{y}^k$
where in the inner minimization σ_i²(t) is chosen from the set Σ.

In some embodiments of the present invention the TOA estimation is performed using a set of signal waveform candidates. In these embodiments, the signal waveform can also be incorporated into the TOA estimation.

For example, let χ={x₀(t), x₁(t), . . . , x_I−1(t)} be a set of I signal wave form candidates x_i(t). A representative and non-limiting example for TOA estimation for this scenario is given by:

$\hat{\tau }=\arg {\min }_{\tau }\sum _{k=0}^{K-1}{\min }_{x_i\in x}\ln \left(R\left(\tau ,x_i\right)\right)+{\left(y^k\right)}^H{\left(R\left(\tau ,x_i\right)\right)}^{-1}{y}^k$

In some embodiments of the present invention a joint estimation of candidate TOA, candidate signal waveform and candidate power delay profile is executed. In these embodiments, a score is optionally and preferably given to each set of candidate signal waveform, candidate TOA and candidate power delay profile. The estimated TOA can then be selected based on the score of the sets. Typically, but not necessarily, the candidate TOA that is associated with the higher score among the sets is selected as the estimated TOA.

For example, let Σ={σ₀²(t), σ₁²(t), . . . , σ_J-1²(t)} be a set of J power delay profile candidates and let χ={x₀(t), x₁(t), . . . , x_I−1(t)} be a set of I signal wave form candidates x_i(t). A representative and non-limiting example for TOA estimation for this scenario is given by:

$\hat{\tau }=\arg {\min }_{\tau }\sum _{k=0}^{K-1}{\min }_{x_i\in x,{\sigma }_i^2\left(t\right)\in \Sigma }\ln \left(R\left(\tau ,{\sigma }_i^2,x_i\right)\right)+{\left(y^k\right)}^H{\left(R\left(\tau ,{\sigma }_i^2,x_i\right)\right)}^{-1}{y}^k$
where the inner minimization is on the set of signal waveform candidates and power delay profile candidates.

When the calculation is employed in the frequency domain, the cost function to be minimized is optionally and preferably expressed as a sum of powers, for example, a polynomial function, and the minimization is executed by means of a transform to the frequency domain, e.g., via FFT or the like. A representative example of such minimization procedure is provided in Example 6 of the Examples section that follows.

Once the TOA or propagation delay is estimated, an indirect position estimation can be executed for estimating the position of the object. This can be done using any technique known in the art. Typically, an optimization procedure is employed. For example, the TOA estimates can be used to calculate a set of distances and a minimization procedure can be executed to minimize the error associated with the correspondence of these distances to the distances obtained from the estimated position. A representative and non-limiting example of a minimization procedure includes minimizing the sum or weighted sum of squared errors between the set of distances and the distances corresponding to the estimated position.

According to an aspect of some embodiments of the present invention there is provided a method of estimating a propagation delay of a signal transmitted between an object and a reference station.

The method comprises receiving a multipath signal having a plurality of components respectively corresponding to a plurality of paths between the object and the reference station. The method further comprises processing the signal using a matrix R featuring a characteristic power delay profile, a characteristic signal waveform and optionally a noise variance to provide a processed signal, as further detailed hereinabove.

The method can comprise receiving more than one received signal. In these embodiments, the propagation delay estimation is performed by utilizing the plurality of received signals. Optionally and preferably one or more of the received signals is a diversity signal as further detailed hereinabove.

The method further comprises calculating an estimate of the propagation delay of the signal based, at least in part, on the processed signal. A representative example of a procedure suitable for estimating the propagation delay between the object and a reference station based on two way signals is provided in Example 7 of the examples section that follows.

The present embodiments contemplate use of a model which describes the relation between the power delay profile and the distance from the object to a respective reference station. In these embodiments, the matrix R is expressed in terms of this model. A representative example of this embodiment is provided in Example 8 of the examples section that follows.

Once the propagation delay is estimated, an indirect position estimation can be executed for estimating the position of the object. This can be done using any technique known in the art. Typically, an optimization procedure is employed. For example, the propagation delays can be used to calculate a set of distances and a minimization procedure can be executed to minimize the error associated with the correspondence of these distances to the distances obtained from the estimate position. A representative and non-limiting example of a minimization procedure includes minimizing the sum or weighted sum of squared errors between the set of distances and the distances corresponding to the estimated position.

As used herein the term “about” refers to ±10%.

The word “exemplary” is used herein to mean “serving as an example, instance or illustration.” Any embodiment described as “exemplary” is not necessarily to be construed as preferred or advantageous over other embodiments and/or to exclude the incorporation of features from other embodiments.

The word “optionally” is used herein to mean “is provided in some embodiments and not provided in other embodiments.” Any particular embodiment of the invention may include a plurality of “optional” features unless such features conflict.

The terms “comprises”, “comprising”, “includes”, “including”, “having” and their conjugates mean “including but not limited to”.

The term “consisting of” means “including and limited to”.

The term “consisting essentially of” means that the composition, method or structure may include additional ingredients, steps and/or parts, but only if the additional ingredients, steps and/or parts do not materially alter the basic and novel characteristics of the claimed composition, method or structure.

As used herein, the singular form “a”, “an” and “the” include plural references unless the context clearly dictates otherwise. For example, the term “a compound” or “at least one compound” may include a plurality of compounds, including mixtures thereof.

It is appreciated that certain features of the invention, which are, for clarity, described in the context of separate embodiments, may also be provided in combination in a single embodiment. Conversely, various features of the invention, which are, for brevity, described in the context of a single embodiment, may also be provided separately or in any suitable subcombination or as suitable in any other described embodiment of the invention. Certain features described in the context of various embodiments are not to be considered essential features of those embodiments, unless the embodiment is inoperative without those elements.

Various embodiments and aspects of the present invention as delineated hereinabove and as claimed in the claims section below find experimental support in the following examples.

EXAMPLES

Reference is now made to the following examples, which together with the above descriptions illustrate some embodiments of the invention in a non limiting fashion. The Examples below contain enumerated mathematical expressions. For clarity of presentation, the enumeration is self-contained within each example. Unless otherwise stated, when a reference is made within a particular example to an equation by its reference numeral, such reference should be understood as a reference to the equation having that reference numeral within that particular example.

Example 1

TOA Estimate

The present Example describes a TOA estimate procedure according to some embodiments of the present invention. The approach taken in this example is to approximate the received multipath signal as a Gaussian process and derive the maximum likelihood estimator. In order to further decrease the computational load of the new algorithm, a low complexity approximation with negligible performance degradation was developed. The algorithm is useful for either single channel realization or multiple channel realizations using diversity either in time, frequency or space. When applying diversity technique a substantial performance gain is attained due to the optimal combining of the channel realizations and thus reliable TOA estimation is attainable even at low SNR. The estimator's performance can be closely predicted by a closed form analytical error expression.

Introduction

The range between a transmitter and a receiver can be determined by measuring the signal propagation time from one to the other. Towards this end, the TOA of the received signal can be used. The time resolution is inversely proportional to the signal bandwidth, hence it is typically desired to use the widest bandwidth available. The large bandwidth frequently imposes a large multipath delay spread relative to the transmitted signal duration. The best example is the Ultra-Wideband (UWB) technology which enables precise ranging due to its remarkable bandwidth.

For ranging estimation only the TOA of the first multipath component needs to be detected, which represent the line-of-sight path. In a typical multipath environment the first arrival path is often attenuated significantly due to destructive superposition with multipath reflections. In the presence of additive thermal noise, robust and accurate TOA estimation is a considerable challenge especially when the signal to noise ratio (SNR) is low. Unfavorable SNR is encountered when the range is relatively large, the transmitter is shadowed or the averaging time is limited due to the need for quick response or a rapidly changing environment. An SNR performance gain can be obtained by averaging multiple observations of the received signal all with the same multipath realization but with different noise realizations. However, in case of destructive multipath and/or low SNR the required averaging time for attaining precise TOA estimation can be very large. By applying diversity techniques the receiver can obtain several independent realizations of the multipath channel, all with the same TOA, and thus by proper processing a substantial performance improvement can be obtained which, for a given number of observations, is significantly larger than obtained by coherently combining the same channel realization (averaging). Diversity can be used if the distance between the transmitter and receiver does not change during all measurements. Examples of diversity techniques are frequency diversity, antenna polarization diversity and space diversity. When space diversity is used, the antennas are optionally separated by half a wavelength for low correlation, and it is assumed that this distance is small such that the difference between the signal time of arrivals at the antennas is negligible. Time diversity can also be used if the measurements are separated (in time) more than the coherence time of the channel.

In general, the optimal combining of the independent realizations for each specific TOA estimation algorithm is not straightforward. The simplest solution is to average the TOA estimates of each individual realization with equal weights. A more sophisticated solution is to use weighted averaging where the weights reflect the quality of each individual TOA estimate. However, it is not clear how to define a good quality criterion for each individual estimate. In many cases the estimation error is not Gaussian and has occasional outliers which must be detected and rejected. This leads, in many cases, to heuristic solutions that do not fully exploit the significant potential gain that can be obtained from diversity.

The statistics of a multipath channel have been investigated and discussed in many articles [1]-[4]. In common base-band multipath channel models for wideband transmission the multipath components are statistically independent, each component is represented by a complex gain and a random delay. The arrival delays are modeled as a Poisson process and the gains have exponential decaying variance. In the classical model [1] the distribution of each complex gain is Gaussian. Other models [3] describe the gain distribution as Nakagami when the signal bandwidth is very large.

Simple approaches for TOA estimation are based on passing the received signal through a matched filter (MF) whose output peak epoch is the TOA estimate [5]-[8]. An alternative solution is to search for a finite number of MF output peak epochs and to estimate the TOA as the earliest peak epoch [9]. Other simple solutions are based on comparing the received signal energy to a threshold and then performing a forward or backwards search with a heuristic TOA estimation criterion [10]-[16]. Also of interest is an algorithm for symbol timing synchronization [17].

Other approaches include the “frequency-domain super-resolution TOA estimation” which includes Multiple Signal Classification (MUSIC) [18], Estimation of Signal Parameters via Rotational Invariance Techniques (ESPRIT) [19] and other algorithms [20],[21]. In case the multipath is sparse then the computational load of these methods can be reduced by sampling at rate lower than Nyquist (the rate of innovation) [22],[23]. Win and Scholtz [24] introduced the generalized maximum likelihood (GML) solution. The method relies on the assumption that the number of multipath components is finite and known and all the multipath coefficients and their arrival times are jointly estimated.

An iterative approximation to the GML was introduced by [25]. The multipath coefficients are estimated one by one and removed from the received signal sequentially until the strongest remaining coefficient is below a preset threshold. The TOA is determined according to the earliest arrival path. In [26] the TOA is estimated as the instant separating the leading segment of samples which is associated with noise only, and the following segment of samples that is associated with the signal multiple arrivals. The separating instant is defined as the point where the energy changes.

The exact maximum likelihood (EML) TOA estimator is developed in Appendix A. The estimator chooses the first arrival time which maximizes the conditional probability density function of the measured channel realizations (one or more) given the first arrival time only. In practice, the a priori statistical parameters of the channel such as the Poisson average arrival rate and the signal exponential decaying variance are not known.

The present example describes a practical TOA estimation algorithm that for one or more independent channel realizations attains the EML performance in the limit when the multipath rays are infinitely dense and independent. The Inventors found that the multipath density for which the estimator essentially closes the performance gap with the EML matches the density of a typical indoor environment. The estimator of the present embodiments is superior to other known practical TOA estimators. In various exemplary embodiments of the invention the estimator does not need to know the channel characteristics accurately, thus, it is robust to various multipaths and very attractive.

System Model

The baseband system model considered in this example is depicted in FIG. 1. A wideband single pulse signal s(t), with bandwidth W_s, is transmitted through the multipath channel h(t) and observed with additive white Gaussian complex noise w(t) which has zero mean and variance σ_w². It is assumed that s(t) is known to the receiver. Let H(ƒ) be the Fourier transform of h(t). The channel is assumed to be a time invariant multipath channel with impulse response given by

$\begin{array}{cc}h\left(t\right)=\sum _{m=0}^{M-1}{\alpha }_m\delta \left(t-{\tau }_m\right)& \left(\mathrm{EQ}.1\right)\end{array}$
where M is the number of multipath components, δ(t) is the Dirac delta, α_m and τ_m are the complex coefficient and the delay of the m-th arrival multipath component, respectively. Both α_m and τ_m are statistically independent. Throughout this Example the notation of the first arrival time, τ₀, will often be replaced by τ for simplifying the exhibition. The received signal can be written as

$\begin{array}{cc}r\left(t\right)=\sum _{m=0}^{M-1}{\alpha }_{m}s\left(t-{\tau }_m\right)+w\left(t\right).& \left(\mathrm{EQ}.2\right)\end{array}$

Typically, in wide band transmission, the time duration of s(t) is shorter than the channel delay spread. At the receiver, the signal r(t) is passed through a pre-sampling low pass filter (LPF) whose frequency response and impulse response are B(ƒ) and b(t), respectively. The filters output y(t) is sampled at rate ƒ_s=1/T_s=2W_b where W_b is the bandwidth of B(ƒ) which satisfies W_b≧W_s. The resulting sampled sequence is denoted by y_n and can be written as
y_n=q_n+v_n  (EQ. 3)

In (3) v_n are the noise samples obtained by filtering w(t) with b(t) and q_n is the sampled received multipath signal. Let x(t) be the result of the convolution between s(t) and b(t) then q_n is obtained by filtering h(t) with x(t) and sampling the result at a rate of 1/T, thus

$\begin{array}{cc}{q}_n={\int }_0^{{\mathrm{nT}}_s}h\left(\lambda \right)x\left({\mathrm{nT}}_s-\lambda \right)d\lambda =\sum _{m=0}^{M-1}{\alpha }_{m}x\left({\mathrm{nT}}_s-{\tau }_m\right)& \left(\mathrm{EQ}.4\right)\end{array}$

Denote by T_D the time duration of x(t). The noise samples v_n are assumed to be independent complex circular Gaussian random variables with zero mean and variance σ_v².

The following notation will be used. ƒ_x(x) is used to represent the probability density function of x. The subscript x is usually omitted to simplify the exhibition. The superscripts *, T and H represent conjugate, transpose and conjugate transpose, respectively.

Derivation of Estimation Algorithms

Most models assume that the multipath components are clustered in time. In the present example, clusters are considered to be separated when the inter cluster gaps are longer than the transmitted signal duration. The preliminary step of the estimation algorithm is a coarse search for the samples containing the first arrival multipath cluster. This can be done by energy detection using a short sliding window over the samples y_n. The beginning of the observation widow is set to be several times the sliding window length prior to the first time that the energy exceeds a rough threshold to insure that the first arrival time is included in the window. The observation window ends when the sliding window energy falls below a threshold level. It was found by the present inventors that the performance of the TOA detector described in this example is not affected significantly if the window includes a portion of noise samples following the first cluster as long as this portion is small relative to the cluster duration.

In this example, it is assumes, without loss of generality, that the beginning of the observation window is at time t=0. Let y=[y₀ y₁ . . . y_L-1]^T be a vector of length L containing the observation window samples which are given by y_n=y(nT_s). Note that τ, the first arrival time with respect to the beginning of the observation window, is a real number. Let N_τ be the integer index of the first sample within the observation vector y that contains the multipath signal, thus, (N_τ−1)T_s<τ≦N_τT_s. In FIG. 2 the first arrival time to be estimated, τ, is pointed out for the case of a single path channel (no multipath) i.e., h(t)=δ(t−τ) and y(t)=x(t−τ). The continuous received signal y(t) and its samples y_n along the observation window are shown in the figure. Note that τ is not a point of high energy in the pulse, just a reference point. For example a detector that looks for the pulse maximum does not find it at τ but at a different point τ+T_D/2. However the constant T_D/2 is typically known and can be subtracted since in the present example the pulse shape is known to the receiver.

The observation vector y can be partitioned into two segments y=[y₀^T y₁^T]^T such that the first segment y₀=[y₀ y₁ . . . y_Nτ−1]^T contains only noise samples, and the second segment y₁=[y_Nτ y_Nτ+1 . . . y_L-1]^T contains the received multipath first cluster plus noise samples as demonstrated in FIG. 3. It is assumed that the segment y₁ terminates at the end of the first cluster. The reason that the vector is truncated at this point is that the rest of the samples do not provide any useful information regarding the TOA.

The estimate of the received signal first arrival time τ within the observation window will now be described. Let v=[v_Nτ v_Nτ+1 . . . v_L-1]^T be the noise samples in y₁ then from (4) we write y₁ as a sum of statistically independent vectors so that

$\begin{array}{cc}{y}_1=\sum _m^{M_y-1}{\eta }_m+v& \left(\mathrm{EQ}.5\right)\end{array}$
where η_m=α_m[x(N_τT_s−τ_m) x((N_τ+1)T_s−τ_m) . . . x((L−1)T_s−τ_m)]^T and M_y is the number of multipath components included in y₁. The vectors η_m are statistically independent since α_m and τ_m are statistically independent. The central limit theorem (CLT) for random vectors [28] states that the sum of independent random vectors is a Gaussian vector when the number of the summed vectors approaches infinity. The number of non zero elements in the summation for each term of y₁ is proportional to the density of the multipath delays (τ_m) within the duration of x(t). It is assumes that the density is large enough to justify the use of the Gaussian approximation. In the following section the Maximum Likelihood (ML) estimator is developed based on the Gaussian model. The simulation results show that in practice only three terms were enough to justify the Gaussian approximation. The complexity of the resulting estimator is much lower than that of the EML (the ML solution for the exact channel model). A simplified solution which has an even further reduced complexity and is practical for implementation is developed hereinunder.
Gaussian Model Maximum Likelihood TOA Estimator (GE)

The TOA ML estimator, assuming the receiver obtains K independent observation vectors y⁰,y¹, . . . , y^K−1 all with the same TOA, is given by

$\begin{array}{cc}\hat{\tau }=\arg \underset{\tau }{\max }f\left(y^0,y^1,\dots ,y^{K-1}|\tau \right)=\arg \underset{\tau }{\max }\prod _{k=0}^{K-1}f\left(y^k|\tau \right)& \left(\mathrm{EQ}.6\right)\end{array}$

Typically, K=1 but diversity can be used to obtain K>1 vectors. For each single observation vector y (superscript is omitted to simplify the exhibition), the segments prior and following time τ are statistically independent hence we can write
ƒ(y|τ)=ƒ(y₀|σ)ƒ(y₁|τ)  (EQ. 7)

The noise samples in y₀ are independent complex Gaussian random variable with zero mean, variance σ_v² and independent real and imaginary components, thus

$\begin{array}{cc}f\left(y_0|\tau \right)=\prod _{n=0}^{N_{\tau }-1}\frac{1}{{\mathrm{\pi \sigma }}_v^2}{e}^{-\frac{{y_n}^2}{{\sigma }_v^2}}& \left(\mathrm{EQ}.8\right)\end{array}$

According to (5) E{y₁}=0 since E{α_m}=0 and E{v_n}=0. The L−N_τ samples in y₁ are assumed to have Gaussian distribution hence their distribution is given by

$\begin{array}{cc}f\left(y_1|\tau \right)=\frac{1}{{\pi }^{\left(L-N_{\tau }\right)}R_1\left(\tau \right)}{e}^{-y_1^H{R}_1^{-1}\left(\tau \right)y_1}& \left(\mathrm{EQ}.9\right)\end{array}$
where R₁(τ) is the autocorrelation matrix of y₁ which is dependent on the value of τ. The autocorrelation of y is derived in Appendix B. Denote by R₁^i,j(τ) the element in the i-th row and j-th column of R₁(τ) which is given in Appendix B by

$\begin{array}{cc}{R}_1^{i,j}\left(\tau \right)=E\left\{y_{N_{\tau }+i}{y}_{N_{\tau }}^{*}+j|\tau \right\}={\int }_{t=\tau }^{\infty }{\sigma }_h^2\left(t\right)x\left(\left(N_{\tau }+i\right)T_s-t\right)x\left(\left(N_{\tau }+j\right)T_s-t\right)dt+{\sigma }_v^2{\delta }_K\left(i-j\right)& \left(\mathrm{EQ}.10\right)\end{array}$
where x(μ) initiates at μ=0, δ_K(n) is the Kronecker delta and σ_h²(t) is the channel power delay profile, which is typically the average power as a function of the delay.

By substituting (8) and (9) into (7) and then into (6), applying ln to the result, and eliminating insignificant terms the following Gaussian Model ML TOA estimator (GE) is found

$\begin{array}{cc}{\hat{\tau }}_{\mathrm{GE}}=\arg \underset{\tau }{\min }\sum _{k=0}^{K-1}{N}_{\tau }\ln \left({\sigma }_v^2\right)+\frac{1}{{\sigma }_v^2}{y_0^k}^2+\ln \left(R_1\left(\tau \right)\right)+{\left(y_1^k\right)}^H{R}_1^{-1}\left(\tau \right)y_1^k& \left(\mathrm{EQ}.11\right)\end{array}$
where ∥y₀^k∥² is the squared norm of y₀^k.

The TOA estimator in (11) was derived according to some embodiments of the present invention based on the assumption that the observation vector y is a Gaussian vector. It was found by the present inventors that the Gaussian approximation is sufficient for any multipath scenario, including dense multipath, a single path channel, and any other density of paths (e.g., only few paths within the pulse duration).

When the channel power delay profile σ_h²(t) is not known or difficult to be found, it can be approximated as constant within the first segment. By truncating the decayed multipath part only a minor part of the signal energy is lost and hence it causes insignificant degradation in performance. Let {tilde over (σ)}_h² be the approximated constant value of σ_h²(t). In cases that the multipath arrivals density is low there is a significant difference between the value of σ_h²(t) at time t=τ and at t>τ and it is preferred to increase the variance at time t=τ such that σ_h²(τ)=Υ{tilde over (σ)}_h², where Υ is a constant factor inversely proportional to the multipath density. Thus Υ≧1 and it approaches one as the multipath arrivals density grows.

If σ_v² and {tilde over (σ)}_h² are not known, they can be estimated as follows. The noise variance can be estimated by averaging the energy of the samples received long before the signal arrival time, and {tilde over (σ)}_h² can be estimated by first evaluating the ratio between the energy of y₁ and its duration and then subtracting the noise variance from the result.

It is noted that unlike conventional techniques, the accuracy of the estimator of the present embodiments is not limited by the sampling rate, provided 1/T_s≧2W, where W_s is the bandwidth of s(t) which is also the bandwidth of x(t). The sampling rate does not affect the estimator resolution since the metric to minimize is a continuous function of τ.

The complexity of the estimator given in (11) will now be evaluated. The cost function to be minimized in (11) is not necessarily convex and may have several local minimum points. In the following, a worst case scenario is taken, according to which the exhaustive search is performed on a certain grid of τ with spacing T_Δ. Let P=T_s/T_Δ be the ratio between the sampling period and the search grid spacing. In total there are LT_s/T_Δ=PL hypotheses tests of τ for an observation vector of L samples. The matrices R₁⁻¹(τ) and the values of ln(|R₁(τ)|) for τ=0,T_Δ,2T_Δ, . . . , (LP−1)T_Δ can be computed in advance for several SNR values and stored in memory. The complexity of the estimator is then governed by the computation of J=(y₁^k)^HR₁⁻¹(τ)y₁^k. For a single vector y₁ of length N the complexity of computing J (assuming R₁⁻¹(τ) is stored in memory) is O(N²+N). Per each different length of y₁ there are P evaluations of J and thus the total computational complexity for all PL hypothesizes of τ is given by

$\begin{array}{cc}O\left(K\sum _{n=1}^{L}P\left(n^2+n\right)\right)\approx O\left(\mathrm{KP}\sum _{n=1}^L{n}^2\right)=O\left(P\frac{L\left(L+1\right)\left(2L+1\right)}{6}\right)\approx O\left(\frac{{\mathrm{KPL}}^3}{3}\right).& \left(\mathrm{EQ}.12\right)\end{array}$

Thus, the complexity is proportional to L³. An approximation of the GE with significant lower complexity is described in the following section.

Approximated Gaussian Model ML TOA Estimator (AGE)

In order to simplify the implementation of the GE given in (11) an approximation to
Ω(τ)=ln(|R₁(τ)|)+y₁^HR₁⁻¹(τ)y₁  (EQ. 13)
is derived. The duration of h(t) from the first arrival until it decays to zero is assumed to be larger than the duration, T_D, of x(t). y₁ is further partitioned into two sub-segments y₁=[y_1,0^T y_1,1^T]^T where y_1,0=[y_Nτ y_Nτ+1 . . . y_Nτd−1]^T, y_1,1=[y_Nτ+d y_Nτ+d+1 . . . y_L-1]^T and d is the number of samples within the interval T_D. It is assumed that the segment y_1,1 terminates once the multipath signal level drops significantly and hence the multipath power delay profile, σ_h²(t), in the segment is approximated as constant. The segment y_1,0 is the transient response of the convolution between x(t) and h(t) where only partial correlation between them takes place and the segment y_1,1 is the steady state response where the full length of x(t) is correlated with h(t). The transient period is where the signal frontier is arriving, and is the most significant part of the entire observation window for TOA estimation. The TOA cannot effect the steady state period y_1,1. Nevertheless, this segment is used to verify that the signal has arrived and the transient is not caused by noise. In other words, the TOA estimator utilizes this segment to reduce false signal detection. The dimension of y_1,0 is fixed for a given x(t) (independent of τ) and is relatively small. The dimension of y_1,1 is L−N_τ−d which is relatively large and not fixed. It was found by the present inventors that such choice of segmentation reduces the complexity while maintaining sufficient level of performance. The length of the transient part can be varied. For a sufficiently large transient segment, the estimator is equivalent to the GE described above.

The autocorrelation matrix of y₁ can be written as

$R_1\left(\tau \right)=\left[\begin{array}{cc}{R}_{1,0}\left(\tau -N_{\tau }\right)& G\\ G& R_{1,1}\left(N_{\tau }\right)\end{array}\right]$
where R_1,0(τ−N_τ)=E{y_1,0y_1,0^H} is the autocorrelation of the transient period which from (10) is a function of τ−N_τ and R_1,1(N_τ)=E{y_1,1y_1,1^H} is the autocorrelation of the steady state period which from (10) is a function of N_τ (assuming constant power delay profile). The matrix G=E{y_1,0y_1,1^H} is the cross-correlation between the transient and steady state periods. The two sub-segments are approximated as independent, thus G=0. By substituting the latter into (13) one obtains:
Ω(τ)=y_1,0^HR_1,0⁻¹(τ−N_τ)y_1,0+ln(|R_1,0(τ−N_τ)|)+y_1,1^HR_1,1⁻¹(N_τ)y_1,1+ln(|R_1,1(N_τ)|)  (EQ. 14)

Let λ_i be the eigenvalues of R_1,1(N_τ). Since the determinant of a matrix is equal to the product of its eigenvalues then

$\begin{array}{cc}\ln \left(R_{1,1}\left(N_{\tau }\right)\right)=\sum _{i=0}^{L-N_{\tau }-d-1}\ln \left({\lambda }_i\right).& \left(\mathrm{EQ}.15\right)\end{array}$

According to the spectral decomposition theorem R_1,1⁻¹(N_τ)=UΛ⁻¹U^H where U is a matrix whose columns are the eigenvectors of R_1,1(N_τ) and Λ is a diagonal matrix with the eigenvalues λ_i of R_1,1(N_τ) along its diagonal. Let Ψ=U^Hy_1,1, then:

$y_{1,1}^H{R}_{1,1}^{-1}\left(N_{\tau }\right)y_{1,1}=y_{1,1}^{H}U{\Lambda }^{-1}{U}^H{y}_{1,1}={\Psi }^H{\Lambda }^{-1}\Psi =\sum _{i=0}^{L-N_{\tau }-d-1}{\frac{1}{\sqrt{{\lambda }_i}}\Psi \left(i\right)}^2$

According to the Parseval's theorem

$\begin{array}{cc}\sum _{i=0}^{L-N_{\tau }-d-1}{\frac{1}{\sqrt{{\lambda }_i}}\Psi \left(i\right)}^2=\sum _{n=0}^{L-N_{\tau }-d-1}{{\omega }_n}^2& \left(\mathrm{EQ}.16\right)\end{array}$
where ω_n is the result of circular convolution between the IDFT of 1/√{square root over (λ_i)} and the IDFT of Ψ. The matrix R_1,1(N_τ) is Toeplitz with dimension L−N_τ−d which is large in relation to the correlation length hence we approximate it as a Circulant matrix. As such, U^H is the discrete Fourier transform (DFT) matrix and the eigenvalues λ_i along the diagonal of Λ are equal to the DFT of the first row of R_1,1(N_τ) [29]. Hence, approximating R_1,1(N_τ) as Circulant means approximating Ψ as the DFT of y_1,1 and λ_i as the frequency response of a low pass filter (since it is the DFT of the autocorrelation of y_1,1). Consequently, 1/√{square root over (λ_i)} is high pass filter (HPF) frequency response and ω_n is approximated as the result of circular convolution between this HPF and y_1,1. Let {tilde over (y)}_1,1 be the result of a linear convolution between y_1,1 and the HPF 1/√{square root over (λ_i)}. Assuming that this filter is short enough we can replace the circular convolution with a linear convolution and then
y_1,1^HR_1,1⁻¹(N_τ)y_1,1≈∥{tilde over (y)}_1,1∥²  (EQ. 17)

The effect of changing the size of the first row of R_1,1(N_τ) according to N_τ is zero padding. In the frequency domain it results in resampling the frequency response λ_i. The summation in (15) is then linearly proportional to the dimension of R_1,1(N_τ) and can be approximated as

$\begin{array}{cc}\ln \left(R_{1,1}\left(\tau \right)\right)\approx \left(\frac{L-N_{\tau }-d}{L}\right)\Gamma & \left(\mathrm{EQ}.18\right)\end{array}$
where Γ=Σ_i=0^L−1 ln({tilde over (λ)}_i), and {tilde over (λ)}_i are the eigenvalues of R_1,1(N_τ=0). The factor

$\left(\frac{L-N_{\tau }-d}{L}\right)$
is the ratio between the length of y_1,1 and the dimension for which {tilde over (λ)}_i is evaluated. Substituting (17) and (18) into (14) yields an approximation to Ω(τ). Substituting this approximation into (11), removing irrelevant terms and reorganizing, the following approximate GE (AGE) solution is obtained

$\begin{array}{cc}{\hat{\tau }}_{\mathrm{AGE}}=\arg \underset{\tau }{\min }\sum _{k=0}^{K-1}\gamma +\frac{1}{{\sigma }_v^2}{y_0^k}^2+{\left(y_{1,0}^k\right)}^H{R}_{1,0}^{-1}\left(\tau -N_{\tau }\right)y_{1,0}^k+{{\tilde{y}}_{1,1}^k}^2& \left(\mathrm{EQ}.19\right)\end{array}$
where

$\gamma =N_{\tau }\left(\ln \left({\sigma }_v^2\right)-\frac{1}{L}\sum _{i=0}^{L-1}\ln \left({\tilde{\lambda }}_i\right)\right).$
Further simplification can be obtained when the sampling rate is slightly higher than the Nyquist rate. In this case the HPF becomes close to an all-pass, and has negligible effect and thus the metric in (19) can be simplified by using the approximation

${{\tilde{y}}_{1,1}^k}^2\approx \frac{1}{{\sigma }_{y1}^2}{y_{1,1}^k}^2$
where σ_y1² is the variance of the samples in y_1,1.

The TOA estimator in (19) is simple to implement. The dimensions of R_1,0(τ−N_τ) is only d×d, where d is the number of samples corresponding to the transmitted pulse duration (typically only a few samples). There are only P=T_s/T_Δ different matrices R_1,0⁻¹(τ−N_τ) and value of ln(|R_1,0(τ−N_τ)|) to be used in (19) for the time differences of τ−N_τ=0,T_Δ,2T_Δ, (P−1)T_Δ. They can be calculated in advance, as well as the HPF coefficients, for several SNR values and stored in memory. The complexity of (19) is governed by the calculation of (y_1,0^k)^HR_1,0⁻¹(τ−N_τ)y_1,0^k which has a complexity of O(d²+d). Thus, the complexity of (19) for K independent observations and all PL hypothesizes of τ is O(KPLd²) which grows only linearly with L. This is a remarkable reduction in complexity compared to the complexity of the GE metric (11) which grows with the third power of L. Furthermore, the complexity of the AGE is only d² times higher than the complexity of the simplest TOA estimation methods. These methods are based on comparing the signal energy to a threshold [14]. For these methods the sampling spacing is equal to the search resolution, T_Δ, and thus there are PL samples in the observation window corresponding to PL hypothesizes of τ. Hence there are PL energy calculations per channel realization thus the complexity is O(KLP). The implementation complexity of SCC is also linearly dependent on L, however, it is more complex than AGE since at each iteration the multipath complex coefficients are estimated by the least squares method which requires an inversion of a matrix with dimension equal to the number of arrivals.

Performance Analysis

The performance of the Gaussian estimator given in (11) will now be evaluated. The analysis is applicable for any multipath density and even for low density where the Gaussian approximation is inaccurate. A cost function (a function to be minimized) in the face of (11) can be written as;

$\begin{array}{cc}C\left(\tau \right)=\sum _{k=0}^{K-1}\left[\ln \left(R\left(\tau \right)\right)+{\left(y^k\right)}^H{R}^{-1}\left(\tau \right)y^k\right],& \left(\mathrm{EQ}.20\right)\end{array}$
where y^k is the k-th full observation vector (including all segments) and R(τ) is the autocorrelation of y^k with elements calculated according to (40). Using the relationship

$\sum _{k=0}^{K-1}{\left(y^k\right)}^H{R}^{-1}\left(\tau \right)y^k=\sum _{k=0}^{K-1}\mathrm{Tr}\left\{R^{-1}\left(\tau \right){y^k\left(y^k\right)}^H\right\}=\mathrm{KTr}\left\{R^{-1}\left(\tau \right)\left[\frac{1}{K}\sum _{k=0}^{K-1}{y^k\left(y^k\right)}^H\right]\right\}$
(20) is written as
C(τ)=K ln(|R(τ)|)+KTr{R⁻¹(τ){circumflex over (R)}}  (EQ. 21)
where Tr{ } denotes the trace of a matrix and

$\hat{R}=\frac{1}{K}\sum _{k=0}^{K-1}{y^k\left(y^k\right)}^H$
is the empirical autocorrelation estimation. From (21) the ML estimator can be interpreted as first estimating the autocorrelation {circumflex over (R)} based on the empirical observations and then choosing τ that yields R(τ) that best “matches” {circumflex over (R)}. The “matching” test is realized by Tr{R⁻¹(τ){circumflex over (R)}}.

The derivative of the cost function C(τ) with respect to τ is given by

$\begin{array}{cc}\begin{array}{c}g\left(\tau ,\hat{R}\right)=\frac{\partial C\left(\tau \right)}{\partial \tau }\\ =\frac{\partial }{\partial \tau }\left[K\ln \left(R\left(\tau \right)\right)+\mathrm{KTr}\left\{R^{-1}\left(\tau \right)\hat{R}\left(\tau \right)\right\}\right]\\ =\mathrm{KTr}\left[R^{-1}\left(\tau \right)\frac{\partial R\left(\tau \right)}{\partial \tau }-R^{-1}\left(\tau \right)\frac{\partial R\left(\tau \right)}{\partial \tau }{R}^{-1}\left(\tau \right)\hat{R}\right]\end{array}& \left(\mathrm{EQ}.22\right)\end{array}$

Let τ_t and R(τ_t) be the true TOA and the true autocorrelation respectively. A first-order Taylor series expansion of g(τ, {circumflex over (R)}) at τ=τ_t and {circumflex over (R)}=R(τ_t) yields
g(τ,{circumflex over (R)})≅g(τ_t,R(τ_t))+Q(τ_t,R(τ_t))[τ−τ_t]+F(τ_t,R(τ_t))[vec({circumflex over (R)})−vec(R(τ_t))]  (EQ. 23)
where vec(A) is the vectorized version of the matrix A,

$Q\left(\tau ,\hat{R}\right)=\frac{\partial }{\partial \tau }g\left(\tau ,\hat{R}\right)\mathrm{and}F\left(\tau ,\hat{R}\right)=\frac{\partial }{\partial \mathrm{vec}\left(\hat{R}\right)}g\left(\tau ,\hat{R}\right).$
For a given empirical autocorrelation matrix {circumflex over (R)}, the value τ={circumflex over (τ)} ({circumflex over (τ)} is the estimated TOA) attains the minimum of the cost function and for {circumflex over (R)}=R(τ_t) (the error free autocorrelation) the value τ=τ_t attains the minimum. Hence g({circumflex over (τ)},{circumflex over (R)})=g(τ_t,R(τ_t))=0. Substituting the latter into (23) yields
Q(τ_t,R(τ_t))[{circumflex over (τ)}−τ_t]≅−F(τ_t,R(τ_t))[vec({circumflex over (R)})−vec(R(τ_t))]  (EQ. 24)

Right multiplication by the conjugate transpose of itself at each side of equation (24) yields
Q²(τ_t,R(τ_t))[{circumflex over (τ)}−τ_t]²=F(τ_t,R(τ_t))ΘF^H(τ_t,R(τ_t))  (EQ. 25)
where Θ=[vec({circumflex over (R)})−vec(R(τ_t))][vec({circumflex over (R)})−vec(R(τ_t))]^H. Taking the expectation of (25) and dividing by Q²(τ_t,R(τ_t)) gives

$\begin{array}{cc}E\left\{{\left(\hat{\tau }-{\tau }_t\right)}^2\right\}\cong \frac{F\left({\tau }_t,R\left({\tau }_t\right)\right)E\left\{\Theta \right\}{F}^H\left({\tau }_t,R\left({\tau }_t\right)\right)}{Q^2\left({\tau }_t,R\left({\tau }_t\right)\right)}& \left(\mathrm{EQ}.26\right)\end{array}$

It has been shown [30] that

$\begin{array}{cc}E\left\{\Theta \right\}=\frac{1}{K}{R\left({\tau }_t\right)}^T\otimes R\left({\tau }_t\right)& \left(\mathrm{EQ}.27\right)\end{array}$
where {circle around (x)} is the Kronecker product.

The denominator of (26) is

$\begin{array}{cc}{Q}^2\left({\tau }_t,R\left({\tau }_t\right)\right)={\left[\mathrm{KTr}\left\{{\left[{R^{-1}\left({\tau }_t\right)\left[\frac{\partial R\left(\tau \right)}{\partial \tau }\right]}_{\tau ={\tau }_t}\right]}^2\right\}\right]}^2& \left(\mathrm{EQ}.28\right)\end{array}$
and in the nominator

$\begin{array}{cc}\begin{array}{c}F\left(\tau ,\hat{R}\right)=\frac{\partial }{\partial \mathrm{vec}\left(\hat{R}\right)}\mathrm{KTr}\left[R^{-1}\left(\tau \right)\frac{\partial R\left(\tau \right)}{\partial \tau }-R^{-1}\left(\tau \right)\frac{\partial R\left(\tau \right)}{\partial \tau }{R}^{-1}\left(\tau \right)\hat{R}\right]\\ =-K\frac{\partial }{\partial \mathrm{vec}\left(\hat{R}\right)}\mathrm{Tr}\left[R^{-1}\left(\tau \right)\frac{\partial R\left(\tau \right)}{\partial \tau }{R}^{-1}\left(\tau \right)\hat{R}\right].\end{array}\mathrm{Let}A\left(\tau \right)=R^{-1}\left(\tau \right)\frac{\partial R\left(\tau \right)}{\partial \tau }{R}^{-1}\left(\tau \right)\mathrm{then}F\left(\tau ,\hat{R}\right)=-K\frac{\partial }{\partial \mathrm{vec}\left(\hat{R}\right)}\sum _{i=0}^{N-1}\sum _{j=0}^{N=1}{A}_{i,j}{\hat{R}}_{j,i}=-\mathrm{Kvec}\left(A\left(\tau \right)\right)& \left(\mathrm{EQ}.29\right)\end{array}$
where A_i,j and {circumflex over (R)}_j,i are the elements in the i-th row and j-th column of A(τ) and {circumflex over (R)}, respectively. By substituting (27), (28) and (29) into (26) we arrive at

$\begin{array}{cc}E\left\{{\left(\hat{\tau }-{\tau }_t\right)}^2\right\}\cong \frac{{\mathrm{vec}\left(A\left({\tau }_t\right)\right)}^H\left\{R^T\left({\tau }_t\right)\otimes R\left({\tau }_t\right)\right\}\mathrm{vec}\left(A\left({\tau }_t\right)\right)}{{K\left[\mathrm{Tr}\left\{{\left[{R^{-1}\left({\tau }_t\right)\left[\frac{\partial R\left(\tau \right)}{\partial \tau }\right]}_{\tau ={\tau }_t}\right]}^2\right\}\right]}^2}& \left(\mathrm{EQ}.30\right)\end{array}$

Using linear algebra equalities [29] the nominator of (30) can be written as s

${\mathrm{vec}\left(A\left({\tau }_t\right)\right)}^H\left\{R^T\left({\tau }_t\right)\otimes R\left({\tau }_t\right)\right\}\mathrm{vec}\left(A\left({\tau }_t\right)\right)={\mathrm{vec}\left(A\left({\tau }_t\right)\right)}^H\mathrm{vec}\left(R\left({\tau }_t\right)A\left({\tau }_t\right)R\left({\tau }_t\right)\right)=\mathrm{vec}\left({A\left({\tau }_t\right)}^{H}R\left({\tau }_t\right)A\left({\tau }_t\right)R\left({\tau }_t\right)\right)=\mathrm{Tr}\left\{{\left[{R^{-1}\left({\tau }_t\right)\left[\frac{\partial R\left(\tau \right)}{\partial \tau }\right]}_{\tau ={\tau }_t}\right]}^2\right\}$
by substituting the latter result back into (30) the MSE approximation is obtained:

$\begin{array}{cc}E\left\{{\left(\hat{\tau }-{\tau }_t\right)}^2\right\}\cong \frac{1}{\mathrm{KTr}\left\{{\left[{R^{-1}\left({\tau }_t\right)\left[\frac{\partial R\left(\tau \right)}{\partial \tau }\right]}_{\tau ={\tau }_t}\right]}^2\right\}}& \left(\mathrm{EQ}.31\right)\end{array}$

In (31) the multipath statistics are realized in R(τ) thus the attained approximated MSE expression is suited for various multipath models. The expression in (31) is also the expression of the Cramer-Rao lower bound (CRLB) on the TOA estimation MSE when the received signal y is a Gaussian vector [31] which is the case when the multipath is dense or when the channel has only a single path. Thus, in these two cases the GE MSE is higher or equal to (31) and it is expected to reach the CRLB as the number of independent observation vectors (K) increases.

Results and Discussion

The GE estimator of the present embodiments was tested against the exact maximum likelihood estimator (EML) using a test channel that approximates the commonly used multipath model for which each reflection coefficient has complex Gaussian distribution and the delays have a Poisson distribution [1]. In the test channel the time axis is discrete and its duration is final, for enabling EML calculation, as described in Appendix A. Since the complexity of the EML is extremely high the observation vector was limited to 20 samples. The transmitted signal is presented in FIG. 4. The sample index of the first multipath coefficient (first arrival) was placed randomly with uniform distribution in the first quarter of the observation window, and for each of the following samples a Bernoulli random variable taking the value of zero or one, according to the average arrival rate, determined the existence of an arrival. All the existing reflection coefficients had Gaussian distribution with zero mean and constant variance σ_α². Herein, T_eff is defined the duration of the section of x(t) that contains 98% of its total energy. Thus, T_eff is considered as the energetic duration of x(t) that has the significant contribution to the result of the convolution between x(t) and the channel h(t). For the transmitted signal presented in FIG. 4, T_eff=5 samples. Denote by λ the average number of arrivals within T_eff not including the first arrival.

FIG. 5 and FIG. 6 present the simulation root mean square error (RMSE) results of EML, GE, and of SCC [25] as a function of λ for Ray SNR (RSNR) of 30 dB and 45 dB, respectively, where RSNR=σ_α²/σ_v². The signal used in these simulations is shown in FIG. 4. In practical UWB systems relatively high RSNR of 45 dB can be attained by transmitting successive pulses and averaging the received signals during the coherence time of the channel. The SCC was chosen for comparison since it is a practical solution that attempts to approximate the ML solution. All methods attain the same RMSE at λ=0. The RMSE of SCC increases with λ until it reaches an RMSE error ceiling of about 0.35 T_eff at λ>2. The RMSE of both GE and EML increases as λ increases until λ=0.25 and for λ>0.25 their RMSE reduces until it reaches an error floor of 0.15 T_eff at λ=4 and 0.06 T_eff at λ=2 for RSNR=30 dB and 45 dB, respectively. The RMSE of GE consolidates with that of EML for λ>2. The GE has lower RMSE than SCC at RSNR of 45 dB for λ>0 and at RSNR of 30 dB for λ>1. These results show that approximating the received signal as a Gaussian vector is effective when on average there are more than two arrivals (excluding the first arrival) within the effective pulse duration (T_eff), or, on the other extreme, a single arrival within the effective pulse duration. For such cases the GE attains the performance of EML and there is a substantial performance advantage with respect to SCC. For lower multipath density the GE outperforms other leading reference methods such as the SCC at high RSNR. In order to apply the above conclusions to real life channels, the ray density of practical channel models was tested. For transmission bandwidth of up to 1 Ghz the value of λ is greater than two for most channel models proposed by the IEEE 802.15.4a study group [3]-[4] as detailed in the sequel.

The behavior of the RMSE as λ increases can be explained as follows. The result of two contradicting effects is observed. The first effect is increasing the probability of decayed first arrival due to destructive superposition with multipath reflections and the second is reducing the distance between the first and following arrivals. The first effect increases the miss-detection probability of the first arrival while the second effect reduces the error when arrivals following the first arrival are detected instead. FIG. 5 and FIG. 6 show that the performance of EML and GE are dominated by the first effect at low values of λ and for λ>0.25 the second effect is the dominant.

The performance of the estimator of the present embodiments has been evaluated for three different channel models. The first tested channel, denoted CH0, is a white complex Gaussian process with exponential decaying variance, thus, it perfectly matches the channel model for which the GE was derived. The multipath duration was 40 nsec and the exponential decay rate is given by σ_h²(t)=e^{−(t−τ)/ξ} where ξ=10 nsec. The other two channels, denoted CH2 and CH3, are channel model 3 (indoor office environment) and 1 (residential indoor environment), respectively from the UWB channel models proposed by the IEEE 802.15.4a study group [4]. The impulse response of the transmitted signal used for testing all three channels is presented in FIG. 7. The total pulse duration is T_D=12 ns, the effective duration is T_eff=2 ns and the side lobe levels absolute peak is 18 dB below the main lobe peak. The main difference between CH1 and CH2 is that the multipath arrivals of CH1 are dense (λ≈6) and those of CH2 are sparse (λ≈0.3) with respect to the T_eff.

The performance of the GE and AGE are compared to SCC [25] and to the Jump Back Search Forward (JBSF) method [14]. In the JBSF the pre-sampling filter was chosen as B(ƒ)=S*(ƒ), which is a matched filter to the transmitted signal and then the TOA was estimated by a search forward from the beginning of the estimation window for the first sample energy that exceeds a threshold level. For all other tested methods B(ƒ) was an anti-aliasing low pass filter. The SCC and JBSF are chosen as reference methods since they represent two ends of solutions considered in the literature. The SCC is a high-end method and the JBSF in contrast, is a widely used low-end method. For all methods a coarse sliding window energy search was conducted as described above in the section entitled “Derivation of Estimation Algorithms.” Note that this algorithm provided the same end results for the various algorithms as genie aided window placement.

The performance evaluation criterion was the RMSE as a function of the SNR defined as SNR=(Σ_n=∞^∞|q_n|²)/σ_v². The algorithms were tested at various SNR levels with Monte Carlo computer simulations where for each channel model 5000 channel realizations were generated. The multipath channel was generated with a bandwidth of 64 W, such that the minimum time between multipath rays was 31.25 psec and the τ search resolution was the same. For SCC and JBSF achieving such resolution requires increasing sampling rate or interpolation to the appropriate sampling rate. For GE and AGE this is not necessary since the metric can be calculated in any resolution as long as the sampling rate is above Nyquist minimum. The sampling rate for GE and AGE was 1 GHz and we have verified that indeed increasing it does not affect the results.

For fair comparison between the estimation methods, each had its own parameters optimized to yield minimum RMSE per channel model (not per each specific channel realization) and SNR value: the threshold level for SCC and JBSF and the ratio between the first arrival variance and noise variance σ_h²(τ)/σ_v² used for GE and AGE. The GE and AGE showed to be more robust against mismatch in the signal and noise variance level compared to SCC and JBSF (data not shown).

The RMSE results for CH0, CH1 and CH2 are presented in FIGS. 8-11. The results for all three channels show that there is an insignificant performance loss between AGE and GE when the TOA estimation is done with a single channel realization i.e. K=1, and a slight degradation at low SNR when K=5. The RMSE of both methods reduces with SNR and by increasing the number of independent realizations.

When the number of realizations was increased from one to five the mean square error (MSE) reduced by a factor larger than 10. For CH0 and CH1 the performance degradation was insignificant when a constant multipath power delay profile was used in calculating (10) for GE and AGE instead of the exact exponential decaying profile. However for CH2 both methods showed a performance gain when the exact power delay profile was used. In fact, the only significant information in the exact power delay profile for CH2 is the differentiation between the line of sight (t=τ) and the rest (t>τ).

The same performance gain was obtained by approximating the exact power delay profile as having only two different levels: one for t=τ and the other for t>τ. The results of SCC reach an RMSE floor of about 1 nsec which does not reduce with SNR, thus as the SNR increases, the RMSE gap between SCC and AGE increases. The RMSE results of AGE and GE are lower than SCC for all three channel models and all SNR ranges except for the low SNR range on CH2 where the SCC has slightly lower RMSE results.

The RMSE of the JBSF method is the highest from all the tested methods at low and medium SNR levels. At high SNR the JBSF has lower RMSE than SCC but still significantly higher than AGE. The developed RMSE analytical expression given by the root of (31) showed to be a loose lower bound for a single channel realization but for five independent realizations (K=5) the lower bound was tight.

The results presented in FIGS. 8-11 show that the performance of the developed low complexity AGE is robust to various multipath channels (both sparse and dense) and has superior performance to JBSF and to SCC which is widely referred to in literature as a potentially high accurate solution which approximates the ML estimation of all channel coefficients and delays. Unlike SCC, the RMSE of AGE reduces monotonically with SNR. The SCC reaches an error floor which is most likely explained by the imprecise successive cancelation of the multipath components that results in a residual multipath signal regardless of the SNR level. The effect of this residual signal on the TOA estimation is similar to the effect of noise. The AGE does not need to know the channel power delay profile accurately. Most of the performance gain is achieved by approximating the channel power delay profile as constant and thus the detector needs to estimate only the ratio between the signal energy and the noise variance.

Using diversity techniques to obtain only a small number of independent channel realizations yields a significant reduction of the MSE. For five independent realizations the MSE was reduced by a factor larger than 10 while averaging the TOA estimation of each individual realization gives a reduction factor of only 5. Hence, by applying diversity techniques, the AGE can achieve reliable estimation even at very low SNR.

The developed analytical RMSE analysis showed to lower bound the simulation results and the bound becomes tight as the number independent realizations increases. For five independent realizations the analysis was accurate. For the Gaussian channel it means that five independent realizations are sufficient to reach the CLRB. It is known that the ML estimator reaches the CRLB as the number of measurements increases. In the present case, the measurements are increased by increasing the number of realizations.

Appendix A

Exact Maximum Likelihood Estimator

In this Appendix, the exact maximum likelihood (EML) TOA estimator for a typical channel model with Poisson arrival times and complex Gaussian distribution coefficients [1] is derived following the notations given in Example 1.

Let {tilde over (τ)} be a vector containing the μ arrival times within the observation window that follow the arrival time of the first multipath component, τ, where 0≦μ≦N_μ and N_μ is the maximal number of arrivals within the observation window. Let Λ(μ) be the set of all the different vectors {tilde over (μ)} of length μ, then the EML estimation of the first arrival time τ based on the observation vector y is given by

$\begin{array}{cc}{\hat{\tau }}_{\mathrm{EML}}=\arg \underset{\tau }{\max }f\left(y|\tau \right)=\arg \underset{\tau }{\max }\sum _{\mu =1}^{N_{\mu }}\sum _{\tilde{\tau }\in \Lambda \left(\mu \right)}f\left(y|\tau ,\tilde{\tau }\right)f\left(\tilde{\tau }|\tau \right)& \left(\mathrm{EQ}.32\right)\end{array}$

The distribution of y given all the arrival times (τ and {tilde over (τ)}) is a Gaussian multivariable distribution thus

$\begin{array}{cc}f\left(y|\tau ,\tilde{\tau }\right)=\frac{1}{{\pi }^LC\left(\tau ,\tilde{\tau }\right)}{e}^{-y^H{C}^{-1}\left(\tau ,\tilde{\tau }\right)y}& \left(\mathrm{EQ}.33\right)\end{array}$
where L is the length of y and C(τ,{tilde over (τ)})=E{yy^H|τ,{tilde over (τ)}}. For practical implementation we suggest the following approximated model. The arrival times are treated as discrete over a finite grid with time spacing T_ε. The arrival times in {tilde over (τ)} have Poisson distribution hence the probability of an arrival within T_ε is given by P_τ=λe^−λTε where λ is the average arrival rate within T_ε which is assumed known. Let T_L be the observation vector duration time, then the maximal number of arrivals within T_L, given τ, is M_μ(τ)=(T_L−τ)/T_ε and thus
ƒ({tilde over (τ)}|τ)=P_τ^μ(1−P_τ)^Mμ(τ)−μ  (EQ. 34)

By substituting (33) and (34) into (32) we arrive at

$\begin{array}{cc}{\hat{\tau }}_{\mathrm{EML}}=\arg \underset{\tau }{\max }\sum _{\mu =1}^{M_{\mu }\left(\tau \right)}\sum _{\tilde{\tau }\in \Lambda \left(\mu \right)}{P_{\tau }^{\mu }\left(1-P_{\tau }\right)}^{M_{\mu }\left(\tau \right)-\mu }\frac{1}{{\pi }^LC\left(\tau ,\tilde{\tau }\right)}{e}^{-y^H{C}^{-1}\left(\tau ,\tilde{\tau }\right)y}& \left(\mathrm{EQ}.35\right)\end{array}$

The total number of summations in (35) is 2^Mμ(τ). For accurate TOA estimation of a wideband signal in a typical indoor environment this value is extremely high and hence the EML TOA estimator is too complex for practical implementation.

Appendix B

Autocorrelation of y

From (3) and (4) of Example 1 the autocorrelation elements of y, given τ, can be written as

$\begin{array}{cc}E\left\{y_i{y}_j^{*}|\tau \right\}=E\left\{y\left({\mathrm{iT}}_s\right)y^{*}\left({\mathrm{jT}}_s\right)|\tau \right\}=E\left\{\sum _m\sum _q{\alpha }_m{\alpha }_q^{*}x\left({\mathrm{iT}}_s-{\tau }_m\right)x\left({\mathrm{jT}}_s-{\tau }_q\right)|\tau \right\}+E\left\{v_i{v}_j^{*}\right)& \left(\mathrm{EQ}.36\right)\end{array}$

Since α_m and α_q are statistically independent for m≠q and so are v_i and v_j for i≠j then we can write (36) as

$\begin{array}{cc}E\left\{y_i{y}_j^{*}|\tau \right\}=E_{{\tau }_1,\dots {\tau }_{M-1}}\left\{\sum _{m=0}^{M-1}\zeta \left(m\right)|\tau \right\}+{\sigma }_v^2{\delta }_K\left(i-j\right)& \left(\mathrm{EQ}.37\right)\end{array}$
where ξ(m)=σ_α²(m)x (iT_s−τ_m)x (jT_s−τ_m), M is the number of multipath components, σ_α²(m)=E_αm{|α_m|²} and δ_K(n) is the Kronecker delta. The multipath components variances σ_α²(m) decay as m increases and hence the value of M is set to be large enough so that σ_α²(M) approaches zero. Let ƒ_τm(τ_m|τ) be the probability density of the m-th multipath component arrival time (m>1), given τ, thus,

$E_{{\tau }_1,\dots {\tau }_{M-1}}\left\{\sum _{m=0}^{M-1}\zeta \left(m\right)\right\}={\int }_{{\tau }_1=\tau }^{\infty }\dots {\int }_{{\tau }_{M-1}}^{\infty }=\tau \sum _{m=0}^{M-1}\zeta \left(m\right)f\left({\tau }_1,\dots ,{\tau }_{M-1}|\tau \right)d{\tau }_1\dots d{\tau }_{M-1}={\sigma }_{\alpha }^2\left(0\right)x\left(i{T}_s-\tau \right)x\left(j{T}_s-\tau \right)+\sum _{m=1}^{M-1}{\sigma }_{\alpha }^2\left(m\right){\int }_{{\tau }_m=\tau }^{\infty }x\left(i{T}_s-{\tau }_m\right)x\left(j{T}_s-{\tau }_m\right)f_{{\tau }_m}\left({\tau }_m|\tau \right)d{\tau }_m$

By changing the integration variable from τ_m to t and replacing the order of the summation and the integral we can write the last expression in (38) as

$\begin{array}{cc}\sum _{m=1}^{M-1}{\sigma }_{\alpha }^2\left(m\right){\int }_{{\tau }_m=\tau }^{\infty }x\left(i{T}_s-{\tau }_m\right)x\left(j{T}_s-{\tau }_m\right)f_{{\tau }_m}\left({\tau }_m|\tau \right)d{\tau }_m={\int }_{t=\tau }^{\infty }\sum _{m=1}^{M-1}{\sigma }_{\alpha }^2\left(m\right)f_{{\tau }_m}\left(t|\tau \right)x\left(i{T}_s-t\right)x\left(j{T}_s-t\right)dt& \left(\mathrm{EQ}.39\right)\end{array}$

By substituting (39) into (38) and then the result into (37) we finally arrive at

$\begin{array}{cc}E\left\{y_i{y}_j^{*}|\tau \right\}={\int }_{t=\tau }^{\infty }{\sigma }_h^2\left(t\right)x\left(i{T}_s-t\right)x\left(j{T}_s-t\right)dt+{\sigma }_v^2{\delta }_K\left(i-j\right)\mathrm{where}& \left(\mathrm{EQ}.40\right)\\ {\sigma }_h^2\left(t\right)={\sigma }_{\alpha }^2\left(0\right)\delta \left(t-\tau \right)+\sum _{m=1}^{M-1}{\sigma }_{\alpha }^2\left(m\right)f_{{\tau }_m}\left(t|\tau \right)& \left(\mathrm{EQ}.41\right)\end{array}$
is the channel power delay profile, i.e. the average power vs delay, as shown next. The power delay profile can be obtained by transmitting narrow pulses with duration going to zero through realizations of a random multipath channel with fixed time of arrival and then measuring the received signal variance over delay. Let p(t) be a transmitted narrow pulse with unit energy given by

$\begin{array}{cc}p\left(t\right)=\underset{\Delta \to \infty }{\lim }\left\{\begin{array}{cc}0& t>\Delta \\ 1/\sqrt{\Delta }& t<\Delta \end{array}\right\}.& \left(\mathrm{EQ}.42\right)\end{array}$

The channel output is given by Σ_mα_mp(t−τ_m) which is the result of convolution between p(t) and h(t). The channel output variance given τ is then

$\begin{array}{cc}E\left\{{\sum _m{\alpha }_{m}p\left(t-{\tau }_m\right)}^2|\tau \right\}=E_{{\tau }_1},\dots {\tau }_{M-1}\left\{\sum _m{\sigma }_{\alpha }^2\left(m\right)p^2\left(t-{\tau }_m\right)|\tau \right\}={\sigma }_{\alpha }^2\left(0\right)p^2\left(t-\tau \right)+\sum _{m=1}^{M-1}{\sigma }_{\alpha }^2\left(m\right){\int }_{{\tau }_m=\tau }^{\infty }{p}^2\left(t-{\tau }_m\right)f_{{\tau }_m}\left({\tau }_m|\tau \right)d{\tau }_m& \left(\mathrm{EQ}.43\right)\end{array}$

Substituting (42) into (43) yields

$\underset{\Delta \to 0}{\lim }{\sigma }_{\alpha }^2\left(0\right)p^2\left(t-\tau \right)+\sum _{m=1}^{M-1}{\sigma }_{\alpha }^2\left(m\right)\underset{\Delta \to 0}{\lim }\frac{1}{\Delta }{\int }_{{\tau }_m=t}^{t+\Delta }{f}_{{\tau }_m}\left({\tau }_m|\tau \right)d{\tau }_m={\sigma }_{\alpha }^2\left(0\right)\delta \left(t-\tau \right)+\underset{\Delta \to 0}{\lim }\sum _{m=1}^{M-1}{\sigma }_{\alpha }^2\left(m\right)\frac{1}{\Delta }{f}_{{\tau }_m}\left(t|\tau \right)\Delta ={\sigma }_{\alpha }^2\left(0\right)\delta \left(t-\tau \right)+\sum _{m=1}^{M-1}{\sigma }_{\alpha }^2\left(m\right)f_{{\tau }_m}\left(t|\tau \right)={\sigma }_h^2\left(t\right)$
thus, the channel power delay profile is equal to σ_h²(t).

References

• [1] A. M. Saleh and R. A. Valenzuela, “A Statistical Model for Indoor Multipath Propagation,” IEEE J. Selected. Areas Commun., vol. SAC-5, no. 2, pp. 128-137, February 1987.
• [2] A. F. Molisch, “Ultrawideband Propagation Channels—Theory, Measurements, and Modeling,” IEEE Trans. on Vehicular Technology, vol. 54, no. 5, pp. 1528-1545, September 2005.
• [3] A. F. Molisch et al., “A Comprehensive Standardized Model for Ultrawideband Propagation Channels,” IEEE Trans. on Antennas and Propagation, vol. 54, no. 11, pp. 3151-3166, November 2006
• [4] A. F. Molisch et al., “IEEE 802.15.4a channel model—Final report,” September 2004.
• [5] W. Chung and D. Ha, “An accurate ultra wideband (UWB) ranging for precision asset location,” in Proc. IEEE Conf. Ultrawideband Syst. Technol. (UWBST), Reston, Va., November 2003, pp. 389-393.
• [6] K. Yu and I. Oppermann, “Performance of UWB position estimation based on time-of-arrival measurements,” in Proc. IEEE Conf. Ultrawideband Syst. Technol. (UWBST), Kyoto, Japan, May 2004, pp. 400-404.
• [7] R. Fleming, C. Kushner, G. Roberts, and U. Nandiwada, “Rapid acquisition for ultra-wideband localizers,” in Proc. IEEE Conf. Ultrawideband Syst. Technol. (UWBST), Baltimore, Md., May 2002, pp. 245-249.
• [8] C. Falsi, D. Dardari, L. Mucchi, and M. Z. Win, “Time of arrival estimation for UWB localizers in realistic environments,” EURASIP J. Appl. Signal Process., vol. 2006, pp. 1-13, 2006.
• [9] V. Lottici, A. D'Andrea, and U. Mengali, “Channel estimation for ultrawideband communications,” IEEE J. Select. Areas Commun., vol. 20, no. 9, pp. 1638-1645, December 2002.
• [10] L. Stoica, A. Rabbachin, and I. Oppermann, “A low-complexity noncoherent IR-UWB transceiver architecture with TOA estimation,” IEEE Trans. Microwave Theory Tech., vol. 54, pp. 1637-1646, June 2006.
• [11] I. Guvenc, Z. Sahinoglu, and P. V. Orlik, “TOA Estimation for IR-UWB Systems With Different Transceiver Types,” IEEE Trans. Microwave Theory and Tech., vol. 54, no. 4, April 2006.
• [12] D. Dardari, C.-C. Chong, and M. Z. Win, “Analysis of threshold-based TOA estimator in UWB channels,” in Proc. Eur. Signal Process. Conf. (EUSIPCO), Florence, Italy, September 2006.
• [13] A. Chehri, P. Fortier and P. M. Tardif, “Time-of-Arrival Estimation For IR-UWB Systems Based on Two Step Energy Detection” Biennial Symposium on Commun., pp. 369-373, June 2008
• [14] I. Guvenc and Z. Sahinoglu, “Threshold-based TOA estimation for impulse radio UWB systems,” in Proc. IEEE in Conf. Ultra-Wideband (ICUWB), Zurich, Switzerland, 2005, pp. 420-425.
• [15] D. Dardari, A. Conti, U. Ferner, A. Giorgetti, M. Z. Win, “Ranging With Ultrawide Bandwidth Signals in Multipath Environments” Proc. of IEEE, Special Issue on UWB Technology and Emerging Applications, vol. 97, no. 2, pp. 404-426, February 2009.
• [16] I. Guvenc and Z. Sahinoglu, “Multiscale Energy Products for TOA Estimation in IR-UWB Systems,” in Proc. IEEE Global Telecommun Conf., St. Louis, Mo., vol. 1, pp. 209-213, December 2005.
• [17] Z. Tian and G. B. Giannakis, “A GLRT approach to Data-Aided Timing Acquisition in UWB Radios—Part I: Algorithms”, IEEE Trans. Wireless Commun., vol. 4, no. 6, pp. 2956-2967, November 2005.
• [18] X. Li, “Super-Resolution TOA Estimation With Diversity for Indoor Geolocation,” IEEE Trans. Wireless Commun., vol. 3, no. 1, January 2004.
• [19] R. Roy, A. Paulraj and T. Kailath, “Estimation of Signal Parameters via Rotational Invariance Techniques—ESPRIT”, IEEE Trans. Acoust., Speech, Signal Process., vol. 37, no. 7, pp. 984-995, July 1989.
• [20] R. J. J. Vidal, M. Najar, “High resolution time-of-arrival detection for wireless positioning systems”, In IEEE Vehicular Technology Conference, Vancouver, Canada, pp. 2283-2287, September 2002.
• [21] M. Navarro and M. Najar, “Joint Estimation of TOA and DOA in IR-UWB”, IEEE Workshop on Signal Processing Advances in Wireless Commun., Helsinki, Finland, pp. 17-20, June 2007.
• [22] K. Gedalyahu, R. Tur, Y. C. Eldar, “Multichannel Sampling of Pulse Streams at the Rate of Innovation”, IEEE Trans. on Signal Processing, vol. 59, no. 4, pp. 1491-1504, April 2011.
• [23] J. Kusuma, I. Maravic, and M. Vetterli. “Sampling with finite rate of innovation: Channel and timing estimation for UWB and GPS”. In IEEE International Conference on Commun (ICC'03), pp. 3540-3544, 2003.
• [24] M. Z. Win and R. A. Scholtz, “Characterization of ultra-wide bandwidth wireless indoor communications channels: A communication theoretic view,” IEEE J. Sel. Areas Commun., vol. 20, pp. 1613-1627, December 2002.
• [25] J. Y. Lee and R. A. Scholtz, “Ranging in a dense multipath environment using an UWB radio link,” IEEE J. Select. Areas Commun., vol. 20, no. 9, pp. 1677-1683, December 2002.
• [26] C. Mazzucco, U. Spagnolini, and G. Mulas, “A ranging technique for UWB indoor channel based on power delay profile analysis,” in Proc. IEEE Vehic. Technol. Conf. (VTC), Los Angeles, Calif., September 2004, pp. 2595-2599.
• [27] Y. Qi, H. Kobayashi, and H. Suda. “On time-of-arrival positioning in a multipath environment,” IEEE Trans. Veh. Technol., vol. 55, no. 5, pp. 1516-1526, September 2006.
• [28] R. C. Mittelhammer, Mathematical Statistics For Economics And Business, Chapter 5.5, Springer-Verlag New York Inc., 1996.
• [29] D. S. Bernstein, Matrix Mathematics Theory Facts and Formulas, Second edition, Princeton University Press, 2009.
• [30] A. J. Weiss and B. F. Friedlander, “Performance Analysis of Spatial Smoothing with Interpolated Arrays,” IEEE Tans. Signal Processing, vol. 41, no. 5, May 1993.
• [31] S. M. Kay, Fundamentals of Statistical Signal Processing: Detection Theory. Upper Saddle River, N.J.: Prentice-Hall, 1998.

Example 2

Direct Position Estimate

The present Example describes a direct position estimate procedure according to some embodiments of the present invention. The data presented below demonstrate that the direct position estimator of the present embodiments has superior performance compared to the conventional indirect method at low SNR over various multipath channel models and the gap narrows as the SNR increases. A closed form analytical expression that tightly approximated the DPE position estimation root mean square error is also presented.

Introduction

Position estimation methods can be partitioned into two important classes: direct and indirect position estimation. In the direct position estimation (DPE) the target position is estimated directly from the received signals (the raw measurements), preferably from all reference stations. In other words the received signals are transferred to the server for joint processing. In indirect position estimation (IPE) a few parameters are first estimated from each received signal, locally for each received signal, and then the position is estimated based on those parameters. Indirect positioning requires passing less data to the location estimation processor. On the other hand, IPE does not utilize all available information in the location estimation. Various solutions for IPE have been developed [7]-[11]. The common approach for the position estimation in the second step is to search for the position that attains the least squares error between the measured distances (from the first step) and the distances corresponding to each position hypothesis. In case a quality measurement is given per each distance measurement then the position can be estimated by weighted least squares [1].

System Model

Consider a wireless network where there are K reference stations at known positions and a target at unknown position. The 3D position of the k-th reference station and of the target are denoted by p_ref^k=[x^k y^k z^k]^T, k=0,1, . . . , K−1, and p=[x y z]^T, respectively. It is assumed that the target is the transmitting side and that the reference stations are the receivers. The system model described in the sequel also applies to the case that the reference stations are the transmitters side and the target is the receiver side since the problem is equivalent for both cases. The target transmits a single pulse s(t) which is assumed known to the reference stations. It is assumed that the reference stations and the target are synchronized so that the pulse transmission time is known to the reference stations. On the way from the target to the k-th reference station the pulse travels through channel h^k(t) which is assumed to be a time invariant multipath channel with impulse response given by

$\begin{array}{cc}{h}^k\left(t\right)=\sum _{m=0}^{M_k-1}{\alpha }_m^k\delta \left(t-{\tau }_m^k\right)& \left(\mathrm{EQ}.1\right)\end{array}$
where k denotes the reference station index, M_k is the number of multipath components in h^k(t), δ(t) is the Dirac delta function, α_m^k and τ_m^k are the complex coefficient and the delay of the m-th multipath component, respectively. All α_m^k and τ_m^k for any k or m are statistically independent. In other words, the multipath reflections are independent within a channel and also the channels of different reference stations are statistically independent. Throughout this Example the notation of the first arrival at the k-th reference station, τ₀^k, is replaced by τ^k for simplifying the exhibition. It is assumed that the first arrival corresponds to the direct path from the transmitter to the receiver. It can be either a line-of-sight component or a direct path through an obstruction. The received signal at the k-th reference station is given by

$\begin{array}{cc}{r}^k\left(t\right)=\sum _{m=0}^{M_k-1}{\alpha }_m^{k}s\left(t-{\tau }_m^k\right)+w^k\left(t\right)& \left(\mathrm{EQ}.2\right)\end{array}$
where w^k(t) is additive white complex Gaussian noise observed at the k-th reference station.

At the receiver, the signal r^k(t) is passed through a pre-sampling low pass filter (LPF) whose impulse response is denoted by b(t) and its bandwidth is denoted by W. It is assumed that W is larger than the bandwidth of s(t). The filter output y^k(t) is sampled at rate ƒ_s=1/T_s=2W. The resulting sampled sequence is denoted by y_n^k and can be written as
y_n^k=q_n^k+v_n^k  (EQ. 3)
where v_n^k are the noise samples obtained by filtering w^k(t) with b(t) and q_n^k are the samples of the received multipath signal. Let x(t) be the result of the convolution between s(t) and b(t) then q_n^k is obtained by convolving h^k(t) with x(t) and sampling the result at a rate of 1/T_s

$\begin{array}{cc}{q}_n^k={\int }_0^{{\mathrm{nT}}_s}{h}^k\left(\lambda \right)x\left({\mathrm{nT}}_s-\lambda \right)d\lambda =\sum _{m=0}^{M_k-1}{\alpha }_m^{k}x\left({\mathrm{nT}}_s-{\tau }_m^k\right).& \left(\mathrm{EQ}.4\right)\end{array}$

The noise samples v_n^k are assumed to be independent complex circular Gaussian random variables with zero mean and variance σ_v². It is assumed that x(t) is causal.

Let y^k=[y₀^k y₁^k . . . y_Lk−1^k]^T be a vector of length L_k containing the k-th reference station observation window samples which are given by y_n^k=y^k(nT_s). Without loss of generality, it is assumed that the beginning of the observation window is at time t=0 which is also the transmission time of s(t). It is assumed that the observation window contains the multipath signal or at least the initial part of it and does not contain noise samples following the multipath signal. A procedure suitable for searching for the observation window is presented in Example 1, above. It is noted that the time-frame definition above is used for convenient mathematical description and is not the physical time frame. The receivers know the time difference between the actual transmission and the observation window start and can compensate for it. The observation vectors from all reference stations are collected by a server which uses them to estimate the position of the target. It is assumed that the signal propagation speed is the speed of light denoted by c.

Direct Position Estimation

In this section DPE is derived for the system model presented above. For given target and reference station positions, the TOA of the first multipath component at the k-th reference station can be expressed as

$\begin{array}{cc}{\tau }^k=\frac{1}{c}\sqrt{{\left(p_{\mathrm{ref}}^k-p\right)}^T\left(p_{\mathrm{ref}}^k-p\right)}.& \left(\mathrm{EQ}.5\right)\end{array}$

There are all together K equations for τ⁰, τ¹, . . . , τ^K−1 thus in principle if K≧3 then the target 3D position can be determined by the time of arrival measurements. The ML DPE is given by

$\begin{array}{cc}\hat{p}=\arg \underset{p}{\max }f\left(y^0,y^1,\dots ,y^{K-1}|p\right).& \left(\mathrm{EQ}.6\right)\end{array}$

Since each position hypothesis p determines the time of arrivals at the K reference stations, according to (5), then the ML estimator in (6) can be written as

$\begin{array}{cc}\hat{p}=\arg \underset{p}{\max }f\left(y^0,y^1,\dots ,y^{K-1}|{\tau }^0,{\tau }^1,\dots ,{\tau }^{K-1}\right)& \left(\mathrm{EQ}.7\right)\end{array}$
where τ^k is a function of p. The multipath channel realization of different reference stations are assumed statistically independent hence the observation vectors y⁰, y¹, . . . , y^K−1 are independent and thus

$\begin{array}{cc}\hat{p}=\arg \underset{p}{\max }\prod _{k=0}^{K-1}f\left(y^k|{\tau }^k\right)& \left(\mathrm{EQ}.8\right)\end{array}$

An approximation to the likelihood function ƒ(y^k|τ^k) was derived in Example 1 above by approximating the received signal as Gaussian. The resulting likelihood function is practical for implementation and approaches the exact likelihood function when the multipath is dense or when there is a single channel component. In the intermediate cases it showed to be a good approximation. In the present Example, the detailed DPE is derived following the approach taken in Example 1.

Each observation vector y^k, given τ^k, is partitioned into two segments y=[(y₀^k)^T (y₁)^T]^T such that the first segment y₀^k contains only noise samples, and the second segment y₁^k contains the received multipath first cluster plus noise samples. Note that τ^k, the first arrival time with respect to the beginning of the observation window, is a real number. Let N_τ^k be the integer index of the first sample within the observation vector y^k that contains the multipath signal, thus, (N_τ^k−1)T_s<τ^k≦N_τ^kT_s. The segment y₁^k is assumed to have a Gaussian distribution. This assumption is strictly true in two cases. The first is when h^k(t)=α₀^kδ(t−τ^k) and α₀^k is a complex Gaussian random variable. The other case is when the multipath density is so high that each sample in y₁^k is a sum of infinite number of random variables. However, it has been shown that even a small number of multipath arrivals within the transmitted pulse duration is enough to justify the Gaussian assumption and that the TOA estimator based on this assumption gives good performance for various multipath models both sparse and dense [14].

For each observation vector y^k, the segments prior and following τ^k are statistically independent hence
ƒ(y^k|τ^k)=ƒ(y₀^k|τ^k)ƒ(y₁^k|τ^k).  (EQ. 9)

The noise samples in y₀^k are independent complex Gaussian random variable with zero mean, variance σ_v² and independent real and imaginary components, thus

$\begin{array}{cc}f\left(y_0^k|{\tau }^k\right)=\prod _{n=0}^{N_{\tau }^k-1}\frac{1}{\pi {\sigma }_v^2}{e}^{-\frac{{y_n^k}^2}{{\sigma }_v^2}}.& \left(\mathrm{EQ}.10\right)\end{array}$

Let L_k be the length of y^k. The L_k−N_τ^k samples in y₁^k are assumed to have Gaussian distribution hence their distribution is given by

$\begin{array}{cc}f\left(y_1^k|{\tau }^k\right)=\frac{1}{{\pi }^{L_i-N_{\tau }^k}R_1^k\left({\tau }^k\right)}{e}^{-{\left(y_1^k\right)}^H{\left(R_1^k\left({\tau }^k\right)\right)}^{-1}{y}_1^k}& \left(\mathrm{EQ}.11\right)\end{array}$
where R₁^k(τ^k) is the autocorrelation matrix of y₁^k which is dependent on the value of τ^k. Denote by ρ_i,j^k(τ^k) the element in the i-th row and j-th column of R₁^k(τ^k) which is given by
ρ_i,j^k(τ^k)=E{y_Nτk+i^k(y_Nτk+j^k)*|τ^k}=∫_t=τ^∞σ_h²(t)x((N_τ^k+i)T_s−t)x((N_τ^k+j)T_s−t)dt+σ_v²δ_K(i−j)  (EQ. 12)
where δ_K(n) is the Kronecker delta and σ_h²(t) is the channel power delay profile, i.e. the average power vs delay [14].

By substituting (10) and (11) into (9) we arrive at the following Gaussian model TOA likelihood function

$\begin{array}{cc}f\left(y^k|{\tau }^k\right)=\left[\prod _{n=0}^{N_{\tau }^k-1}\frac{1}{\pi {\sigma }_v^2}{e}^{-\frac{{y_n^k}^2}{{\sigma }_v^2}}\right]\frac{1}{{\pi }^NR_1^k\left({\tau }^k\right)}{e}^{-{\left(y_1^k\right)}^H{\left(R_1^k\left({\tau }^k\right)\right)}^{-1}{y}_1^k}.& \left(\mathrm{EQ}.13\right)\end{array}$

By substituting (13) into (8), applying ln to the result, and eliminating insignificant terms we arrive at the following DPE

$\begin{array}{cc}{\hat{p}}_{\mathrm{DPE}}=\arg \underset{p}{\min }\left\{\sum _{k=0}^{K-1}{N}_{\tau }^k\ln \left({\sigma }_v^2\right)+\frac{1}{{\sigma }_v^2}·y_0^k{·}^2+\ln \left(R_1^k\left({\tau }^k\right)\right)+{\left(y_1^k\right)}^H{\left(R_1^k\left({\tau }^k\right)\right)}^{-1}{y}_1^k\right\}.& \left(\mathrm{EQ}.14\right)\end{array}$

Recall that τ^k is a function of p. Note that further complexity reduction can be obtained by applying an approximation to (11) [14]. The search for the position that minimizes (14) can be done in two stages. In a first stage (14) is calculated over a discrete grid with relatively wide spacing and the U position hypothesis that attain minimum metric are set as initial estimations. In a second stage a fine search is conducted in the vicinity of each initial estimation to find a single position that minimizes (14). The fine search can be done by calculating (14) over a fine grid or by applying numerical minimization methods such as steepest descent or Gaussian-Newton [15].

Performance Analysis

The performance of the Gaussian estimator given in (14) will now be evaluated. An expression which approximates the position estimation mean square error (MSE) under the assumption that the error is small has been derived. This condition holds as the SNR increases and also when the number of measurements grow, which in the present Example can be realized by increasing the number of reference stations. The results presented below show that with four reference stations the small error conditions holds for a wide SNR range. The analysis is applicable for any multipath density and even for low density where the Gaussian approximation is inaccurate.

A cost function (a function to be minimized) in the face of (14) can be written as:

$\begin{array}{cc}C\left(p\right)=\sum _{k=0}^{K-1}\left[\ln \left(R^k\left({\tau }^k\right)\right)+{\left(y^k\right)}^H{\left(R^k\left({\tau }^k\right)\right)}^{-1}{y}^k\right]& \left(\mathrm{EQ}.15\right)\end{array}$
where y^k is the k-th full observation vector (including the segment before and after the TOA), R^k(τ^k)=E{y^k(y^k)^H} and the elements E{y_i^k(y_j^k)^H} are given in (14). For the analysis it is assumed, without loss of generality, that the observation vectors y^k, k=0,1, . . . K−1, all have the same length L. Using the relation (y^k)^H(R^k(τ^k))⁻¹ y^k=Tr{(R^k(τ^k))⁻¹ y^k(y^k)^H} (15) can be written as:

$\begin{array}{cc}C\left(p\right)=\sum _{k=0}^{K-1}\left[\ln \left(R^k\left({\tau }^k\right)\right)+\mathrm{Tr}\left\{{\left(R^k\left({\tau }^k\right)\right)}^{-1}{\hat{R}}^k\right\}\right]& \left(\mathrm{EQ}.16\right)\end{array}$
where Tr{ } denotes the trace of a matrix and {circumflex over (R)}^k=y^k(y^k)^H is the empirical autocorrelation estimation. Let p_t and {circumflex over (p)} be the true and estimated positions, respectively. Denote by vec(X) the vectorized version of some matrix X. Let v_t=└vec(R⁰(τ⁰))^T vec(R¹(τ¹))^T . . . vec(R^K−1(τ^K−1))^T┘ and {circumflex over (v)}=[vec({circumflex over (R)}⁰)^T vec({circumflex over (R)}¹)^T . . . vec({circumflex over (R)}^K−1)^T] be a vectorized representation of the true and empirical autocorrelation matrices, respectively. Let

$G\left(p,\hat{v}\right)=\frac{\partial }{\partial p}C\left(p\right).$
A first order Taylor expansion of G at p=p_t and {circumflex over (v)}=v_t gives
G(p,{circumflex over (v)})≅G(p_t,v_t)+Q[p−p_t]+F[{circumflex over (v)}−v_t]  (EQ. 17)
where

$Q={\left[\frac{\partial }{\partial p}G\right]}_{p=p_t,\hat{v}=v_t}\mathrm{and}F={\left[\frac{\partial }{\partial \hat{v}}G\right]}_{p=p_t,\hat{v}=v_t}.$
For a given realization {circumflex over (v)} the cost function C attains maximum at p={circumflex over (p)} and for {circumflex over (v)}=v_t the cost function is maximized for p=p_t hence G({circumflex over (p)},{circumflex over (v)})=G(p_t,v_t)=0. Substituting the later into (17) gives
Q[{circumflex over (p)}−p_t]≅−F[{circumflex over (v)}−v_t]  (EQ. 18)

Right multiplying (18) with the conjugate transpose of itself yields
Q[{circumflex over (p)}−p_t][{circumflex over (p)}−p_t]^HQ^H≅F[{circumflex over (v)}−v_t][{circumflex over (v)}−v_t]^HF^H  (EQ. 19)

Left and right multiplying (19) by Q⁻¹ and (Q⁻¹)^H respectively and then taking the expectation results in
E{[{circumflex over (p)}−p_t][{circumflex over (p)}−p_t]^H}≅Q⁻¹FE{[{circumflex over (v)}−v_t][{circumflex over (v)}−v_t]^H}F^H(Q⁻¹)^H   (EQ. 20)

Let Ω_i,j=E{[vec({circumflex over (R)}ⁱ)−vec(Rⁱ(τⁱ))][vec({circumflex over (R)}^j)−vec(R^j(τ^j))]^H} then the expectation on the right hand side of (20) can be written as

$\begin{array}{cc}E\left\{{\left[\hat{v}-v_t\right]\left[\hat{v}-v_t\right]}^H\right\}=\left[\begin{array}{cccc}{\Omega }_{0,0}& {\Omega }_{0,1}& \dots & {\Omega }_{0,K-1}\\ {\Omega }_{1,0}& {\Omega }_{1,1}& \dots & {\Omega }_{1,K-1}\\ \dots & \dots & \dots & \dots \\ {\Omega }_{K-1,0}& {\Omega }_{K-1,1}& \dots & {\Omega }_{K-1,K-1}\end{array}\right]& \left(\mathrm{EQ}.21\right)\end{array}$
since E{vec({circumflex over (R)}ⁱ)}=vec(Rⁱ(τⁱ)) and {circumflex over (R)}ⁱ,{circumflex over (R)}^j are statistically independent for i≠j then Ω_i,j=0 for i≠j. It has been shown [16] that Ω_i,j=Rⁱ(τⁱ){circle around (x)}Rⁱ(τⁱ) where {circle around (x)} is the Kronecker product.
Hence (20) can be written as:

$\begin{array}{cc}E\left\{{\left[\hat{p}-p_t\right]\left[\hat{p}-p_t\right]}^H\right\}\cong Q^{-1}F\left[\begin{array}{cccc}{R}^0\left({\tau }^0\right)\otimes R^0\left({\tau }^0\right)& 0& \dots & 0\\ 0& R^1\left({\tau }^1\right)\otimes R^1\left({\tau }^1\right)& \dots & 0\\ \dots & \dots & \dots & \dots \\ 0& 0& \dots & R^{K-1}\left({\tau }^{K-1}\right)\otimes R^{K-1}\left({\tau }^{K-1}\right)\end{array}\right]{F^H\left(Q^{-1}\right)}^H& \left(\mathrm{EQ}.22\right)\end{array}$

Q is calculated as follows. Let τ=[τ⁰ τ¹ . . . τ^K−1]^T then G is given by

$\begin{array}{cc}G=\frac{\partial C\left(p\right)}{\partial p}={\left(\frac{d\tau }{dp}\right)}^T\frac{\partial C\left(p\right)}{\partial \tau }={{\left(\frac{d\tau }{dp}\right)}^T\left[\begin{array}{cccc}\frac{\partial C\left(p\right)}{\partial {\tau }^0}& \frac{\partial C\left(p\right)}{\partial {\tau }^1}& \dots & \frac{\partial C\left(p\right)}{\partial {\tau }^{K-1}}\end{array}\right]}^T& \left(\mathrm{EQ}.23\right)\end{array}$
and the derivative of G with respect to p is given by

$\begin{array}{cc}\frac{\partial }{\partial p}G=\frac{\partial }{\partial p}\left({\left(\frac{d\tau }{dp}\right)}^T\frac{\partial }{\partial \tau }C\left(p\right)\right)={\left(\frac{d\tau }{dp}\right)}^T\left[\frac{{\partial }^2}{{\partial }^2\tau }C\left(p\right)\right]\left(\frac{d\tau }{dp}\right)={\left(\frac{d\tau }{dp}\right)}^T\left[\begin{array}{ccc}\frac{{\partial }^2}{{\partial }^2{\tau }^0}C\left(p\right)& \frac{\partial \partial }{\partial {\tau }^0\partial {\tau }^1}C\left(p\right)& \frac{\partial \partial }{\partial {\tau }^0\partial {\tau }^{K-1}}C\left(p\right)\\ \frac{\partial \partial }{\partial {\tau }^1\partial {\tau }^0}C\left(p\right)& \dots & \dots \\ \frac{\partial \partial }{\partial {\tau }^{K-1}\partial {\tau }^0}C\left(p\right)& \dots & \frac{{\partial }^2}{{\partial }^2{\tau }^{K-1}}C\left(p\right)\end{array}\right]\left(\frac{d\tau }{dp}\right)& \left(\mathrm{EQ}.24\right)\end{array}$
Since

$\frac{\partial \partial }{\partial {\tau }^j\partial {\tau }^i}C\left(p\right)=0$
for i≠j then from (24) Q can be written as

$\begin{array}{cc}Q={\left[\frac{\partial }{\partial p}G\right]}_{p=p_t,\hat{v}=v_t}=J^T\left[\begin{array}{ccc}{\gamma }_0& 0& 0\\ 0& \dots & 0\\ 0& 0& {\gamma }_{K-1}\end{array}\right]J\mathrm{where}{\gamma }_i={\left[\frac{{\partial }^2}{{\partial }^2{\tau }^i}C\left(p\right)\right]}_{p=p_t,\hat{v}=v_t}\mathrm{and}& \left(\mathrm{EQ}.25\right)\\ J={\left[\frac{d\tau }{dp}\right]}_{p=p_t}={\left[\begin{array}{ccc}\frac{\partial {\tau }^0}{\partial x}& \frac{\partial {\tau }^0}{\partial y}& \frac{\partial {\tau }^0}{\partial z}\\ \frac{\partial {\tau }^1}{\partial x}& \frac{\partial {\tau }^1}{\partial y}& \frac{\partial {\tau }^1}{\partial z}\\ ⋮& ⋮& ⋮\\ \frac{\partial {\tau }^{K-1}}{\partial x}& \frac{\partial {\tau }^{K-1}}{\partial y}& \frac{\partial {\tau }^{K-1}}{\partial z}\end{array}\right]}_{p=p_t}.& \left(\mathrm{EQ}.26\right)\end{array}$

Note that by substituting (5) into (26) it can be shown that the i-th row of J is the normalized vector connecting the i-th reference station position and the target position.

Next we calculate F which is given by

$\begin{array}{cc}F={\left[\frac{\partial }{\partial v}G\right]}_{p=p_t,\hat{v}=v_t}={J^T\left[\begin{array}{c}\frac{\partial \partial }{\partial \hat{v}\partial {\tau }^0}C\left(p\right)\\ \frac{\partial \partial }{\partial \hat{v}\partial {\tau }^1}C\left(p\right)\\ ⋮\\ \frac{\partial \partial }{\partial \hat{v}\partial {\tau }^{K-1}}C\left(p\right)\end{array}\right]}_{p=p_t,\hat{v}=v_t}& \left(\mathrm{EQ}.27\right)\end{array}$
where

$\frac{\partial \partial }{\partial \hat{v}\partial {\tau }^i}C\left(p\right)$
is a row vector given by

$\begin{array}{cc}\frac{\partial \partial }{\partial \hat{v}\partial {\tau }^i}C\left(p\right)=[\frac{\partial \partial }{\partial \mathrm{vec}\left({\hat{R}}^0\right)\partial {\tau }^i}C\left(p\right)\hspace{1em}\frac{\partial \partial }{\partial \mathrm{vec}\left({\hat{R}}^1\right)\partial {\tau }^i}C\left(p\right)\dots \frac{\partial \partial }{\partial \mathrm{vec}\left({\hat{R}}^{K-1}\right)\partial {\tau }^i}C\left(p\right)]& \left(\mathrm{EQ}.28\right)\end{array}$
and

$\frac{\partial \partial }{\partial \mathrm{vec}\left({\hat{R}}^i\right)\partial {\tau }^j}C\left(p\right)$
is a row vector with the length of vec({circumflex over (R)}ⁱ) which is L². It can be shown that

$\frac{\partial \partial }{\partial \mathrm{vec}\left({\hat{R}}^i\right)\partial {\tau }^j}C\left(p\right)=\left[\begin{array}{cccc}0& 0& \dots & 0\end{array}\right]$
for i≠j hence (27) can be written as

$\begin{array}{cc}F=J^T\left[\begin{array}{ccc}{\left[\frac{\partial \partial }{\partial \mathrm{vec}\left({\hat{R}}^0\right)\partial {\tau }^0}C\left(p\right)\right]}_{p=p_t,\hat{v}=v_t}& Z& Z\\ Z& \dots & Z\\ Z& \dots & {\left[\frac{\partial \partial }{\partial \mathrm{vec}\left({\hat{R}}^{K-1}\right)\partial {\tau }^{K-1}}C\left(p\right)\right]}_{p=p_t,\hat{v}=v_t}\end{array}\right]& \left(\mathrm{EQ}.29\right)\end{array}$
where Z is a row vector with zero elements and length L². (29) and (25) are substituted into (22) to provide

$\begin{array}{cc}E\left\{{\left[\hat{p}-p_t\right]\left[\hat{p}-p_t\right]}^H\right\}\cong {\left[J^T\left[\begin{array}{ccc}{\gamma }_0& 0& 0\\ 0& \dots & 0\\ 0& 0& {\gamma }_{K-1}\end{array}\right]J\right]}^{-1}{J}^T\left[\begin{array}{cccc}{\psi }_0& 0& \dots & 0\\ 0& {\psi }_1& \dots & 0\\ \dots & \dots & \dots & \dots \\ 0& 0& \dots & {\psi }_{K-1}\end{array}\right]{J\left[J^T\left[\begin{array}{ccc}{\gamma }_0& 0& 0\\ 0& \dots & 0\\ 0& 0& {\gamma }_{K-1}\end{array}\right]J\right]}^{-1}\mathrm{where}& \left(\mathrm{EQ}.30\right)\\ {\psi }_i={\left[\frac{\partial \partial }{\partial \mathrm{vec}\left({\hat{R}}^i\right)\partial {\tau }^i}C\left(p\right)\right]}_{p=p_t,\hat{v}=v_t}{R}^i\left({\tau }^i\right)\otimes {R^i\left({\tau }^i\right)\left[\frac{\partial \partial }{\partial \mathrm{vec}\left({\hat{R}}^i\right)\partial {\tau }^i}C\left(p\right)\right]}_{p=p_t,\hat{v}=v_t}^T& \left(\mathrm{EQ}.31\right)\end{array}$

In Appendix A it is shown that

$\begin{array}{cc}{\psi }_i={\gamma }_i=\mathrm{Tr}\left\{{\left[{{\left(R^i\left({\tau }^i\right)\right)}^{-1}\left[\frac{\partial R^i\left(\tau \right)}{\partial \tau }\right]}_{\tau ={\tau }^i}\right]}^2\right\}.& \left(\mathrm{EQ}.32\right)\end{array}$
(32) is substituted into (30) to obtain the DPE covariance given by

$\begin{array}{cc}E\left\{{\left[\hat{p}-p_t\right]\left[\hat{p}-p_t\right]}^H\right\}\cong {\left[J^T\left[\begin{array}{ccc}{\psi }_0& 0& 0\\ 0& \dots & 0\\ 0& 0& {\psi }_{K-1}\end{array}\right]J\right]}^{-1}& \left(\mathrm{EQ}.33\right)\end{array}$
and the DPE MSE given by

$\begin{array}{cc}E\left\{{\hat{p}-p_t}^2\right\}\cong \mathrm{Tr}\left\{{\left[J^T\left[\begin{array}{ccc}{\psi }_0& 0& 0\\ 0& \dots & 0\\ 0& 0& {\psi }_{K-1}\end{array}\right]J\right]}^{-1}\right\}& \left(\mathrm{EQ}.34\right)\end{array}$

Under small error conditions the MSE of the ML TOA estimation that maximizes the TOA likelihood function in (13) is given by [14]

$\begin{array}{cc}E\left\{{\left({\hat{\tau }}^k-{\tau }^k\right)}^2\right\}\cong \frac{1}{{\psi }_k}& \left(\mathrm{EQ}.35\right)\end{array}$

Let {circumflex over (p)}_ie be the IPE. It has been shown [18] that under the assumption of small error in the TOA estimation, the covariance of the IPE error is given by
E{[{circumflex over (p)}_ie−p_t][{circumflex over (p)}_ie−p_t]^H}≅Tr{[J^TW¹J]⁻¹}  (EQ. 36)
where

$\begin{array}{cc}W=E\left\{\left[\begin{array}{c}\left({\hat{\tau }}^0-{\tau }^0\right)\\ \left({\hat{\tau }}^1-{\tau }^1\right)\\ ⋮\\ \left({\hat{\tau }}^{K-1}-{\tau }^{K-1}\right)\end{array}\right]\left[\left({\hat{\tau }}^0-{\tau }^0\right)\left({\hat{\tau }}^1-{\tau }^1\right):\left({\hat{\tau }}^{K-1}-{\tau }^{K-1}\right)\right]\right\}& \left(\mathrm{EQ}.37\right)\end{array}$

Since the received signals at different reference stations are statistically independent and the TOA error is assumed unbiased then E{({circumflex over (τ)}ⁱ−τⁱ)({circumflex over (τ)}^j−τ^j)}=0 for i≠j. Thus the IPE covariance can be written as

$\begin{array}{cc}E\left\{{\left[{\hat{p}}_{\mathrm{ie}}-p_t\right]\left[{\hat{p}}_{\mathrm{ie}}-p_t\right]}^H\right\}\cong \mathrm{Tr}\left\{{\left[J^T\left[\begin{array}{ccc}1/E\left\{{\left({\hat{\tau }}^0-{\tau }^0\right)}^2\right\}& 0& 0\\ 0& \dots & 0\\ 0& 0& 1/E\left\{{\left({\hat{\tau }}^{K-1}-{\tau }^{K-1}\right)}^2\right\}\end{array}\right]J\right]}^{-1}\right\}& \left(\mathrm{EQ}.38\right)\end{array}$

(33), which is the covariance of the DPE, is obtained by substituting (35) into (38).

The above proves that when the small error conditions holds for both DPE and IPE then there performance are equal.

It is noted that there is a difference between the small error conditions for which the covariance expressions of DPE and IPE given in (33) and (38) respectively are accurate. The SNR for which the small error is valid is far lower for the DPE than for the IPE. Thus, the performance of the DPE and IPE consolidate only at high SNR. It is also noted that even at low SNR the small error approximation is valid for the DPE and the analysis is expected to give a good approximation.

In case the TOA estimation MSE is equal to all reference stations then (34) can be written as
E{∥{circumflex over (p)}−p_t∥²}≅E{({circumflex over (τ)}^k−τ^k)²}×GDOP²  (EQ. 39)
where GDOP is the geometric dilution of precision factor given by GDOP=√{square root over (Tr{[J^TJ]⁻¹})}. The GDOP is the ratio between the root mean square errors of the location and the TOA estimations. It gives a relation between the position estimation error gain and the spatial geometry of the reference stations and target positions [18].
Results and Discussion

The performance of the developed DPE was tested with Monte Carlo simulation. The reference stations and target where placed in a 50×50×50 meters cubic volume and each received signal was obtained by filtering the transmitted pulse with a multipath realization taken from channel model 3 of the UWB channel models proposed by the IEEE 802.15.4a study group [6]. This channel models an indoor office environment with relatively dense multipath. The impulse response of the target transmitted signal is presented in FIG. 12. The total pulse duration is 12 ns. The width of its effective energetic duration is 2 ns, which means that the signal energy in this width contains most of the signal energy. The received signal was sampled at 1 Ghz. The performance criterion was the position error root mean square error RMSE vs SNR where RMSE=√{square root over (E{∥{circumflex over (p)}−p∥²})} and SNR=(Σ_n=−∞^∞|q_n|²)/σ_v². The RMSE of the DPE given in (14) was compared to the traditional least squares IPE. In the first step of the IPE the TOA of each reference station was estimated by maximizing the approximate TOA likelihood function given in (11). This TOA estimation method showed to give good results and outperform other conventional TOA estimation methods in multipath conditions [14]. In the second step the position was estimated by

$\begin{array}{cc}{\hat{p}}_{\mathrm{IPE}}=\arg \underset{p}{\min }\sum _{k=0}^{K-1}{\left({\hat{\tau }}^k-{\tau }^k\right)}^2& \left(\mathrm{EQ}.40\right)\end{array}$
where {circumflex over (τ)}^k are the estimated time of arrivals and τ^k are the time of arrivals corresponding to the position hypothesis. For both direct and indirect estimators the cost functions were minimized by first performing a coarse search for three initial estimations (per each method) that attain minimum cost function. The coarse search was over the entire cube dimension of 50×50×50 meter with 1 meter spacing. Then a fine grid search with 10 cm spacing was performed over cubes with dimension 1×1×1 meter centered at each initial estimation. We have verified that the coarse grid search spacing is small enough so that it does not degrade the estimation performance.

FIGS. 13 and 14 show RMSE vs SNR results obtained when the reference stations were located at the four upper corners of the cubic volume which are (0,0,50)m, (50,0,50)m, (0,50,50)m, (50,50,50)m and two different target locations (25,25,25)m and (25,25,47.5)m, respectively. The target at (25,25,25)m has GDOP=1.5 which is low compared to the other target position with GDOP=7.15.

As shown, for any SNR or method, the point with the higher GDOP has larger RMSE. The RMSE is monotonically reducing as the SNR increases in both methods. However, the RMSE of DPE is significantly lower than IPE at low SNR and the gap diminishes as the SNR increases. At SNR values between 17.5 and 22.5 dB the RMSE of DPE is about 2.5 times lower than IPE. For SNR values between 25 dB and 35 dB the RMSE reduction factor is about 1.5 and for SNR higher than 40 dB the RMSE of both methods consolidate. A comparison of the SNR for the same RMSE, reveals about 5 dB or 3 dB advantage for the DPE for SNR in the range 17.5-35 dB, for GDOP=7.15 and GDOP=3 dB, respectively.

The analytical RMSE expression given by the square root of the expression in (34) lower bounds the simulation results. The bound becomes tighter as the SNR increases, but a small constant gap remains. Both the DPE and the IPE converge close to the small error approximation at high enough SNR starting from some SNR threshold, but the SNR threshold of DPE is far less in the DPE than the IPE. The analysis gives a good approximation to the DPE simulation results for SNR>25 dB in FIG. 13 and SNR>20 dB in FIG. 14 and to the IPE results only for SNR>40 dB in both figures.

FIGS. 15 and 16 show RMSE vs SNR results for the same target and reference stations positions as FIGS. 13 and 14, respectively, with two additional reference stations located at (25,50,50)m and (25,0,50)m. Thus there were all-together six reference stations. The GDOP values for target locations (25,25,25) and (25,25,47.5) are 1.26 and 5.09, respectively. Compared to the results in FIGS. 13 and 14 the RMSE reduction of DPE is larger than of IPE at low SNR. At low SNR the RMSE of DPE is 3 times smaller than IPE. For SNR<30 dB IPE achieves the RMSE of DPE with SNR larger by about 7.5 dB. For SNR>35 dB the performance of DPE and IPE consolidate. The RMSE analysis gives a good approximation to simulation results for SNR>20 dB and SNR>17.5 dB in FIGS. 15 and 16 respectively. These SNR values are smaller by 2.5 and 5 dB compared to the results in FIGS. 13 and 14, respectively. At high SNR there is no noticeable gap between the analysis and simulation results of the target with GDOP=1.26 and for the target with GDOP=5.09 a small gap remains similarly to the results in FIG. 14. It is important to note that a substantial performance gap between DPE and IPE was also observed over other multipath channel models from the IEEE 802.15.4a study group [6] as well as other target and reference stations positions. These results are not presented here due to space limitations.

The results in FIGS. 13-16 show that the DPE attains significant improved position estimation performance compared to IPE at wide SNR range, and that the improvement diminishes with increasing SNR. The performance gap is increased with the increase of the number of reference stations from 4 to 6. The analytical RMSE expression in (34) closely lower bounds the DPE RMSE for most of the interest SNR range at which the position estimation error is smaller than a few meters.

The effect of the number of reference stations will now be discussed. As the number of reference stations increase, there are two contradicting effects that contribute to the performance gap between DPE and IPE. The first is the increase in the probability of a large TOA estimation error (in the first step of IPE) since there are more TOA estimations. The second effect is that the influence of each individual TOA estimation error on the indirect positioning cost function reduces as the number of reference stations increase. The first effect increases the gap and the second effect reduces it. The results presented in this Example show that when increasing the number of reference stations from four to six the performance gap between the direct and indirect methods increases hence the first effect is more dominant.

Without wishing to be bound to any theory, it is assumed that the high GDOP value of the target at (25,25,47.5)m compared to the target at (25,25,25)m is the reason for the analysis small gap from simulation results at high SNR reduced when increasing the number of reference stations from four to six only for the target at (25,25,25)m. As the GDOP increases the effect of the reference stations on the position error is less equal. In case of high GDOP, some reference stations contribute to the position error with larger weight than others. The fact that the increase in the number of reference stations from four to six did not reduce significantly the GDOP of the target at (25,25,47.5)m indicates that the weights of the additional two readers are low compared to the original four reference stations. Thus, when increasing the number of reference stations from four to six the effective number of reference stations for the target at (25,25,47.5)m did not grow as for the target at (25,25,25) which has low GDOP.

Appendix A

Proof of (32)

The derivative of the cost function in (16) with respect to τ^k is given by

$\begin{array}{cc}\frac{\partial C\left(p\right)}{\partial {\tau }^k}=\frac{\partial }{\partial {\tau }^k}\left[\ln \left(R^k\left({\tau }^k\right)\right)+\mathrm{Tr}\left\{{\left(R^k\left({\tau }^k\right)\right)}^{-1}{\hat{R}}^k\right\}\right]=\mathrm{Tr}\left[{\left(R^k\left({\tau }^k\right)\right)}^{-1}\frac{\partial R^k\left({\tau }^k\right)}{\partial {\tau }^k}-{\left(R^k\left({\tau }^k\right)\right)}^{-1}\frac{\partial R^k\left({\tau }^k\right)}{\partial {\tau }^k}{\left(R^k\left({\tau }^k\right)\right)}^{-1}{\hat{R}}^k\right]& \left(\mathrm{EQ}.41\right)\end{array}$

The derivative of (41) with respect to τ^k is given by

$\begin{array}{cc}\frac{{\partial }^{2}C\left(p\right)}{{\partial }^2{\tau }^k}=\mathrm{Tr}\left\{\left[{\left(R^k\left({\tau }^k\right)\right)}^{-1}\frac{\partial R^k\left({\tau }^k\right)}{\partial {\tau }^k}\right]\left[2{\left(R^k\left({\tau }^k\right)\right)}^{-1}{\hat{R}}^k-I\right]+\Gamma \left[I-{\left(R^k\left({\tau }^k\right)\right)}^{-1}{\hat{R}}^k\right]\right\}& \left(\mathrm{EQ}.42\right)\end{array}$
where

$\Gamma =\left[{\left(R^k\left({\tau }^k\right)\right)}^{-1}\frac{{\partial }^2{R}^k\left({\tau }^k\right)}{{\partial }^2{\tau }^k}\right].$
By substituting {circumflex over (R)}^k=R^k (τ^k) into (42) we obtain γ_k given by

$\begin{array}{cc}{\gamma }_k={\left[\frac{{\partial }^2}{{\partial }^2{\tau }^k}C\left(p\right)\right]}_{p=p_t,\hat{v}=v_t}=\mathrm{Tr}\left\{{\left[{{\left(R^k\left({\tau }^k\right)\right)}^{-1}\left[\frac{\partial R^k\left(\tau \right)}{\partial \tau }\right]}_{\tau ={\tau }^k}\right]}^2\right\}& \left(\mathrm{EQ}.43\right)\end{array}$

ψ_k of (31) is derived as follows. From (41), the following can be derived

$\begin{array}{cc}\frac{\partial }{\partial \mathrm{vec}\left({\hat{R}}^i\right)}\left[\frac{\partial C\left(p\right)}{\partial {\tau }^i}\right]=\frac{\partial }{\partial \mathrm{vec}\left({\hat{R}}^i\right)}\mathrm{Tr}\left[-{\left(R^k\left({\tau }^k\right)\right)}^{-1}\frac{\partial R^k\left({\tau }^k\right)}{\partial {\tau }^k}{\left(R^k\left({\tau }^k\right)\right)}^{-1}{\hat{R}}^k\right]& \left(\mathrm{EQ}.44\right)\\ A={\left(R^k\left({\tau }^k\right)\right)}^{-1}\frac{\partial R^k\left({\tau }^k\right)}{\partial {\tau }^k}{\left(R^k\left({\tau }^k\right)\right)}^{-1}& \left(\mathrm{EQ}.45\right)\end{array}$
then

$\begin{array}{cc}\frac{\partial }{\partial \mathrm{vec}\left({\hat{R}}^i\right)}\mathrm{Tr}\left[-{\left(R^k\left({\tau }^k\right)\right)}^{-1}\frac{\partial R^k\left({\tau }^k\right)}{\partial {\tau }^k}{\left(R^k\left({\tau }^k\right)\right)}^{-1}{\hat{R}}^k\right]=\frac{\partial }{\partial \mathrm{vec}\left({\hat{R}}^k\right)}\sum _{i=0}^{L-1}\sum _{j=0}^{L=1}{A}_{i,j}{\hat{R}}_{j,i}^k=-{\mathrm{vec}\left(A\right)}^T& \left(\mathrm{EQ}.46\right)\end{array}$
where A_i,j and {circumflex over (R)}_j,i^k are the elements in the i-th row and j-th column of A and {circumflex over (R)}^k, respectively. (46) is substituted into (31) to provide
ψ_k=vec(A)^TR^k(τ^k){circle around (x)}R^k(τ^k)vec(A)  (EQ. 47)

Using linear algebra equalities [17] (47) can be written as
ψ_k=vec(A)^Tvec(R^k(τ^k)AR^k(τ^k))=Tr(A^TR^k(τ^k)AR^k(τ^k))  (EQ. 48)

(32) is then obtained by substituting (45) into (48) and comparing the result to (43).

References

• [1] S. Gezici and H. V. Poor “Position Estimation via Ultra-Wide-Band Signals”, Proceeding of the IEEE, vol. 97, no. 2, pp. 386-403, February 2009.
• [2] S. Gezici et al. “Localization via Ultra-Wideband Radios”, IEEE Signal Processing Magazine, pp. 70-84, July 2005.
• [3] A. M. Saleh and R. A. Valenzuela, “A Statistical Model for Indoor Multipath Propagation,” IEEE J. Selected. Areas Commun., vol. SAC-5, no. 2, pp. 128-137, February 1987.
• [4] A. F. Molisch, “Ultrawideband Propagation Channels—Theory, Measurements, and Modeling,” IEEE Trans. on Vehicular Technology, vol. 54, no. 5, pp. 1528-1545, September 2005.
• [5] A. F. Molisch et al., “A Comprehensive Standardized Model for Ultrawideband Propagation Channels,” IEEE Trans. on Antennas and Propagation, vol. 54, no. 11, pp. 3151-3166, November 2006
• [6] A. F. Molisch et al., “IEEE 802.15.4a channel model-Final report,” September 2004.
• [7] I. Guvenc and C. C. Chong, “A Survey on TOA Based Wireless Location and NLOS Mitigation Techniques”, IEEE Comm Surveys and Tutorials, vol. 11, no. 3, pp. 107-124, 2009.
• [8] K. W. K. Lui, H. C. So and W. K. Ma, “Maximum A Posteriori Approach to Time-of-Arrival-Based Localization in Non-Line-of-Sight Environment”, IEEE Trans. on Vehicular Technology, vol. 53, no. 3, pp. 1517-1523, Match 2010.
• [9] P. C. Chen, “A non-line-of-sight error mitigation algorithm in location estimation”, in Proc. IEEE Int. Conf. Wireless Commun. Netw. (WCNC), vol. 1, pp. 316-320, September 1999.
• [10] A. Al-Jazzar and J. J. Caffery, “ML and Bayesian TOA location estimators for NLOS environments”, in Proc. IEEE Veh. Technol. Conf. (VTS), vol. 2, pp. 1178-1181, September 2002.
• [11] Y. T. Chan, H. Y. Chin Hang and Pak-chung Ching, “Exact and Approximate Maximum Likelihood Localization Algorithms”, IEEE Trans. on Vehicular Technology, vol. 55, no. 1, pp. 10-16, January 2006.
• [12] A. J. Weiss, “Direct Position Determination of Narrowband Radio Frequency Transmitters”, IEEE Signal Processing Letters, vol. 11, no. 5, May 2004.
• [13] M. Rydstrom, L. Reggiani, E. G. Strom, and A. Svensson, “Suboptimal Soft Ranging Estimators With Applications in UWB Sensor Networks”, IEEE Trans. on Signal Processing, vol. 56, no. 10, pp. 4856-4866, October 2008.
• [14] See Example 1
• [15] J. E. Dennis, Jr. and B. Schnabel, Numerical Methods for Unconstrained Optimization and Non-linear Equations, Englewood Cliffs, N.J.: Prentice-Hall, 1983.
• [16] A. J. Weiss and B. F. Friedlander, “Performance Analysis of Spatial Smoothing with Interpolated Arrays,” IEEE Tans. Signal Processing, vol. 41, no. 5, May 1993.
• [17] D. S. Bernstein, Matrix Mathematics Theory Facts and Formulas, Second edition, Princeton University Press, 2009.
• [18] R. B. Langley, “Dilution of Precision”, GPS World, pp. 52-59, May 1999.

Example 3

Use of Correlated Diversity

TOA Estimation

Following is an exemplified derivation of the TOA estimator for the case in which the diversity information comprises data received from a plurality of correlated channels of a respective reference station

Let y⁰, y¹, . . . , y^K−1 be K received signal column vectors, each being a different multipath vector. Two or more of the vectors can correspond to the same reference station (for example, the same reference station can transmit in K different frequencies). Collectively, the vectors constitute a diversity scheme. Herein, the vectors y⁰, . . . , y^K−1 are referred to as realization vectors.

In the present Example, the K realizations are correlated but all are assumed to have the same time of arrival. Given the TOA each vector y^k is assumed to have a Gaussian distribution. Let {tilde over (y)}=[(y⁰)^T (y¹)^T . . . (y^K−1)^T]^T be a concatenation of all the received vectors and {tilde over (R)}=E{{tilde over (y)}{tilde over (y)}^H} be the autocorrelation of {tilde over (y)}. The vector {tilde over (y)} is assumed to be a Gaussian vector hence the Maximum Likelihood TOA estimation for this case is given by

$\hat{\tau }=\arg {\max }_{\tau }f\left(y^0,y^1,\dots ,y^{K-1}❘\tau \right)=\frac{1}{{\pi }^N\tilde{R}\left(\tau \right)}{e}^{-{\tilde{y}}^H{\tilde{R}}^{-1}\left(\tau \right)\tilde{y}}$

Applying a logarithm and eliminating insignificant terms the following expression is obtained:
{circumflex over (τ)}=arg min_τ{ ln(|{tilde over (R)}(τ)|)+{tilde over (y)}^H{tilde over (R)}⁻¹(τ){tilde over (y)}}  (1)

In cases in which the K realizations are all with the same multipath but different noise measurements, the TOA estimation in (1) results in
{circumflex over (τ)}=arg min_τ{K ln(|R(τ)|)+z^HR⁻¹(τ)z}  (2)
where

$z=\sum _{k=0}^{K-1}\left({\tilde{y}}^k\right)$
and R(τ)=E{y⁰(y⁰)^H}.
Position Estimation

Following is an exemplified derivation of the direct position estimator for a situation in which there is more than one received signal per each link between a reference station and the object. Let there be Q reference stations, wherein the qth station features a rank—K_q diversity scheme (for example, the station can transmit in K_q different frequencies).

Let Ω^q=└y₁⁰ y_q¹ . . . y_q^Kq−1┘ be the set of vectors obtained per link between the target and the q-th reference station. These vectors are referred to as realization vectors. Thus subscript q of y is the reference station index (denoting the link between the target and the q-th reference station) and the superscript of y is the realization index. The vectors in Ω^q may be all received at the reference station or all at the object or part at the target and part at the reference station. Let p be the position vector. The direct position maximum likelihood estimator is given by
arg max_p=ƒ(Ω⁰,Ω¹, . . . ,Ω^Q−1|p)  (3)

It is assumed that Ωⁱ and Ω^j are independent for i≠j hence (3) can be written as

$\begin{array}{cc}p=\arg {\max }_P=\prod _{q=0}^{Q-1}f\left({\Omega }^q❘p\right)& \left(4\right)\end{array}$

Applying a logarithm to (4) the following expression is obtained:

$\begin{array}{cc}p=\arg {\max }_P=\prod _{q=0}^{Q-1}\ln \left(f\left({\Omega }^q❘p\right)\right)& \left(4a\right)\end{array}$

Following is a derivation of ln(ƒ(Ω^q|p)) in (4a) for the general case that the K_q vectors in Ω^q have joint Gaussian distribution.

Let {tilde over (y)}_q=[(y_q⁰)^T (y_q¹)^T . . . (y_q^Kq−1)^T]^T be a concatenation of all vectors in the set Ω^q and {tilde over (R)}_q=E{{tilde over (y)}_q{tilde over (y)}_q^H} be the autocorrelation of {tilde over (y)}_q. Then
ln(ƒ(Ω^q|p))=−N_q ln(π)−ln(|{tilde over (R)}_q(τ)|)−{tilde over (y)}_q^H{tilde over (R)}_q⁻¹(τ){tilde over (y)}_q  (5)
where N_q is the length of {tilde over (y)}_q.

Then the Maximum Likelihood direct position estimator is given by

$\begin{array}{cc}\hat{p}=\arg {\min }_P=\sum _{q=0}^{Q-1}\ln \left({\tilde{R}}_q\left(\tau \right)\right)+{\tilde{y}}_q^H{\tilde{R}}_q^{-1}\left(\tau \right){\hat{y}}_q& \left(6\right)\end{array}$

When the realizations of Ω^q are independent then (5) can be written as:

$\ln \left(f\left({\Omega }^q❘p\right)\right)=\sum _{k=0}^{K_q}\ln \left(f\left(y_q^k|p\right)\right)=-N_q\ln \left(\pi \right)-K_q\ln \left(R_q\left(\tau \right)\right)-\sum _{k=0}^{K_q}{\left(y_q^k\right)}^H{R}_q^{-1}\left(\tau \right)y_q^k$
where R_q(τ)=E{y_q^k(y_q^k)^H}.

When the realizations of the reference stations have identical multipath but different noise measurements then
ln(ƒ(Ω^q|p))=−N_q ln(π)−K_q ln(|R_q(τ)|)−z_q^HR_q⁻¹(τ)z_q
where

$z_q=\sum _{k=0}^{K_q-1}\left(y_q^k\right)$
and R_q(τ)=E{y_q⁰(y_q⁰)^H}.

Example 4

Unsynchronized Signals

In this example a general expression for the TOA likelihood is derived and then used for solving direct position estimation for various synchronization conditions.

Solution of the TOA Likelihood Function

A wideband pulse signal, which is assumed known to the receiver, is transmitted through a multipath channel and observed with additive white Gaussian complex noise. The samples of the received signal are given by

$\begin{array}{cc}{y}_n=\sum _{m=0}^{M-1}{\alpha }_{m}x\left({\mathrm{nT}}_s-{\tau }_m\right)+v_n& \left(\mathrm{EQ}.1\right)\end{array}$

where v_n are statistically independent complex Gaussian noise samples with zero mean and variance σ_v², x(t) is the convolution between the transmitted pulse and the receiver's filters, M is the number of multipath components, α_m and τ_m are the complex coefficient and the delay of the m-th arrival multipath component, respectively. Both α_m and τ_m are statistically independent. The notation of the first arrival time, τ₀, will often be replaced by τ for simplifying the exhibition.

Let y=[y₀ y₁ . . . y_L-1]^T be a vector of length L containing the receiver's observation window samples. Without loss of generality, it is assumed that the beginning of the observation window is at time t=0. Note that τ, the first arrival time with respect to the beginning of the observation window, is a real number. Let N_τ be the integer index of the first sample within the observation vector y that contains the multipath signal, thus, (N_τ−1)T_s<τ≦N_τT_s. The observation vector y can be partitioned into two segments y=[y₀^T y₁^T]^T such that the first segment y₀=[y₀ y₁ . . . y_Nτ−1]^T contains only noise samples, and the second segment y₁=[y_Nτ y_Nτ+1 . . . y_L-1]^T contains the received multipath first cluster plus noise samples as demonstrated in FIG. 3. It is assumed that the segment y₁ terminates at the end of the first multipath cluster. Let v=[v_Nξ v_Nτ+1 . . . v_L-1]^T be the noise samples in y₁ then from (1) y₁ can be written as a sum of statistically independent vectors so that

$\begin{array}{cc}{y}_1=\sum _m^{M_y-1}{\eta }_m+v& \left(\mathrm{EQ}.2\right)\end{array}$
where η_m=α_m[x(N_τT_s−τ_m) x((N_τ+1)T_s−τ_m) . . . x((L−1)T_s−τ_m)]^T and M_y is the number of multipath components included in y₁. The vectors η_m are statistically independent since α_m and τ_m are statistically independent. The central limit theorem (CLT) for random vectors states that the sum of independent random vectors is a Gaussian vector when the number of the summed vectors approaches infinity. The number of non zero elements in the summation for each term of y₁ is proportional to the density of the multipath delays (τ_m) within the duration of x(t). The density is assumed to be large enough to justify the development of the TOA likelihood function based on the Gaussian approximation. Simulation performed by the present inventors show that in practice only three terms were enough to justify the Gaussian approximation.

For each single observation vector y the segments prior and following time τ are statistically independent hence
ƒ(y|τ)=ƒ(y₀|τ)ƒ(y₁|τ)  (EQ. 3)

The noise samples in y₀ are independent complex Gaussian random variable with zero mean, variance σ_v² and independent real and imaginary components and the L−N_τ samples in y₁ are assumed to have Gaussian distribution hence the TOA likelihood function in (3) is given by

$\begin{array}{cc}f\left(y❘\tau \right)=\left[\prod _{n=0}^{N_{\tau }-1}\frac{1}{{\mathrm{\pi \sigma }}_v^2}{e}^{-\frac{{y_n}^2}{{\sigma }_v^2}}\right]\frac{1}{{\pi }^{\left(L-N_{\tau }\right)}R_1\left(\tau \right)}{e}^{-y_1^H{R}_1^{-1}\left(\tau \right)y_1}& \left(\mathrm{EQ}.4\right)\end{array}$
where R₁(τ) is the autocorrelation matrix of y₁ which is dependent on the value of τ. Denote by R₁^i,j(τ) the element in the i-th row and j-th column of R₁(τ) which is given by
R₁^i,j(τ)=E{y_Nτ+iy*_Nτ+j|τ}=∫_t=τ^∞σ_h²(t)x((N_τ+i)T_s−t)x((N_τ+j)T_s−t)dt+σ_v²δ_K(i−j)  (EQ. 5)
where x(μ) initiates at μ=0, δ_K(n) is the Kronecker delta and σ_h²(t) is the channel power delay profile, i.e. the average power vs delay.
TOA Estimation Algorithm

The TOA ML estimator, assuming the receiver obtains K independent observation vectors y⁰, y¹, . . . , y^K−1, all with the same TOA, is given by

$\begin{array}{cc}\hat{\tau }=\arg \underset{\tau }{\max }f\left(y^0,y^1,\dots ,y^{K-1}❘\tau \right)=\arg \underset{\tau }{\max }\prod _{k=0}^{K-1}f\left(y^k❘\tau \right)& \left(\mathrm{EQ}.6\right)\end{array}$

Typically, K=1 but diversity can be used to obtain K>1 vectors. Substituting (4) into (6), applying a logarithm to the result, and eliminating insignificant terms the following Gaussian Model ML TOA estimator (GE) is obtained:

$\begin{array}{cc}{\hat{\tau }}_{\mathrm{GE}}=\arg \underset{\tau }{\min }\left\{\sum _{k=0}^{K-1}{N}_{\tau }\ln \left({\sigma }_v^2\right)+\frac{1}{{\sigma }_v^2}{y_0^k}^2+\ln \left(R_1\left(\tau \right)\right)+{\left(y_1^k\right)}^H{R}_1^{-1}\left(\tau \right)y_1^k\right\}& \left(\mathrm{EQ}.7\right)\end{array}$
where ∥Y∥² denotes the squared norm of the vector Y. The optimal combination of K independent observations is thus realized as the sum of their metrics.

The TOA estimator in (7) was derived based on the assumption that the observation vector y is a Gaussian vector. In simulations performed by the preset inventors it was found that the results are very close to the EML even for average density of only few paths within the pulse duration (x(t)).

The noise variance can be estimated by averaging the energy of the samples received long before the signal arrival time and σ_h²(t) can be approximate as constant within y₁ and then estimated by first evaluating the ratio between the energy of y₁ and its duration and then subtracting the noise variance from the result.

To further simplify the implementation of the GE given in (7) y₁ is further partitioned into two sub-segments y₁=[_1,0^T y_1,1^T]^T where y_1,0=[y_Nτ y_Nτ+1 . . . y_Nτ+d−1]^T, y_1,1=[y_Nτ+d y_Nτ+d+1 . . . y_L-1]^T and d is the number of samples within the duration of x(t). An approximation to the GE estimator is derived by approximating the segments y_1,0 and y_1,1 as statistically independent and obtained the approximate GE (AGE) is given by

$\begin{array}{cc}{\hat{\tau }}_{\mathrm{AGE}}=\arg \underset{\tau }{\min }\left\{\sum _{k=0}^{K-1}\gamma +\frac{1}{{\sigma }_v^2}{y_0^k}^2+{\left(y_{1,0}^k\right)}^H{R}_{1,0}^{-1}\left(\tau -N_{\tau }\right)y_{1,0}^k+{{\tilde{y}}_{1,1}^k}^2\right\}& \left(\mathrm{EQ}.8\right)\end{array}$
where {tilde over (y)}_1,1 is the result of linear convolution between y_1,1 and a high pass filter, R_1,0=E{y_1,0y_1,0^H}, γ=N_τ(ln(σ_v²)−ρ) and ρ is a constant.

The TOA estimator in (8) is simple to implement and has relatively low implementation complexity. The dimension of R_1,0 is d×d which is fixed and small compared to the dimension (L−N_τ)×(L−N_τ) of R₁ and hence the complexity of the TOA estimator in (8) is linearly dependent on L (the length of y) while the complexity of (7) grows with the third power of L. The complexity of the AGE is on the order of the complexity of the simplest TOA estimation methods. These methods are based on comparing the signal energy to a threshold.

The present inventors have developed an analytical expression that approximates the mean square error of the GE estimator in (7) under the assumption that the error is small. This condition holds as the number of independent measurements of the TOA grows which in our case means increasing the number of independent observation vectors, K, by means of diversity. The analytical expression is given by:

$\begin{array}{cc}E\left\{{\left(\hat{\tau }-{\tau }_t\right)}^2\right\}\cong \frac{1}{\mathrm{KTr}\left\{{\left[{R^{-1}\left({\tau }_t\right)\left[\frac{\partial R\left(\tau \right)}{\partial \tau }\right]}_{\tau ={\tau }_t}\right]}^2\right\}}& \left(\mathrm{EQ}.9\right)\end{array}$

The expression in (9) is also the expression of the Cramer-Rao lower bound (CRLB) on the TOA estimation MSE when the received signal y is a Gaussian vector which is the case when the multipath is dense or when the channel has only a single path. Thus, in these two cases the GE MSE is higher or equal to (9) and it is expected to reach the CRLB as the number of independent observation vectors (K) increases.

Direct Position Estimation

The ML DPE will now be formulated for various cases of synchronization between the reference stations and object.

Consider a wireless network where there are Q reference stations at known positions and a target at unknown position. The 3D position of the q-th reference station and of the target are denoted by p_ref^q=[x^q y^q z^q]^T and p=[x y z]^T, respectively where q=0,1, . . . , Q−1. Let ζ^q(p) be the signal propagation time between the target and the q-th reference station which is given by

$\begin{array}{cc}{\zeta }^q\left(p\right)=\frac{1}{c}\sqrt{{\left(p_{\mathrm{ref}}^q-p\right)}^T\left(p_{\mathrm{ref}}^q-p\right)}.& \left(\mathrm{EQ}.10\right)\end{array}$
where c is the speed of light.

Denote by Ω^q the set of all observation vectors related to the link between the target and the q-th reference station. The ML DPE is given by

$\begin{array}{cc}\hat{p}=\arg \underset{p}{\max }f\left({\Omega }^0,{\Omega }^1,\dots ,{\Omega }^{Q-1}|p\right)& \left(\mathrm{EQ}.11\right)\end{array}$

For each target position p there are Q corresponding propagation times ζ⁰(p), ζ¹(p), . . . , ζ^Q−1(p). It is assumed that the channels between the target and different reference stations are statistically independent thus Ωⁱ and Ω^j are statistically independent for i≠j hence the ML estimator in (11) can be written as

$\begin{array}{cc}\hat{p}=\arg \underset{p}{\max }\prod _{q=0}^{Q-1}f\left({\Omega }^q|{\zeta }^q\left(p\right)\right)& \left(\mathrm{EQ}.12\right)\end{array}$

The solution of (12) depends on the clocks synchronization conditions of the reference stations and target. When the clocks of the reference stations are synchronized then the target position can be estimated by transmission from one side only (target or reference stations). It is assumed, without loss of generality, that the target is the transmitting side. Assuming, for now, only one received signal per reference station thus Ω^q=y^q. Let t₀ be the transmission time and τ_y^q be the TOA at the q-th reference station which is given by τ_y^q=t₀+ζ^q. In case the reference stations clocks are synchronized to the target clock then the time t₀ is known to the reference stations and the ML DPE is given by

$\begin{array}{cc}{\hat{p}}_s=\arg \underset{p}{\max }\prod _{q=0}^{Q-1}f\left(y^q|{\tau }_y^q=t_0+{\zeta }^q\left(p\right)\right)& \left(\mathrm{EQ}.13\right)\end{array}$

When the reference stations are synchronized among themselves but not synchronized to the target then for the reference stations receivers the target transmission time (t₀) is a random variable with uniform distribution. In this case the position estimator could be obtained by either averaging the product in (13) over t₀ or by maximizing it over t₀. The later is the Generalized ML (GML) approach which is simpler and not expected to be degraded since t₀ has uniform distribution. The GML DPE for synchronized reference stations which are not synchronized to the target is given by

$\begin{array}{cc}{\hat{p}}_{\mathrm{sr}}=\arg \underset{p}{\max }\left[\underset{t_0}{\max }\prod _{q=0}^{Q-1}f\left(y^q|{\tau }_y^q=t_0+{\zeta }^q\left(p\right)\right)\right]& \left(\mathrm{EQ}.14\right)\end{array}$

In the next considered scenario, the reference stations clocks are not synchronized one to the other nor to the target's clock. In this case measuring the distance between the target and reference stations requires transmission by both target and reference stations. Consider a single link between the q-th reference station and target. Let t₁ and t₂ be the transmission time of the reference station and target respectively. Let Δ^q=t₂−t₁ be the unknown time difference between the transmission times. Let τ_y^q be the time difference between the TOA at the reference station and the reference station transmission time t₁, and τ_x^q be the time difference between the TOA at the target and the target transmission time t₂ i.e.
τ_y^q=t₂+ζ^q(p)−t₁=ζ^q(p)+Δ^q  (EQ. 15)
τ_x^q=t₁+ζ^q(p)−t₂=ζ^q(p)−Δ^q  (EQ. 16)

Let x^q be the received samples at the target obtained from the q-th reference station transmission. In the asynchronous case one has Ω^q={x^q, y^q}. By substituting (15) and (16) into (12), assuming Δ^q(p) has uniform distribution and applying the GML approach the following asynchronous DPE is obtained:

$\begin{array}{cc}{\hat{p}}_{\mathrm{as}}=\arg \underset{P}{\max }\prod _{q=0}^{Q-1}\underset{{\Delta }_q}{\max }f\left(y^q,x^q|{\tau }_y^q={\zeta }^q\left(p\right)+{\Delta }^q,{\tau }_x^q={\zeta }^q\left(p\right)-{\Delta }^q\right)& \left(\mathrm{EQ}.17\right)\end{array}$

Solving the joint likelihood function ƒ(y^q, x^q|τ_y^q=ζ^q(p)+Δ^q,τ_x^q=ζ^q(p)−Δ^q) depends on the correlation assumption between the channel from target to reference station and vice versa. In case they are statistically independent (for example when the target and reference stations have different carrier frequencies) then the position estimator in (17) can be written as

$\begin{array}{cc}{\hat{p}}_{\mathrm{as}}=\arg \underset{P}{\max }\prod _{q=0}^{Q-1}\underset{{\Delta }_q}{\max }f\left(y^q|{\tau }_y^q={\zeta }^q\left(p\right)+{\Delta }^q\right)f\left(x^q|{\tau }_x^q={\zeta }^q\left(p\right)-{\Delta }^q\right)& \left(\mathrm{EQ}.18\right)\end{array}$

Another case is two way location (transmission from target and also base station sides) where the base stations are synchronized among themselves but not synchornized to the target. In this case the unknown time difference between the transmittions of the target and base station is the same for all base stations, i.e., Δ⁰=Δ¹= . . . =Δ^Q−1=Δ. The direct position estimation for this case can be written as:

$\begin{array}{cc}{\hat{p}}_{\mathrm{as}}=\arg \underset{P}{\max }\prod _{q=0}^{Q-1}\underset{{\Delta }_q}{\max }f\left(y^q|{\tau }_y^q={\zeta }^q\left(p\right)+{\Delta }^q\right)f\left(x^q|{\tau }_x^q={\zeta }^q\left(p\right)-{\Delta }^q\right)& \left(\mathrm{EQ}.19\right)\end{array}$

For the various cases of synchronization the present inventors have succeeded to formulate the ML direct position estimators as a product of the TOA likelihood function ƒ(y|τ). This likelihood function was approximated and derived above. By substituting the result given in (4) into (13), (14) and (18) novel direct position estimators are obtained. The estimators complexity can be reduced by applying the approximations that yield (8) from (7).

Example 5

General Formulation of the Maximum Likelihood Direct Position Estimation

Let k={0,1, . . . , K−1} be the index of the reference node. Let there be Q_k observation vectors corresponding to the link between the target and the k-th base station. These vectors can be obtained with use of diversity technique such as space, time, frequency or antenna polarization. The Q_k vectors can be statistically independent or also correlated. Let y_q^k be the q-th observation column vector from the Q_k vectors that correspond to the link between the target and the k-th base station. Let {tilde over (y)}^k=[(y₀^k)^T (y₁^k)^T . . . (y_Qk−1^k)^T]^T be a column vector which is a concatenation of all received column vectors per link between the target and the k-th reference. The vector {tilde over (y)}^k may include vectors which were received at the reference node or target or both. It is assumed the {tilde over (y)}^k has a Gaussian distribution. Let R^k(p)=E{{tilde over (y)}^k({tilde over (y)}^k)^H|p} be the autocorrelation of {tilde over (y)}^k given the position p and μ^k(p)=E{{tilde over (y)}^k|p} be the mean of {tilde over (y)}^k given p. It is assumed that {tilde over (y)}^k and {tilde over (y)}^m are statistically independent for k≠m. The Maximum Likelihood direct position estimation (ML DPE) is given by

$\begin{array}{cc}\hat{p}=\arg \underset{p}{\max }\prod _{k=0}^{K-1}f\left({\tilde{y}}^k|p\right)& \left(\mathrm{EQ}.1\right)\end{array}$
where ƒ({tilde over (y)}^k|p) is the probability distribution function which is assumed Gaussian with autocorrelation R^k(p) and mean μ^k(p). Then the ML DPE can be written as

$\begin{array}{cc}\hat{p}=\arg \underset{p}{\min }\sum _{k=0}^{K-1}\ln \left(R^k\left(p\right)\right)+{\left({\tilde{y}}^k-{\mu }^k\right)}^H{\left(R^k\left(p\right)\right)}^{-1}\left({\tilde{y}}^k-{\mu }^k\right)& \left(\mathrm{EQ}.2\right)\end{array}$

Let ξ^k(p) be the signal propagation time between the target at position p and the k-th reference node. Let τ^k be the time of arrival at the k-th reference node which is given by τ^k=t₀+ξ^k(p) where t₀ is the transmission time.

If the transmission time t₀ is known then (2) can be written as

$\begin{array}{cc}\hat{p}=\arg \underset{p}{\min }\sum _{k=0}^{K-1}\ln \left(R^k\left({\tau }^k=t_0+{\xi }^k\left(p\right)\right)\right)+{\left({\tilde{y}}^k-{\mu }^k\right)}^H{\left(R^k\left({\tau }^k=t_0+{\xi }^k\left(p\right)\right)\right)}^{-1}\left({\tilde{y}}^k-{\mu }^k\right)& \left(\mathrm{EQ}.3\right)\end{array}$

If the transmission time t₀ is unknown but the all reference nodes are tune synchronized then the position can be obtained by

$\begin{array}{cc}\hat{p}=\arg \underset{p}{\min }\left\{{\min }_{t_0}\sum _{k=0}^{K-1}\ln \left(R^k\left({\tau }^k=t_0+{\xi }^k\left(p\right)\right)\right)+{\left({\tilde{y}}^k-{\mu }^k\right)}^H{\left(R^k\left({\tau }^k=t_0+{\xi }^k\left(p\right)\right)\right)}^{-1}\left({\tilde{y}}^k-{\mu }^k\right)\right\}& \left(\mathrm{EQ}.4\right)\end{array}$

Example 6

Frequency Domain

In this example a TOA Maximum Likelihood estimation in the frequency domain, according to some embodiments of the present invention is described.

Let y be a vector of length N which is the DFT of the received vector and let F be the autocorrelation of y, so that F=E{yy^H}. The Maximum Likelihood estimation of the TOA, {circumflex over (τ)}, can be written as
{circumflex over (τ)}=arg min_τy^HG(τ)FG^H(τ)y  (EQ. 1)
where G(τ)=diag{1,z, . . . , z^N−1}; z=e^j2πτ/T and T is the time duration of vector y.
The cost function to be minimized in Equation (1) can be written as,
Q(z)=y^HG(τ)FG^H(τ)y=g^H(z)Y^HFYg(z)  (EQ. 2)
where Y=diag{y} and g(z)=[1,z,z², . . . , z^N−1]^T. Note that Q(z) is a polynomial in z given by
Q(z)=a_—(N−1)z^−(N−1)+a_—(N−2)z^−(N−2)+ . . . a_N−2z^N−2+a_N−1z^N−1  (EQ. 3)

The coefficients of the polynomial terms are the sum of the diagonals of the matrix Y^HFY. Since this matrix is Hermitian the coefficients have the property: a_i=a*_−i. Hence, the cost function in Equation (3) can be written as
Q(z)=Re{a*₀+2a*₁z⁻¹ . . . +2a*_N−2z^−(N−2)+2a*_N−1z^−(N−1)}.  (EQ. 4)
Let x=[a*₀,2a*₁, . . . , 2a*_N−1] then the ML estimator in (1) can be implemented by the following stages:

• • (a) Applying a Discrete Fourier Transform (DFT) or Fast Fourier Transform (FFT) on the vector x with any desired resolution using zero padding prior to the DFT/FFT operation or first applying decimation to vector x followed by zero padding and then DFT/FFT. Note that the result is a complex expression having a real part and an imaginary part. For the next stage, only the real part can be used.
  • (b) Searching for the frequency that attains minimum real DFT/FFT. From the definition of z=e^−j2πτ/T let ƒ=τ/T be the frequency and ƒ_M be the frequency that attains minimum real value then, according to some embodiments of the present invention the TOA estimation is extracted from ƒ_M by {circumflex over (τ)}=ƒ_M×T.

Alternatively, the ML estimator in (1) can be implemented by calculating DFT at relatively large frequency spacing (coarse resolution) that span the range of minimum and maximum TOA values. Followed by searching for the frequency that attains minimum real value and then calculating DFT at the vicinity of the minimal frequency with small frequency spacing (fine resolution) and conducting another search for the minimal frequency. From the definition of z=e^−j2πτ/T let ƒ=τ/T be the frequency and ƒ_M be the frequency that attains minimum real value then, according to some embodiments of the present invention the TOA estimation is extracted from ƒ_M by {circumflex over (τ)}=ƒ_M×T.

Example 7

Propagation Delay Estimation Based on Two Way Signals

In this example, two way signals are used according to some embodiments of the present invention to estimate the propagation delay between the object and a reference station. This technique is particularly useful for estimating the range to the object in an indirect position estimation.

Let y_ref and y_tar be column vectors of the received signals at the reference and target respectively. Let {tilde over (y)}=[y_ref^T y_tar^T]^T be column vector which is the concatenation of y_ref and y_tar. Let ζ be the signal propagation time between the object and the reference station. Let Δ be the difference between the transmission times of the object and the reference station. Let τ_ref be the difference between the time of arrival at the reference station and the transmission time of the reference station. Let τ_tar be the difference between the time of arrival at the target and the transmission time of the target.

It is assumed, without loss of generality that y_ref starts from the transmission time of the reference station and y_tar starts from the transmission time of the target. The ML estimator of the propagation delay (using the notations of Example 4 above) is given by
{circumflex over (ζ)}=arg min_ζ max_Δ ƒ(y_ref,y_tar|τ_ref=ζ+Δ,τ_tar=ζ−Δ)=arg min_ξ min_Δ ln(|{tilde over (R)}(ζ,Δ)|)+{tilde over (y)}^H{tilde over (R)}⁻¹(ζ,Δ){tilde over (y)}
where {tilde over (R)}(ζ,Δ)=E{{tilde over (y)}{tilde over (y)}^H|ζ,Δ}.

Example 8

Position or Range Estimation by Propagation Delay Modeling

The power delay profile relates to the distance between the object and the reference station. For example, the channel delay spread (the effective duration of the power delay) is shorter as the distance becomes smaller and also the ratio between the first arrival energy and the following multipath arrivals energies is larger as the distance is smaller.

The relation between the distance and the power delay profile can according to some embodiments of the present invention be modeled, and incorporated into the direct position estimation or into the propagation delay estimation. For each candidate position or candidate propagation delay the distance between the object and reference station is given and the power delay profile used in the matrix R is according to the given distance.

A representing example for the propagation delay estimation is given by
{circumflex over (ζ)}=arg min_ξ min_Δ ln(|{tilde over (R)}(ζ,Δ,σ_h²(ζ))|)+{tilde over (y)}^H{tilde over (R)}⁻¹(ζ,Δ,σ_h²(ζ)){tilde over (y)}
where σ_h²(ζ) is the power delay profile as a function of the propagation delay ζ between the object and the respective reference station and Δ is the difference between the transmission times of the reference station and target.

A representing example for the position estimation where the object and reference stations are assumed to be synchronized (which is not to be considered as limiting) is given by

$\hat{p}=\arg \underset{p}{\min }\sum _{k=0}^{K-1}{\ln \left(R^k({\tau }^k,{\sigma }_h^2\left({\tau }^k\right))+{\left({\tilde{y}}^k\right)}^H\left(R^k{\tau }^k,{\sigma }_h^2\left({\tau }^k\right)\right)\right)}^{-1}{y}^k$
where τ^k is the TOA at the kth reference station (corresponding to the position hypothesis p), which is proportional to the distance, and σ_h²(τ^k) is the power delay profile as a function of τ^k.

Although the invention has been described in conjunction with specific embodiments thereof, it is evident that many alternatives, modifications and variations will be apparent to those skilled in the art. Accordingly, it is intended to embrace all such alternatives, modifications and variations that fall within the spirit and broad scope of the appended claims.

All publications, patents and patent applications mentioned in this specification are herein incorporated in their entirety by reference into the specification, to the same extent as if each individual publication, patent or patent application was specifically and individually indicated to be incorporated herein by reference. In addition, citation or identification of any reference in this application shall not be construed as an admission that such reference is available as prior art to the present invention. To the extent that section headings are used, they should not be construed as necessarily limiting.

Claims

1. A method of estimating a position of an object, comprising:

receiving a plurality of signals transmitted between each of a respective plurality of reference stations of known position and the object;

for each received signal: processing said signal using a matrix featuring a characteristic power delay profile and a characteristic signal waveform to provide a processed signal, and partitioning said signal into at least a first segment containing noise and a second segment containing at least a multipath signal corresponding to a plurality of paths between the object and a respective reference station; and

calculating a direct estimate of the position of the object based, at least in part, on said processed signals and on said segments.

2. The method according to claim 1, wherein said matrix is based on a noise variance.

3. The method according to claim 1, wherein said matrix is based on an autocorrelation matrix of said received signal and is a function of the position of the object.

4. The method according to claim 1, wherein said processing comprising modeling said signal using a function characterized by said matrix, said function being selected from the group consisting of a Gaussian function, a Lorentzian function, a hyperbolic secant function, a logistic distribution, a Fourier transform and a wavelet transform.

5. The method according to claim 1, wherein at least one of the signals is a diversity signal having diversity measurements, and the method further comprising processing said diversity measurements, wherein said estimate is based on said processed diversity measurements.

6. The method according to claim 1, wherein at least two of the reference stations are synchronized, and the object is unsynchronized with any reference station.

7. The method according to claim 1, wherein at least two of the reference stations are unsynchronized, wherein the object is unsynchronized with any reference station, and wherein a signal transmitted between the object and each of said at least two reference stations is a two-way signal.

8. The method according to claim 1, wherein at least two of the reference stations are synchronized, wherein the object is unsynchronized with any reference station, and wherein a signal transmitted between the object and each of said at least two reference stations is a two-way signal.

9. The method according to claim 1, wherein said reference stations are located at indoor positions, and the method is utilized for estimating an indoor position of the object.

10. The method according to claim 1, wherein the object is selected from the group consisting of a moving object, an inventory item, a human, an animal, a vehicle and a mobile phone.

11. The method according to claim 1, wherein said reference stations are located at outdoor positions, and the method is utilized for estimating an outdoor position of the object.

12. The method according to claim 1, wherein the reference stations are Access Points, transmitting and receiving over a wireless network according to the LAN protocol, and wherein the object is a client of said wireless network.

13. The method according to claim 1, wherein the reference stations system are satellites are transmitting according to a Global Navigation Satellite System (GNSS) protocol, and the object comprises a receiver for GNSS signals, where said GNSS is selected from the group consisting of a Global Positioning System (GPS), a Global Orbiting Navigation Satellite System (GLONASS) and a Galileo navigation system.

14. The method according to claim 1, wherein said characteristic power delay profile is calculated based on an estimate of a signal-to-noise ratio.

15. The method according to claim 1, further comprising jointly estimating a plurality of pairs of candidate power delay profiles and respective candidate direct estimations of the position of the object, scoring each pair of said plurality of pairs, and selecting the direct estimation of the position of the object based on said score.

16. The method according to claim 1, further comprising jointly estimating a plurality of pairs of candidate signal waveforms and respective candidate direct estimations of the position of the object, scoring each pair of said plurality of pairs, and selecting the direct estimation of the position of the object based on said score.

17. A system for estimating a position of an object, comprising:

a plurality of reference stations of known positions; and

a data processor, configured for executing the method according to claim 1.

18. A non-transitory computer software product, comprising a computer-readable medium in which program instructions are stored, which instructions, when read by a data processor, cause the data processor to receive a plurality of signals transmitted between each of a respective plurality of reference stations of known position and an object and execute the method according to claim 1.

19. A method of estimating a time-of-arrival (TOA) of a signal transmitted between an object and a reference station, comprising:

receiving a multipath signal having a plurality of components respectively corresponding to a plurality of paths between the object and the reference station;

processing said signal using a matrix featuring a characteristic power delay profile and a characteristic signal waveform to provide a processed signal; and

calculating an estimate of the TOA of the signal based, at least in part, on said processed signal;

wherein a signal transmitted between the object and each of said at least two reference stations is a two-way signal, and wherein said matrix features a difference between a signal transmission time of the object and a signal transmission time of a respective reference station.

20. The method according to claim 19, wherein said reference stations are located at indoor positions, and the method is utilized for estimating a distance between a reference station and an object.

21. The method according to claim 19, wherein said matrix is the same for all received signals.

22. The method according to claim 19, wherein said matrix is different for at least two of the received signals.

23. The method according to claim 19, wherein said characteristic power delay profile is calculated based on an estimate of a signal-to-noise ratio.

24. The method according to claim 19, further comprising updating said characteristic power delay profile using data collected from signal transmissions between said reference stations and other objects.

25. The method according to claim 19, further comprising, obtaining at least one additional power delay profile, wherein said processing said signal is executed separately for each power delay profile.

26. The method according to claim 19, further comprising jointly estimating a plurality of pairs of candidate power delay profiles and respective candidate TOAs, scoring each pair of said plurality of pairs, and selecting the TOA estimate based on said score.

27. The method according to claim 19, further comprising, obtaining at least one additional signal waveform, wherein said processing said signal is executed separately for each signal waveform.

28. The method according to claim 19, further comprising:

jointly estimating a plurality of sets of: candidate power delay profiles, candidate signal waveforms, and respective candidate TOAs,

scoring each set of said plurality of sets, and

selecting the TOA estimate based on said score.

29. A system for estimating a position of an object, comprising:

a plurality of reference stations of known positions; and

a data processor, configured for executing the method according to claim 19 to provide a TOA of signals transmitted between the object and each of said reference stations and for calculating the position based on said TOA.

30. A non-transitory computer software product, comprising a computer-readable medium in which program instructions are stored, which instructions, when read by a data processor, cause the data processor to receive a plurality of signals transmitted between each of a respective plurality of reference stations of known position and an object and execute the method according to claim 19.

31. A method of estimating a position of an object, comprising:

receiving a plurality of signals transmitted between each of a respective plurality of reference stations of known position and the object;

for each signal, estimating a time-of-arrival (TOA) of said signal using the method according to claim 19, thereby providing a plurality of TOA estimates; and

calculating an indirect estimate of the position of the object based, at least in part, on said plurality of TOA estimates.

32. The method according to claim 31, wherein said calculation of said indirect estimate comprises utilizing at least one procedure selected from the group consisting of a weighted least-squares procedure, a non-weighted least-squares procedure and a Kalman filter procedure.

33. A method of estimating a position of an object, comprising:

receiving a plurality of signals transmitted between each of a respective plurality of reference stations of known position and the object;

for each received signal: processing said signal using a matrix featuring a characteristic power delay profile and a characteristic signal waveform to provide a processed signal;

calculating a direct estimate of the position of the object based, at least in part, on said processed signals;

jointly estimating a plurality of pairs of candidate power delay profiles and respective candidate direct estimations of the position of the object;

scoring each pair of said plurality of pairs; and

selecting the direct estimation of the position of the object based on said score.

34. A method of estimating a position of an object, comprising:

receiving a plurality of signals transmitted between each of a respective plurality of reference stations of known position and the object;

for each received signal: processing said signal using a matrix featuring a characteristic power delay profile and a characteristic signal waveform to provide a processed signal; and

calculating a direct estimate of the position of the object based, at least in part, on said processed signals;

wherein the method comprises:

jointly estimating a plurality of pairs of candidate signal waveforms and respective candidate direct estimations of the position of the object;

scoring each pair of said plurality of pairs; and

selecting the direct estimation of the position of the object based on said score.

35. A method of estimating a time-of-arrival (TOA) of a signal transmitted between an object and a reference station, comprising:

receiving a multipath signal having a plurality of components respectively corresponding to a plurality of paths between the object and the reference station;

processing said signal using a matrix featuring a characteristic power delay profile and a characteristic signal waveform to provide a processed signal; and

calculating an estimate of the TOA of the signal based, at least in part, on said processed signal;

wherein the method comprises:

jointly estimating a plurality of pairs of candidate power delay profiles and respective candidate TOAs;

scoring each pair of said plurality of pairs; and

selecting the TOA estimate based on said score.

36. A method of estimating a time-of-arrival (TOA) of a signal transmitted between an object and a reference station, comprising:

receiving a multipath signal having a plurality of components respectively corresponding to a plurality of paths between the object and the reference station;

processing said signal using a matrix featuring a characteristic power delay profile and a characteristic signal waveform to provide a processed signal; and

calculating an estimate of the TOA of the signal based, at least in part, on said processed signal;

wherein the method comprises:

jointly estimating a plurality of sets of: candidate power delay profiles, candidate signal waveforms, and respective candidate TOAs;

scoring each set of said plurality of sets; and

selecting the TOA estimate based on said score.

37. A method of estimating a position of an object, comprising:

receiving a plurality of signals transmitted between each of a respective plurality of reference stations of known position and the object;

for each signal, estimating a time-of-arrival (TOA) of said signal executing a method of estimating a TOA of said signal to provide a plurality of TOA estimates; and

calculating an indirect estimate of the position of the object based, at least in part, on said plurality of TOA estimates;

wherein said method of estimating a TOA of said signal comprises: receiving a multipath signal having a plurality of components respectively corresponding to a plurality of paths between the object and the reference station, processing said signal using a matrix featuring a characteristic power delay profile and a characteristic signal waveform to provide a processed signal, and calculating an estimate of said TOA of said signal based, at least in part, on said processed signal; and

wherein said calculation of said indirect estimate comprises utilizing at least one procedure selected from the group consisting of a weighted least-squares procedure, a non-weighted least-squares procedure and a Kalman filter procedure.