Ranging method and apparatus

Abstract

The present invention relates to a method of determining a distance or location of a remote device or reflector, the method comprising: receiving a signal from a remote signal transmitter associated with or contained in or attached to the remote device; estimating a first propagation delay associated with the received signal, wherein the first propagation delay represents a first candidate for a correct propagation delay; deriving a relationship between the first candidate and one or more other candidates for a correct propagation delay from the received signal; determining a plurality of other candidates for a correct propagation delay based on said relationship; generating a likelihood histogram based on said candidates for a correct propagation delay; selecting a propagation delay from said candidates based on the likelihood histogram; and determining a distance or location of the remote device or reflector using the selected propagation delay.

Description

CROSS-REFERENCE TO RELATED APPLICATIONS

This application is a National Stage Application under 35 U.S.C. § 371 of PCT Application No. PCT/EP2016/055380, filed Mar. 11, 2016, which claims priority to German Patent Application No. 102015103605.3, filed Mar. 11, 2015, the entire disclosures of which are hereby incorporated by reference.

BACKGROUND OF THE INVENTION

Ranging at high accuracy is a core problem in many technical fields. For example, high precision ranging can increase efficiency and productivity in agriculture by using robots for pruning, weeding and crop-spraying. In marine navigation, ranging with high accuracy can help to enter a port with big ships. Also the landing phase of an airplane can be automated by reliable high precision ranging techniques.

Ranging with radio systems is usually performed by measuring the propagation delay of an electro-magnetic wave with known structure. If the radio signal is transmitted on a high carrier frequency, it is understood that the carrier phase conveys significant information about the delay parameter. However, as the mapping between phase and delay parameter is ambiguous, it is believed that the carrier phase can only be exploited by combining measurements attained with different signal sources.

In the application of satellite-based synchronization and navigation (GPS, GLONASS, Galileo, etc.), high precision is therefore achieved by performing three independent steps, as illustrated in FIG. 1. An acquisition algorithm delivers some initial knowledge about the range between transmitter 1 and the receiver for all available transmitters l=1, . . . , L. Here this initial knowledge is characterized by a Gaussian random variable with mean μ_init^(l) and variance σ_init². For each transmitter an individual tracking module then measures and tracks the baseband delay and the carrier phase of the radio signal as independent parameters. In practice, this is done with two control loops, the delay-locked loop (DLL) and the phase-locked loop (PLL). With the DLL only a coarse ranging solution can be obtained. The carrier phase ζ^(l) can be measured with much higher precision. However, the carrier phase is periodic with 2π and the measurement ζ^(l), is, thus, only given by some fraction of a cycle, i.e. the integer number Ψ^(l) of whole cycles is not known at the receiver

${\widehat{ϛ}}^{\left(l\right)}={\omega }_c{\tau }^{\left(l\right)}+2{\mathrm{\pi \Psi }}^{\left(l\right)}+e_{\varsigma }^{\left(l\right)},{\widehat{ϛ}}^{\left(l\right)}={\omega }_c{\tau }^{\left(l\right)}+2{\mathrm{\pi \Psi }}^{\left(l\right)}+e_{\varsigma }^{\left(l\right)},$ where

$w_c=2\pi f_c$ is the carrier frequency,

${\tau }^{\left(l\right)}=\frac{r^{\left(l\right)}}{c}$ is the propagation-delay, c velocity of light and e_ζ^(l) the measurement error. Resolving the integer Ψ^(l) precisely with measurements from one transmitter is not possible. Measurements from multiple transmitters and multiple time instances must be combined in order to resolve the ambiguity problem and obtain a ranging solution with high accuracy.

Approaches pertaining to the technological background of the present invention are described in the following references:

• [1] G. Seco-Granados, J. A. Lopez-Salcedo, D. Jumenez-Baňos and G. Lopez-Risueňo, “Challenges in Indoor Global Navigation Satellite Systems,” IEEE Signal Processing Magazine, vol. 29, no. 2, pp. 108-131, 2012.
• [2] P. Misra and P. Enge, “Global Positioning System—Signals, Measurements, and Performance”, Second Edition, Ganga-Jamuna Press, 2006.
• [3] G. Blewitt, “Carrier Phase Ambiguity Resolution for the Global Positioning System Applied to Geodetic Baselines up to 2000 km,” Journal of Geophysical Research, vol. 94, no. B8, pp. 10187-10203, 1989.
• [4] P. Teunissen, “Least-Squares Estimation of the Integer GPS Ambiguities”, Invited lecture, Section IV “Theory and Methodology,” Proc. Of Gen. Meet. of the Int. Assoc. of Geodesy, Bejing, China, pp. 1-16, 1993.
• [5] P. Teunissen, “A new method for fast carrier phase ambiguity estimation,” Proc. of IEEE Pos., Loc. and Symp. (PLANS), Las Vegas, USA, pp. 562-573, 1994.
• [6] P. Teunissen, “The least-squares ambiguity decorrelation adjustment: a method for fast GPS ambiguity estimation,” Journal of Geodesy, vol. 70, pp. 65-82, 1995.
• [7] P. Teunissen, “Statistical GNSS carrier phase ambiguity resolution: A review,” Proc. of the l-th IEEE Workshop of Statistical Signal Processing (SSP), pp. 4-12, 2001.
• [8] C. Günther, P. Henkel, “Integer Ambiguity Estimation for Satellite Navigation,” IEEE Transactions on Signal Processing, vol. 60, no. 7, pp. 3387-3393, 2012.
• [9] S. M. Kay, “Fundamentals of Statistical Signal Processing: Estimation Theory,” Pretice Hall, 1993.
• [10] B. Ristic, S. Arulampalam and N. Gordon, “Beyond the Kalman Filter—Particle Filters for Tracking Applications,” Artech House Inc., 2004.
• [11] A. Doucet and A. Johansen, “A Tutorial on Particle Filtering and Smoothing: Fifteen Years later,” Oxford Handbook of Nonlinear Filtering, Oxford University Press, 2011.
• [12] A. Kong, J. Liu and W. Wong, “Sequential imputations and Bayesian missing data problems,” Journal of the American Statistical Association, vol. 89, no. 425, pp. 278-288, 1994.
• [13] A. Doucet, S. Godsill and C. Andrieu, “On sequential Monte Carlo sampling methods for Bayesian filtering,” Statistics and Computing, vo. 10, pp. 197-208, 2000.
• [14] N. J. Gordon, D. J. Salmond and A. F. M. Smith, “Novel approach to nonlinear/non-Gaussian Bayesian state estimation,” IEE Proc. F Radar Signal Process., vol. 140, no. 2, pp. 107-113, 1993.
• [15] T. Li, T. P. Sattar, Q. Han and S. Sun, “Roughening Methods to Prevent Sample Impoverishment in the particle PHD filter,” 16th International Conference on Information Fusion (FUSION), 2013.
• [16] K. Borre, D. M. Akos, N. Bertelsen, P. Rinder, S. H. Jensen, A Software-Defined GPS and Galileo Receiver: A Single-Frequency Approach, Birkhäuser Boston, 2007.
• [17] M. Simandl, J. Kralovec, P. Tichavsky, “Filtering, Predictive, and Smoothing Cram'er-Rao Bounds for Discrete-time Nonlinear Dynamic Systems”, Automatica, vol. 37, no. 11, pp. 1703-1716, 2001.
• [18] H. L. Van Trees, K. L. Bell, “Bayesian Bounds for Parameter Estimation and Nonlinear Filtering/Tracking,” lth ed. Piscataway, N.J.: Wiley-IEEE Press, 2007.
• [19] P. Misra and P. Enge, “Global Positioning System—Signals, Measurements, and Performance”, Second Edition, Ganga-Jamuna Press, 2006.
• [20] B. Ristic, S. Arulampalam and N. Gordon, “Beyond the Kalman Filter—Particle Filters for Tracking Applications,” Artech House Inc., 2004.

The present invention aims to solve the above mentioned problems. In particular, the present invention aims to provide a method and apparatus that provides robust ranging at high accuracy and high processing rate.

SUMMARY OF THE INVENTION

The present invention is recited in the independent claims. Optional features are recited in the sub-claims.

In an embodiment, the present invention enables to exploit the carrier phase information directly, i.e. to resolve the ambiguity problem directly and independently for each signal source in a tracking module. To this end, the noise-free part of the receive signal is modelled as an exact function of the propagation delay. In particular, the dependency between the carrier phase and the delay parameter is taken into account in an explicit way. Additionally, a statistical model for the temporal evolution of the delay parameter is used to align a tracking algorithm to the grid of possible delay solutions. Over subsequent observation blocks then a long integration ambiguity histogram (LIAH) is constructed with the likelihood function evaluated on the ambiguity grid. This allows to resolve the ambiguity issue and to output a precise unbiased measurement of the delay parameter. The potential of this approach can be demonstrated within a Global Navigation Satellite System (GNSS) scenario where each millisecond it becomes possible to measure the range between a fast moving satellite and a GPS receiver with millimeter accuracy. The present invention is not limited to a GPS application and may also be implemented, for example, in FMCW radar positioning.

In another embodiment, the present invention demonstrates that it is possible to combine carrier and baseband information at the physical layer in order to obtain an ultra-precise and unambiguous time-delay estimate with a single narrowband transmit signal. To this end, a precise observation model is derived and side-information about the temporal evolution of the delay process is taken into account. By applying a Bayesian error bound, the present invention explores the fundamental performance limit for delay estimation with the established models. In order to make this theoretic accuracy accessible at moderate complexity, a two-step approximation of the optimum tracking algorithm is formulated. The presented estimation theoretic analysis and the Monte-Carlo simulations of the developed algorithm according to the present invention show that it is possible to increase the ranging and synchronization accuracy of wireless systems by a factor of 100.

BRIEF DESCRIPTION OF THE DRAWINGS

Embodiments of the invention are now described with reference to the accompanying drawings in which:

FIG. 1 is an illustration of a convention positioning system;

FIG. 2 is a flow chart illustrating the steps performed by a tracking module with LIAH in accordance with an embodiment of the present invention;

FIG. 3 is a flow chart illustrating the steps performed by a tracking module with DLL/PLL in accordance with an embodiment of the present invention;

FIG. 4 a RF front-end for use with an embodiment of the present invention;

FIG. 5 illustrates a maximum likelihood (ML) function generated in a method in accordance with an embodiment of the present invention;

FIG. 6 illustrates a two-dimensional GNSS scenario used for the simulation of a method in accordance with an embodiment of the present invention;

FIG. 7 illustrates the range estimation result obtained in a method in accordance with an embodiment of the present invention;

FIG. 8 illustrates the range estimation result obtained in a conventional DLL/PLL approach.

FIG. 9 illustrates Joint Carrier-Baseband Ranging and SNR Gain in accordance with an embodiment of the present invention;

FIG. 10 illustrates the normalized expected value of the likelihood function for baseband-oriented ranging in accordance with an embodiment of the present invention;

FIG. 11 illustrates the normalized expected value of the likelihood function for carrier-baseband ranging;

FIG. 12 illustrates a traffic control scenario in accordance with an embodiment of the present invention;

FIGS. 13a-d respectively illustrate range estimation results obtained in a method in accordance with an embodiment of the present invention

DESCRIPTION OF EMBODIMENTS OF THE INVENTION

An embodiment is described which allows a receive device to measure precisely and unambiguously the time-delay of a radio wave signal with known structure s(t) and carrier frequency ω_c. The procedure is performed in the digital domain, i.e. the receive device preprocesses the analog receive sensor signal by a radio front-end consisting of amplifiers, filters, two mixers and two samplers.

Brief Description of an Embodiment

Signal Model

The two sampling devices operate at a frequency

$f_s=\frac{1}{T_c}..$ Within a time period of T=NT_s, which is referred to as the duration of one processing block, a buffering device collects N=T/T_s samples from each digital output of the radio front-end. Therefore in the k-th processing block the digital data

$\begin{array}{c}{y}_k=\left[\begin{array}{c}{y}_{k,1}\\ ⋮\\ y_{k,N}\end{array}\right]\\ ={\gamma }_k\left[\begin{array}{c}{x}_{k,i}\left({\theta }_k\right)\\ ⋮\\ x_{k,N}\left({\theta }_k\right)\end{array}\right]+\left[\begin{array}{c}{n}_{k,1}\\ ⋮\\ n_{k,N}\end{array}\right]\\ ={\gamma }_k{x}_k\left({\theta }_k\right)+n_k\end{array}$ =γ_kx_k(θ_k)+n_k

is available, whereγ_k,n=γ_kx_k,n(θ_k)+n_k,n are the two samples at each sampling instance, in the following referred to as snapshot. Each snapshot contains the sampled version of the received radio wavex_k,n(θ_k)=T(τ_k)d_n(ν_k)s_k,n(τ_k,ν_k),scaled by the signal attenuation γ_k, with the matrix

$T\left({\tau }_k\right)=\lfloor \begin{array}{cc}\cos \left({\omega }_c{\tau }_k\right)& \sin \left({\omega }_c{\tau }_k\right)\\ -\sin \left({\omega }_c{\tau }_k\right)& \cos \left({\omega }_c{\tau }_k\right)\end{array}\rfloor ,$ due to the propagation delay τ_k and the carrier frequency ω_c of the radio wave from the transmitter to the receiver and a vector

$d_n\left(v_k\right)=\left[\begin{array}{c}\cos \left({\omega }_c{v}_k\left(n-1\right)T_s\right)\\ -\sin \left({\omega }_c{v}_k\left(n-1\right)T_s\right)\end{array}\right],$ due to the rate of change ν_k of the propagation delay (relative velocity between transmitter and receiver). The signal s^{sk,n(τk,νk)} is the sampled version of the transmitted baseband wave s(t). The snapshot component n_k,n is independent additive random noise with covarianceE[n_k,nn_k,n^T]=I₂.

For brevity of notation, the signal parameters (delay and velocity) are summarized in vector notationθ_k=[τ_kν_k]^T.

The temporal evolution of the two parameters (delay and velocity) over subsequent processing blocks follows an autoregressive model of first order

${\theta }_{k+1}=\left[\begin{array}{c}{\tau }_{k+1}\\ v_{k+1}\end{array}\right]=F{\theta }_k+z_k,$ where the process matrix is

$F=\left[\begin{array}{cc}1& T\\ 0& 1\end{array}\right]$ and z_k is additive process noise with the covariance matrix

$E\left[z_k{z}_k^T\right]={\sigma }_z^2\left[\begin{array}{cc}\frac{T^2}{3}& \frac{T^2}{2}\\ \frac{T^2}{2}& T\end{array}\right].$ Precise Time-Delay Measurement with LIAH

The aim of the receiver is to output accurate measurements (estimates) of the signal parameters {circumflex over (γ)}_k, {circumflex over (τ)}_k and {circumflex over (ν)}_k for a large number K of subsequent processing blocks. As accurate estimates of {circumflex over (τ)}_k exhibit a strong ambiguity the procedure contains a tracking module (ambiguity alignment) and an ambiguity resolution module (LIAH). The two modules perform individual steps and exchange information within each processing block.

Ambiguity Alignment Filter

The alignment filter provides the estimates {circumflex over (γ)}_k and {circumflex over (ν)}_k and produces an alignment center estimate {circumflex over (τ)}_A,k, which allows to deduce the grid of possible ambiguous solutions for {circumflex over (τ)}_k. In each block a signal generator produces replicas x_k(θ_k¹), x_k(θ_k²), . . . , x_k (θ_k^j) of the transmit signal for J different versions of the parameters θ_k. The signal generator forwards the replicas to J calculation units. The signal generator also generates a replica x_k({tilde over (θ)}_k) and forwards it to the ambiguity alignment unit. The j-th unit performs the calculation

${\hat{\gamma }}_k^j=\frac{y_k^T{x}_k\left({\theta }_k^j\right)}{x_k^T\left({\theta }_k^j\right)x_k\left({\theta }_k^j\right)}\mathrm{and}$ $f\left(y_k;x_k\left({\theta }_k^j\right)\right)=\frac{1}{2}{\left(y_k-{\hat{\gamma }}_k^j{x}_k\left({\theta }_k^j\right)\right)}^T\left(y_k-{\hat{\gamma }}_k^j{x}_k\left({\theta }_k^j\right)\right)=f_k^j.$

The calculation result ƒ_k^j is forwarded to the alignment unit. The alignment unit holds J weighting _k^j factors. With the calculations ƒ_k^j the weightings are updated_k^j=_k−1^jexp(−ƒ_k^j)and subsequently normalized such that

$\sum _{j=1}^j{w}_k^j=1.$

The tracking estimates are calculated

${\hat{\tau }}_{𝒜,k}=\sum _{j=1}^J{w}_k^j{\tau }_k^j$ ${\hat{v}}_k=\sum _{j=1}^J{w}_k^j{v}_k^j.$

The effective filter size

$J_{\mathrm{eff}}=\frac{1}{\sum _{j=1}^J{w}_k^j}$ is calculated. If J_eff<βJ, a resampling step is performed, which recalculates the parameter points θ_k¹, θ_k², . . . , θ_k^J by systematic resampling, adds roughening noiseθ_k^j=θ_k^j+r_k^j and sets

$w_k^j=\frac{1}{J}$ for all j=1, 2, . . . , J. Finally, the parameter points θ_k¹, θ_k², . . . , θ_k^J are propagated by the process modelθ_k+1^j=Fθ_k^j+z_k^j for all j=1, 2, . . . , J. A prediction for the next block is produced by

${\stackrel{\_}{\theta }}_{k+1}=\sum _{j=1}^J{w}_k^j{\theta }_{k+1}^j.$

The estimates {circumflex over (τ)}_A,k and {circumflex over (ν)}_k are forwarded to the LIAH module. The attenuation estimate is calculated

${\hat{\gamma }}_k=\frac{y_k^T{x}_k\left({\tilde{\theta }}_k\right)}{x_k^T\left({\tilde{\theta }}_k\right)x_k\left({\tilde{\theta }}_k\right)}.$

Initialization: For the initialization of the ambiguity alignment filter, i.e. for the first processing block k=0 the parameter points θ₀¹, θ₀², . . . , θ₀^J are initialized byτ₀^j˜[μ_init,τ−0.5Δ,μ_init,τ+0.5Δ]ν₀^j˜(μ_init,ν,σ_init,ν²),where μ_init,τ and μ_init,ν are the estimates of an acquisition algorithm and σ_init,ν is the standard deviation of the algorithm with respect to the velocity. Δ is proportional to half the carrier wave-length

$\Delta =\frac{\pi }{{\omega }_c}.$ The weighting factors are initialized

${\omega }_0^j=\frac{1}{J}$ for all j=1, 2, . . . , J and the prediction for the first block is produced by

${\tilde{\theta }}_1=\frac{1}{J}\sum _{j=1}^J{\theta }_0^j.$

Shift Event: If the LIAH filter triggers the shift event and provides a shift l_shift, all parameter versions are shifted such thatτ_k^j←τ_k^j+l_shiftΔ.LIAH Filter

The long integration ambiguity histogram (LIAH) filter provides an likelihood histogram p_k on the grid of possible ambiguous solutions for {circumflex over (τ)}_k. In each block a signal generator produces replicas of the transmit signal

$x_k\left({\theta }_k^1\right),x_k\left({\theta }_k^2\right),\dots ,x_k\left({\theta }_k^A\right)$ with A (A is an odd number) different versions of the parameters

${\theta }_k^a={\left[{\alpha }_{k,a}\Delta +{\hat{\tau }}_{𝒜,k}{\hat{v}}_k\right]}^T$ where α_k∈^A×1 is a vector containing the number of intervals Δ between the ambiguity grid point a and the alignment center point {circumflex over (τ)}_A,k. The signal generator forwards the replicas to A calculation units. The α-th unit performs the calculation

${\hat{\gamma }}_k^a=\frac{y_k^T{x}_k\left({\theta }_k^a\right)}{x_k^T\left({\theta }_k^a\right)x_k\left({\theta }_k^a\right)}$ $\mathrm{and}$ $\begin{array}{c}f\left(y_k;x_k\left({\theta }_k^a\right)\right)=\frac{1}{2}{\left(y_k-{\hat{\gamma }}_k^a{x}_k\left({\theta }_k^a\right)\right)}^T\left(y_k-{\hat{\gamma }}_k^a{x}_k\left({\theta }_k^a\right)\right)\\ =f_k^a.\end{array}$

The calculation result ƒ_k^α is forwarded to the LIAH unit. The LIAH unit holds A histogram values p_k^α. With the calculations ƒ_k^α the histogram values are updated

$p_k^a=p_{k-1}^a\exp \left(-f_k^a\right)$ and subsequently normalized such that

$\sum _{j=1}^A{p}_k^a=1.$

The delay estimate {circumflex over (τ)}_k is determined by

${\hat{\tau }}_k=a_{k,a}*\Delta +{\hat{\tau }}_{A,k},\mathrm{where}$ $a^{*}=\arg \underset{a\in \left\{1,\dots ,A\right\}}{\max }{p}_k^a.$

The LIAH unit holds A counters c_k^α. If a histogram value exceeds p_k^α>ρ, where

$\frac{1}{A}<\rho <1,c_k^a$ is incremented by one. If c_k^α>C the filter triggers a shift event, sets all counters c_k^α to zero, resets the histogram

$p_k^a=\frac{1}{A}$ and feedbacks the mitt l_shift (number of intervals Δ between old and new alignment center) to the alignment filter.

Select Checked Ambiguities: The LIAH unit checks a set of A ambiguities within the delay interval[{circumflex over (τ)}_A,k−A_maxΔ,{circumflex over (τ)}_A,k+A_maxΔ].

Therefore, the checked positions on the ambiguity grid are

${\alpha }_{k,a}=\mathrm{round}\left(\left(-1+\left(a-1\right)\frac{2}{A-1}\right)A_{\max }\right).$

Every time the shift event occurs the interval width is updated by

$A_{\max }=\max \left(\chi ,\frac{A-1}{2}\right)$ $\mathrm{with}$ $\chi =\underset{u,v\in \left\{1,\dots ,A\right\}}{\max }\left|{\alpha }_{k,u}-{\alpha }_{k,v}\right|$ under the constraint |u−ν|=1.

Initialize Checked Ambiguities: In the beginning of the procedure, i.e. k=0 the initialization is

$A_{\max }=\max \left(\lceil \frac{\partial {\sigma }_{\mathrm{init},\tau }}{\Delta }\rceil ,\frac{A-1}{2}\right)$ with ϵ being a design parameter σ_init,τ and the standard deviation of the acquisition algorithm with respect to the delay parameter.Integration into a Conventional Receiver

A possible approach to integrate the presented technique into a conventional receiver is to use a conventional DLL/PLL receiver as alignment filter. The DLL/PLL provides a carrier-phase measurement {circumflex over (ψ)}k and a coarse time-delay estimate {circumflex over (τ)}_k from which an alignment center point can de deduced by solving

${\hat{i}}_k=\arg \underset{l\in Z}{\min }{\left(\frac{{\hat{\phi }}_k+l\pi }{{\omega }_c}-{\hat{\tau }}_k\right)}^2$ and setting

${\hat{\tau }}_{A,k}=\frac{{\hat{\phi }}_k+{\tilde{l}}_k\pi }{{\omega }_c}.$

Detailed Description of an Embodiment

Observation Model

Consider a scenario with one radio transmitter and a receiver. The transmitter emits an electromagnetic wave of known periodic structurex′(t)=s′(t)cos(ω_ct),where s′(t)∈ is a periodic baseband signal and ω_c is the carrier frequency. The power of the transmitted signal x′(t) is assumed to be normalized∫_−∞^∞|X′(ω)|²dω=1,where |X′(ω)|² is the Fourier transform of the autocorrelation function of x′(t). The signal at the receive sensor

$\begin{array}{c}{y}^{\prime }\left(t\right)={\gamma }^{\prime }\left(t\right)x^{\prime }\left(t-\tau \left(t\right)\right)+n^{\prime }\left(t\right)\\ ={\gamma }^{\prime }\left(t\right)s^{\prime }\left(t-\tau \left(t\right)\right)\cos \left({\omega }_c\left(t-\tau \left(t\right)\right)\right)+n^{\prime }\left(t\right).\end{array}$ is characterized by a time-dependent propagation delay τ(t)∈ and an attenuation γ′(t)∈, while the additive random noise n′(t) is assumed to have flat power spectral density (PSD) Φ′(ω)=N₀. The receive signal ′(t) is demodulated with two orthogonal functionsd₁(t)=cos(ω_ct)d₂(t)=−sin(ω_ct)oscillating at carrier frequency (see FIG. 4). The signals in the two demodulation channels can, hence, be written as

$\begin{array}{cc}{y}_1^{\prime }\left(t\right)=& y^{\prime }\left(t\right)d_1\left(t\right)\\ =& \gamma \left(t\right)s^{\prime }\left(t-\tau \left(t\right)\right)(\cos \left({\omega }_c\tau \left(t\right)\right)+\\ & \cos \left(2{\omega }_{c}t-{\omega }_c\tau \left(t\right)\right))+n_1^{\prime }\left(t\right)\end{array}$ $\mathrm{and}$ $\begin{array}{cc}{y}_2^{\prime }\left(t\right)=& y^{\prime }\left(t\right)d_2\left(t\right)\\ =& \gamma \left(t\right)s^{\prime }\left(t-\tau \left(t\right)\right)(-\sin \left({\omega }_c\tau \left(t\right)\right)-\\ & \sin \left(2{\omega }_{c}t-{\omega }_c\tau \left(t\right)\right))+n_2^{\prime }\left(t\right),\end{array}$ $\mathrm{where}$ $n_i^{\prime }\left(t\right)=d_i\left(t\right)n^{\prime }\left(t\right),i\in \left\{1,2\right\},\mathrm{and}$ $\gamma \left(t\right)=\frac{{\gamma }^{\prime }\left(t\right)}{2}.$

Note that n₁′(t) and n′₂(t) are uncorrelated, i.e. E[n′₁(t)n′₂(t)]=0, ∀t. Moreover, the PSD of the additive white Gaussian noise components n′_i(t) is given by

${\Phi }_i\left(\omega \right)=\frac{N_0}{2},i\in \left\{1,2\right\}.$ The ideal low-pass filters h(t) of the two channels are assumed to have one-sided bandwidth B. The filtered analog signals can, thus, be written as

$\begin{array}{c}{y}_1\left(t\right)=y_1^{\prime }\left(t\right)*h\left(t\right)\\ =\gamma \left(t\right)s\left(t-\tau \left(t\right)\right)\cos \left({\omega }_c\tau \left(t\right)\right)+{\eta }_1\left(t\right)\\ y_2\left(t\right)=y_2^{\prime }\left(t\right)*h\left(t\right)\\ =-\gamma \left(t\right)s\left(t-\tau \left(t\right)\right)\sin \left({\omega }_c\tau \left(t\right)\right)+n_2\left(t\right),\end{array}$ $\mathrm{where}$ $\begin{array}{c}s\left(t\right)=s^{\prime }\left(t\right)*h\left(t\right)\\ n_i\left(t\right)=n_i^{\prime }\left(t\right)*h\left(t\right),i\in \left\{1,2\right\},\end{array}$ with * being the convolution operator. The signals of the two channels can be written in compact matrix-vector representation

$y\left(t\right)=\left[\begin{array}{c}{y}_1\left(t\right)\\ y_2\left(t\right)\end{array}\right]=\gamma \left(t\right)b\left(\tau \left(t\right)\right)s\left(t;\tau \left(t\right)\right)+n\left(t\right)$ $\mathrm{with}$ $\begin{array}{c}s\left(t;\tau \left(t\right)\right)=s\left(t-\tau \left(t\right)\right)\\ n\left(t\right)={\left[n_1\left(t\right)n_2\left(t\right)\right]}^T\end{array}$ $\mathrm{and}$ $b\left(\tau \left(t\right)\right)={\left[\cos \left({\omega }_c\tau \left(t\right)\right)-\sin \left({\omega }_c\tau \left(t\right)\right)\right]}^T.$

After filtering, the analog signals are sampled at a rate of

$f_s=\frac{1}{T_s}.$ In the following one observation block consists of N samples from each channel, i.e. _k∈^2N is the observation in block k. The time-delay process is assumed to be approximately linear within one blockτ(t)≈τ_k+ν_k(t−t_k),t∈[t_k;t_k+NT_s),where τ_k is the time-delay of the first sample in block k and ν_k is the relative velocity (normalized by speed of light c) between receiver and transmitter in block k. The signal strength γ(t) is assumed to be constant over one blockγ(t)=γ_k,t∈[t_k;t_k+NT_s).

For brevity of notation, the two parameter vectorsθ′_k=[τ_kν_kγ_k]^T θ_k=[τ_kν_k]^T are introduced. The n-th sample in the k-th block is given by

$\begin{array}{c}{y}_{k,n}=\left[\begin{array}{c}{y}_{k,n}^{\left(1\right)}\\ y_{k,n}^{\left(2\right)}\end{array}\right]=\left[\begin{array}{c}{x}_{k,n}^{\left(1\right)}\left({\theta }_k^{\prime }\right)\\ x_{k,n}^2\left({\theta }_k^{\prime }\right)\end{array}\right]+\left[\begin{array}{c}{n}_{k,n}^{\left(1\right)}\\ n_{k,n}^{\left(2\right)}\end{array}\right]\\ =x_{k,n}\left({\theta }_k^{\prime }\right)+n_{k,n}\\ ={\gamma }_k{b}_n\left({\theta }_k\right)s_{k,n}\left({\theta }_k\right)+n_{k,n}\end{array}$ $\mathrm{with}$ $\begin{array}{c}{s}_{k,n}\left({\theta }_k\right)=s\left(t_k+\left(n-1\right)T_s-{\tau }_{k,n}\left({\theta }_k\right)\right)\\ n_{k,n}=\left[\begin{array}{c}{n}_1\left(t_k+\left(n-1\right)T_s\right)\\ n_2\left(t_k+\left(n-1\right)T_s\right)\end{array}\right]\end{array}$ $\mathrm{and}$ $b_n\left({\theta }_k\right)={\left[\cos \left({\omega }_c{\tau }_{k,n}\left({\theta }_k\right)\right)-\sin \left({\omega }_c{\tau }_{k,n}\left({\theta }_k\right)\right)\right]}^T,$ where τ_k,n is the delay of the n-th sample in the k-th blockτ_k,n(θ_k)=τ_k+ν_k(n−1)T_s.

The vector b_n(θ_k) can be decomposedb_n(θ_k)=T(τ_k)d_n(ν_k)into a matrix

$T\left({\tau }_k\right)=\left[\begin{array}{cc}\cos \left({\omega }_c{\tau }_k\right)& \sin \left({\omega }_c{\tau }_k\right)\\ -\sin \left({\omega }_c{\tau }_k\right)& \cos \left({\omega }_c{\tau }_k\right)\end{array}\right],$ which only depends on the delay parameter τ_k and a vector

$d_n\left(v_k\right)=\left[\begin{array}{c}\cos \left({\omega }_c{v}_k\left(n-1\right)T_s\right)\\ -\sin \left({\omega }_c{v}_k\left(n-1\right)T_s\right)\end{array}\right],$ which only depends on the relative velocity ν_k. The receive signal can, thus, be modeled asx_k,n(θ′_k)=γ_kT(τ_k)d_n(ν_k)s_k,n(θ_k).

The vector of one observation block _k is defined as

$\begin{array}{c}{y}_k=\left[\begin{array}{c}{y}_{k,1}\\ ⋮\\ y_{k,N}\end{array}\right]=\left[\begin{array}{c}{x}_{k,1}\left({\theta }_k^{\prime }\right)\\ ⋮\\ x_{k,N}\left({\theta }_k^{\prime }\right)\end{array}\right]+\left[\begin{array}{c}{n}_{k,1}\\ ⋮\\ n_{k,N}\end{array}\right]\\ =x_k\left({\theta }_k^{\prime }\right)+n_k.\end{array}$

The noise n_k is assumed to be uncorrelated. Hence, the noise covariance matrix is given asR=E[n_kn_k^T]=BN₀I_2N,where I_2N is the identity matrix of dimension 2N.ML Estimation with a Single Block

The receiver is interested in estimating the parameters of the receive signal in order to gain information about the propagation channel between transmitter and receiver. Using only one observation block for estimation, the maximum likelihood (ML) estimator is the best unbiased estimator. The ML estimator is attained by solving the optimization

${\hat{\theta }}_k^{\prime }=\arg \underset{{\theta }^{\prime }\in {\Theta }^{\prime }}{\max }{f}_{\mathrm{ML}}^{\prime }\left(y_k;{\theta }_k^{\prime }\right),\mathrm{where}$ $f_{\mathrm{ML}}^{\prime }\left(y_k;{\theta }_k^{\prime }\right)=p_y\left(y_k|{\theta }_k^{\prime }\right)$ $\mathrm{and}$ $p_y\left(y_k|{\theta }_k^{\prime }\right)=\frac{e^{\left(-\frac{1}{2{\mathrm{BN}}_0}{\left(y_k-x_k\left({\theta }_k^{\prime }\right)\right)}^T\left(y_k-x_k\left({\theta }_k^{\prime }\right)\right)\right)}}{{\left(2\pi {\mathrm{BN}}_0\right)}^N}.$

The ML estimate for the signal γ_k strength can be computed in closed form as a function of θ_k and is given by

${\hat{\gamma }}_k\left({\theta }_k\right)=\frac{{\sum }_{n=1}^N{y}_{k,n}^{T}T\left({\tau }_k\right)d_n\left(v_k\right)s_{k,n}\left({\theta }_k\right)}{{\sum }_{n=1}^N{\left(s_{k,n}\left({\theta }_k\right)\right)}^2}.$ Substituting γ_k in the ML function results in a compact version

$f_{\mathrm{ML}}\left(y_k;{\theta }_k\right)=\frac{{\left({\sum }_{n=1}^N{y}_{k,n}^{T}T\left({\tau }_k\right)d_n\left(v_k\right)s_{k,n}\left({\theta }_k\right)\right)}^2}{{\sum }_{n=1}^N{\left(s_{k,n}\left({\theta }_k\right)\right)}^2}$ and the problem can be reformulated

${\hat{\theta }}_k=\arg \underset{\theta \in \Theta }{\max }{f}_{\mathrm{ML}}\left(y_k;{\theta }_k\right).$

Note that only the maximization with respect to τ_k and ν_k is required, while the ML estimate for γ_k can be computed in closed form with the solution {circumflex over (θ)}_k. In FIG. 5, the normalized noise-free ML function ƒ_ML(x_k({tilde over (θ)}′_k); θ_k) with {tilde over (θ)}′_k=[0 0 1]^T is plotted. A GPS-signal (C/A L1, Sat. 1)

$s^{\prime }\left(t\right)=\sum _{m=-\infty }^{\infty }{\left[b\right]}_{\mathrm{mod}\left(m,M\right)}ℊ\left(t-{\mathrm{mT}}_C\right)$ with block length N=2046 and chip duration Ct=977.53 ns is used, where b∈{−1, +1} is a sequence of M=1023 binary symbols, each of duration T_C, mod^(⋅) is the modulo operator and g(t) is a bandlimited transmit pulse. The carrier frequency is given as ƒ_c=1575.42 MHz. The one-sided bandwidth of the ideal low-pass filter at the receiver is equal to B=T_C⁻¹=1.023 MHz: The sampling frequency is chosen according to the sampling theorem ƒ_s=2B=2.046 MHz. The ML function is plotted in a range from 0 m to 0.4 m in τ_k-direction. There is no clear global maximum within this range, however, there are many local maxima. The distance between two neighboring maxima is half the wavelength

$\Delta =\frac{z}{2f_c}=0.0951$ m as the sign or the signal amplitude is not known at the receiver.

As the height of the local maxima decays slowly in the direction of τ_k, the multi-modal shape of the ML function makes estimation with one observation block impossible. The idea of the following sections is to resolve this ambiguity issue with the help of an likelihood histogram which is constructed over a long integration time.

Nearly Constant Velocity Model

In order to realize a long integration time within a dynamic scenario, the temporal evolution of the channel parameters has to be modeled precisely. Here an autoregressive model of first order is used

${\theta }_{k+1}=\left[\begin{array}{c}{\tau }_{k+1}\\ v_{k+1}\end{array}\right]=F{\theta }_k+{\omega }_k,$ where the matrix F∈^2×2 is the process matrix and _k is additive process noise with the covariance matrix

$E\left[{\omega }_k{\omega }_k^T\right]=Q\in {ℝ}^{2\times 2}.$

This simple model turns out to be quiet accurate for practical GNSS scenarios. A meaningful assumption for practical scenarios is that the first derivative {dot over (τ)}(t) of the continuous time-delay process τ(t) given in (19) is nearly constant over the duration of one block. Consequently, higher order derivatives are almost equal to zero and the second order derivative {umlaut over (τ)}(t) can be modeled as a zero mean white noise process {umlaut over (τ)}(t)=(t) with

$E\left[\omega \left(t\right)\omega \left(t^{\prime }\right)\right]={\sigma }_w^2\delta \left(t-t^{\prime }\right).$

The process matrix F is then given as

$F=\left[\begin{array}{cc}1& T\\ 0& 1\end{array}\right].$ where T=t_k+1−t_k is the duration of one block, i.e. T=NT_s and the covariance matrix Q as

$Q={\sigma }_{\omega }^2\left[\begin{array}{cc}\frac{T^3}{3}& \frac{T^2}{2}\\ \frac{T^2}{2}& T\end{array}\right].$ Precise Delay Estimation with Ambiguity Resolution

Apart from the observation model and the process model, the prior knowledge from the acquisition algorithm is assumed to be Gaussian, i.e. τ₁˜N(μ_init,τ,σ_init,τ²) and ν₁˜N(μ_init,ν,σ_init,τ²). Combining all available information, it is possible to estimate and track the time-delay process τ(t) with very high accuracy. The proposed low complexity estimation process consists of two steps for each block. The first step estimates and tracks one arbitrary time-delay ambiguity {circumflex over (τ)}_A,k (i.e., a first candidate for a correct propagation delay) with a particle filter (PF). Moreover, an estimate for the relative velocity {circumflex over (ν)}_k is provided. The second step exploits the structural information of the likelihood function by updating a likelihood histogram formed on a subset of points on the ambiguity grid A_k. Based on this long integration ambiguity histogram (LIAH) the algorithm finally decides for the most probable time-delay solution {circumflex over (τ)}_k from a plurality of candidates for a correct propagation delay.

Ambiguity Grid Alignment with a Particle Filter

The optimal estimator for the considered estimation problem is the conditional mean estimator (CME). Since the observation model shows severe non-linearities, the CME cannot be stated in closed form. Hence, suboptimal approaches need to be used. An estimation method, which approximates the CME and is able to handle strong non-linearities, is particle filtering. Note that the PF is identical to the CME only for an infinite number of particles. However, a large number of particles results in high computational complexity. In order to guarantee a correct and precise delay-estimation with a small number of particles, step 1 only focuses on finding one arbitrary ambiguity. Therefore, the particles are initializedτ₁^j˜[μ_init,τ−0.5Δ,μ_init,τ+0.5Δ]ν₁^j˜N(μ_init,ν,σ_init,ν²)for j=1, . . . , J, where J is the number of used particles. Note that the time-delay particles are initialized uniformly in the range of one ambiguity. With this initialization, it is possible to estimate and track one ambiguity with high accuracy. The weight _k^j of particle j is updated_k^j∝_k−1^jp(_k|θ_k^j,{circumflex over (γ)}_k(θ_k^j))exploiting the observation _k of block k. Note that the weights are initialized uniformly, i.e.

${\omega }_0^j=\frac{1}{j}$ and are normalized in every step such that Σ_j=1^J_k^j=1. The tracking estimates are

${\hat{\tau }}_{v^{A,k}}=\sum _{j=1}^J{\omega }_k^j{\tau }_k^j$ ${\hat{v}}_k=\sum _{j=1}^J{\omega }_k^j{v}_k^j.$

In every block the effective sample size

$J_{\mathrm{eff}}=\frac{1}{{\sum }_{j=1}^J{\omega }_k^j}$ is computed. If J_eff<0.5*J, a resampling and roughening step is needed to guarantee the stability of the PF. Here systematic resampling is used. The process model is used to update the particlesθ_k+1^j=Fθ_k^j+_k^j,j∈{1, . . . , J}. Ambiguity Resolution with LIAH

In the second step the structural information of the likelihood function is exploited. As the distance Δ between two neighboring ambiguities is known, the positions of all ambiguities, referred to as the ambiguity grid A_k, can be estimated from {circumflex over (τ)}_A,k. In the following the algorithm only considers the ambiguities within the interval[{circumflex over (τ)}_A,k−ϵσ_init,τ,{circumflex over (τ)}_A,k+ϵσ_init,τ].

Note that there is a trade-off between complexity and reliability which needs to be taken into account when choosing the design parameter ϵ. The probability that the true ambiguity, i.e. the time-delay τ_k is within the interval[μ_init,τ−ϵσ_init,τ,μ_init,τ+ϵσ_init,τ]can be computed with the initial knowledge of the acquisition. The core part of the presented delay estimation process consists of different stages. In every stage the likelihood of a fixed number A≥5 of ambiguities out of the grid A_k is tested. Without loss of generality, A is assumed to be odd. Within one stage the algorithm decides for one of the A ambiguities. The search is refined in the next stage. The A ambiguities for the first stage are

$a_k^i=\mathrm{round}\left(\left(-1+\left(i-1\right)\frac{2}{A-1}\right)a_{\max }\right)$ with round^(⋅) being the rounding operator, i∈{1, . . . , A} and

$a_{\max }=\max \left(\lceil \frac{\partial {\sigma }_{\mathrm{init},\tau }}{\Delta }\rceil ,\frac{A-1}{2}\right),$ where ┌ ┐ is the ceiling operator, i.e. the time-delay valuesα_kⁱΔ+{circumflex over (τ)}_A,k, i∈{1, . . . ,A}are checked. In order to decide for one of these ambiguities, a probability measure p_kⁱ is introduced and assigned to each of the tested ambiguities. The likelihood histogram p_kⁱ is initialized

$p_0^i=\frac{1}{A},i\in \left\{1,\dots ,A\right\}$ and updated in each blockp_kⁱ∝p_k−1ⁱp(_k|{circumflex over (θ)}_kⁱ,{circumflex over (γ)}_k({circumflex over (θ)}_kⁱ)).with{circumflex over (θ)}_kⁱ=[α_kⁱΔ+{circumflex over (τ)}_A,k{circumflex over (ν)}_k]^T,where {circumflex over (τ)}_A,k and {circumflex over (ν)}_k are the estimates of the aligning filter. The histogram is then normalized such that Σ_i=1^Ap_kⁱ=1. The delay parameter can be determined with the histogram

${\hat{\tau }}_k=a_k^i\Delta +{\hat{\tau }}_{𝒜,k},$ where i*=arg max_{i∈{1, . . . ,A}}p_kⁱ. Clearly, the true ambiguity needs not to be among the tested ambiguities in the first stage if the initial search interval is wider than A ambiguities. In order to find the true ambiguity, the search needs to be refined. Thus, a counter c_kⁱ, i=1, . . . , A, is introduced and initialized with zero, c_nⁱ=0,∀i. Every time p_kⁱ>ρ, where

$\frac{1}{A}<\rho <1,$ c_kⁱ by one. If a block is reached where c_kⁱ exceeds the design parameter C∈, the algorithm decides for ambiguity α_kⁱ and refines the search on the ambiguity grid. The particles τ_kⁱ of the grid aligning filter are shifted

${\tau }_k^j\leftarrow {\tau }_k^j+a_k^i\Delta ,$ i.e. also the estimate of the ambiguity τ_A,k is shifted

${\hat{\tau }}_{𝒜,k}\leftarrow {\hat{\tau }}_{𝒜,k}+a_k^i\Delta ,$ and α_max is updated as

$a_{\max }=\max \left(\chi ,\frac{A-1}{2}\right)$

with

$\chi =\underset{u,v\in \left\{1,\dots ,A\right\}}{\max }a^u-a^v$ under the constraint |u−ν|=1. The interval which is considered in the following blocks is a [{circumflex over (τ)}_A,k−α_maxϕ, {circumflex over (τ)}_A,k+α_maxΔ]. The tested ambiguities are computed as stated above. Note that if a

$a_{\max }=\frac{\stackrel{⋓}{A}-1}{2}$ no further refinement is necessary since all ambiguities in the considered interval are tested.Simulations

For the simulations a simple two-dimensional GNSS scenario depicted in FIG. 6 is considered. The receiver is positioned on a circle with radius R_E and center M. The receiver Rx is assumed to be static. The transmitter Tx moves on a circular orbit around the center M. If the transmitter is in the zenith, the distance between transmitter and receiver is h. The range r(t) between transmitter and receiver for this scenario depends on the angle α(t)=≮(R_x, M, T_x) which is given as

$\alpha \left(t\right)={\alpha }_0-2\pi \frac{t}{T_0},$ where T₀ is the circulation time of the transmitter Tx and α₀=α(0). Applying the law of cosine results in

$r\left(t\right)=\sqrt{R_E^2+R^2-2R_{E}R\cos \left(\alpha \left(t\right)\right)}$ for the range between transmitter and receiver, where R=R_E+h. The velocity is given by

$\stackrel{.}{\tau }\left(t\right)=\frac{R_{E}R\sin \left(\alpha \left(t\right)\right)\stackrel{.}{\alpha }\left(t\right)}{r\left(t\right)}\mathrm{with}\stackrel{.}{\alpha }\left(t\right)=-\frac{2\pi }{T_b}.$ For the simulations, the geometry is chosen according to a GPS scenario. R_E=6371*10³ m is equal to the radius of the earth. T₀=11 h 58 min and h=20200*10³ m are chosen according to the satellites of GPS

${\alpha }_0=\frac{\pi }{4}.$ The parameter σ² of the nearly constant velocity model is determined with a least squares approach. For the scenario this results in σ²=2.6279*10⁻¹⁴ The GPS-signal s′(t), the bandwidth B, the carrier frequency ƒ_c and the sampling frequency are chosen as described above. The signal strength is assumed to be 55 dB-Hz. The initial uncertainty of the acquisition is

${\sigma }_{\mathrm{init},\tau }=75m,{\sigma }_{\mathrm{init},v}=50\frac{m}{s}.$

The design parameters for the algorithm are ϵ=3.5, A=9, J=100 and C=10. In FIG. 7, the absolute value of the bias

${\mathrm{BIAS}}_k=E\left[{\tau }_k-{\hat{\tau }}_k\right],$ the mean error (MSE)

${\mathrm{MSE}}_k=E\left[{\left({\tau }_k-{\hat{\tau }}_k\right)}^2\right]$ and the variance

${\mathrm{VAR}}_k={\mathrm{MSE}}_k-{\mathrm{BIAS}}_k^2$ are measured for the estimation via LIAH with 250 realizations. It is observed, that the RMSE decreases to millimeter level within 1000 observation blocks

$\left(\sqrt{{\mathrm{MSE}}_{1000}}\approx {10}^{-3}m\right).$

Apart from that, the estimation result is unbiased, i.e.

${\mathrm{BIAS}}_k^2⪡{\mathrm{MSE}}_k$ for k sufficiently large. As a reference, the range estimation result attained with a standard DLL/PLL approach, which is described and implemented in [16], is plotted in FIG. 8. The DLL/PLL damping ratio is 0.7 and the DLL bandwidth is 2 Hz, respectively. The PLL bandwidth is 25 Hz and the DLL correlator spacing is 0.5 Ct. 2500 realizations are used to measure the MSE, bias and variance. It can be observed that it is possible to estimate the range with the accuracy of a few meters with this standard approach. Interestingly, the mean square error (MSE)

${\mathrm{MSE}}_k={\mathrm{BIAS}}_k^2+{\mathrm{VAR}}_k$ of this prior art method is dominated by the bias, i.e. a systematic estimation error. Ignoring this error and considering the variance, it is possible to estimate the range on meter-level

$\left(\sqrt{{\mathrm{VAR}}_{1000}}\approx 2.0m\right)$ with this conventional approach.

Detailed Description of a Further Embodiment

In a further embodiment of the invention, the ranging problem is addressed by means of a frequency-modulated continuous-wave (FMCW) radar example with a single transmit and a single receive antenna. A complex-valued representation of a wireless transmit signal x(t)∈C can in general be factorized

$x\left(t\right)=x_c\left(t\right)x_b\left(t\right),$ where the carrier signal x_c(t)∈C is a complex sinusoid

$x_c\left(t\right)=e^{j{\omega }_{c}t},$ oscillating at carrier frequency ω_C∈R. The signal x_b(t)∈C denotes an application-dependent modulating baseband component. For FMCW radar the baseband signal is given by a complex-valued chirp

$x_b\left(t\right)=e^{j\mu t-2t},0\le t<T_o,$ where μ denotes the rate of change in rad per s² and T_o the duration of the chirp signal. The bandwidth of the baseband signal is

$B=\frac{\mu T_o}{4\pi }.$

The receiver observes a delayed and attenuated version of the transmit signal

$z\left(t\right)=\gamma x\left(t-\tau \right)+\eta z\left(t\right),$ distorted by an additive white noise process ηz(t)∈C. The attenuation parameter γ∈R models the effect of the path loss and the antenna gain while τ∈R represents the time-delay due to the propagation of the electro-magnetic radio wave between the transmitter and the receiver. For FMCW applications the receiver usually performs compression in the analog domain by multiplying the receive signal with the complex-conjugate transmit signal

$\begin{array}{cc}\begin{array}{cc}\mathrm{ }& y\left(t\right)=z\left(t\right)x^{*}\left(t\right)\\ f_s=& =\gamma e^{j{\omega }_c\left(t-\tau \right)}{e}^{j\mu \left(t-\tau -\frac{{\tau }_o}{2}\right)\left(t-\tau \right)}{e}^{-j{\omega }_{c}t}{e}^{-j\mu \left(t-\frac{T_0}{2}\right)t}+{\eta }_y\left(t\right)\\ \mathrm{ }& =\gamma e^{-j{\omega }_c\tau }{e}^{j\left(-2\mu \mathrm{\tau t}+\mathrm{\mu \tau }\frac{T_o}{2}+{\mathrm{\mu \tau }}^2\right)}+{\eta }_y\left(t\right).\end{array}& \left(6\right)\end{array}$

The compressed signal is low-pass filtered to a bandwidth of B_s>μτ/π and sampled at a rate

$Fs=\frac{1}{T_s}=2B$ such that we finally obtain N digital receive samples

$y=\gamma e^{-j{\omega }_c\tau }b\left(\tau \right)+\eta$ with the vector entries

$b_n\left(T\right)=e^{j\left(-2\mathrm{\mu \tau }\left(n-1\right)T_s+\mathrm{\mu \tau }\frac{T_0}{2}+{\mathrm{\mu \tau }}^2\right)}$ ${\eta }_n=\eta \left(\left(n-1\right)T_s\right).$ wherein η(t) is the low-pass filtered version of η_y(t) such that the digital noise samples η exhibit a white temporal covariance matrix

$E_{\eta }\left[{\mathrm{\eta \eta }}^H\right]=I.$

The ranging problem deals with the question how to extract a high-resolution estimate {umlaut over (τ)}() for the time-delay parameter r given the receive vector y.

Estimation Theoretic Performance

The achievable mean squared error (MSE) with an unbiased estimator is lower bounded by the Cramér-Rao lower bound (CRLB).

$\begin{array}{c}\mathrm{MSE}\left(\tau \right)=E_{y\text{:}\tau }\left[{\left(\hat{\tau }\left(y\right)-\tau \right)}^2\right]\\ \ge \frac{1}{F\left(\tau \right)}\text{,}\end{array}$ where the Fisher information with respect to the delay parameter is defined

$F\left(\tau \right)={\int }_Y{\left(\frac{\partial \ln p\left(y;\tau \right)}{\partial \tau }\right)}^2\mathrm{d\gamma }.$

The optimum estimator which asymptotically achieves the theoretic limit is the maximum-likelihood estimator (MLE)

$\hat{\tau }\left(\gamma \right)=\arg \underset{\tau \in R}{\max }\ln p\left(\gamma ;\tau \right).$

Based on the above three formulas, it becomes clear that the performance and complexity of the ranging problem depends on the statistical receive model p(y;τ).

Conventional Baseband-Oriented Ranging:

The conventional approach to the ranging problem is based on the simplification that the effect of the propagation delay τ onto the carrier component x_c(t) can be modeled as a constant independent phase shift ϕ, i.e.,

$\begin{array}{c}\gamma =\gamma e^{-j{\omega }_c\tau }b\left(\tau \right)+\eta \\ =\gamma e^{j\varphi }b\left(\tau \right)+\eta ,\end{array}$ such that the parametric probability density function of the received signal is modeled as

$p\left(y;\tau \right)\sim \exp \left(-{y-\gamma e^{j\phi }b\left(\tau \right)}_2^2\right).$

The Fisher information with respect to the delay parameter is given by

$\begin{array}{c}{F}_b\left(\tau \right)=2{\gamma }^2\mathrm{Re}\left\{{\left(\frac{\partial b\left(\tau \right)}{\partial \tau }\right)}^H\frac{\partial b\left(\tau \right)}{\partial \tau }\right\}\\ =2{\gamma }^2\sum _{n=1}^N{\left(\frac{\partial b_n\left(\tau \right)}{\partial \tau }\right)}^2\frac{\partial b_n\left(\tau \right)}{\partial \tau }\\ =2{\gamma }^2\sum _{n=1}^N{\left(-2\mu \left(n-1\right)T_s+\mu \frac{T_o}{2}+2\mathrm{\mu \tau }\right)}^2,\end{array}$ where it has been used that

$\frac{\partial b_n\left(\tau \right)}{\partial \tau }=j\left(-2\mu \left(n-1\right)T_s+\mu \frac{T_o}{2}+2\mathrm{\mu \tau }\right)·e^{j\left(-2\mathrm{\mu \tau }\left(n-1\right)T_s+\mathrm{\mu \tau }\frac{T_0}{2}+{\mathrm{\mu \tau }}^2\right)}.$

Joint Carrier-Baseband Ranging:

In the case where the exact model and the parametric probability density function (the carrier-oriented model)

$p\left(y;\tau \right)\sim \exp \left(-{y-\gamma e^{-j{\omega }_c\tau }b\left(\tau \right)}_2^2\right)$ is used. Then the time-delay Fisher information is

$\begin{array}{cc}\begin{array}{c}{F}_c\left(\tau \right)=2{\gamma }^2\mathrm{Re}\left\{{\left(\frac{\partial e^{-j{\omega }_c\tau }b\left(\tau \right)}{\partial \tau }\right)}^H\frac{\partial e^{-j{\omega }_c\tau }b\left(\tau \right)}{\partial \tau }\right\}\\ =2{\gamma }^2\sum _{n=1}^N{\left(-2\mu \left(n-1\right)T_s+\mu \frac{T_o}{2}+2\mathrm{\mu \tau }-{\omega }_c\right)}^2.\end{array}& \left(19\right)\end{array}$

In FIG. 9 the information ratio

$\chi \left(\tau \right)=\frac{F_c\left(\tau \right)}{F_b\left(\tau \right)}$ between both Fisher information measures is depicted for a system with carrier frequency f_c=24 GHz and with ƒ_c=77 GHz. The values T_o=10 μs, τ=40 ns are set and χ(t) calculated as a function of the bandwidth B in MHz. It can be observed that switching from the model applied in the conventional approach to the carrier-oriented model is associated with a significant Fisher information gain. Measured in equivalent signal-to-noise ratio (SNR) the carrier-oriented model provides an information gain of 40 dB for a system with ƒ_c=24 GHz and B=200 MHz which, due to the relation between the MSE and the Fisher information, is equivalent to diminishing the root-mean squared error √{square root over (MSE(τ))} by a factor of 100. This approximately corresponds to the ratio

$\frac{f_c}{B}.$ Consequently, the performance gain which is associated with the configuration ƒ_c=77 GHz is higher than for the system with ƒ_c=24 GHz.Ranging Complexity

While the carrier-oriented model seems highly attractive with respect to its estimation theoretic performance, it has a serious drawback when considering the complexity associated with practical estimation algorithms. In order to envision this aspect, the maximum-likelihood estimator for the model applied in the conventional approach is considered. Given a specific ranging parameter τ the complex attenuation can be substituted

$\gamma e^{j\varphi }=\frac{b^H\left(\tau \right)\gamma }{b^H\left(\tau \right)b\left(\tau \right)}=\frac{b^H\left(\tau \right)y}{N},$ such that the time-delay maximum-likelihood estimate is found by solving (maximization task)

$\begin{array}{c}\stackrel{.}{\tau }=\arg \underset{\tau \in R}{\max }\frac{|b^H\left(\tau \right)y{|}^2}{b^H\left(\tau \right)b\left(\tau \right)}\\ =\arg \underset{\tau \in R}{\max }{f}_b\left(y;\tau \right).\end{array}$

Depicted in FIG. 10 is the normalized expected value of the likelihood function

$\begin{array}{c}\frac{E_{y\text{:}{\tau }_0}\left[f_c\left(y;\tau \right)\right]}{N}=\frac{|b^H\left(\tau \right)b\left({\tau }_0\right){|}^2}{N^2}\\ ={\stackrel{\_}{f}}_b\left(\tau ;{\tau }_0\right)\end{array}$ under the true parameter τ₀ for the exemplary system with ƒ_c=24 GHz and B=200 MHz. It can be observed that the likelihood function decays smoothly and that the maximization task can be solved at low complexity as the objective function has a one pronounced global extremum. Using the carrier-oriented model we obtain the substitution

$f=\frac{\mathrm{Re}\left\{{\left(e^{-j{\omega }_c\tau }b\left(\tau \right)\right)}^{H}y\right\}}{N},$ for the real-valued attenuation γ, such that the likelihood function obtains the form

$f_c\left(y;\tau \right)=\frac{\mathrm{Re}{\left\{{\left(e^{-j{\omega }_c\tau }b\left(\tau \right)\right)}^{H}y\right\}}^2}{N}.$

In FIG. 11 the normalized expected likelihood function is depicted

$\begin{array}{cc}\begin{array}{c}\frac{E_{y\text{:}{\tau }_a}\left[f_c\left(y;\tau \right)\right]}{N}=\frac{\mathrm{Re}{\left\{{\left(e^{-j{\omega }_c\tau }b\left(\tau \right)\right)}^H\left(e^{-j{\omega }_c{\tau }_0}b\left({\tau }_0\right)\right)\right\}}^2}{N^2}\\ ={\stackrel{\_}{f}}_c\left(\tau ;{\tau }_0\right).\end{array}& \left(26\right)\end{array}$

It becomes obvious that the expected likelihood function has a significantly sharper extremum at τ₀ than the baseband-oriented counterpart normalized expected value of the likelihood function. However, maximization of the likelihood function ƒ_c(γ;τ) under an individual realization of the observation data γ is difficult as the likelihood function exhibits pronounced side lobes around the true parameter. This also causes the problem that for certain noise realizations a side lobe attains a higher likelihood value than the main lobe which is close to the true parameter, therefore causing a high estimation error. However, the discussed FMCW radar example shows that using the joint carrier-baseband model is in theory associated with a significant accuracy gain but might in practice require high computational complexity and high integration time in order to benefit from this estimation theoretic fact.

Example of Implementation

As will be described in more detail below, the present invention discloses a method wherein it is possible to derive a reliable low-complexity algorithm which enables the full ranging precision provided by the exact measurement (carrier-oriented) model of the line-of-sight propagation delay after a short transient phase to be obtained. According to this exemplary implementation of the present invention, a joint carrier-baseband time-delay estimation is formulated as a state-space estimation problem wherein access to side information about the temporal evolution of the propagation delay process is utilized and signal processing may be performed over subsequent observation blocks. Based on the state-space model, a theoretic Bayesian performance bound and a nonlinear tracking algorithm is formulated which allows to accurately align to the moving set of likelihood extrema. In addition, a LIAH is constructed by evaluating and accumulating the likelihood function in an efficient way on a subset of possible range estimation solutions. After a short transient phase the histogram allows to detect the most probable ambiguity out of the solution set and to output an ultra-precise and unambiguous delay estimate for every observation block. With the applications of FMCW radar and GPS satellite-based positioning the efficiency of this approach is demonstrated by the present invention and the possible performance level highlighted. For a traffic control scenario, a sub-millimeter accuracy with a 24 GHz FMCW radar (240 MHz bandwidth) is obtained by performing ranging with the proposed approach of the present invention. For a GPS setup with ƒ_c=1.57 GHz (6.138 MHz bandwidth) it is shown that LIAH allows unambiguous low SNR tracking of the range parameter with respect to a fast moving narrowband transmit satellite at a RMSE level of a few millimeters.

System Model

Analog Transmit and Receive Signal

A mobile sender x(t_x) and a mobile receive sensor z(t) both comprise individual clocks t_x and t. The transmitter emits an analog radio wavex(t_x)=e^jωctxx′_b(t_x)on the carrier frequency ƒ_c with a baseband component x′_b(t_x)∈. At the receive sensor z(t) the transmit signal arrives attenuated by γ(t)∈R and delayed by τ(t)∈R

$\begin{array}{c}z\left(t\right)=\gamma \left(t\right)x\left(t-\tau \left(t\right)\right)+{\eta }^{\prime }\left(t\right)\\ =\gamma \left(t\right)e^{j{\omega }_c\left(t-\tau \left(t\right)\right)}{x}_b^{\prime }\left(t-\tau \left(t\right)\right)+{\eta }^{\prime }\left(t\right).\end{array}$

After demodulation and low-pass filtering to a bandwidth of B_s, the analog output (receive signal) is obtained

$\begin{array}{c}y\left(t\right)=h\left(t\right)*\left(e^{-j{\omega }_{c}t}z\left(t\right)\right)\\ =\gamma \left(t\right)e^{-j{\omega }_c\tau \left(t\right)}{x}_b\left(t-\tau \left(t\right)\right)+\eta \left(t\right).\end{array}$

This receive signal is further processed in blocks. In the k-th block, where it is assumed that the rate of change of the delay process and the attenuation are approximately(k−1)T₀≤t<kT₀,constant, i.e.,τ(t)≈τ_k+ν_k(t−(k−1)T₀)γ(t)≈γ_k.Digital Receive Signal

Sampling the analog signal at a rate of ƒ=2B_s the k-th receive block is obtained

$\begin{array}{c}{y}_k={\gamma }_k{e}^{-j{\omega }_c{\tau }_k}v\left(v_k\right)\odot x\left({\tau }_k,v_k\right)+\eta \\ ={\gamma }_{k}s\left({\theta }_k\right)+\eta \end{array}$ $\mathrm{with}$ $\begin{array}{c}{v}_n\left(v_k\right)=e^{-j{\omega }_c{v}_k\left(n-1\right)T_s}\\ x_n\left({\tau }_k,v_k\right)=x_b\left(\left(1-v_k\right)\left(n-1\right)T_s-{\tau }_k\right)\\ {\eta }_a=\eta \left(\left(n-1\right)T_s\right)\end{array}$ $\mathrm{and}$ $s\left({\theta }_k\right)=e^{-j{\omega }_c{\tau }_k}v\left(v_k\right)\odot x\left({\tau }_k,v_k\right).$ where all parameters of the delay process have been summarized into a state-space vectorθ_k=[τ_kν_k]^T.

The measurement model is therefore fully characterized by the parametric probability density function

$p\left(y_k|{\gamma }_k,{\theta }_k\right)=\frac{1}{{\pi }^N}\exp \left(-||y_k-{\gamma }_{k}s\left({\theta }_k\right){||}_2^2\right).$

Substituting the attenuation parameter by its maximum-likelihood estimate

${\gamma }_k=\frac{\mathrm{Re}\left\{s^H\left({\theta }_k\right)y\right\}}{s^H\left({\theta }_k\right)s\left({\theta }_k\right)}.$ a compressed measurement model can be stated in the form

$\begin{array}{cc}\tilde{p}\left(y_k|{\theta }_k\right)=\frac{1}{{\pi }^N}\exp \left(-||y_k-\frac{\mathrm{Re}\left\{s^H\left({\theta }_k\right)y_k\right\}s\left({\theta }_k\right)}{s^H\left({\theta }_k\right)s\left({\theta }_k\right)}{||}_2^2\right).& \left(39\right)\end{array}$ State-Space Model

Due to the assumption that during one observation block the rate of change of the delay process ν_k stays nearly constant, additionally to the measurement equation a state-space model can be formulated

${\theta }_k=Q{\theta }_{k-1}+{\omega }_k$ $\mathrm{with}$ $Q=\left[\begin{array}{cc}1& T_o\\ 0& 1\end{array}\right]$ $\mathrm{and}$ $\begin{array}{c}{E}_{{\omega }_k}\left[{\omega }_k{\omega }_k^T\right]={\sigma }_{\omega }^2\left[\begin{array}{cc}\frac{T_o^3}{3}& \frac{T_o^2}{2}\\ \frac{T_o^2}{2}& T_o\end{array}\right]\\ =W,\end{array}$ where σ_w is an application-specific variable, such that the state-space model is fully characterized by the process model

$\begin{array}{cc}p\left({\theta }_k|{\theta }_{k-1}\right)=\frac{\exp \left(-\frac{1}{2}{\left({\theta }_k-Q{\theta }_{k-1}\right)}^T{W}^{-1}\left({\theta }_k-Q{\theta }_{k-1}\right)\right)}{2\pi \sqrt{\det W}}.& \left(43\right)\end{array}$

For the initialization of the state-space model it is assumed that prior knowledge

$\rho \left({\theta }_0\right)=\frac{\exp \left(-\frac{1}{2}{\left({\mu }_0-{\theta }_0\right)}^T{\Sigma }_0^{-1}\left({\mu }_0-{\theta }_0\right)\right)}{2\pi \sqrt{\det {\Sigma }_0}}$ is available. According to the process model, the marginal p(θ_k) is a Gaussian random variable with mean

$\begin{array}{c}{\mu }_k=Q{\mu }_{k-1}\\ =Q^k{\mu }_0\end{array}$ and covarianceΣ_k=^Q(Σ_k−1+W)^QT.Nonlinear Filtering Problem

The resulting signal processing task of deducing an accurate estimate for the delay process z(t) is identified as a nonlinear Bayesian tracking problem. The solution to such a problem is to combine the measurement model and the process (state-space) model in order to calculate the optimum filter (optimum estimate){circumflex over (θ)}_k(Y_k)∫_Θkθ_kp(θ_k|Y_k)dθ_k of the block-wise state parameter vector θ_k, where the receive matrix

$Y_k=\left[\begin{array}{cccc}{y}_1& y_2& \dots & y_k\end{array}\right]$ contains the receive signals of the k-th block and all k−1 blocks received in the past. The average tracking error in the k-th block is given by the MSE matrix

$R_k=E_{{\theta }_k,Y_k}\left[\left(\hat{\theta }\left(Y_k\right)-{\theta }_k\right){\left(\hat{\theta }\left(Y_k\right)-{\theta }_k\right)}^T\right].$ Tracking with a Low-Complexity Filter

In order to calculate the optimum estimate for the problem at hand it is required to find an appropriate representation of the a posteriori probability density

$p\left({\theta }_k|Y_k\right)=\frac{\tilde{p}\left(y_k|{\theta }_k\right)p\left({\theta }_k|Y_{k-1}\right)}{p\left(y_k|Y_{k-1}\right)}$ by using the measurement model while, with the process (state-space) model, the a priori probability density is given by

$p\left({\theta }_k|Y_{k-1}\right)={\int }_{{\Theta }_{k-1}}p\left({\theta }_k|{\theta }_{k-1}\right)p\left({\theta }_{k-1}|Y_{k-1}\right)d{\theta }_{k-1}.$

Achieving a performance close to the theoretic limit through a low-complexity approximation of the optimum estimate is challenging. This is due to the fact that for a problem with nonlinear measurement model no closed-form solutions for the required calculation of these aforementioned formulas exist.

Resolving the Ambiguity Problem

Additionally, here the probability density function of the measurement model has a highly multimodal shape within the delay parameter τ_k. In the direction of the velocity parameter ν_k, the probability density function decays smoothly. This problem is known as ambiguity problem and is visualized by the expected likelihood function ƒ_c(τ; τ0) depicted in FIG. 3. However, it is also observed that the baseband component x_b(t) produces a slight decay in the absolute values of the extrema of ⁻ƒ_c(τ; τ₀). This indicates that over time one of the modes of p(θ_k|Y_k) will become dominant such that precise and unambiguous estimation of the time-delay parameter is possible. Using a standard particle filter for tracking the a posterior density requires a extremely high number of particles in order to accurately approximate its multimodal shape. For the problem at hand this makes state estimation with conventional methods computationally prohibitively complex. Following the concepts developed in the first embodiment of the present invention, it is proposed to use a low-complexity approximation of the a posteriori probability density, where the knowledge about the multimodal structure of the measurement model is exploited. The approach according to this embodiment of the present invention comprises two main steps. First, a particle filter is used in order to precisely align to a single mode of the a posteriori probability density. Then, advantage is taken of the fact that, due to the measurement model, the distance of the multiple modes in the direction of the delay parameter τ_k are placed at a regular distance of

$\Delta =\frac{\lambda }{2c}=\frac{\pi }{{\omega }_c},$ where λ is the wave-length and the speed of light. Therefore, the values of the a posteriori probability density are tracked exclusively at points which are located lΔ, l∈Z around the aligning mode. This is achieved by constructing a LIAH over subsequent blocks. Combining the information of the aligning filter and the ambiguity histogram enables deduction of an unambiguous and precise solution for the ranging problem.Ambiguity Aligning with a Gaussian Particle Filter

For the adaptive alignment and tracking of the ambiguity grid a Gaussian particle filter is used. To this end, a single mode of the a posteriori density in block k is approximated by a multivariate normal distribution N (μ_k,Σ_k) with mean μ_k and covariance Σ_k and the single-mode a priori density is approximated by N (^˜μ_k,^˜Σ_k) with mean ^˜μ_k and covariance ^˜Σ_k. In order to update these representations, in the k-th block J particles are drawn from the approximation of the a priori densityθ_k^j˜N({tilde over (μ)}_k,{tilde over (Σ)}_k)and the measurement model is used to calculate the weights_k^j˜{tilde over (p)}(_k|θ_k^j)with the received signal y_k. After normalization of the particle weights

${\omega }_k^j=\frac{{\stackrel{\_}{\omega }}_k^j}{{\sum }_{j=1}^j{\stackrel{\_}{\omega }}_k^j}$ the mean and covariance of the a posteriori density are estimated by

${\hat{\mu }}_k=\sum _{j=1}^J{\omega }_k^j{\theta }_k^j,{\hat{\Sigma }}_k=\sum _{j=1}^J{\omega }_k^j\left({\theta }_k^j-{\mu }_k\right){\left({\theta }_k^j-{\mu }_k\right)}^T.$

Then a correction is applied in order to form the a posteriori meanμ_k={circumflex over (μ)}_k+μ_Δ(l_k),whereμ_Δ(l)=[l_Δ0]^T is a shift of the delay estimate by l∈Z modes. The shift information l_k for each block is provided by the ambiguity resolution algorithm described in the next section. In a prediction step the a priori for the (k+1)-th block is formed by{tilde over (μ)}_k+1=Q_μk {tilde over (Σ)}_k+1=^QΣ_k^QT+Z,

In order to focus on a single mode in the direction of the delay parameter, in the first block k=1 the alignment filter is initialized with the a priori knowledge

${\tilde{\Sigma }}_l=\left[\begin{array}{cc}{\left({\rho }_P\frac{\Delta }{6}\right)}^2& 0\\ 0& {\delta }_{v,1}^2\end{array}\right],$ where ρ_P<1. As the measurement model is highly informative with respect to the delay parameter, the estimate of the a posterior covariance matrix degenerates within the first tracking blocks, i.e., the entries of the a posterior covariance matrix become very small. Further, with the update by the precise state-space model, the degeneracy propagates into the next block. This forms a severe problem for the convergence of the velocity parameter estimates ν_k which is essential for the correct behavior of the LIAH algorithm. In order to counteract the effect of fast degeneration of the covariance matrices, while preserving the ability of the algorithm to take advantage of the precise measurement and state-space model, the a posterior covariance estimate is smoothed by the a prior covariance matrix with a smoothing factor ρ_S<1. Further, in order to prevent divergence of the algorithm, the individual entries of the a prior covariance are restricted to the values of the initialization matrix

${\left[{\tilde{E}}_{k+1}\right]}_{11}\le {{\left[{\tilde{\Sigma }}_1\right]}_{11}\left[{\tilde{\Sigma }}_{k+1}\right]}_{22}\le {\left[{\tilde{\Sigma }}_1\right]}_{22}.$ while preserving positive semi-definiteness by

$\left|{\left[{\tilde{\Sigma }}_{k+1}\right]}_{12}\right|=\left|{\left[{\tilde{\Sigma }}_{k+1}\right]}_{21}\right|\le \sqrt{{{\left[{\tilde{\Sigma }}_{k+1}\right]}_{11}\left[{\tilde{\Sigma }}_{k+1}\right]}_{22}}.$

These steps result in a smooth convergence of the covariance matrices during the alignment process and an increasing tracking accuracy over time while avoiding degeneracy.

Ambiguity Resolution—Likelihood Histogram (LIAH)

While the alignment filter provides local information about the a posteriori density around the aligning mode, the ambiguity resolution algorithm aims at the likelihood characterization of a set A of surrounding modes. In order to select the best mode for the solution for the ranging problem a likelihood ambiguity histogram (LIAH) with |A|=A elements is adaptively updated in each block by calculating the histogram weights_k^α˜_k−1^α{tilde over (p)}(_k|θ_k^α)withθ_k^α={circumflex over (μ)}_k+μ_Δ(l_k^α)and performing the normalization

${\omega }_k^a=\frac{{\tilde{\omega }}_k^a}{{\sum }_{a=1}^A{\tilde{\omega }}_k^a}.$

The sequence l_a^k characterizes the set A of modes which are evaluated around the aligning mode μ_k. In order to reduce and fix the complexity of the algorithm the sequence l_a^k is designed such that the A bins of the histogram cover valuesl^α∈[−σ_H;σ_H],i.e., not all modes are checked when σH>A−12. Note, that in order to obtain a symmetric histogram A=|A| is restricted to odd numbers. After each update, the histogram is monitored by calculating its mean and standard deviation

$l_{\mathrm{LIAH},k}=\sum _{a=1}^A{\omega }_k^a{l}_k^a,{\sigma }_{\mathrm{LIAB},k}=\sqrt{\sum _{a=1}^A{{\omega }_k^a\left(l_k^a-l_{\mathrm{LIAH},k}\right)}^2}.$

In order to adjust the aligning algorithm to the dominant mode, a counter c^a* is incremented if the histogram uncertainty falls below

$\frac{{\sigma }_{\mathrm{LIAH},k}}{{\sigma }_H}<{\rho }_U,$ where the most likely mode is selected from the histogram by

$a^{*}=\arg \underset{a\in 𝒜}{\max }{\omega }_k^a.$

If the uncertainty is higher than the aforementioned threshold, all counters c^a are set to zero. In the case that c^a exceeds the threshold c^a>C, a histogram shift is triggered. Therefore, the ambiguity resolution algorithm informs the aligning filter to shift the a priori and a posteriori estimates by l_k=l^a*_k modes. In the case of no shift event l_k=0 is set. After each shift event the histogram size is adjusted

${\sigma }_H=\max \left(\frac{A-1}{2},{\rho }_H{\sigma }_H\right)$ with ρH<1, a new sequence l_a^k is formed according to the A bins of the histogram values and the histogram weights are initialized by

${\omega }_k^a=\frac{1}{A}.$

The computational complexity of the LIAH algorithm is restricted by checking only a subset of modes within a certain delay interval. In the initial phase of the tracking process, the correct ambiguity is not necessarily part of the ambiguity set A. The approach of diminishing the size of the histogram range σH assures that the set of tested ambiguities is refined with each shift event in a conservative way until A neighboring ambiguities, centered by the aligning mode, are checked. When the histogram size has reached the smallest range, i.e.,

${\sigma }_H=\frac{A-1}{2},$ with each subsequent shift event the shift threshold C is increased and the histogram initialized by_k^α=N(l=l_k^α;0,ρ_Wσ_H)and a subsequent normalization, in order to diminish the probability of further shift events and favor bins in the middle of the histogram.Applications and Simulation Results

In order to demonstrate the possible ranging performance under a precise observation model and exploitation of the available side information provided by the state-space model two technical applications are considered by the present invention i.e. ranging with a FMCW radar and synchronization with GPS. For both examples, a typical scenario and a strict reference is defined in order to measure the estimation error with Monte-Carlo simulations.

Ranging with a FMCW Radar

In addressing the initial radar problem, a velocity control scenario on a highway with a FMCW radar system is considered. A car is driving at a constant velocity ν along a road. A FMCW transceiver with a single transmit and a single collocated receive sensor is placed at a distance h from the middle of the traffic lane (see FIG. 12). The car reflects the FMCW radio wave with baseband structure

$x_b\left(t\right)=e^{j\frac{\mu }{2}\left(l-\frac{T_0}{2}\right)t},0\le t\le T_o.$

By measuring the propagation delay process τ(t) from the receive blocks γ_k the radar system is able to deduce the velocity v of the car and check for compliance with the traffic rules. In this example, the position of the car is parameterizedβ(t)=β₀+νt.

This allows calculation of an exact reference for the distance

$r\left(t\right)=\sqrt{h^2+{\left({\beta }_0+vt\right)}^2}$ and the relative velocity

$\stackrel{.}{r}\left(t\right)=\frac{\left({\beta }_0+\mathrm{vt}\right)v}{\sqrt{h^2+{\left({\beta }_0+\mathrm{vt}\right)}^2}},$ between car and the FMCW radar transceiver, such that it is possible to use

$\tau \left(t\right)=\frac{r\left(t\right)}{c}$ in order to evaluate the range estimation error and

$v\left(t\right)=\frac{\stackrel{.}{r}\left(t\right)}{c}$ for the estimation error of the relative velocity. For this example a car velocity of ν=120 km/h is assumed and a measurement distance of h=10 m. For the FMCW transceiver the carrier frequency

$\frac{{\omega }_c}{2\pi }=24$ GHz is set, the observation time for each block to To=1 μs and the bandwidth to B=240 MHz. Measuring the error of the autoregressive model (43) by

${\sigma }_{\omega }=\sqrt{\frac{1}{2\left(K_0+1\right)}\sum _{k=-K_0}^{K_0}{\left(\frac{\stackrel{.}{r}\left({\mathrm{kT}}_o\right)}{c}-\frac{\stackrel{.}{r}\left(\left(k-1\right)T_o\right)}{c}\right)}^2}$ for K₀=20000 blocks around the starting point

${\beta }_0=\beta \left(0\right)=-33.33m,$ a process error level of

${\sigma }_w=0.00264.$ is obtained. For each run the algorithm is initialized randomly by

${\tau }_0\sim 𝒩\left(\tau \left(0\right),\frac{1m}{c}\right).v_0\sim 𝒩\left(v\left(0\right),\frac{50m\text{/}s}{c}\right).$

For the alignment algorithm J=100 particles and ρ_P=¼. are used. The histogram has size A=11, a shift threshold of C=10, an initial width of σ_H=9*1m/Δc modes, a resizing factor of ρ_H=0.5 and a low-level initialization factor of ρ_W=0.33. For covariance smoothing ρ_S=0.99 is used. FIGS. 13a and 13b show the estimation performance of the proposed LIAH algorithm averaged over 100 noise realizations for each of the K=2000 blocks in a scenario with a receive power level of SNR=0 dB. It can be observed in FIG. 13a, that after approximately 0.5 ms the LIAH algorithm provides an unambiguous FMCW ranging solution with a root-mean squared error (RMSE) below 0.1 mm. The achievable performance with a classical baseband oriented approach is within a RMSE of a few cm. Therefore, LIAH allows to improve the FMCW ranging accuracy by approximately a factor of 100. FIG. 6 shows the accuracy of the relative velocity estimate, where after 2 ms the LIAH algorithm reaches a RMSE level of a few cms.

Satellite-based Positioning (GNSS)

As a second technical application, the present invention contemplates satellite-based position with the GPS system. For the example scenario a single in-view transmitting satellite moving on a circular orbit around the earth is considered. The height of the satellite is h and the radius of the earth is R_E (see FIG. 6).

The scenario is parameterized by a time-dependent angle

$\alpha \left(t\right)={\alpha }_0-\frac{2\pi }{T_R}t,$ where α₀ is the starting angle at t=0 and T_R is the duration of one full orbit. This simple model allows the generation of a rigorous referencer(t)=√{square root over (R_E²+R²−2R_ER cos(α(t)))}for the distance between satellite and receiver and its rate of change

$\stackrel{.}{r}\left(t\right)=-\frac{2R_{\varepsilon }R\pi \sin \left(\alpha \left(t\right)\right)}{T_R\sqrt{R_E^2+R^2-2R_{E}R\cos \left(\alpha \left(t\right)\right)}}.$

A reasonable setup according to a GPS application is R_E=6371·10³ m, T_R=11 h 58 min and h=20200·10³ m. The GPS baseband signal has the structure

$x_b\left(t\right)=\sum _{m=-\infty }^{m=\infty }{d}_mℊ\left(t-{\mathrm{mT}}_c\right),$ where d_m is a satellite-dependent chip sequence of length M=1023.

$T_c=\frac{1}{1.023\mathrm{MHz}}$ is the duration of each chip and g(t) is a bandlimited transmit pulse. Due to the signal structure the duration of one signal period is T_o=1 ms. A typical GPS carrier frequency is ω_c/2π=1.57 GHz and a bandwidth of B=6·1.023=6.138 MHz is assumed. Measuring the error of the auto-regressive model for K₀=20000 blocks around the starting angle

${\alpha }_0=\alpha \left(0\right)=\frac{45}{180}\pi ,$ a process error level of

${\sigma }_w=0.0003925.$ is obtained. For each run the algorithm is initialized randomly by

${\tau }_0\sim 𝒩\left(\tau \left(0\right),\frac{10m}{c}\right),v_0\sim 𝒩\left(v\left(0\right),\frac{50m\text{/}s}{c}\right).$

For the alignment algorithm J=25 particles and ρ_P=4 are used. The histogram has size A=11, a shift threshold of C=10, an initial width of

${\sigma }_H=9\frac{10m}{\Delta c}$ modes, a resizing factor of ρ_H=0.5 and a low-level initialization factor of ρ_W=0.33. For covariance smoothing ρ_S=0.9975 is used. FIGS. 13c and 13d visualize the GPS synchronization accuracy of the LIAH algorithm averaged over 100 noise realizations for each of the K=5000 blocks in a scenario with a receive power level of SNR=−12 dB. For the considered GPS setup FIG. 13c shows, that after approximately 2 s the LIAH algorithm provides an unambiguous synchronization solution with a root-mean squared error (RMSE) of around 1 mm. The achievable performance with a classical baseband oriented approach is within a RMSE of a few tenth cm. FIG. 13d shows the accuracy of the estimate of the relative velocity of the receiver and the GPS satellite. The LIAH algorithm achieves a velocity RMSE level of a few

$\frac{cm}{s}$ after approximately 5 s.

CONCLUSION

As described in respect of embodiments of the present invention, unbiased ranging with high precision with one single transmitter is possible in the tracking module of a receiver at moderate complexity. This is achieved by modeling the carrier phase as an exact function of the propagation delay parameter in the statistical model of the receive signal. The ambiguity issue is resolved by means of a tracking based alignment filter and a long integration histogram which assigns probabilities to each ambiguity. In a satellite-based synchronization and positioning application (GPS), the presented approach can outperform prior art ranging methods (DLL/PLL) with respect to the RMSE.

Claims

1. A method of determining a distance or location of a remote device or reflector, the method comprising: receiving a signal from a remote signal transmitter associated with or contained in or attached to the remote device;estimating and tracking a first propagation delay associated with the received signal, wherein the first propagation delay represents a first candidate for a correct propagation delay;deriving a relationship between the first candidate and one or more other candidates for a correct propagation delay from the received signal;determining a plurality of other candidates for a correct propagation delay based on said relationship;generating a likelihood histogram based on said candidates for a correct propagation delay;selecting a propagation delay from said candidates based on the likelihood histogram; anddetermining a distance or location of the remote device or reflector using the selected propagation delay.

2. The method of claim 1, wherein said relationship is a distance between neighboring candidates for a correct propagation delay in a grid of said candidates.

3. The method of claim 1, wherein said candidates correspond to local maxima in a likelihood function or distribution.

4. The method of claim 1, further comprising applying a maximum likelihood estimation to the received signal, thereby to determine said first propagation delay.

5. The method of claim 1, further comprising: tracking said candidates for a correct propagation delay; andupdating the likelihood histogram in a predetermined number of iterations or over a predetermined period of time.

6. The method of claim 1, further comprising: selecting the first propagation delay by selecting one of a plurality of candidates for a correct propagation delay using particle filtering.

7. The method of claim 6, further comprising uniformly initializing propagation delay particles within a predetermined range of said selected one of said candidates for a correct propagation delay.

8. The method of claim 1, further comprising: processing the received signal to generate a plurality of blocks of signal samples; andgenerating the likelihood histogram over subsequent blocks of signal samples.

9. The method of claim 1, wherein the first propagation delay is a time-delay propagation delay.

10. A signal and data processing device comprising: a receiver to receive a signal from a remote signal transmitter associated with or contained in or attached to a remote device;a signal processing module to convert the received signal into digital data; anda data processing module to: determine and track a first propagation delay associated with the received signal, wherein the first propagation delay represents a first candidate for a correct propagation delay;determine a plurality of other candidates for a correct propagation delay based on a relationship between the first candidate and one or more other candidates for a correct propagation delay;generate a likelihood histogram based on said candidates for a correct propagation delay;select a propagation delay from said candidates based on the likelihood histogram; anddetermine a distance or location of the remote device using the selected propagation delay.

11. The signal and data processing device of claim 10, comprising a DLL/PLL receiver to perform a propagation delay estimate to determine said first propagation delay.

12. The signal and data processing device of claim 10, wherein the first propagation delay is a time-delay propagation delay.