2.1. Fundamentals of pulsar navigation

The fundamentals of pulsar navigation can be expressed by the following equation

(1)$${{\bi r}_{{\bi 12}}} \cdot {\bi n} = c \cdot ({t_e} - \Delta t)$$

where r₁₂ is the relative position vector between two spacecraft, n is the unit vector in the direction of the pulsar, c is the speed of light, t _e is the time delay of pulsar arrival between two spacecraft, and Δt is the clock error. In the navigation process, the value of time delay is obtained by estimating the phase delay of pulsar signals.

2.2. Epoch folding procedure

Epoch folding is usually the first step to estimate the phase delay. In this subsection, a brief review of the epoch folding procedure is given.

To describe the rate of photons arrival, an overall rate function λ(t) is defined by (Golshan and Sheikh, Reference Golshan and Sheikh2007)

(2)$$\lambda (t) = {\lambda _b} + {\lambda _s}h(\varphi (t))$$

where λ _b is a background arrival rate, λ _s is a source arrival rate, h(ϕ) is the periodic pulsar profile, and ϕ(t) is the phase of the pulsar signal.

Epoch folding is the procedure that recovers the X-ray pulsar intensity density function from the measured photons' Times Of Arrival (TOAs). This procedure was formulated mathematically by Emadzadeh et al. (Reference Emadzadeh, Golshan and Speyer2009).

Assume that the spacecraft continually observes a pulsar for N _P pulsar periods, and a pulsar period is divided into N _b bins equally. So the length of each bin is T _b = P/N _b, where P is the period of the pulsar. Let c(t _i) represent the number of photons in the ith bin whose centre is t _i, then the empirical rate function can be written as (Emadzadeh and Speyer, Reference Emadzadeh, Golshan and Speyer2009)

(3)$$\hat \lambda ({t_i}) = \displaystyle{1 \over {{N_P}{T_b}}}\sum\limits_{\,j = 1}^{{N_P}} {{c_j}({t_i})} $$

and it can be proved that for a long observation time T _obs,where $\hat n({t_i})$ is an uncorrelated noise that satisfies

(5)$$E[\hat n({t_i})] = 0$$

(6)$${\mathop{\rm var}} [\hat n({t_i})] = \displaystyle{{{N_b}} \over {{T_{obs}}}}\lambda ({t_i})$$

The overall rate function λ(t) can be regarded as a probability density function of the arrival photons. The purpose of epoch folding is to recover the probability density function using the frequency distribution histogram of the samples (i.e., TOAs of photons).

Figure 1 shows a typical frequency distribution histogram, where a photon (black spot) is represented by a rectangle with width of T _b and a unit height. By adding up all these rectangles, we can obtain the empirical rate function $\hat \lambda ({t_i})$.

Figure 1. A frequency distribution histogram.

Assume that the number of photons arriving during the observation time is N. For every photon, it is necessary to find the exact bin it falls in and update the corresponding ${c_j}({t_i})$. So, the computation complexity of the epoch folding procedure is $O(N)$. In Section 3.2, it will be shown that representing a photon by a delta-function instead of a rectangle can diminish the computational complexity.

2.3. Phase delay estimators

In this subsection, three methods to estimate the phase delay are presented: maximum-likelihood estimator, non-linear least squares estimator and cross-correlation estimator.

The fundamental of maximum-likelihood estimator is to find a phase delay that maximises the joint probability of the observed TOAs (Sala et al., Reference Sala, Urruela, Villares, Jordi and Sebastia2008). In practice, log-likelihood is commonly used

(7)$$\Psi (\varphi ) = \sum\limits_{i = 1}^N {\ln (\lambda ({t_i};\varphi ))} $$

And the estimation of phase delay is

(8)$$\hat \varphi = \arg \max \Psi (\varphi )$$

which requires optimising a non-convex function, and the grid search is a general choice. Assume that we search N _b candidates of ϕ, and for every ϕ, ψ(ϕ) should be computed, which requires a computational complexity of $O(N{N_b})$. Besides, finding the maximum of all the ψ(ϕ)s needs $O({N_b})$ complexity. Therefore, the total computational complexity is O(NN _b + N _b) for the maximum-likelihood estimator.

The non-linear least squares estimator minimises a so-called cost function (Emadzadeh et al., Reference Emadzadeh, Golshan and Speyer2009)

(9)$$J(\varphi ) = \sum\limits_{n = 1}^{{N_b}} {{{(\hat \lambda (n;\varphi ) - \lambda (n))}^2}} $$

where λ(n) and $\hat \lambda (n;\varphi )$ are the discretizations of the overall and empirical rate functions, respectively. The estimation of phase delay is

(10)$$\hat \varphi = \arg \min J(\varphi )$$

This process needs epoch folding with the complexity of O(N), and for every ϕ we should compute J(ϕ), which requires a complexity of O(N _b²). To find the minimal of J(ϕ), O(N _b) a computational cost is generated. Thus, the total computational complexity of the non-linear least squares estimator is O(N + N _b² + N _b).

The cross-correlation estimator maximises the cross correlation function of the overall rate function and the empirical rate function (Emadzadeh and Speyer, Reference Emadzadeh and Speyer2011a), which can be expressed as

(11)$$R(\varphi ) = \int_{ - \infty }^\infty {{\lambda ^*}(t)\hat \lambda (t;\varphi )dt} $$

And the estimation of phase delay is

(12)$$\hat \varphi = - \arg \max R(\varphi )$$

In practice, DFT is usually employed to compute the cross correlation function

(13)$$\eqalign{{\Lambda _1}(k) =\, & DFT[\lambda (n)],\;{\Lambda _2}(k) = DFT[\hat \lambda (n)], \cr & \quad R(n ) = IDFT[{\Lambda _1}(k){\ast} {{\Lambda _2}^{\ast}}(k)]} $$

This process also needs epoch folding, and the DFT process requires a computational cost of $O({N_b}\log {N_b})$ with the fast Fourier transformation (FFT) algorithm (Cooley and Tukey, Reference Cooley and Tukey1965). Then the maximum of R(n) can be found with a complexity of O(N _b). So the total complexity is $O(N + {N_b}\log {N_b} + {N_b})$ for the cross-correlation estimator.

Table 1 summarises the computational complexities of the above three estimators. It can be seen that when the number of photons N and the number of bins N _b are large, the cross-correlation estimator has the least computational cost.

Table 1. Computational complexities of estimators.

3.1. Proposal of the DFT based method

In this section, a method for phase delay estimation based on the time-shift property of DFT is proposed.

According to the time-shift property of DFT, a cyclic shift of a sequence in time domain is equivalent to a phase shift in frequency domain. That is, if

(14)$$\left\{ {\matrix{ {X(k) = DFT[x(n)]} \cr {y(n) = x{{(n - m)}_N}} \cr}} \right.$$

then,

(15)$$DFT[\,y(n)] = X(k)\exp \left( { - j\displaystyle{{2\pi} \over N}mk} \right)$$

where j ² = −1.

As a consequence, the time shift in time domain m can be obtained by computing the phase shift $ - \displaystyle{{2\pi} \over N}mk$ in frequency domain.

The intention of phase delay estimation is to find a phase delay $\hat \varphi $ in time domain such that

(16)$$\hat \lambda (n) = \lambda {(n - \hat \varphi )_{{N_b}}} + \hat n(n)$$

where $\hat n(n)$ is an uncorrelated noise.

Assume that

(17)$$\left\{ {\matrix{ {{\Lambda _1}(k) = DFT[\lambda (n)]} \cr {{\Lambda _2}(k) = DFT[\hat \lambda (n)]} \cr}} \right.$$

and according to the time-shift property of DFT,

(18)$$DFT[\lambda (n - \hat \varphi )] = {\Lambda _1}(k)\exp \left( { - j\displaystyle{{2\pi} \over {{N_b}}}\hat \varphi k} \right)$$

It can then be obtained by performing a DFT on Equation (16) that

(19)$${\Lambda _2}(k) = {\Lambda _1}(k)\exp \left( { - j\displaystyle{{2\pi} \over {{N_b}}}\hat \varphi k} \right) + N(k),\quad 0 \le k \le {N_b} - 1$$

where $N(k) = DFT[\hat n(n)]$ is the DFT of the noise, and Equation (19) is valid for any $k \in [0,{N_b} - 1]$.

So, the estimation of phase delay is

(20)$$\hat \varphi = - \displaystyle{{{N_b}} \over {2\pi k}}ang\displaystyle{{{\Lambda _2}(k)} \over {{\Lambda _1}(k)}}$$

where $ang\displaystyle{{{\Lambda _2}(k)} \over {{\Lambda _1}(k)}}$ is the phase angle of the complex $\displaystyle{{{\Lambda _2}(k)} \over {{\Lambda _1}(k)}}$, and Equation (20) is also valid for any $k \in [0,{N_b} - 1]$.

Instead of computing the whole DFT sequence, this method only computes Λ₁(k) and Λ₂(k) of a certain $k \in [0,{N_b} - 1]$ to obtain the estimation of phase delay. Using the definition of DFT, it can be obtained that

(21)$${\Lambda _2}(k) = \sum\limits_{n = 1}^{{N_b}} {\hat \lambda (n)\exp \left( { - j\displaystyle{{2\pi} \over {{N_b}}}nk} \right)} $$

With the Horner Algorithm (Horner, Reference Horner1819), the computations of Λ₁(k) and Λ₂(k) by Equation (21) only require O(N _b) complexity. This process also needs epoch folding with the complexity of O(N). As a result, the total computational complexity is O(N + N _b), smaller than that of estimators presented in Section 2.

As stated above that for any $k \in [0,{N_b} - 1]$, the phase delay can be estimated by Equation (20). The best estimation of phase delay can be obtained by choosing the ‘best’ k which minimises the effect of the noise $\hat n(n)$.

Equation (19) leads to

(22)$$\displaystyle{{{\Lambda _2}(k)} \over {{\Lambda _1}(k)}} = \exp \left( { - j\displaystyle{{2\pi} \over {{N_b}}}\hat \varphi k} \right) + \displaystyle{{N(k)} \over {{\Lambda _1}(k)}}$$

Since the amplitude and phase angle of N(k) are random, the k to be chosen is the one that maximises the amplitude of Λ₁(k).

(23)$$k = \arg \max \vert{\Lambda _1}(k)\vert$$

In this way, the effect of $\displaystyle{{N(k)} \over {{\Lambda _1}(k)}}$ can be minimised.

Searching for the ‘best’ k requires the computation of the whole DFT sequence of λ(n), but the overall rate function λ(n) is ‘known’ before the spacecraft launched, so the ‘best’ k and corresponding Λ₁(k) can be found ‘offline’, which means the computations of k and corresponding Λ₁(k) impose no complexity for on board computation.

3.2. Delta-Function Approximation

Since k and Λ₁(k) are obtained ‘offline’, during the process of navigation, only Λ₂(k) is computed using Equation (21), where the empirical rate function $\hat \lambda (n)$ can be obtained through the epoch folding process. In this subsection, a delta-function approximation is put forward to replace epoch folding. With this approximation, one can directly obtain Λ₂(k) in a lower complexity.

As discussed before, the epoch folding procedure regards photons as rectangles and uses a histogram to estimate the overall probability density function. Similarly, we can regard the photons as delta functions

(24)$$\delta (t - \mu ) = \left\{ {\matrix{ {\infty, t - \mu = 0} \cr {0,t - \mu \ne 0} \cr}} \right.$$

whose centres are the ‘fractions’ of TOAs. Assume that the TOA of a photon contains ${N_P}$ pulsar periods, then the ‘fraction’ is defined as

(25)$$\mu = TOA - {N_P}P$$

The empirical rate function $\hat \lambda (t)$ can be redefined by

(26)$$\hat \lambda (t) = \sum\limits_{i = 1}^N {\delta (t - {\mu _i})} $$

where μ _i is the ‘fraction’ of the ith photon.

By noting the time-shift property of DFT and DFT [δ(t)] = 1, it can be derived that

(27)$$DFT[\delta (t - \mu )] = \exp ( - jk\mu )$$

With the use of DFT on Equation (26), it can be obtained that

(28)$${\Lambda _2}(k) = \sum\limits_{i = 1}^N {{e^{ - jk{\mu _i}}}} $$

where μ _i is the ‘fraction’ of the ith photon.

Equation (28) means that, with the delta-function approximation, the computation of empirical rate function $\hat \lambda (t)$ can be skipped and Λ₂(k) can be obtained directly from TOAs of photons. As a consequence, the total computational cost can be reduced to O(N), significantly less than that of the proposed method in Section 3.1. On the other hand, employing the delta-function approximation will not significantly lose accuracy, which will be shown in Section 4.

Furthermore, with the delta-function approximation, the memory cost is also decreased, because there is no need to save the empirical rate function $\hat \lambda (n)$, whose length is N _b.

3.3. Scheme of the proposed method

Based on the above description, the scheme of the proposed method is clarified in this subsection.

When k > 1, the phase angle of the complex $\displaystyle{{{\Lambda _2}(k)} \over {{\Lambda _1}(k)}}$ may exceed 2π, meaning that the 2Mπ part of $ang\displaystyle{{{\Lambda _2}(k)} \over {{\Lambda _1}(k)}}$ cannot be found directly by using trigonometric calculation. Consequently, the estimation of phase delay cannot be obtained directly through Equation (20) and it should be rewritten as

(29)$$\hat \varphi = - \displaystyle{{{N_b}} \over {2\pi k}}\left( {ang\displaystyle{{{\Lambda _2}(k)} \over {{\Lambda _1}(k)}} + 2M\pi} \right)$$

where M is an integer.

To solve this problem, let us set k = 1 at first, and make a ‘rough’ estimation of $\hat \varphi $ by

(30)$${\hat \varphi _1} = - \displaystyle{{{N_b}} \over {2\pi}} ang\displaystyle{{{\Lambda _2}(1)} \over {{\Lambda _1}(1)}}$$

Then, substitute the ‘rough’ estimation into Equation (29) to solve M

(31)$$M = round\left( { - \displaystyle{{k{{\hat \varphi} _1}} \over {{N_b}}} - \displaystyle{1 \over {2\pi}} ang\displaystyle{{{\Lambda _2}(1)} \over {{\Lambda _1}(1)}}} \right)$$

where round(x) is the arithmetical operation to find the closest integer of x. Using the solution of M, the estimation of phase delay $\hat \varphi $ can be finally determined by Equation (29).

The scheme of the DFT-based method with the delta-function approximation can be summarised as follows:

1. (1) Offline procedure:

   1. a) Compute Λ₁(k) = FFT [λ(n)]

   2. b) Find k = arg max[Λ₁(k)]

2. (2) Online procedure:

   1. a) Compute ${\Lambda _2}(1) = \sum\limits_{i = 1}^N {{e^{ - j{\mu _i}}}} $, where μ = TOA−N _PP

   2. b) Compute ${\Lambda _2}(k) = \sum\limits_{i = 1}^N {{e^{ - jk{\mu _i}}}} $

   3. c) ‘Rough’ estimation ${\hat \varphi _1} = - \displaystyle{{{N_b}} \over {2\pi}} ang\displaystyle{{{\Lambda _2}(1)} \over {{\Lambda _1}(1)}}$

   4. d) Find $M = round\left( { - \displaystyle{{k{{\hat \varphi} _1}} \over {{N_b}}} - \displaystyle{1 \over {2\pi}} ang\displaystyle{{{\Lambda _2}(1)} \over {{\Lambda _1}(1)}}} \right)$

   5. e) $\hat \varphi = - \displaystyle{{{N_b}} \over {2\pi k}}\left( {ang\displaystyle{{{\Lambda _2}(k)} \over {{\Lambda _1}(k)}} + 2M\pi } \right)$