光子学报  2017, Vol. 46 Issue (5): 532002-  DOI: 10.3788/gzxb20174605.0532002
0

引用本文  

宋佳凝, 徐国栋, 董立珉, 徐振东, 张兆祥. 基于TOA信息的脉冲星信号周期估计方法[J]. 光子学报, 2017, 46(5): 532002-. DOI: 10.3788/gzxb20174605.0532002.
SONG Jia-ning, XU Guo-dong, DONG Li-min, XU Zhen-dong, ZHANG Zhao-xiang. Pulse Period Estimation Method Based on TOA Information[J]. Acta Photonica Sinica, 2017, 46(5): 532002-. DOI: 10.3788/gzxb20174605.0532002.

Foundation item

The National High Technology Research and Development Program of China (No. 2008AA8051602)

First author

SONG Jia-ning (1991-), female, Ph.D candidate, mainly focuses on pulsar based navigation and signal processing.Email:hitsjn@163.com

Corresponding author

XU Guo-dong (1961-), male, professor, mainly focuses on communication and information engineering, spacecraft design and pulsar based navigation. Email:xgdong_61@163.com

Article History

Received: Nov. 24, 2016
Accepted: Feb. 15, 2017
基于TOA信息的脉冲星信号周期估计方法
宋佳凝, 徐国栋, 董立珉, 徐振东, 张兆祥    
(哈尔滨工业大学 卫星技术研究所, 哈尔滨 150001)
摘要:针对脉冲星信号周期估计需逐步搜索、计算量大等问题,提出脉冲星信号周期直接估计算法.根据周期不确定误差对时域上脉冲到达时间(Pulse Time of Arrival,TOA)估计方法的影响,推导了基于TOA信息的脉冲星信号周期估计方法的数学模型.该方法将一组观测到的脉冲光子到达时间序列(Photon Time of Arrivals,TOAs)进行等时间间隔分段,在时域上对每段TOAs进行相位估计获得相应的TOA信息,并根据建立的TOA信息与周期误差的关系采用最小二乘原理估计脉冲星信号周期.对周期准确性的评价不同于传统的依赖折叠轮廓的好坏,而是以TOA与时间图像的斜率作为依据,这一依据更直观、易于理解.理论分析与物理、计算机仿真数据结果表明,本文所提的周期估计方法能够根据短时间脉冲光子到达时间序列获得高精度和高分辨率的脉冲信号周期,有利于X射线脉冲星导航的工程应用.
关键词脉冲星    光子    导航    脉冲到达时间    周期估计    
中图分类号:O434.19;V446+.9      文献标识码:A      文章编号:1004-4213(2017)05-0532002-10
Pulse Period Estimation Method Based on TOA Information
SONG Jia-ning, XU Guo-dong, DONG Li-min, XU Zhen-dong, ZHANG Zhao-xiang    
(Harbin Institute of Technology, Research Center of Satellite, Harbin 150001, China)
Foundation item: The National High Technology Research and Development Program of China (No. 2008AA8051602)
First author: SONG Jia-ning (1991-), female, Ph.D candidate, mainly focuses on pulsar based navigation and signal processing.Email:hitsjn@163.com
Corresponding author: XU Guo-dong (1961-), male, professor, mainly focuses on communication and information engineering, spacecraft design and pulsar based navigation. Email:xgdong_61@163.com
Received: Nov. 24, 2016; Accepted: Feb. 15, 2017
Abstract: Aiming at the problem of huge computation in pulsar period estimation methods, a direct pulsar period estimation algorithm was proposed. Based on the discussion of the influence of period errors on pulse time of arrival estimation methods, the mathematical model of the said period estimation method using TOA information is derived. In this method, a set of photon time of arrivals(TOAs) is divided into several segments in an equal time interval, and TOA information corresponding to each segment is calculated using time domain methods. According to the formula of the TOA information and period error, the Least Square method is adopted to estimate the period. A new intuitionistic criterion for the precision of pulsar period is developed by the value of the slope of the φ-t graph, which is different from the scheme that searches the perfect period based on maximum peak principle. Theoretical analysis and the results of experiments utilizing physical and numerical data are demonstrated that the presented pulsar period estimation method can achieve a precise and high-resolution period from a short observation of photon time of arrivals, which can help to realize the engineering application of X-ray pulsar navigation.
Key words: Pulsar    Photon    Navigation    Pulse time of arrival    Period estimation    
OCIS Codes: 320.5550;030.5290;000.4430
0 Introduction

With full autonomy and uninterrupted advantages, X-ray pulsar-based navigation(XNAV) provides spacecraft plentiful navigation information, such as position, velocity and time, and has become one of the prevailing in autonomous navigation[1-3]. In the XNAV system, spacecraft receives photons emitted from pulsars and records the photon time of arrivals(TOAs) which is the fundamental measurement of the navigation system. After the pulse time of arrival(TOA) is estimated from a set of observed TOAs, the navigation measurement equation can be established for an autonomous navigation system[4-6]. The accuracy of estimated TOA information mainly determines the positioning accuracy of XNAV system. Acquiring high precise TOA in short observation time is the key to XNAV. Nowadays, TOA estimation algorithms in the time domain are comparatively mature and can be classified into two categories, including the Maximum Likelihood Estimation(MLE) method using TOAs[7] and the phase estimation methods based on epoch folding[8]. The time-domain TOA estimation methods demand a stringent signal period. For epoch-folded methods, the quality of folded profile is affected by the period accuracy. For the MLE method, in order to reduce the enormous computation, the TOA is estimated using the given period. Thus an exact period is necessary to guarantee a good accuracy of the estimated TOA. Nevertheless, the period of pulsar signals observed by spacecraft is shifting due to the doppler effect caused by spacecraft motion. Besides, the unpredicted glitches of pulsars also have an impact on the period of pulsar signals[9]. Therefore, spacecraft is expected to have a certain ability to estimate and correct the given period of pulsar signals according to the received TOAs onboard.

Many researchers have studied the period estimation methods of pulsar signals. Aiming at the radio signals measured by radio astronomy observatories, researchers proposed fast Fourier transform algorithm[10], cyclo-period search method based on maximum correlation[11], quick search method of pulsar period using Lomb algorithm[12] and so forth. In frequency domain methods of pulsar period estimation, sampling frequency should be selected carefully to reduce the spectrum leakage and fence effect. What's more, enormous computation is needed to estimate an accurate period. Aiming at pulse time of arrival array in X-ray band, Zhang[13] compared the rebuilt profile set generated by a fast folding algorithm with the standard profile, and use the cross-correlation method to estimate pulsar period. While the performance of Zhang′s method is affected by the sampling rate and amplitudes.

In this paper, a direct pulsar period estimation algorithm is studied to reduce the complexity of the traditional methods. We firstly discuss the influence of period errors on the estimated TOA, and then present mathematical models of the proposed period correction algorithm. The accuracy and computational complexity of the proposed period estimation algorithm are also analyzed by comparing with the estimation method. Finally, the feasibility and the performance of the presented algorithm are illustrated by experiments based on physical and numerical simulation data.

1 Period estimation method using TOA information 1.1 The proposed method using TOA information from epoch-folding TOA estimation algorithms

All the time tags during an observation time are folded into a single given pulse period to derive the empirical pulsar profile in pulse delay estimation algorithms based on epoch folding. TOA is estimated by comparing the difference between the empirical and the standard pulsar profiles. It can be seen clearly that pulsar period errors have an influence on the process of epoch folding and the result of the estimated TOA.

Assuming a given pulsar period is P, the difference between the given pulsar period and the true pulsar period is dP, the length of the observation time is Tobs which consists of Np pulsar periods. If the given pulsar period equals to the true one exactly, that is dP=0, each cycle is time aligned and the folded profile has a sharp peak without any broadening. While if the given pulsar period differs from the true one (dP≠0), the phase difference between two adjacent cycle profiles can be given by

${\rm{d}}\varphi = \frac{{{\rm{d}}P}}{{P + {\rm{d}}P}}$ (1)

It should be pointed out that empirical profiles and dφ can hardly be obtained in a single cycle due to the low energy flow density of pulsars. Additionally, dP is very small, usually in the order of nanosecond. Thus, in order to estimate TOA and correct periods, we divide a long observation time into many pieces which contain considerable and reasonable cycles. Subsequently, the TOA of each piece is calculated by folding the correspondent piece into a profile. The phase difference between two pieces within time T can be described as

${\rm{d}}\varphi \left( T \right) = \frac{{{\rm{d}}P}}{{P + {\rm{d}}P}}\cdot\frac{T}{{P + {\rm{d}}P}}$ (2)

where, dφ(T)=φ(t+T)-φ(t), φ(t) is the pulse phase(TOA) at time t corresponding to a piece TOAs of the whole observation time. Next, the calculated pulse phases are plotted on the φ-t graph, and the slope kφ-t can be estimated by the least square method. The equation of kφ-t and period errors is

${k_{\varphi - t}} = \frac{{{\rm{d}}P}}{{{{\left( {P + {\rm{d}}P} \right)}^2}}}$ (3)

where dP is the positive real root of the above quadratic equation. On account of dPP, the denominator on the right side of Eq.(3) can be approximated to P2 and the approximate solution can be calculated easily.

1.2 The proposed method using TOA information from a maximum likelihood estimator

Employing the probability density function associated with the detected TOAs, a maximum-likelihood estimator can be formulated to estimate the initial phase φ0 and the unknown frequency f. The unknowns can be found by solving the following[14]

$\left( {{{\hat \varphi }_0},\hat f} \right) = {\rm{arg}}\;{\rm{max}}\sum\limits_{i = 1}^N {{\rm{ln}}\left[ {\beta + \alpha h\left( {{{\hat \varphi }_0} + \left( {{t_i} - {t_0}} \right)\hat f} \right)} \right]} $ (4)

where β, α are the known effective background rate and source arrival rate respectively. h(·) is the standard periodic pulsar profile and usually defined on the phase interval φ∈[0, 1] (cycle). Furthermore, the function h(φ) is non-negative and normalized according to $\int_0^1 {h\left( \varphi \right){\rm{d}}\varphi = 1} $, and minφh(φ)=0. As for the periodical pulsar signal, there is h(φ)=h(φ+n), where n is a positive integer.

Using the given pulsar period to calculate the frequency in Eq.(4), the two-dimensional parameter estimation problem degenerates into a one-dimensional phase estimation problem. Period errors influence the result of estimated phases with time continuing, showing as follows

$\begin{array}{l} {{\hat \varphi }_0} = {\rm{arg}}\;{\rm{max}}\sum\limits_{i = 1}^N {{\rm{ln}}\left[ {\beta + \alpha h\left( {{{\hat \varphi }_0} + \left( {{t_i} - {t_0}} \right)\left( {{f_{{\rm{ref}}}} + {\rm{d}}f} \right)} \right)} \right]} = \\ \quad {\rm{arg}}\;{\rm{max}}\sum\limits_{i = 1}^N {{\rm{ln}}\left[ {\beta + \alpha h\left( {\left[ {{{\hat \varphi }_0} + \left( {{t_i} - {t_0}} \right){\rm{d}}f} \right] + \left( {{t_i} - {t_0}} \right){f_{{\rm{ref}}}}} \right)} \right]} \end{array}$ (5)

where, fref is the reference frequency corresponding to the given pulsar period P, df is the error of the fref relative to the true frequency of the pulsar signal. Comparing Eq.(5) to Eq.(4), it is clear to see that the df has a linear impact on the estimated phase with time t increasing

$\varphi \left( t \right) = {{\hat \varphi }_0} + \left( {t - {t_0}} \right){\rm{d}}f$ (6)

Similarly, by plotting the φ-t graph and calculating the slope kφ-t, where kφ-t=df. The period error can be expressed as

${\rm{d}}P = P - \frac{1}{{{P^{ - 1}} + {k_{\varphi - t}}}}$ (7)
2 Accuracy and computational complexity analysis

The proposed period estimation method can be implemented as follows: firstly, a set of photon arrival time series is segmented, illustrating in Fig. 1, where, Tp is the time length of each segment, Tm is the time interval of the neighboring segment, and n is the total number of all segments. Next, the TOA of each segment is estimated in time domain by using the given pulsar period P. Finally, the period is corrected by calculating the slope of the φ-t graph and Eq. (3) or (7).

Fig.1 Schematic diagram of data grouping
2.1 Accuracy analysis

From the above process, the accuracy of the proposed method is affected by TOA estimation methods, parameters Tp and Tm. The parameter Tp reflects the signal accumulative time of the process of the TOA estimating and has a direct influence on the accuracy of TOA information according to Ref.[14]. There is an inversely proportional function of the accumulative time Tp and the theoretical standard deviation δφ(Tp) based on different TOA estimators, including a Nonlinear Least-Squares (NLS) estimator, a Cross Correlation (CC) estimator, and a Maximum Likelihood (ML) estimator.

We now take a CC estimator as an example to discuss the precision of the proposed method in theory. According to Eq.(2), it follows that if the parameter Tm is too small, there will be dφ(Tm) < δφ(Tp). In that case, dφ(Tm) is submerged in δφ(Tp) and the phase difference between two neighboring segments cannot reflect the phase bias caused by period errors. Hence we have the following restriction: parameters Tm and Tp should meet dφ(Tm) > δφ(tp) to make sure dφ(Tm) can be detected. While if we use many segments to estimate phases and calculate period errors, the above restriction can be loosened by

${\rm{d}}\varphi \left( {{T_m}} \right) \cdot n > \delta \varphi \left( {{T_p}} \right)$ (8)

where n is the number of segments. With the help of Eq.(2), we can show Eq.(8) clearly. The left side of Eq.(8) can be expressed as

${\rm{d}}\varphi \left( {{T_m}} \right) \cdot n = \frac{{{\rm{d}}P}}{{P + {\rm{d}}P}} \cdot \frac{{{T_m}}}{{P + {\rm{d}}P}} \cdot n = {\rm{d}}\varphi \left( {{T_m} \cdot n} \right)$ (9)

Eq.(9) shows that increasing the number of segments is equal to enlarge the parameter Tm in some sense. It can also be explained by the simulating results of Fig. 2. Where three or four points correspond to the same TOA value(y axis), looking like 'footsteps'. These 'footsteps' are caused by dφ(Tm) < δφ(Tp). While the number of segments is large enough and the parameters are satisfied to Eq.(8), the slope kφ-t still can be obtained.

Fig.2 The φ-t graph

Combining Eq.(8) and (2), the resolution of the estimated period in theory δP is defined as

$\delta P \approx \frac{{\delta \varphi \left( {{T_p}} \right) \cdot {P^2}}}{{{T_{{\rm{obs}}}}}}$ (10)

Next, we turn to study the precision of the proposed method in the aspect of the numerical solution method. Let Nb be the number of folding 'bins' in CC and NLS estimators. If we expect to distinguish the phase bias among two neighboring segments from epoch-folded estimation methods, there should be dφ(Tm)>Nb-1. Actually, the above restriction is too strong and can be loosened by

${\rm{d}}\varphi \left( {{T_m}} \right) \cdot n > 1/{N_b}$ (11)

The Eq.(11) illustrates that the phase bias between the last and the first segment of an observation time Tobs should be greater than the estimated TOA resolution. Combining Eqs.(11) and (2), we have another theoretical resolution of the estimated period δP,

$\delta P \approx \frac{{{P^2}}}{{{T_{{\rm{obs}}}}}} \cdot \frac{1}{{{N_b}}}$ (12)

The same conclusion can also be given when a ML estimator is utilized to estimate TOA. The phase resolution of a ML estimator is determined by the search step noted as 'Sstep'. Therefore substituting 'Sstep'' for '1/Nb' in Eq.(12), we have a δP based on a ML phase estimator

$\delta P \approx \frac{{{P^2}}}{{{T_{{\rm{obs}}}}}} \cdot {S_{{\rm{step}}}}$ (13)

By Eqs.(10), (12) and (13), the theoretical resolution of the estimated period δP of the proposed method can be presented as

$\delta P = {\rm{max}}\left( {\frac{{\delta \varphi \left( {{T_p}} \right) \cdot {P^2}}}{{{T_{{\rm{obs}}}}}},\frac{{{P^2}}}{{{T_{{\rm{obs}}}}}} \cdot \frac{1}{{{N_b}}}\quad {\rm{or}}\quad \frac{{{P^2}}}{{{T_{{\rm{obs}}}}}} \cdot {S_{{\rm{step}}}}} \right)$ (14)

From Eq.(14), the estimated period δP is determined by the max of dφ(Tm) and 1/Nb (or Sstep). If we let 1/Nb be 0.001 and dφ(Tm) < 0.001(cycle), the δP of a millisecond pulsar during 1000s observation time will be in the order of nanosecond. For most TOA estimation methods, the accuracy of the estimated TOA is better than 0.001(cycle). Thus the proposed method can have a good period estimation with less demand on the accuracy of TOA information from the above evaluation. It also should be pointed out that the period error is closer to the theoretical resolution of the estimated period, the accuracy of the proposed method is worse. This is because during an observation time Tobs, the phase bias caused by period errors is limited and the effective number of estimated phases corresponding to different segments becomes small.

Noting that the time derivative of pulsar barycentric rotation frequency is not taken into account in our method. While in the actual situations, the periods of pulsar signals decrease with the increasing time. For example, the time derivative of Crab pulsar frequency is about 1×10-10 Hz/s and the frequency will have a 1×10-7 Hz shift during 1 000 s observation time. In that case, if the period error dp is over 1ns, the changed frequency caused by period error will be -9×10-7 Hz and the period error can be obtained by our method. Otherwise, our method will not work well because the changed frequency caused by the period error is submerged in the frequency shift. Fortunately, for most pulsars, the time derivative of pulsar frequency is in the order of 1×10-13 Hz/s or smaller[15]. Hence, the time derivative of pulsar frequency has a very small impact on the proposed method.

2.2 Discussion about selecting parameters Tp, Tm and n

Now considering the influence of the parameters Tp and n on the estimated kφ-t, the estimated error variance of kφ-t from the least square method can be given[16]

${\rm{var}}\left( {{k_{\varphi - t}}} \right) = \frac{{{\sigma ^2}}}{{\sum\limits_{i = 1}^n {{{\left( {{x_i} - E\left( {{x_i}} \right)} \right)}^2}} }}$ (15)

where, φi~N(b+kφ-txi, σ2), σ=δφ(Tp) and xi=i·(Tm). The relationship between parameters Tp, Tm, n, and Tobs can be approximated by

${T_p} + n \cdot {T_m} \approx {T_{{\rm{obs}}}}$ (16)

The inversely proportional function of δφ(Tp) and Tp can be expressed as[14]

$\delta \varphi \left( {{T_p}} \right) = \frac{1}{{{T_p}}} \cdot \frac{{\int_0^1 {\lambda \left( \varphi \right){{\dot \lambda }^2}\left( \varphi \right){\rm{d}}\varphi } }}{{{{\left( {\int_0^1 {{{\dot \lambda }^2}\left( \varphi \right){\rm{d}}\varphi } } \right)}^2}}} = \frac{1}{{{T_p}}} \cdot A$ (17)

where $A = \frac{{\int_0^1 {\lambda \left( \varphi \right){{\dot \lambda }^2}\left( \varphi \right){\rm{d}}\varphi } }}{{{{\left( {\int_0^1 {{{\dot \lambda }^2}\left( \varphi \right){\rm{d}}\varphi } } \right)}^2}}}$. By substituting Eqs.(16) and (17) into (15), we have

$\begin{array}{l} {\rm{var}}\left( {{k_{\varphi - t}}} \right) = \frac{{{{\left[ {\delta \varphi \left( {{T_p}} \right)} \right]}^2}}}{{\sum\limits_{i = 1}^n {\left( {i{T_m} - \frac{{\sum\nolimits_{i = 1}^n {i{T_m}} }}{{{n^2}}}} \right)} }} = \\ \frac{{12{{\left( {\frac{A}{{{T_p}}}} \right)}^2}}}{{T_m^2\left( {n - 1} \right)n\left( {n + 1} \right)}} \approx \frac{{12{A^2}}}{{n{{\left[ {{T_p}\left( {{T_{{\rm{obs}}}} - {T_p}} \right)} \right]}^2}}} \end{array}$ (18)

In order to acquire an accurate kφ-t, we choose the parameters Tp, Tm and n to make the Eq.(18) be as small as possible. Usually, if period errors and the chosen Tp are so large that the folded profile has a gross distortion, the peak of the folded profile gradually degenerates, which has an adverse impact on TOA estimation. Therefore, we may try different parameter combinations before finding a proper one sometimes. Some proper parameter combinations are given in the following experiments.

2.3 Computational complexity analysis

The computational complexity of the proposed period estimation is mainly determined by the TOA estimators. The computational costs for epoch folding, CC estimator, NLS estimator and ML estimator have been studied by Emadzadeh Amir Abbas and Speyer Jason in Ref.[14]. Here we list their results in Table 1, as well as the computational cost of the Least Square method.

In Table 1, Ng is the number of search grids in the interval (0, 1) cycle and M is the total photon number. For the proposed method, n times TOA calculation are needed to obtain the TOA information. From Eq.(3), we can see two multiplications and one division are needed for computing the slope kφ-t. From Eq.(7), one addition, one subtraction and two divisions are needed.

Tab.1 Computational costs for TOA estimation and epoch folding

Now we construct the cost function of the traditional χ2 period estimation method to make a comparison. In χ2 estimation method, different test periods are utilized to fold a profile and the index function ${\chi ^2} = \sum\nolimits_{j = 1}^{{N_b}} {{{\left( {{m_j} - \bar m} \right)}^2}} $[17]. Where mj is the number of photons in the jth bin and m is the mean number of photons in the folded profile. Hence, the total computational cost of one time χ2 estimation method is the sum of the epoch folding cost function and (Nb-1) additions, Nb subtractions, Nb multiplications and one division.

For ease of comparison and notations, let Ng be equal to Nb. Substitute [Tp/P] and Tp(α+β) for Np and M respectively, the computation cost results are shown in Table 2.

Tab.2 Computational cost for TOA estimation and epoch folding

Where nx is the number of χ2 estimation method times for searching the accurate period, and nx is determined by the searching interval and searching step. Usually, nx is in the order of hundred or thousand. n is relatively small and could be about one hundred. Thus, the computational cost of the proposed method using CC or NLS estimator is less than the computational cost of χ2 estimation method. Since the amount of calculations for the proposed method using ML estimator significantly increases as the number of photons Tp (α+β) becomes longer. Hence, when we using a ML estimator to calculate TOA information for the proposed method Tp is set to be not very large, which is consist with the conclusion in Section 2.2. And compared to the χ2 estimation method case, it is hard to say which one need more calculations. Because the amount of calculation for the proposed method using ML estimators and χ2 estimation method is affected by different parameters.

3 Experiments and performance

In this section, the physical simulation data of the ground X-ray source and MCP detectors and the numerical simulation data based on Crab pulsar parameters are utilized to verify the feasibility and performance of the proposed method.

3.1 Experiment one

Experiments using the physical simulation data can reflect the actual stimulating situation and help to improve the physical simulation system. The data from the ground X-ray source and MCP detectors is given in column, showing photon time of arrivals recorded by MCP detectors. The file also provides the period of pulsar signals used to experiments, with 33.4 ms. As experimenters illustrate that the X-ray source and MCP detectors don't move during tests, the estimated phase of each segment picked up from the entire TOAs will fluctuate around a constant if the given pulsar period is exact. Otherwise, estimated phases change with time increasing and the slope kφ-t of the φ-t graph is not equal to zero.

In accordance with the above discuss, we first divide TOAs into segments as indicated in Fig. 1. Then the given pulsar period P=33.4 ms is adopted to estimate TOA of each segment based on a CC estimator and a ML estimator. In order to acquire a suitable TOA precision, let Tp and Tm be 200 s and 7 s respectively. The Nb of the CC estimator and Sstep of the ML estimator are set to be 1 024 and 0.001 respectively. Tobs, the total length of the TOAs is 997 s. The analysis results see in Fig. 2. The amount operations for the proposed method using CC and ML estimators are 8.23×108 and 4.26×1010 respectively. If we use the above parameters to search period by χ2 estimation method, 3.06×107nx calculations are needed. It is clear to see that only when nx is equal to 27, the computational cost of χ2 estimation method is same with the proposed method based on CC estimator. However, it is quite difficult for nx to be so small in most cases. And there is the same computational cost with the proposed method based on ML estimator when nx equals to 1 395.

By Fig. 2, the slope kφ-t is a non-zero constant, which means the given pulsar period has period errors. Slopes for the CC estimator and the ML estimator are calculated using Least Square method, with -3.3959×10-5 s-1 and -3.3855×10-5 s-1, respectively. Correspondent period errors according to Eqs.(3) and (7) are -37.883 6 ns and -37.766 8 ns. As we do not have the exact pulsar period from the physical simulation system, we fold all the TOAs data into a profile using the given period and the calibrated period for comparison. It is apparent from Fig. 3 that folded profiles using the corrected period have a higher peak than the folded profile using the given period. Hence the corrected period is regarded as a more accurate one in common sense. It is also noticed that a constant phase difference exists between the CC estimator and the ML estimator in Fig. 2. We guess this is because of different data processing procedures. Here we only care about the slopes of the graph and the phase difference doesn′t impact the slope estimation. So we do not analyze this phenomenon in detail.

Fig.3 Contrast figure of folded profiles

Similarly, analysis with another group of physical simulation data in low signal-to-radio (SNR) is executed. The total length of this group is 2 042 s. Let Tp and Tm be 200 s and 7 s respectively. Slopes of φ-t graphs in Fig. 4 are -1.5856×10-5 s-1 and -1.5889×10-5 s-1 corresponding to the CC estimator and the ML estimator, respectively. And period errors from Eq. (3) and (7) are -17.688 0 ns(CC estimator) and -17.725 4 ns(ML estimator) respectively. Fig. 5 also illustrates the comparison of folded profiles using different pulsar periods. From this group analysis, it can be concluded that the proposed period estimation method can calculate period errors even though the observation has a low SNR.

Fig.4 The φ-t graph in low SNR
Fig.5 Contrast figure of folded profiles in low SNR

The reason why there are period errors in physical simulation experiments is mainly because of the system noise when X-ray source generates pulsar signals and clock errors when MCP detectors record the arrival time of photons. The period error is not expected in the physical simulation system. Since it has a negative influence on testing TOA estimation algorithms and other experiments based on the generated data. Therefore, it is necessary to correct the period error to guarantee a reliable and accurate period for further use.

3.2 Experiment two

In this experiment, the RXTE observation data[18] is utilized to simulate photon Time Of Arrivals(TOAs) in order to investigate the accuracy of the proposed method. By analyzing the data package FS_b36a20f-b36acab of Crab pulsar observed by RXTE, we have the ratio of the effective background rate and source arrival rate, about 9:1. And the event of TOAs can be demonstrated by a non-homogeneous Poisson process. By setting the exact pulsar period be 33.4 ms, a two-thousand-second TOAs observation is generated, with the effective background and source arrival rates at 450 photons/s and 50 photons/s respectively.

The experiment is carried out as follows: let tp=50 s and Tm=19 s. Firstly, the TOA of each segment is calculated over 100 Monte Carlo trails. The Root Mean Squared error(RMS) of TOA by CC and ML estimators are evaluated as 9.817×10-4 cycle and 8.652×10-4 cycle respectively. Secondly, the above RMS is used to calculate the theoretical resolution of the estimated period by Eq.(9), with the results of 0.55 ns (CC estimator) and 0.48 ns (ML estimator). While taking into account the Nb and Sstep that are chosen to be 1 024 and 0.001 in these two estimators, theoretical resolutions of the estimated period corresponding to CC and ML estimators should be 0.55 ns and 0.56 ns according to Eq. (14), respectively. Thirdly, different period errors, including 0.8 ns, 1 ns, 10 ns, 0.1 μs, and 1 μs, are added to the exact period(33.4 ms) to re-estimate the TOA of all segments. In order to reduce the computation, we only use CC estimator to estimate TOA. Finally, slopes of φ-t graphs are computed using Least Square method, and the corrected periods are calculated according to Equation (3).

Fig. 6 shows different φ-t graphs of various period errors. Considering the integer part of the estimated phase, we redefine φ on the phase interval φ∈[0, +∞)(cycle). With dP increasing, the linear relationship becomes apparent and the absolute value of the slope kφ-t increases gradually. In Fig. 6(a) and (b), as dP is small and δφ(Tp) is relatively large, the linear function between φ and t is not obvious, which agrees with the analysis in Section 2. Fig. 7 compares the folded profiles using the different periods. Table 3 also shows the statics of period errors over 100 Monte Carlo trials. The RMS of different periods displayed in Table 3 is in the same order with about 0.35 ns and does not affect by the value of period errors. The results of this experiment show that though the period has a big error, making the folded profile a serious distortion, the proposed period estimation method is able to calculate periods in a quick and direct manner.

Fig.6 φ-t graphs in different period errors
Fig.7 Folded profiles in different periods
Tab.3 Statics results of the period error estimation

It should be commented that an iteration algorithm can be developed to acquire a more precise estimated period in the further studies since there are some approximations in the data process. The iteration algorithm can be designed as follows: the said method is first to execute to calculate a new period, then using the new period we do the same procedure again based on the said method. When the difference between the new period and the last one reaches the threshold that we set, the iteration algorithm stops.

4 Conclusion

In this paper, we demonstrate a pulsar period estimation method based on TOA information. The proposed method established the formula between period errors and estimated phases of each segment from a set of TOAs. The main contribution of the proposed method is to provide a direct method to estimate periods in the time domain with a new criterion for the precision of pulsar period. It is totally different from other period estimation methods in frequency domain or methods that search the perfect period based on maximum peak principle. The proposed method can estimate period in the order of nanosecond from a hundred seconds of observation time. Using the simulating data with a photon arrival rate 500 photons/s and total observation time 2 000 s, the proposed method has a theoretical resolution of the estimated period at 0.6ns, estimating the period at 0.8 ns successfully. The proposed pulsar period estimation method can help to improve the ground physical simulation system and be applied in engineering.

References
[1] MITCHELL J W, WOOD K S, LITCHFORD R J, et al. SEXTANT-station explorer for X-ray timing and navigation technology[C]. AIAA Guidance, Navigation, and Control Conference, 2015.
[2] LUI Xiu-ping, JING Jun-feng, SUN Hai-feng, et al. Denoising of X-ray pulsar signal based on wavelet-Fisz transformation[J]. Acta Photonica Sinica, 2014, 43(12): 1204001. DOI:10.3788/gzxb
[3] WANG Yi-di, ZHENG Wei, SUN Shou-ming, et al. X-ray pulsar-based navigation using time-differenced measurement[J]. Aerospace Science and Technology, 2014, 36: 27-35. DOI:10.1016/j.ast.2014.03.007
[4] SHEIKH S I. The use of variable celestial X-ray source for spacecraft navigation[D]. Department of Aerospace Engineering University of Maryland, 2005.
[5] WINTERNITZ L M B, MITCHELL J W, HASSOUNEH M A, et al. SEXTANT X-ray pulsar navigation demonstration: flight system and test results[C]. IEEE Aerospace Conference, 2016.
[6] FENG Dong-zhu, GUO He-he, WANG Xin, et al. Autonomous orbit determination and its error analysis for deep space using X-ray pulsar[J]. Aerospace Science and Technology, 2014, 32(1): 35-41. DOI:10.1016/j.ast.2013.12.008
[7] ANDERSON K D, PINES D J, SHEIKH S I. Validation of pulsar phase tracking for spacecraft navigation[J]. Journal of Guidance, Control, and Dynamics, 2015, 38(10): 1885-1897. DOI:10.2514/1.G000789
[8] EMADZADEH A A, SPEYER J L. On modeling and pulse phase estimation of X-ray pulsars[J]. IEEE Transactions on Signal Processing, 2010, 58(9): 4484-4495. DOI:10.1109/TSP.2010.2050479
[9] SUN Hai-feng, BAO Wei-min, FANG Hai-yan, et al. Effect of stability of X-ray pulsar profiles on range measurement accuracy in X-ray pulsar navigation[J]. Acta Physica Sinica, 2014, 63(06): 441-448.
[10] BURNS W R, CLARK B G. Pulsar search technology[J]. Astronomy and Astrophysics, 1969(02): 280-287.
[11] LI Jian-xun, KE Xi-zheng, ZHAO Bao-sheng. A new time-domain estimation method for period of pulsars[J]. Acta Physica Sinica, 2012, 61(06): 537-543.
[12] ZHOU Qing-yong, JI Jian-feng, REN Hong-fei. Quick search algorithm of X-ray pulsar period based on unevenly spaced timing data[J]. Acta Physica Sinica, 2013, 62(01): 533-540.
[13] ZHANG Xin-yuan, SHUAI Ping, HUANG Liang-wei. Profile folding distortion and period estimation for pulsar navigation[J]. Journal of Astronautics, 2015, 36(09): 1056-1060.
[14] EMADZADEH A A, SPEYER J L. Navigation in Space by X-ray Pulsars[M]. New York: Springer, 2011.
[15] ATNF Pulsar Catalogue Version 1. 55[DB/OL]. 2016-11-14[2017-01-20]. http://www.atnf.csiro.au/research/pulsar/psrcat/
[16] XU Guo-dong, SONG Jia-ning, LI Peng-fei. Pulsar navigation adaptive filtering algorithm based on information quality[J]. Optics and Precision Engineering, 2015, 23(3): 827-837. DOI:10.3788/OPE.
[17] MAO Yue. Research on X-ray pulsar navigation algorithms[D]. PLA Information Engineering University, 2009.
[18] NASA HEASARC FTP. 40805-01-05-000[DB/OL]. 2012[2016-08-01]. http://heasarc.gsfc.nasa.gov/FTP/rxte/data/archive/AO4/P40805/40805-01-05-000.