Fast Butterfly Epoch Folding-Based X-Ray Pulsar Period Estimation With a Few Distorted Profiles

The traditional pulsar period estimation method using a single distorted profile averts multiple epoch folding processes, but still has space for accuracy to improve. To give consideration to both accuracy and computational load, we combine the fast butterfly epoch folding (FBEF) with the single distorted profile (SDP)-based pulsar period estimation, and propose the FBEF-based pulsar period estimation with a few distorted profiles (FBEF-FDP), which merely uses a few epoch folding processes to enhance accuracy sharply. Firstly, the fast butterfly epoch folding is used to produce a few distorted pulsar profiles with different periods. And then the single distorted profile-based pulsar period estimation explores each distorted prolife to provide a corresponding distortion degree. Finally, using these distortion degree values, the nonlinear search method for multiple minimums can get an optimal distortion degree estimation. Besides, we prove that distortion degree estimation errors of different accumulative distorted profiles from the same pulsar photon sequence are incompletely correlated. This is the theorical basis that the distortion degree estimation values from the same pulsar photon sequence can be fitted to obtain more highly accurate distortion degree estimation. The computational complexity demonstrates that the computational load increases with the number of groups. Thus, to ensure the real-time of the FBEF-FDP, the number of groups should not be too large. Simulation results demonstrate that the accuracy of the fast butterfly epoch folding-based pulsar period estimation using a few distorted profiles enhances at the cost of a small computational load.


I. INTRODUCTION
In the field of autonomous navigation for deep space exploration [1], [2], X-ray pulsar navigation has great potential for development. Its whole frame can be summarized as follows: X-ray pulsar navigation system uses X-ray sensor to receive the X-ray signals from X-ray pulsars, and then extracts the pulse Time-of-Arrival (TOA) from these X-ray signals. Finally, using the TOAs as measurements, the extended Kalman filter estimates the position and velocity of spacecraft [3], [4]. Time-of-Arrival is a basic observation of X-ray pulsar navigation [5], [6]. The X-ray detector at the spacecraft collects the photons radiated by the X-ray pulsar [7], and forms the accumulative pulsar profiles [8] by epoch folding. And then it is combined with the standard The associate editor coordinating the review of this manuscript and approving it for publication was Ali Zemouche. pulsar profile, pulsar timing model and other information to estimate the pulsar TOA. However, the spacecraft motion and the pulsar glitches [9], [10] lead to a change in the period of the received pulsar signals. The change of the pulsar period gives rise to the pulsar profile accumulated in the inherent period to be distorted, which causes the pulsar TOA to drift. The pulsar period estimation is of great significance for highprecision estimation of the pulsar TOA, and has become a research hotspot.
The pulsar period estimation methods can be divided into two categories: one uses the drift of the pulsar TOA to estimate the period error, and the other one explores the pulsar profile distortion to inverse the pulsar period. The latter is the current mainstream.
The basic idea of the period estimation methods based on pulsar TOA drift is as follows: they establish the relationship model between the period error and the drift of the VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see http://creativecommons.org/licenses/by/4.0/ pulsar TOA. On this basis, the least square method is used to estimate the pulsar TOA and the period error [11], [12]. In order to achieve the purpose of suppressing the impact of the period error, other tools can be utilized, such as: the EKF-CMBSEE (Extended Kalman Filter-Correlated Measurement Bias and State Estimation Error) [13], and the tracking filter [14]. The basic idea of the pulsar estimation methods based on the profile distortion is as follows: They try to fold the pulsar signals according to different periods to obtain a series of the distorted pulsar profiles, and then find the accumulative pulsar profile with the smallest distortion. The corresponding period is the inherent period. Currently, typical methods is as follows: the time domain χ 2 method, the fast butterfly epoch folding (FBEF) [15], the Fourier spectrum method [16], the coherent statistic of cyclostationary signal [17], the search method based on maximum correlation variance [18], the Lomb-based χ 2 algorithm [19], the pulsar profile entropy-based method [20], the time-delay compensation method [21] and the profile feature method [22]. All of the above methods need to try to fold the pulsar signals in different periods. However, the amount of pulsar radiation signal data is huge, which leads to problems such as large calculation amount, and is not suitable for on-board computer. For instance, when the area of X-ray sensor and the observation time are 1000cm 2 and 5000s, respectively, the amount of pulsar radiation signal data is on the order of 10 7 . For the χ 2 methods, due to the lack of exploration of the pulsar profile, their accuracy is low.
In order to solve this problem, Liu et al. have proposed the single distorted profile (SDP)-based pulsar period estimation [23]. In this method, the compressive sensing [24] is used to compress the pulsar radiation signal. And then the pulsar inherent period is estimated directly from the distortion of the individual pulsar profile. As the SDP explores the pulsar profile, its accuracy is higher than those of the χ 2 methods. This method avoids the multiple folding of pulsar signals and greatly reduces the amount of calculation, which makes it possible to estimate the pulsar TOA in orbit. However, there is room for improvement in its calculation accuracy. The reason is that the single distorted profile-based pulsar period estimation uses only one distortion profile for estimation. The information utilization rate is limited.
The profile distortion-based pulsar period estimation methods require to fold pulsar signals by multiple periods, which leads to the problem of a large computational load. Fortunately, using the butterfly summation calculation method similar to the fast Fourier transform, the redundant summation operation of the folding process in different periods is reduced. The calculation speed of the fast butterfly epoch folding [15] is improved.
In this paper, we combine the advantages of both the single distorted profile-based pulsar period estimation and the FBEF, and propose the FBEF-based X-ray pulsar period estimation with a few distorted profiles (FBEF-FDP). The FBEF folds the pulsar photon sequence by multiple periods to obtain a series of accumulative pulsar profiles. Then, the single distorted profile-based pulsar period estimation is used to estimate the distortion degree of each accumulative pulsar profile. We intend to prove that the noises of these accumulative pulsar profiles are incompletely correlated. Therefore, the estimation errors of these distortion degrees are also not completely correlated. This is a prerequisite for combining these distortion errors to obtain higher accuracy of these pulsar period estimates.
The paper is organized as follows: After the introduction, the fast butterfly epoch folding is reviewed and simplified in Section 2. In Section 3, the FBEF-based X-ray pulsar period estimation with a few distorted profiles is developed. The correlation of pulsar noises and the computational complexity are analyzed in Section 4 and 5, respectively. In Section 6, simulations demonstrate the superiority and feasibility of the FBEF-FDP. Finally, conclusions are drawn in Section 7.

II. REVIEW AND SIMPLIFICATION OF FAST BUTTERFLY EPOCH FOLDING
The tentative epoch folding adopts multiple pulsar period to fold the pulsar photon sequence. This leads to a high computational load. Fortunately, similar to the fast Fourier transform, the fast butterfly epoch folding reduces the redundant summation operation in the epoch folding processes with different periods. And its computational load declines sharply. The FBEF [15] includes two steps.
Step 1: Pulsar intensity matrix. Assume that the number of bins in a pre-estimate period is N . Divide the pulsar photon sequence into N g pre-estimate periods. The pulsar intensity matrix (N g * N ) is constructed. N g satisfies inequation (1): (1) Step 2: Butterfly accumulation. N g pre-estimate periods is divided into N g /2k subsets (k = 1, 2, . . . ,log 2 N g ). And each subset is further divided into the first half and the last half regions. The pulsar signals of the j-th period in the first half region is shifted to the left by j bins and j+1 bins. And then they are accumulated by the pulsar signals of the j-th period in the last half region. A new pulsar intensity matrix is obtained by above operations. Thus, we have N g /2k matrices from the N g /2k subsets. Repeat Step 2. Until k = log 2 N g , we obtain the final pulsar intensity matrix. From inequation (1), we can see that N g is very large, as N is on the order of ten thousand or more. This results in an extremely high computational load of Step 2. To reduce the computational load sharply, we simplify the fast butterfly epoch folding. Its diagram is presented as Fig. 1. The details are given as follows: Step 1: Sub-profile. Assume that the total observation time is T . The pulsar signals in the observation interval [0, T ] are equally divided into N g groups, where N g is power of 2. Fold the pulsar signals in each group to get each sub-profile. In the simplified FBEF, N g is not too large. Commonly, its value is several or several tens. Thus, the computational load of the simplified FBEF declines sharply. Step 2: Butterfly accumulation. The process is the same as the traditional one.

III. X-RAY PULSAR PERIOD ESTIMATION WITH A FEW DISTORTED PROFILES
The number of the distorted pulsar profiles produced by the simplified FBEF is very small. This results in the fact that the pulsar profile without distortion may not belong to the set of these distorted pulsar profiles. The SDP operates normally with a single distorted profile, and avoids this problem essentially. Theoretically, it still has space for accuracy to improve when multiple distorted profiles are explored.
In this paper, we combine the FBEF with the SDP, and propose the FBEF-based pulsar period estimation with a few distorted profiles. In this method, the simplified FBEF produces a few distorted profiles. And then the single distorted profile-based pulsar period estimation explores each distorted prolife to provide a distortion degree. Finally, using these distortion degree values, the nonlinear search method for multiple minimums gets an optimal distortion degree estimation.
The diagram of this method includes the following three steps: (1) the fast butterfly epoch folding, (2) the single distorted profile-based matching, and (3) the search method for multiple minimums.
The fast butterfly epoch folding and its simplification have been introduced in Section 2, so no longer elaborated.

A. SINGLE DISTORTED PROFILE-BASED MATCHING
Traditional SDP method utilizes fast Fourier transform and compressive sensing to estimate the X-ray pulsar period. However, fast Fourier transform introduces imaginary number. Thus, in this paper, we directly adopt the crosscorrelation to achieve the single distorted profile-based matching. The diagram of this method includes the following three steps: Step 1: Distorted profile dictionary. The distorted pulsar profile dictionary is constructed by the circular correlation of the standard pulsar profile and the spread vector.
The distorted pulsar profile dictionary, (N × M ), can be designed to where m represents the distortion spread degree. M is the maximum distortion spread degree. The atom, ϕ m (N ×1), stands for the m-th atom in the dictionary , which can be expressed by where h (N × 1) is the standard pulsar profile, % denotes for the circular correlation. G m with size N × 1 represents the spread vector, which can be expressed as follows, Step 2: Matching matrix. Using the correlation between the dictionary, (N × M ), and the accumulated pulsar profile matrix,H, we can construct the matching matrix, S as follows: whereH is the set of the shifted pulsar profiles.
where H ±p denotes for the left shift or the right shift of the accumulated pulsar profile, and P is the range of the phase shift.
Step 3: Super-resolution estimation with maximum elements. Using the row and column subscripts of maximum value in the matching matrix, the super-resolution estimation is obtained.
To further improve the estimation accuracy, we use the super-resolution estimation with maximum elements, which is presented as follows: where n and m are the pre-estimated row and column subscripts of the maximal value of S, respectively. The elements max S m+1 and max S m−1 are the maximal elements in the (m+1) th and (m-1) th columns, respectively. The process of the SDP is shown in Fig. 3.

B. NONLINEAR SEARCH METHOD FOR MULTIPLE MINIMUMS
We can get multiple distortion degree values from those distorted pulsar profiles in the previous section. In this section, we use the search method to obtain the optimal distortion degree from those distortion degree values. As we known, the core of the search method is to construct the objective function. Next, we introduce the objective function. After correlating the accumulative pulsar profile H j with the dictionary atom, the distortion estimation m j is obtained by using the super-resolution estimation with maximum elements, where j = 0, 1, 2, . . . , N g -1. Assume that the real value of m j is d + j, where d denotes for the distortion bias caused by the period bias. The difference between m j and d +j is the estimation error, and its absolute value is: N g accumulative pulsar profiles correspond to N g estimation errors, respectively. The objective function is the sum of the absolute values of N g estimated errors, which is shown as: When d s is minimum, the difference between those estimated distortion values and the true ones is the smallest. At this time, the corresponding distortion bias, d, is the solve.
Obviously, the above objective function is nonlinear. Thus, we use the search method to get the optimal estimation. If there are multiple minimum values, the mean of the qualified d is taken as the optimal distortion degree.
The traditional super-resolution estimation using the largest element can merely estimate the error's amplitude, and not the plus-minus sign. Fortunately, our search method overcomes this problem.
Assume that f m is the frequency offset corresponding to the atom, which can be described by T is the corresponding pulsar period estimation error, which can be expressed by where T 0 represents the inherent pulsar period.
Therefore, the pulsar period estimation can be given as follows.

IV. CORRELATION OF THE PULSAR NOISES
We utilize the same pulsar photon sequence to accumulate multiple distorted pulsar profiles, and then estimate the distortion degrees of these accumulative pulsar profiles. Only when these distortion degree errors are incompletely correlated, the nonlinear search method for multiple minimums can obtain more highly accurate distortion degree estimation. Next, we investigate the correlation of the pulsar noises.
Let the pulsar profiles be h i (n), i = 0, 1, 2, . . . N − 1, the standard pulsar profile be h(n). The relationship between the two is as follows: where δ b i (n) is the background noise, whose variance is σ 2 . The accumulative pulsar profiles H j (n), j = 0, 1, 2, . . . N g − 1, can be expressed as follows: For the two accumulative pulsar profiles from the same pulsar photon sequence, H j (n) and H j+b (n), the correlation coefficient of their noise, r n , is as follows: As the noise correlation coefficient of different pulsar profiles is zero, the above equation can be simplified as: If Eq.(18) holds, then: Conversely, if Eq.(18) does not hold, then: For the sake of analysis, this article only discusses the necessary conditions for the establishment of Eq.(18): In Eq.(21), when the difference between two distortion degree values, b, increases, i decreases. Assume that I is the maximum value at which Ieq. (21) holds. We can get Ieq. (22).
The correlation coefficient of noises satisfies: That is, the noises of two accumulative pulsar profiles are not completely correlated. Besides, from above equations, we can see that the increase of the difference between two distortion degree values leads to the decline of the correlation coefficient of noises.

V. ANALYSIS OF COMPUTATIONAL COMPLEXITY
In this section, we perform a complexity analysis of the FBEF-FDP. As the dictionary construction can be performed on the ground and does not occupy the computing resources of spacecraft, it is not considered. The search method for multiple minimums is computationally insignificant, and can be ignored. Computing resources are mainly consumed in the fast butterfly epoch folding and the single distorted profilebased pulsar period estimation.
(1) FBEF. The FBEF requires log 2 N g steps of the butterfly calculation. Each step of the butterfly calculation requires N g additions for the pulsar profile. The number of the pulsar profile bins is N . That is, the calculation amount of addition operation for two pulsar profile is N . Therefore, the amount of calculation of the FBEF is N · N g · log 2 N g .
(2) SDP. The essence of the single distorted profile-based pulsar period estimation is the matching of a single pulsar profile to a dictionary. Assume that the number of atoms in the dictionary is M . To ensure that the pulsar period estimation algorithm is immune to the phase, each distorted pulsar profile in the directory is phase shifted to form a series of the pulsar profiles with different phases. Assume that the range of the phase shift is [-(P-1)/2, (P-1)/2]. As the number of the pulsar profile bins is N , the computational amount of the accumulative profile with a certain phase multiplied by a single atom is 2N . Therefore, the calculation amount of one SDP is 2N · P · M . For N g distorted profiles, the calculation amount for the matching estimation is 2N · P · M · N g .
That is, the total calculation amount of the FBEF-FDP is: For instance, suppose that N , P, M and N g are 33000, 161, 80, 32, respectively. The total calculation amount of the FBEF-FDP is 2.7 × 10 10 MAC.
From the above equation, we can see that the computational load of the FBEF-FDP is proportional to the number of the pulsar profile bins, the range of the phase shift, the number of distorted profiles and the number of groups. Obviously, to reduce the computational load of the FBEF-FDP, we should reduce the values of the above four variables. The number of the pulsar profile bins is determined by the time resolution of X-ray detector, and cannot be changed arbitrarily. Under the priori conditions, we minimize the number of distorted profiles and the range of the phase shift. In order to ensure a moderate amount of calculation, the number of groups is not too large.

VI. SIMULATION
In order to manifest the feasibility and superiority of the FBEF-FDP, we compare it with the SDP-based pulsar period estimation. Firstly, we search the optimal number of groups. And then, we study the correlation coefficient of the distortion degree from the same pulsar photon sequence. Finally, the influence of both the X-ray detector's effective area and the observation time on the estimation error of the pulsar period is analyzed. Besides, we also investigate the influence of both the background noise and the pulsar photon flux.

A. SIMULATION CONDITIONS
The simulation conditions are presented in Table 1. The standard profile of PSR B0531+21 can be found in the European Pulsar Network [25]. The experiments were carried out on a Dell notebook with an Intel Core i7-7700HQ CPU @ 2.80 GHz and 8 GB RAM.

B. OPTIMAL NUMBER OF GROUPS
We investigate the performance of the FBEF-FDP vs. the number of groups as the number of groups is a vital factor of the FBEF-FDP. In Fig.4, the pulsar period estimation error declines as the number of groups, whilst the running time increases. The reason is that these distortion degree errors from one pulsar photon sequence are incompletely correlated, which is proved in Section 4. Namely, the available information increases as the number of groups. On the other hand, the computational load is approximately proportional to the number of groups, which is proved in Section 5. Therefore, we must make a compromise between the computational load and the estimation error of the pulsar period. In this paper, we choose 32 as the optimal number of groups. Besides, according to actual situation, readers can select other values as the optimal number of groups.   We compare the FBEF-FDP with the SDP and the χ 2 algorithm. Table 2 shows the simulation results. The accuracy of the FBEF-FDP is higher than the SDP at the cost of the running time. In practice, the users can select the optimal number of groups to find a compromise between accuracy and computational load. Besides, the performance of the FBEF-FDP excels that of the χ 2 algorithm.

C. CORRELATION COEFFICIENT
We utilize the same pulsar photon sequence to accumulate multiple distorted pulsar profiles, and then estimate the distortion degrees of these accumulative pulsar profiles. Only when these distortion degree errors are incompletely correlated, the nonlinear search method for multiple minimums can obtain more highly accurate distortion degree estimation. Fig.5 shows two correlation coefficient curves. One is the correlation coefficients between the errors of zero distortion degree and the errors of others. The other one is the correlation coefficients between the errors of 63 distortion degree and the errors of others. In Fig.5, with the increase of the difference between two distortion degrees, the correlation coefficient declines first and then stabilizes. These results demonstrate that the estimation errors of these distorted pulsar profiles are incompletely correlated. This is consistent with the theoretical analysis in Section 5. Therefore, we fit incompletely correlated estimation errors to enhance the estimation accuracy.
Besides, the increase of the number of groups means that more incompletely correlated estimation values. And more estimation information is beneficial to improve the final estimation accuracy. However, due to the limit of computational load, we should select an optimal number of groups.

D. EFFECTIVE AREA AND OBSERVATION TIME
Both the X-ray detector's effective area and the observation time affect the estimation accuracy of the pulsar period. In this section, we investigate the relationship between these two factors and the estimation error of the X-ray pulsar period. The simulation results are shown in Fig.6.
It can be seen from Fig. 6(a) that the accuracy of the X-ray pulsar period estimation is inversely proportional to the X-ray detector's effective area and the observation time. Therefore, increasing the observation time and the detector area can reduce the pulsar estimation error.
For the convenience of analysis, the estimation error is studied in three typical observation periods of 100s, 1000s and 10000s. As shown in Fig.6 (b), it can be seen that when the observation time is extended by two orders of magnitude, the accuracy of the period estimation is improved by more than two orders of magnitude. And when the effective area of the X-ray detector is increased by two orders of magnitude, the accuracy of the period estimation is improved by less than an order of magnitude. Therefore, under the condition that the accuracy of the period estimation is required to be improved, it is more advantageous to improve the estimation accuracy by increasing the observation time than increasing the effective area of X-ray detector.

E. PHOTON FLUX AND BACKGROUND NOISE
The pulsar radiation photon flux and the background noise are two important factors influencing the estimation accuracy of the pulsar period. We study the effect of these two factors on the accuracy of the pulsar period estimation. The results are shown in Fig.7.
It can be seen from Fig. 7 that as the background noise increases, the period estimation error increases. As the pulsar radiation flux increases, the period estimation error decreases. Besides, when the pulsar photon flux is smaller than 1ph/cm 2 /s, the accuracy of the pulsar period estimation is significantly affected by the background noise and the pulsar photon flux. Otherwise, the pulsar period estimation accuracy is high and does not change much. The reason can be given that when the pulsar radiation photon flux is high, the background noise is not the major noise source, and the quantization error caused by the search step becomes the major error source. At this time, the influence of noise on the estimation accuracy of the pulsar period is small.

F. SEARCH STEP SIZE
In this section, we study the relationship between different search step sizes and the accuracy of pulsar estimation, and verify that the relationship between the search step size and the period accuracy is almost consistent when the search step size is small. The simulation conditions of four different situations are shown in Table 3. We take Situation A as a benchmark. Situations B, C and D adjust the area of X-ray sensor, the observation and the photon flux, respectively. Figure 8 shows the simulation results. It can be seen from Figure 8 that with the decrease of the search step size, VOLUME 8, 2020   the accuracy of the period error first decreases and then presents a steady trend. This is because when the search step size is reduced to a certain value, the quantization error caused by the search step size is ignorable. In Figure 8, when the search step size is equal to 0.01, it is the critical point of accuracy change for four situations shown in Table 3. Therefore, we choose 0.01 as the optimal search step size.

G. INITIAL VALUES OF THE PULSAR PERIOD
In this section, we investigate the relationship between the initial pulsar period and the estimation accuracy. The simulation results are as Table 4.
As can be seen from the above table, as the pulsar period increases, the accuracy of the pulsar period estimation increases or decreases irregularly. Therefore, the initial pulsar period does not affect the accuracy of the period estimation.

H. SIZES OF THE DISTORTED PROFILE DICTIONARY
In this section, we investigate the relationship between the distorted profile dictionary's size and the performance of the FBEF-FDP. The simulation results are as Table 5. The size of the distorted profile dictionary is greater than the value range of the estimation period.
As can be seen from the above table, as the size of the distorted profile dictionary increases, the accuracy of the pulsar period estimation is almost unchanged, and the running time increases. Therefore, the small distorted profile dictionary is benefit under the premise that the size of the distorted profile dictionary is greater than the value range of the estimation period.

VII. CONCLUSION
Considering the fact that the single distorted profile-based pulsar period estimation still has space for accuracy to improve, we combine it with the fast butterfly epoch folding, and propose the FBEF-based pulsar period estimation with a few distorted profiles. The FBEF-FDP method can adjust the number of groups to enhance accuracy with the computational limit of on-board computer. Besides, theoretically, we analyze the correlation of the pulsar noises and the computational complexity. Simulation results manifest the feasibility and superiority of the FBEF-FDP method.
The FBEF-FDP method has following virtues: 1. High accuracy. We prove that although different accumulative distorted profiles stem from the same pulsar photon sequence, distortion degree estimation errors are incompletely correlated. This is the theorical basis that the distortion degree estimation values from the same pulsar photon sequence is fitted to obtain more highly accurate distortion degree estimation. The single distorted profile-based pulsar period estimation merely explores the information of one distorted profile, whilst the FBEF-FDP method synthesizes multiple distorted profiles. Namely, compared with the single distorted profile-based one, the FBEF-FDP method more deeply explores the received pulsar photon sequence, and get higher accuracy.
2. Adjustable computational load. The number of groups is an important factor for both accuracy and computational load. Adjustable number of groups is beneficial to find a compromise between accuracy and computational load. The computational complexity analysis results manifest that the computational load is approximately proportional to the number of groups. In fact, to enhance accuracy, we can select the maximal number of groups with the computational limit of on-board computer.
3. Pure real computation. The FBEF-FDP method excludes imaginary number. Thus, its computation process is simple, which favors on-board computer.
The FBEF-FDP method requires a pulsar profile with high SNR. However, in fact, the obtainment of the pulsar profile with high SNR requires a great mass of observation data. Our next work plan is to research a novel algorithm without pulsar profile. Besides, the pulsar period evaluation without real period is a problem deserving of study. ZHI-WEI KANG was born in Xiangtan, China in November 1962. He received the Ph.D. degree from Hunan University, in 2007. He is currently a Professor with Hunan University. He is mainly working on signal and information processing. VOLUME 8, 2020