METHOD AND APPARATUS FOR PERFORMING JOINT CHANNEL AND TIME ESTIMATION IN A GNSS RECEIVER

Abstract

A receiver is provided for use with a global navigation satellite system (GNSS) comprising multiple satellites. The receiver comprises a receiver clock and at least one antenna for receiving multiple signals over multiple respective channels, each channel being defined by a transmitting satellite and a receiving antenna at opposing ends of the channel. The receiver further comprises at least one correlator for calculating cross-correlation functions between (i) the signals received over the multiple channels and (ii) reference versions of the navigation signals provided by the receiver. The receiver is configured to use the calculated cross-correlation functions to perform a joint estimate of (i) a clock bias of the receiver clock relative to the time reference maintained by the GNSS, and (ii) a composite channel comprising the combined contribution of the multiple channels as a function of time-delay.

Claims

1. A receiver for use with a global navigation satellite system (GNSS) comprising multiple satellites, wherein each satellite transmits a respective navigation signal containing a spreading code, and wherein the GNSS maintains a time reference which is encoded into the navigation signals, the receiver comprising: a receiver clock; at least one antenna for receiving multiple signals over multiple respective channels, each channel being defined by a transmitting satellite and a receiving antenna at opposing ends of the channel; and at least one correlator for calculating cross-correlation functions between (i) the signals received over the multiple channels and (ii) reference versions of the navigation signals provided by the receiver; wherein the receiver is configured to use the calculated cross-correlation functions to perform a joint estimate of (i) a clock bias of the receiver clock relative to the time reference maintained by the GNSS, and (ii) a composite channel comprising the combined contribution of the multiple channels as a function of time-delay.

2. The receiver of claim 1, wherein the signals received over the respective channels include contributions from one or more of the navigation signals.

3. The receiver of claim 2, wherein the receiver is configured to detect a given navigation signal included in a signal received over a given channel based on the calculated cross-correlation function between the signal received over the given channel and the reference version of the given navigation signal.

4. The receiver of claim 3, wherein the receiver is further configured to use the calculated cross-correlation function between the signal received over the given channel and the reference version of the given navigation signal to determine a code-phase alignment between the received given navigation signal and the provided reference version of the given navigation signal.

5. The receiver of claim 4, wherein the code-phase alignment is dependent upon: (i) the propagation time of the given navigation signal from the satellite to the receiver and (ii) the clock bias between the receiver clock and the time reference maintained by the GNSS.

6. The receiver of claim 4, wherein the calculated cross-correlation functions used for performing the joint estimate are centred on the determined code-phase alignment for each satellite.

7. The receiver of any of claim 2, wherein the contributions for the one or more navigation signals include multiple non-line-of-sight components.

8. The receiver of claim 1, wherein the time-delay is referenced to arrival of a line-of-sight component.

9. The receiver of claim 1, wherein the clock bias of the receiver clock is estimated as an offset from an approximate clock bias.

10. The receiver of claim 9, wherein the approximate clock bias is equal to an initial estimated clock bias having a typical accuracy (1−sigma) in the range 0.1-1 microseconds, or in the range 0.2-0.8 microseconds, or in the range 0.25-0.6 microseconds.

11. The receiver of claim 9, wherein the receiver is configured to estimate the approximate clock bias, prior to performing the joint estimation, by using the calculated cross-correlation functions to determine a maximum of a composite cross-correlation function for the composite channel.

12. The receiver of claim 11, wherein the composite cross-correlation function is defined by $\hat{\Lambda }\left({\Phi }_{0,n}\right)=\underset{j=1}{\overset{N}{.Math.}}{x}_j^{H}B\left({\Phi }_{0,n}-\Delta {\Phi }_j\right){\left(B^H\left({\Phi }_{0,n}-\Delta {\Phi }_j\right)B\left({\Phi }_{0,n}-\Delta {\Phi }_j\right)\right)}^{-1}{B}^H\left({\Phi }_{0,n}-\Delta {\Phi }_j\right)x_j.$

13. The receiver of claim 1, wherein performing the joint estimate of the clock bias includes searching a predefined range of possible clock bias values.

14. The receiver of claim 1, wherein the composite channel is specified by a set of coefficients, each coefficient corresponding to a particular time-delay, and representing the aggregate level of power transmission through the composite channel for that time-delay.

15. The receiver of claim 14, wherein the receiver includes at least one probability distribution function for each coefficient in the set of coefficients.

16. The receiver of claim 15, wherein the provided probability distribution functions are based on an estimate of the composite channel obtained using received signals from a previous epoch.

17. The receiver of claim 16, wherein the receiver is configured to update the probability distribution functions based on the estimated composite channel.

18. The receiver of claim 15, wherein receiver performs the joint estimate using a Monte Carlo search based on trial values of the composite channel, generated in accordance with the probability density functions, to find the trial values and clock bias that best fit the calculated cross-correlation functions.

19. The receiver of claim 1, wherein performing the joint estimate is defined by equation (3).

20. The receiver of any of claim 1, wherein the receiver includes multiple antennas, and the joint estimate takes advantage of spatial correlation between the signals received by different antennas.

21. The receiver of claim 20, wherein the performing the joint estimate is defined by equation (4).

22. The receiver of claim 1, wherein the estimated clock bias has an accuracy of less than 100 nanoseconds, or less than 80 nanoseconds, or in the range 25 to 75 nanoseconds, or of the order of 50 nanoseconds.

23. A method of operating a receiver for use with a global navigation satellite system (GNSS) comprising multiple satellites, wherein each satellite transmits a respective navigation signal containing a spreading code, and wherein the GNSS maintains a time reference which is encoded into the navigation signals, the receiver comprising a receiver clock and at least one antenna for receiving multiple signals over multiple respective channels, each channel being defined by a transmitting satellite and a receiving antenna at opposing ends of the channel, the method comprising: calculating cross-correlation functions between (i) the signals received over the multiple channels and (ii) reference versions of the navigation signals provided by the receiver; performing a joint estimate of (i) a clock bias of the receiver clock relative to the time reference maintained by the GNSS, and (ii) a composite channel comprising the combined contribution of the multiple channels as a function of time-delay, using the calculated cross-correlation functions.

24. The method of claim 23, wherein the receiver is located in an indoor environment.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0015] Various implementations of the disclosure will now be described in detail by way of example only with reference to the following drawings:

[0016] FIG. 1 is a schematic diagram of a receiver for receiving a GNSS navigation signal, according to one or more embodiments shown and described herein;

[0017] FIG. 2 is a plot showing an example of a distortion introduced into a correlation function by an indoor propagation channel for a given antenna and satellite;

[0018] FIG. 3 is a plot showing examples of distortions introduced into a correlation function by an indoor propagation channel for different antennas and a given satellite;

[0019] FIG. 4 shows examples of the power delay profile for simulated indoor conditions, according to one or more embodiments shown and described herein;

[0020] FIG. 4A shows a subset of the data of FIG. 4, showing the power delay profile for a single satellite, according to one or more embodiments shown and described herein; and

[0021] FIG. 4B show another subset of the data from FIG. 4, showing the power delay profile for a single, according to one or more embodiments shown and described herein.

DETAILED DESCRIPTION

[0022] FIG. 1 is a high-level schematic block diagram of a receiver 701 suitable for use with a GNSS for receiving navigation signals from respective satellites. The receiver 701 may be provided as a stand-alone unit, or may be incorporated into some larger device, for example a mobile (cellular) telephone, a cellular radio access node (e.g., a 5G small cell), a computer, an automobile or other form of vehicle, a hospital bed, an aircraft or ship, a freight container, and so on. In operation, receiver 701 includes at least one antenna 715 for receiving GNSS signals. The antenna 715 links to an RF conditioning, down-conversion and A/D conversion unit 720, which in turn passes the base-band complex envelope of the signal received to channel acquisition and tracking unit 725.

[0023] The receiver 701 also includes a code generator 712, which can generate (or store) the spreading codes (pseudo-random noise sequences) for the satellite constellation. The channel acquisition and tracking unit 725 is provided by the code generator 712 with the relevant spreading codes so that the incoming channel(s) can be acquired from the base-band complex envelope of the received signal. The channel acquisition unit 725 utilises the appropriate number of spreading codes to support the sequential or parallel operation of the channel acquisition unit as appropriate, e.g., in accordance with the number of receiver channels. In the illustrated embodiment, the channel acquisition and tracking unit 725 includes R correlators 728.

[0024] The receiver 701 is configured to perform a joint estimate of (i) the clock bias and (ii) a composite channel, as described in more detail below. This processing may be integrated into the channel acquisition unit 725, or may be performed in a separate unit within the receiver 701 (not shown in FIG. 1), for example, in a processor, a digital signal processor (DSP), or any other suitable component.

[0025] Some navigation systems may include local elements (sometimes referred to as pseudolites). These are local positioning systems, for example at airports, which may supplement the positioning signals from a satellite navigation system to provide additional location information. It will be appreciated that a receiver 701 as described herein is intended to receive and process any suitable navigation signal (whether from a satellite, pseudolite, or any other relevant source), and to perform the joint estimation processing accordingly.

[0026] FIG. 2 is a plot showing an example of the distortion introduced into the correlation function by the indoor propagation channel for a given antenna and satellite (as obtained with the channel model described in T. Jost et al., “A Wideband Satellit-to-Indoor Channel Model for Navigation Applications,” IEEE Transactions on Antennas and Propagation, vol. 62, no. 10, October 2014, pp. 5307-5320, and Report ITU-R P.2145-2, “Model Parameters for the Physical-Statistical Wideband Model in Recommendation ITU-R P.681,” September 2017). The reference correlation function (dashed line) is included for comparison, and has a peak aligned with zero. The observed correlation function (solid line) is distorted and slightly delayed (shifted to the right) due to the presence of (delayed) NLOS components, which also participate in the correlation function and hence bias the outcome (in contrast the reference correlation function can be considered as the correlation function that would be produced by receiving only the LOS component). Note in particular that the maximum of the distorted correlation function for the indoor channel in FIG. 2 is biased (delayed) with respect to the reference correlation function, which will then introduce a bias into any position and/or timing estimations.

[0027] FIG. 3 is a plot showing examples of distortion introduced in the correlation function by the indoor propagation channel for different antennas and a given satellite (again as obtained with the channel model in T. Jost et al., “A Wideband Satellit-to-Indoor Channel Model for Navigation Applications,” IEEE Transactions on Antennas and Propagation, vol. 62, no. 10, October 2014, pp. 5307-5320, and Report ITU-R P.2145-2, “Model Parameters for the Physical-Statistical Wideband Model in Recommendation ITU-R P.681,” September 2017). Note that the correlation functions are normalized such that the antenna having the highest correlation has a value of unity. (The different antennas typically have different locations and orientations, and hence experience a slightly different channel from each other). It can be seen that the antenna having the strongest correlation value also has the best shape in terms of a well-defined peak for the correlation function. In contrast, the correlation functions for the other antennas are generally somewhat distorted, having lower and less well-defined peaks (due to distortion from NLOS components and at least partial loss of the LOS component). Therefore any position and/or timing estimations based directly on these distorted correlation functions are susceptible to bias or error.

[0028] The approach described herein helps to overcome the bias/distortion such as illustrated in FIGS. 2 and 3. We consider a receiver having an array of N antennas (with N≥1) and receiving M navigation signals (with M≥1, all with equivalent correlation properties and transmitted in the same frequency band) from a given GNSS constellation (or any other satellite-based navigation system). The received signals from the M satellites are considered to be impacted by an indoor propagation channel. Based on this defined multiple-input multiple-output (MIMO) system, and considering a composite channel comprising the combined contribution of the M×N propagation channels of the MIMO system, a method is provided for the joint estimation of the composite channel and time synchronization (between the receiver clock and the time reference of the GNSS), based on the received satellite navigation signals.

[0029] The following pieces of information are considered to be known a priori: [0030] a) a baseline (phase offset) between each of the antennas and the phase centre of the receiver.
b) the position of the receiver;
c) a coarse clock bias for the receiver, ;
Regarding (a), this is a property of the receiver itself, and may be determined using known measuring equipment. Regarding (b) and (c), these may be estimated using conventional techniques, such as implemented in existing GNSS receivers. For example, a receiver may initially perform a conventional acquisition of the navigation signals to obtain this position and timing (clock bias) information before performing the additional processing described herein. Alternatively (or additionally), the position and/or timing information may be known from other sources e.g., the receiver may be located in a building whose position is already known, e.g., from mapping information, likewise timing information may be available from an external source (such as a telecommunications network or similar).

[0031] The composite channel is defined as the aggregate contribution of the M×N propagation channels—corresponding to the N received versions (one per antenna) of the M satellite signals—as a function of the time-delay with respect to the LOS signal (even though such a LOS signal may not be visible to the receiver). The composite channel gives a measure of the power delay profile for the resulting aggregated MIMO channel—i.e., the total aggregated power of all the received signal contributions as a function of the time-delay. Based on an estimation of the composite channel, an estimation of the clock bias of the receiver can also be derived—i.e., a joint estimation process is performed.

[0032] In general, each of the coefficients of the composite channel (one coefficient per time-delay, each coefficient reflecting the level of power transmission through the channel for that time-delay) will follow a probability density function (PDF) which is a priori unknown, and highly dependent on the propagation conditions of the receiver environment. As described below, the estimation of the composite channel is performed as a multi-hypothesis estimation problem in which different PDF assumptions for each coefficient of the composite channel are assessed.

[0033] The estimation of the composite channel is based on the composite cross-correlation function Λ(δt), which is dependent on the clock bias. Accordingly, the estimation of the composite cross-correlation function at different discrete clock bias values δt.sub.n is of interest. In practice, the composite cross-correlation function is also dependent on other parameters, like the clock drift {dot over (δ)}.sub.t and the phase center position of the receiver relative to the N antennas (among others).

[0034] Thus Δ(δt.sub.n)Δ(δt.sub.n,{dot over (δ)}t,p.sub.0.sup.T)=Λ(Φ.sub.0,n), where the receiver's state vector Φ.sub.0,n=[δt.sub.n,{dot over (δ)}t,p.sub.0.sup.T].sup.T and the sub-index 0 refers to the phase center of the array of antennas (i.e., the phase centre position of the receiver). The composite cross-correlation function for the receiver's state vector Φ.sub.0,n is derived, following J. A. Garcia-Molina and J. A. Fernandez-Rubio, “Exploiting Spatial Diversity in Low-cost SDR Platforms: the MIMO-GNSS approach,” 6th International Colloquium on Scientific and Fundamental Aspects of GNSS/Galileo, Valencia, Spain, October 2017, as:

[00001]$\begin{array}{cc}\hat{\Lambda }\left({\Phi }_{0,n}\right)=\underset{j=1}{\overset{N}{.Math.}}{x}_j^{H}B\left({\Phi }_{0,n}-\Delta {\Phi }_j\right){\left(B^H\left({\Phi }_{0,n}-\Delta {\Phi }_j\right)B\left({\Phi }_{0,n}-\Delta {\Phi }_j\right)\right)}^{-1}{B}^H\left({\Phi }_{0,n}-\Delta {\Phi }_j\right)x_j,& \left(\mathrm{equation}1\right)\end{array}$

where the following definitions apply: [0035] a) Φ.sub.0,n−ΔΦ.sub.j is equivalent to the state vector per antenna Φ.sub.j,n=[δt.sub.n,{dot over (δ)}t,p.sub.j.sup.T].sup.T=[δt.sub.n,{dot over (δ)}t,p.sub.0.sup.T].sup.T−[0,0,Δp.sub.j.sup.T].sup.T. Note that p.sub.j=p.sub.0−Δp.sub.j, with Δp.sub.j the vector between the receiver's phase center and the j-th antenna. In practice: Φ.sub.j,nΦ.sub.j,n(t)=[δt.sub.n(t),{dot over (δ)}t(t),p.sub.j.sup.T(t)].sup.T, [0036] b) x.sub.j=[x.sub.j(t.sub.0) . . . x(t.sub.K-1)].sup.T∈.sup.K×1 is the snapshot of K samples (with sampling period T.sub.δ=t.sub.K−t.sub.K-1) for the j-th antenna corresponding to the integration period being considered in the composite channel estimation (x.sub.j therefore represents the received signal at respective antenna j), [0037] c) B(Φ.sub.j)=[b.sub.1(Φ.sub.j) . . . b.sub.M(Φ.sub.j)]∈.sup.K×M is the basis function matrix for the state vector Φ.sub.j, which is composed by the basis function vectors for each satellite of the j-th antenna b.sub.i(Φ.sub.j)=[b.sub.i(t.sub.=0, (Φ.sub.j) . . . b.sub.i(t.sub.k-1,Φ.sub.j)].sup.T∈.sup.K×1 [0038] d) the basis function b.sub.i(t,Φ.sub.j) for the i-th satellite and the j-th antenna is defined at time instant t as b.sub.i(t,Φ.sub.j)=d.sub.i(t−τ.sub.i,j(t,Φ.sub.j))exp{j2πf.sub.i,j(t,Φ.sub.j)t}, with d.sub.i the complex baseband model of the signal transmitted by the i-th satellite.

[0039] The derivation is shown below:

[0040] A signal vector gathering all the signal snapshots observed by the N antennas is defined as:

y=S(Φ.sub.0)c+e,  (i)

where y=[x.sub.1.sup.T . . . x.sub.N.sup.T].sup.T∈.sup.nK×1, the vector c=[a.sub.1.sup.T . . . a.sub.N.sup.T]∈.sup.MN×1 gathers all the complex amplitudes for the MN paths in the MIMO-GNSS system, e=[n.sub.1.sup.T . . . n.sub.N.sup.T]∈.sup.MK×1 gathers the noise of all the antennas (with e˜CN(0,σ.sup.2I)), and the multi-antenna basis function matrix is defined based on the single-antenna basis functions B(Φ.sub.0−ΔΦ.sub.j) as:

[00002]$\begin{array}{cc}S\left({\Phi }_0\right)=\left[\begin{array}{cccc}B\left({\Phi }_0-\Delta {\Phi }_1\right)& 0& \mathrm{.Math.}& 0\\ 0& B\left({\Phi }_0-\Delta {\Phi }_2\right)& \mathrm{.Math.}& 0\\ \mathrm{.Math.}& \mathrm{.Math.}& \ddots & \mathrm{.Math.}\\ 0& 0& \mathrm{.Math.}& B\left({\Phi }_0-\Delta {\Phi }_N\right)\end{array}\right].& \left(\mathrm{ii}\right)\end{array}$

[0041] The cost function is defined as Λ.sub.1(Φ.sub.0,c)=∥y−S(Φ.sub.0)c∥.sup.2, such that the MIMO maximum likelihood estimate (MLE) of the receiver's state vector Φ.sub.0 is found by minimizing Λ.sub.1(Φ.sub.0, c):

[00003]$\begin{array}{cc}{\hat{\Phi }}_0,\hat{c}=\arg \underset{{\Phi }_0,c}{\min }{\Lambda }_1\left({\Phi }_0,c\right).& \left(\mathrm{iii}\right)\end{array}$

[0042] To find an efficient implementation of the MLE of the receiver's state vector, the cross-correlations are defined as {circumflex over (r)}.sub.yy=y.sup.Hy, {circumflex over (r)}.sub.ys(Φ.sub.0)=S.sup.H(Φ.sub.0)y, and {circumflex over (R)}.sub.ss(Φ.sub.0)=S.sup.H((Φ.sub.0)S(Φ.sub.0). Based on the previous definitions, and taking into account that the least-squares estimator of the complex amplitudes vector c is:

ĉ={circumflex over (R)}.sub.ss.sup.−1(Φ.sub.0){circumflex over (r)}.sub.ys(Φ.sub.0),  (iv)

it can be readily shown that the cost function can be re-written as:

Λ.sub.2(Φ.sub.0)={circumflex over (r)}.sub.yy−{circumflex over (r)}.sub.ys.sup.H(Φ.sub.0){circumflex over (R)}.sub.ss.sup.−1(Φ.sub.0){circumflex over (r)}.sub.ys(Φ.sub.0)  (v)

Therefore, the MLE of the receiver's state vector can be re-defined as:

[00004]$\begin{array}{cc}{\hat{\Phi }}_0=\arg \underset{{\Phi }_0}{\max }\left\{{\hat{r}}_{ys}^H\left({\Phi }_0\right){\hat{R}}_{ss}^{-1}\left({\Phi }_0\right){\hat{r}}_{ys}\left({\Phi }_0\right)\right\}.& \left(\mathrm{vi}\right)\end{array}$

To apply this MLE of the state vector of the receiver to a real-time receiver, it is convenient to reformulate the maximization problem in (vi) based on the cross-correlations for each antenna of the receiver: {circumflex over (r)}.sub.x.sub.j.sub.b(Φ.sub.j)={circumflex over (r)}.sub.x.sub.j.sub.b(Φ.sub.0−ΔΦ.sub.j)=B.sup.H(Φ.sub.0−ΔΦ.sub.j)x.sub.j, and {circumflex over (R)}.sub.bb(Φ.sub.j)={circumflex over (R)}.sub.bb(Φ.sub.0−ΔΦ.sub.j)=B.sup.H(Φ.sub.0−ΔΦ.sub.j)B(Φ.sub.0−ΔΦ.sub.j), resulting in

[00005]$\begin{array}{cc}{\hat{\Phi }}_0=\arg \underset{{\Phi }_0}{\max }\underset{j=1}{\overset{N}{.Math.}}{\hat{r}}_{x_{j}b}^H\left({\Phi }_0-\Delta {\Phi }_j\right){\hat{R}}_{bb}^{-1}\left({\Phi }_0-\Delta {\Phi }_j\right){\hat{r}}_{x_{j}b}\left({\Phi }_0-\Delta {\Phi }_j\right)& \left(\mathrm{vii}\right)\end{array}$

where the cost function is referred to as the “composite cross-correlation function” to be maximized w.r.t. the state vector Φ.sub.Q, which can be re-written as:

[00006]$\begin{array}{cc}{\Lambda }_3\left({\Phi }_0\right)=\underset{j=1}{\overset{N}{.Math.}}{x}_j^{H}B\left({\Phi }_0-\Delta {\Phi }_j\right){\left(B^H\left({\Phi }_0-\Delta {\Phi }_j\right)B\left({\Phi }_0-\Delta {\Phi }_j\right)\right)}^{-1}{B}^H\left({\Phi }_0-\Delta {\Phi }_j\right)x_j& \left(\mathrm{viii}\right)\end{array}$

(corresponding to Equation 1 above).

[0043] The coarse clock bias .sub.c can be estimated, as noted above, using any suitable existing technique, e.g., as implemented within current receivers. In order to help decrease the computational complexity of the joint composite channel and time estimation process, the clock bias error typically obtained by state-of-the-art techniques in indoor conditions (usually of the order of hundreds of ns) may be reduced further to help to constrain the search space in the joint estimation problem. For example, an intermediate clock bias estimation may be performed to refine the initial coarse estimation .sub.c by exploiting the composite cross-correlation function Λ(δt). The intermediate clock bias estimation .sub.int can then be derived as:

[00007]$\begin{array}{cc}{\mathrm{int}}_{}=\arg \underset{\delta t}{\max }\hat{\Lambda }\left(\delta t\right)& \left(\mathrm{equation}2\right)\\ s.t..Math.\delta t-c_{}.Math.<\alpha ,& \end{array}$

where α is a threshold dependent on the quality of the coarse clock bias estimation .sub.e at the current epoch/integration period.

[0044] The received signal at each antenna can be regarded as the sum of multiple separate components, each component corresponding to the transmitted navigation signal but subject to a relative delay and change in amplitude. The line-of-sight (LOS) may be regarded as having zero delay in this respect, and all other components as having a positive delay. In some situations, the LOS component may not be visible.

[0045] To formalize this, the composite channel vector h∈.sup.D×1 may contain the coefficients of the composite channel for D discrete time-delays of interest (uniformly distributed covering the channel delay spread), such that [h].sub.mh.sub.m∈ is the m-th coefficient of the composite channel vector; the coefficient h.sub.1 corresponding to the LOS contributions marks the actual clock bias of the receiver, i.e., relative to the time reference of the GNSS, (the actual clock bias is to be estimated by the procedure described herein). In effect therefore, for a given component, h.sub.m, the value of m specifies the time-delay (relative to the LOS component), while the value of h.sub.m determines the amplitude (or other measure of intensity) of the component having a time-delay m.

[0046] The coefficients of the components of the composite channel {[h].sub.m}.sub.m=1.sup.D are considered to follow an arbitrary set of PDFs {p([h].sub.m)}.sub.m=1.sup.D{p.sub.l,m}.sub.m=1.sup.D, with p.sub.l,m the PDF adopted for the l-th hypothesis and the m-th coefficient, and with a total of N.sub.l hypotheses for the PDF followed by each coefficient. We define the fine clock bias estimator .sub.f based on the intermediate clock bias .sub.int and the (positive or negative) integer offset β∈ (to be estimated) as .sub.f=.sub.int−βT.sub.δ, where T.sub.δ is the time-delay resolution of the composite channel vector. Based on this, the joint estimator of h and β for the l-th hypothesis can be defined as the solution to the constraint optimization problem, as per Equation (3):

[00008]$\begin{array}{cc}{\hat{h}}_l,{\hat{\beta }}_l=\arg \underset{h,\beta }{\min }{.Math.\hat{\Lambda }-\underset{m=1}{\overset{D}{.Math.}}{\left[h\right]}_m{.Math.s\left({\mathrm{int}}_{}-\left(\beta -m+1\right)T_{\delta }\right).Math.}^2.Math.}^2& \left(\mathrm{equation}3\right)\\ s.t.{\left[h\right]}_m\sim p_{l,m},{\forall }_m\mathrm{and}.Math.\beta .Math.<W,& \end{array}$

where the following definitions apply: [0047] a. ĥ.sub.1 is the composite channel vector estimated under the l-th hypothesis for the channel PDFs{p([h].sub.m)}.sub.m=1.sup.D{p.sub.l,m}.sub.m=1.sup.D, with ĥ.sub.l=[ĥ.sub.l,1 . . . ĥ.sub.l,D].sup.T∈.sup.D×1, [0048] b. {circumflex over (β)}.sub.l is the integer offset estimated under the l-th hypothesis, [0049] c. W∈ constrains the values evaluated for the integer offset β (in total, N.sub.w=2W−1 values evaluated) based on the expected error of the intermediate clock bias .sub.int, [0050] d. {circumflex over (Λ)}=[{circumflex over (Λ)}(δt.sub.1) . . . {circumflex over (Λ)}(δt.sub.Q)].sup.T∈.sup.Q×1 is the estimated composite cross-correlation vector for the MIMO system, which is uniformly sampling (with a time-delay resolution T.sub.δ) the composite cross-correlation function Λ(δt) at Q clock offset values, being centered at the intermediate clock bias estimation .sub.int, such that δt.sub.int|.sub.n=(Q+1)/2=.sub.int, [0051] e. s(.sub.int−(β−m+1)T.sub.δ)∈.sup.Q×1 is a vector containing the samples of the reference correlation function (for the GM signal being considered) centered at the clock bias .sub.int−(β−m+1)T.sub.δ.
Note that it is the value of β which is ultimately of interest to define the time at the receiver with reference to the time within the GNSS.

[0052] In high-level terms, equation (3) defines a sum of squared differences between (i) the observed signal values (represented by {circumflex over (Λ)}, as per (d) above), and (ii) the predicted values for this hypothesis, which are determined based on (a) the original navigation signal (represented by s, as per (e)) above, (b) the hypothesized pattern of reflections for the signal received at each antenna, according to the predicted channel ĥ.sub.l, as per (a) above, and the hypothesized time offset of the receiver clock, represented by {circumflex over (β)}.sub.1, as per (b) above. The better the hypothesized values of ĥ.sub.l and {circumflex over (β)}.sub.l match the received data, the smaller the sum of squared differences.

[0053] The optimization is subject to two conditions. Firstly, a limited range of values (±W) for {circumflex over (β)}.sub.l are defined, which reduces the search space about the intermediate clock bias. Note that this is primarily a computational issue—with more computing resources available a wider search range could be explored for {circumflex over (β)}.sub.1, to the extent that in some cases, the search might be performed using the coarse clock bias, without having to first derive an intermediate clock bias .sub.int. Secondly, the values of ĥ.sub.l are constrained by a probability density function, represented by p.sub.l,m, as discussed in more detail below. One way of looking at p.sub.l,m, is as a form of regularization, which ensures the different values of ĥ.sub.l are physically plausible (rather than allowing any arbitrary values of ĥ.sub.l); this also helps to reduce the risk of over-fitting.

[0054] The optimization problem set out in equation (3) can be solved, for example, using Monte Carlo methods, whereby the coefficients of the composite channel vector h are treated as random variables drawn from the corresponding arbitrary PDFs p.sub.l,m (with N.sub.c Monte Carlo runs assessed for each coefficient in each hypothesis). In practice, the actual distribution of the coefficients depends heavily on the indoor channel propagation conditions. Thus multiple PDF hypotheses for the coefficients of the composite channel are evaluated, and following this approach, the estimated composite channel vector ĥ is the one out of the N.sub.l estimated channel vectors ĥ.sub.l that results in the minimum mean square error between (i) the estimated composite cross-correlation and (ii) the reference model based on the composite channel vector being assessed.

[0055] Algorithm 1, as set out below, summarizes the composite channel and fine clock bias estimation (i.e., ĥ and .sub.f) based on equation (3) when considering N.sub.1 PDF hypotheses and N.sub.c Monte Carlo runs for each of the coefficients of the composite channel.

TABLE-US-00001 Algorithm 1 Joint composite channel and fine clock bias estimation Input:  [00009]${\left\{x_j\right\}}_{j=1}^N,{\left\{{\left\{b_i\left(t,{\Phi }_j\right)\right\}}_{j=1}^N\right\}}_{i=1}^M,{\widehat{\delta t}}_c,\alpha ,N_l,N_c,W,{\left\{{\left\{p_{l,m}\right\}}_{m=1}^D\right\}}_{l=1}^{N_l},Q,D,s,T_{\delta }$ Output: ĥ,    .sub.f   1: Initialize the minimum mean square error ξ.sub.min = ∞   2: Estimate the intermediate clock bias    .sub.int based on equation (2)   3: Compute {circumflex over (Λ)} = [{circumflex over (Λ)}(δt.sub.1) . . . {circumflex over (Λ)}(δt.sub.Q)].sup.T centered at    .sub.int based on     equation (1)   4: for l = 1 to N.sub.l do   5:  for c = 1 to N.sub.c do   6:   Draw the normalized channel coefficients [h.sub.l,c.sup.norm].sub.m ~ p.sub.l,m, ∀.sub.m   7:   h.sub.l,c = h.sub.l,c.sup.norm∥{circumflex over (Λ)}∥/∥Σ.sub.m=1.sup.D[h.sub.l,c.sup.norm].sub.m|s(   .sub.int + (m − 1)T.sub.δ)|.sup.2∥   8:   for β = −W + 1 to W − 1 do   9:    Compute ξ.sub.l,c,β = ∥{circumflex over (Λ)} − Σ.sub.m=1.sup.D[h.sub.l,c].sub.m|s(   .sub.int −       (β − m + 1)T.sub.δ)|.sup.2∥.sup.2  10:    if ξ.sub.l,c,β < ξ.sub.min then  11:     ξ.sub.min = ξ.sub.l,c,β  12:     ĥ = h.sub.l,c  13:     {circumflex over (β)} = β  14:    end if  15:   end for  16:  end for  17: end for  18:    .sub.f =    .sub.int − {circumflex over (β)}T.sub.δ

[0056] Algorithm 1A, as set out below, summarizes the composite channel and fine clock bias estimation based on equation (3) for the particular case in which N=1 (a single antenna is used by the receiver) and N.sub.l=1 (a single hypothesis for each coefficient is considered in the estimation process), such that p.sub.1,m is the PDF for the m-th coefficient, Φ.sub.1 is the state vector of the receiver's antenna, and x.sub.1 is the snapshot of K samples gathered by the receiver's antenna.

TABLE-US-00002 Algorithm 1A Joint composite channel and fine clock bias estimation for the particular case in which N = 1 and N.sub.l = 1. Input: x.sub.1, {b.sub.i(t, Φ.sub.1)}.sub.i=1.sup.M,    .sub.c, α, N.sub.c, W, {p.sub.1,m}.sub.m=1.sup.D, Q, D, s, T.sub.δ Output: ĥ,    .sub.f   1: Initialize the minimum mean square error ξ.sub.min = ∞   2: Estimate the intermediate clock bias    .sub.int based on equation (2)   3: Compute {circumflex over (Λ)} = [{circumflex over (Λ)}(δt.sub.1) . . . {circumflex over (Λ)}(δt.sub.Q)].sup.T centered at    .sub.int based on     equation (1)   4: for c = 1 to N.sub.c do   5:  Draw the normalized channel coefficients [h.sub.c.sup.norm].sub.m ~ p.sub.1,m, ∀.sub.m   6:  h.sub.c = h.sub.c.sup.norm∥{circumflex over (Λ)}∥/∥Σ.sub.m=1.sup.D[h.sub.c.sup.norm].sub.m|s(   .sub.int + (m − 1)T.sub.δ)|.sup.2∥   7:  for β = −W + 1 to W − 1 do   8:  [00010]$\mathrm{Compute}{\xi }_{c,\beta }={.Math.\hat{\Lambda }-{{\Sigma }_{m=1}^D\left[h_c\right]}_m|s\left({\widehat{\delta t}}_{\mathrm{int}}-\left(\beta -m+1\right)T_{\delta }\right){|}^2.Math.}^2$   9:   if ξ.sub.c,β < ξ.sub.min then  10:    ξ.sub.min = ξ.sub.c,β  11:    ĥ = h.sub.c  12:    {circumflex over (β)} = β  13:   end if  14:  end for  16: end for  17:    .sub.f =    .sub.int − {circumflex over (β)}T.sub.δ

[0057] In general, the PDFs {p.sub.m}.sub.m=1.sup.D can be set to any desired distribution. One option is to define a set of default PDFs representative of different types of environments, from light-indoor to deep-indoor conditions, as derived from theoretical models, ray-tracing models, and/or from real data for those environments—e.g., using a database of representative indoor conditions. For static receivers, another possibility is to estimate the PDFs {p.sub.m}.sub.m=1.sup.D for the D coefficients of the composite channel based on the estimated coefficients {[ĥ].sub.m}.sub.m=1.sup.D of the previous epochs available—i.e., building or refining the PDFs based on the coefficients estimated so far. For periods of time which are sufficiently long (e.g., several days) such that all possible geometries between the GNSS satellites and the receiver antennas are encountered, the resulting {{circumflex over (p)}.sub.m}.sub.m=1.sup.D should characterize the propagation conditions of the receiver for the indoor conditions in which the receiver is located. In this situation, {{circumflex over (p)}.sub.m}.sub.m=1.sup.D can be regarded as a signature of the environment surrounding the receiver. Such estimation may be performed in a sequential (iterative) way for each epoch, updating the estimated PDFs {{circumflex over (p)}.sub.m}.sub.m=1.sup.D of the previous epoch to provide the new composite channel vector estimation ĥ for the current epoch.

[0058] In an alternative approach, a joint estimator of h and β is defined by taking advantage of the spatial correlation of the NLOS multipath components between antennas. In other words, the NLOS components “seen” by the different antennas (located in slightly different positions) may have some level of mutual correlation. The level of this correlation depends on the distance between the antennas and the indoor propagation conditions.

[0059] The joint estimator can be defined as the solution to the constraint optimization problem

[00011]$\begin{array}{cc}{\hat{h}}_l,{\hat{\beta }}_l=\arg \underset{h,\beta }{\min }\gamma \left(h,\beta \right)& \left(\mathrm{equation}4\right)\\ s.t.{\left[h\right]}_m\sim p_{l,m},{\forall }_m\mathrm{and}.Math.\beta .Math.<W,& \end{array}$

where the following definitions apply: [0060] a)

[00012]$\gamma \left(h,\beta \right)=\frac{\frac{1}{N}\underset{j=1}{\overset{N}{.Math.}}{\hat{\lambda }}_j\left(h,\beta \right)}{{\left(\underset{j=1}{\overset{N}{.Math.}}{\hat{\lambda }}_j\left(h,\beta \right)\right)}^{1/N}},$

where {circumflex over (λ)}.sub.j(h,β) is the j-th eigenvalue associated to the spatial covariance matrix of the multipath mis-modelling errors, which is estimated as Ĉ.sub.ε=(h,β)=Ê(h,β)Ê(h,β).sup.T, [0061] b) Ê(h,β)=[{circumflex over (ε)}.sub.1(h,β) . . . {circumflex over (ε)}.sub.N(h,β)].sup.Tε.sup.Q×N is a matrix containing the vectors with the estimated multipath mis-modelling errors {circumflex over (ε)}.sub.j for each antenna, which are defined as:

[00013]$\begin{array}{cc}{\widehat{\mathrm{.Math.}}}_j\left(h,\beta \right)={\hat{\theta }}_j-{\kappa }_j\left(h\right)\underset{m=1}{\overset{D}{.Math.}}{\left[h\right]}_m{.Math.s\left({\mathrm{int}}_{}-\left(\beta -m+1\right)T_{\delta }\right).Math.}^2,& \left(\mathrm{equation}5\right)\end{array}$

[0062] with κ.sub.j(h)=∥{circumflex over (θ)}.sub.j∥/∥Σ.sub.m=1.sup.D[h].sub.m|s(.sub.int+(m−1)T.sub.δ)|.sup.2∥, [0063] c) {circumflex over (θ)}.sub.j=[{circumflex over (θ)}.sub.j(δt.sub.1) . . . {circumflex over (θ)}.sub.j(δt.sub.Q)].sup.T∈.sup.Q×1 is the estimated composite cross-correlation vector for the j-th antenna, which is uniformly sampling (with a time-delay resolution T.sub.δ) the composite cross-correlation function θ.sub.j(δt) at Q clock offset values, being centred at the intermediate clock bias estimation .sub.int, such that δt.sub.n|.sub.n=(Q+1)/2=.sub.int, [0064] d) {circumflex over (θ)}.sub.j(δt.sub.n)θ.sub.j(δt.sub.n,{dot over (δ)}t,p.sub.0.sup.T)={circumflex over (θ)}.sub.j(Φ.sub.0,n) is estimated as:

{circumflex over (θ)}.sub.j(Φ.sub.0,n)=x.sub.j.sup.HB(Φ.sub.0,n−ΔΦ.sub.j)(B.sup.H(Φ.sub.0,n−ΔΦ.sub.j)B(Φ.sub.0,n−ΔΦ.sub.j)).sup.−1B.sup.H(Φ.sub.0,n−ΔΦ.sub.j)x.sub.j.   (equation 6)

[0065] Algorithm 2 summarizes the alternative composite channel and fine clock bias estimation based on equation (4) when considering N.sub.l PDF hypotheses and N.sub.c Monte Carlo runs for each of the coefficients of the composite channel.

TABLE-US-00003 Algorithm 2 Alternative joint composite channel and fine clock bias estimation Input: {x.sub.j}.sub.j=1.sup.N, {{b.sub.i(t, Φ.sub.j)}.sub.j=1.sup.N}.sub.i=1.sup.M, .sub.c, α, N.sub.l, N.sub.c, W, {{p.sub.l,m}.sub.m=1.sup.D}.sub.l=1.sup.N.sup.l, Q, D, s, T.sub.δ Output: ĥ, .sub.f  1: Initialize γ.sub.min = ∞  2: Estimate the intermediate clock bias .sub.int based on equation (2)  3: Compute {circumflex over (θ)}.sub.j = [{circumflex over (θ)}.sub.j(δt.sub.1) . . . {circumflex over (θ)}.sub.j(δt.sub.Q)].sup.T for each antenna centered at .sub.int based on equation (6)  4: Compute {circumflex over (Λ)} = [{circumflex over (Λ)}(δt.sub.1) . . . {circumflex over (Λ)}(δt.sub.Q)].sup.T centered at .sub.int based on equation (1)  5: for l = 1 to N.sub.l do  6:  for c = 1 to N.sub.c do  7:   Draw the normalized channel coefficients [h.sub.l,c.sup.norm].sub.m ~p.sub.l,m, ∀.sub.m  8:   for j = 1 to N do  9:    h.sub.l,c,j = h.sub.l,c.sup.norm ∥{circumflex over (θ)}.sub.j∥/∥Σ.sub.m=1.sup.D[h.sub.l,c.sup.norm].sub.m|s( .sub.int + (m − 1)T.sub.δ)|.sup.2∥ 10:   end for 11:   for β = −W + 1 to W − 1 do 12:    for j = 1 to N do 13:     ε.sub.l,c,β,j = {circumflex over (θ)}.sub.j − Σ.sub.m=1.sup.D[h.sub.l,c,j].sub.m|s( .sub.int − (β − m + 1)T.sub.δ)|.sup.2 14:    end for 15:    Form E.sub.l,c,β = [ε.sub.l,c,β,1 . . . ε.sub.l,c,β,N].sup.T 16:    C.sub.l,c,β = E.sub.l,c,βE.sub.l,c,β.sup.T 17:    Estimate the eigenvalues {λ.sub.j}.sub.j=1.sup.N of C.sub.l,c,β 18:    [00014]${\gamma }_{l,c,\beta }=\frac{\frac{1}{N}{.Math.}_{j=1}^N{\lambda }_j}{{\left({.Math.}_{j=1}^N{\lambda }_j\right)}^{1\text{/}N}}$ 19:    if γ.sub.l,c,β < γ.sub.min then 20:     γ.sub.min = γ.sub.l,c,β 21:     ĥ = h.sub.l,c.sup.norm ∥{circumflex over (Λ)}∥/∥Σ.sub.m=1.sup.D[h.sub.l,c.sup.norm].sub.m|s( .sub.int + (m − 1)T.sub.δ)|.sup.2∥ 22:     {circumflex over (β)} = β 23:    end if 24:   end for 25:  end for 26: end for 27: .sub.f = .sub.int − {circumflex over (β)}T.sub.δ

[0066] Note that the inputs and the outputs in Algorithm 2 are the same as for Algorithm 1. Accordingly, Algorithm 2 can be seen as an alternative approach to performing the same overall analysis as performed in Algorithm 1. Furthermore, the skilled person will understand that these approaches are provided by way of example, and other approaches may be adopted according to the circumstances of any given implementation. For example, rather than using a Monte Carlo approach, one could perform a systematic search across a predetermined grid covering the space defining the possible (expected) values of ĥ.sub.l and {circumflex over (β)}.sub.l. Furthermore, rather than performing a Monte Carlo or grid search, some implementations might use an optimization (feedback) scheme, which seeks to converge on some minimum for Equation (3) (analogous to a hill-climbing routine, only seeking a minimum rather than a maximum). In all cases, it will be understood that finding a minimum configuration will generally imply finding an approximate minimum (rather than an absolute minimum), e.g., due to limitations in the analysis performed (such as considering only N.sub.w values for {circumflex over (β)}.sub.l in equation (3)).

[0067] The method described herein has been simulated using the implementation of Algorithm 1 via a semi-analytical approach at post-correlation level by considering the realistic indoor propagation channel from T. Jost et al., “A Wideband Satellit-to-Indoor Channel Model for Navigation Applications,” IEEE Transactions on Antennas and Propagation, vol. 62, no. 10, October 2014, pp. 5307-5320, and Report ITU-R P.2145-2, “Model Parameters for the Physical-Statistical Wideband Model in Recommendation ITU-R P.681,” September 2017 and a receiver having a linear array of four antennas (with a separation of 10 cm between each of them). In the simulation N.sub.1=3 (i.e., three hypotheses on the PDF of each of the composite channel coefficients were considered) and N.sub.e=10. The spatial correlations between the signals received for each antenna are considered in the model. As a reference, the timing solution obtained based on the received signals from one of the antennas of the array and a conventional weighted least squares (WLS)-based solution is used. A single integration period is considered in both cases.

[0068] FIG. 4 shows an example of the power delay profile for simulated indoor conditions (based on the channel model of T. Jost et al., “A Wideband Satellit-to-Indoor Channel Model for Navigation Applications,” IEEE Transactions on Antennas and Propagation, vol. 62, no. 10, October 2014, pp. 5307-5320), with attenuations above 20 dB and important NLOS components for all the satellites. In particular, FIG. 4 shows the simulation results for a set of 11 visible satellites (SV1, SV2, . . . SV11), where the results for each satellite are plotted with a symbol having a particular shape/color for that respective satellite. FIGS. 4A and 4B show a subset of the data from FIG. 4, with FIG. 4A showing results for one satellite (denoted SV2) and FIG. 4B showing results for another satellite (denoted SV8).

[0069] Each symbol having the particular shape/color for a given satellite represents a signal component from that satellite. It can be seen that the number of signal components resulting from a single navigation signal (i.e., from one satellite) is relatively high—generally more than 10 or 20, typically in the range 30-120, e.g., 50-100, but potentially higher still. Each signal component is plotted based on:

i) x-axis—the delay for that component in nanoseconds, relative to the ideal LOS component. By way of scaling, a delay of 10 nanoseconds corresponds to an additional travel distance of around 3 meters.
ii) y-axis—normalized power (dB), where normalization here indicates that the total received power (across all components) for each satellite is set to a consistent level (i.e., the same for all satellites).
Note that there is considerable variation in the signal components received from different satellites, For example, FIG. 4B for satellite SV8 shows significantly more components with a higher delay (>100 ns) than FIG. 4A for satellite SV2). In addition, it will be appreciated that the desired accuracy in the delay estimation (of around 50 ns) is significantly less than the delay of some NLOS components (even for satellite SV2).

[0070] The simulation results are summarized in Table 1, which shows the timing errors (RMS and 95-percentile) for the WLS-based reference method and the approach described herein, specifically adopting algorithm 1 (note that additional calibration and satellite system dependent errors are not considered in the simulations). In particular, these Figures relate to the timing errors associated with determining {circumflex over (β)}.sub.l, the clock bias of the receiver compared with the time reference maintained by the GNSS. The results obtained show a clear improvement in the time estimation of the clock bias performed with the method described herein compared with a conventional time estimation approach. This improvement is achieved by the approach described herein making better use of the NLOS signals (mitigating the bias introduced in the clock bias estimation), and by the exploitation of the signal diversity introduced by having an array of multiple antennas (although the approach described herein can be used even if only a single antenna is available, such as by utilizing Algorithm 1A as presented above).

TABLE-US-00004 TABLE 1 SUMMARY OF RESULTS Method RMS (ns) 95% (ns) WLS-based 130 249 Algorithm 1  27  51

[0071] Although the method described herein has a higher computational complexity than conventional timing architectures, it is well-suited to a parallel implementation in order to reduce the computational burden. Furthermore, the computation of the composite channel may be performed over a longer period, taking advantage of the relatively slow decorrelation of the indoor channel for static receivers (typically of the order of several seconds).

[0072] The approach described herein allows satellite navigation signals (including, but not limited to, current GNSS) to help achieve timing accuracies below around 50 ns (1−sigma) in indoor conditions. Such timing accuracies may help to support the use of current GNSS signals in indoor timing applications (without having to rely on other commercial services), such as the proposed 5G network operations mentioned above. In the approach described herein, a joint composite channel and time estimation method is provided to address the low carrier-to-noise (C/No) conditions and the high impact of NLOS multipath components which are generally experienced in such an indoor environment. As described above, this approach has been assessed via software simulations considering a realistic indoor propagation channel and shows a clear improvement of the resulting time synchronization performance with respect to a conventional time estimation approach.

[0073] Although the approach described herein is primarily intended to allow an accurate determination of the timing offset of a receiver clock for the case of indoor reception of GNSS signals, it can also be used more generally to determine the multipath components (such as shown in FIG. 4) of the received signals. This approach may be used to help remove bias (distortion) from the cross-correlation functions, thereby providing a more accurate estimate of the location (instead of or as well as the clock offset) of the receiver in an indoor environment. In addition, the approach may be used to obtain information about the indoor environment. For example, if the observed multipath components change, this may indicate some change in conditions within the environment (such as a door opening or closing), and hence might be utilized, for example, for some form of monitoring system.

[0074] Although the approach described herein is primarily intended for use with existing or future GNSS such as GPS and Galileo, the term GNSS should be understood to include any system which provides suitable satellite (or pseudolite) signals for performing the methods described herein, even if such signals are not necessarily provided for navigation purposes (or at least, this is not their main purpose). For example, the satellite signals utilized herewith might be primarily intended for time synchronisation and/or for telecommunications (rather than for position location).

[0075] The approach described herein provides a receiver for use with a global navigation satellite system (GNSS) comprising multiple satellites. Each satellite transmits a respective navigation signal containing a spreading code. The GNSS maintains a time reference which is encoded into the navigation signals. The receiver comprises a receiver clock and at least one antenna for receiving multiple signals over multiple respective channels, each channel being defined by a transmitting satellite and a receiving antenna at opposing ends of the channel. The receiver further comprises at least one correlator for calculating cross-correlation functions between (i) the signals received over the multiple channels and (ii) reference versions of the navigation signals provided by the receiver. The receiver is configured to use the calculated cross-correlation functions to perform a joint estimate of (i) a clock bias of the receiver clock relative to the time reference maintained by the GNSS, and (ii) a composite channel comprising the combined contribution of the multiple channels as a function of time-delay.

[0076] The signals received over the respective channels typically include contributions from multiple navigation signals; these are contributions are usually below the noise level, but can be discriminated by cross-correlating with the known codes for the different navigation signals. Each cross-correlation function is generally calculated based on (i) a corresponding received signal (for a given channel) and (ii) a corresponding navigation signal (for a given satellite). The receiver will normally only detect a subset of the navigation signals—the remainder may correspond, for example, to satellites that are below the horizon. (In some cases the receiver may have almanac information available and know upfront which satellites are currently below the horizon) It is also possible that a given navigation signal may be detected only on a subset of the channels, for example, due to strong attenuation or fading on one or more other channels. If a given cross-correlation function for a corresponding channel and satellite does not locate the navigation signal (i.e., no correlation peak is found), this cross-correlation function may be discounted from the joint estimate.

[0077] In general, the receiver is configured to use the calculated cross-correlation function between the signal received over a given channel and the reference version of the given navigation signal to determine a code-phase alignment between the received navigation signal and the reference version of this navigation signal. This code-phase alignment is dependent upon: (i) the propagation time of the given navigation signal from the satellite to the receiver and (ii) the clock bias between the receiver clock and the time reference maintained by the GNSS. For a receiver of known position, propagation time can be determined, thereby exposing the clock bias to measurement.

[0078] The receiver described herein is well-suited for use in indoor conditions, which generally suffer from (i) fading or sometimes complete absence of a direct (line-of-sight) component, and (ii) the presence of usually multiple non-line-of-sight components. The joint estimation utilizes information received over all channels (including information from multipath components) to help obtain a higher accuracy for the clock bias than can generally be achieved with existing techniques.

[0079] In some implementations, the clock bias of the receiver clock is estimated as an offset from an approximate clock bias. The approximate clock bias may be derived from conventional techniques and have an uncertainty which is greater than the desired accuracy for the clock bias provided herein. For example, the approximate clock bias may have a typical accuracy (1−sigma) in the range 0.1-1 microseconds, or in the range 0.2-0.8 microseconds, or in the range 0.25-0.6 microseconds. In some implementations, the receiver may be configured to estimate the approximate clock bias, prior to performing the joint estimation, by using the calculated cross-correlation functions to determine a maximum of a composite cross-correlation function for the composite channel. (This latter approach will generally provide a more accurate preliminary estimate for the approximate clock bias than conventional techniques, and so reduce the search space for the joint estimate, as described below).

[0080] The composite channel may represent a multiple-input multiple-output (MIMO) channel, where the multiple inputs correspond to the various satellites, and the multiple outputs correspond to different antennas. However, in some cases the receiver may have just a single antenna, in which case there is only a single output. Note that in cases where there are multiple antennas (i.e., a MIMO channel), the joint estimate may take advantage of spatial correlation between the signals received by different antennas to enhance the accuracy of the joint estimate.

[0081] The joint estimate of the clock bias and composite channel can be considered as an optimization problem, i.e., finding a clock bias and composite channel that would most closely reproduce the calculated cross-correlation functions. This optimization can be performed using any suitable search technique, for example, a grid search, an iterative technique (such as hill-climbing), and so on. For the clock bias, the search may be performed within a window centered (at least approximately) on an overall code-phase alignment of the calculated cross-correlation functions, or by using any other suitable method. As noted above, having a prior known estimate for the clock bias reduces the search space, the more accurate the prior known estimate, the smaller the search space can be.

[0082] The composite channel represents a more complex parameter space. For example, the composite channel may be specified by a set of coefficients, each coefficient corresponding to a particular time-delay, and representing the aggregate level of power transmission through the composite channel for that time-delay (e.g., as a fraction of the total power received into the channel). In some implementations, the receiver includes at least one probability distribution function for each coefficient in the set of coefficients, and this probability distribution function is used to generate trial values for the coefficients. In some cases, the search space may encompass the use of multiple different probability distribution functions for the coefficients. There are various ways to define or provide such probability distribution functions. For example, the probability distribution functions may be determined from simulations of typical indoor propagation channels, such as shown in FIGS. 4, 4A and 4B. Another possibility is that the provided probability distribution functions are based on an estimate of the composite channel obtained using received signals from a previous epoch (if available). For example, the first epoch of measurement may use probability distribution functions obtained from simulations, however, subsequent epochs may use probability distribution functions determined from earlier epochs. In some cases, the receiver may be configured to update the existing probability distribution functions after each new joint estimate of the composite channel becomes available, hence allowing a continuous and ongoing (iterative) enhancement of the probability distribution functions used for the coefficients.

[0083] Also provided herein is a method of operating a receiver for use with a global navigation satellite system (GNSS) comprising multiple satellites, wherein each satellite transmits a respective navigation signal containing a spreading code, and wherein the GNSS maintains a time reference which is encoded into the navigation signals. The receiver comprises a receiver clock and at least one antenna for receiving multiple signals over multiple respective channels, each channel being defined by a transmitting satellite and a receiving antenna at opposing ends of the channel. The method comprises calculating cross-correlation functions between (i) the signals received over the multiple channels and (ii) reference versions of the navigation signals provided by the receiver; and performing a joint estimate of (i) a clock bias of the receiver clock relative to the time reference maintained by the GNSS, and (ii) a composite channel comprising the combined contribution of the multiple channels as a function of time-delay, using the calculated cross-correlation functions. This approach is particularly suitable for use in an indoor environment, which is typically subject to multipath (non-line-of-sight) components and/or fading of the direct (line-of-sight) components. This method can generally benefit from the same features and enhancements discussed above.

[0084] The approach described herein may be implemented into some form of receiver, whether provided as a standalone receiver or incorporated into some other type of device e.g., a mobile telephone, a cellular radio access node (e.g., a 5G small cell), a tablet, a vehicle GPS (car, plane, ship, bicycle, etc.), camera, or any other appropriate location monitoring/tracking system. The receiver may implement such a method at least in part by executing program instructions on one or more processors (such as DSPs, GPUs, etc.), and/or the receiver may implement such a method at least in part by utilising special purpose hardware e.g., one or more chips specially designed to support the method, e.g., by providing multiple hardware correlators, etc. A computer program may also be provided for implementing the approach described herein. Such a computer program typically comprises computer program instructions which may be provided on a non-transitory computer readable storage medium, such as an optical disk (CD or DVD), hard disk, or flash memory. The computer program may be loaded into a computer memory from such a storage medium, or may be downloaded into the computer memory over a network, such as the Internet. The apparatus, e.g., receiver, loading the computer program may comprise one or more processors for executing the computer program, which comprises instructions that cause the apparatus to implement a method such as described above.

[0085] In conclusion, a variety of implementations have been described herein, but these are provided by way of example only, and many variations and modifications on such implementations will be apparent to the skilled person and fall within the scope of the present disclosure, which is defined by the appended claims and their equivalents.