- Email: [email protected]
A fast pulse time-delay estimation method for X-ray pulsars based on wavelet-bispectrum
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Contents lists available at ScienceDirect
Optik journal homepage: www.elsevier.com/locate/ijleo
Original research article
A fast pulse time-delay estimation method for X-ray pulsars based on wavelet-bispectrum Yangkun Wua, Zhiwei Kanga,*, Jin Liub,c a
College of Computer Science and Electronic Engineering, Hunan University, Changsha 410082, China College of Information Science and Engineering, Wuhan University of Science and Technology, Wuhan 430081, China Engineering Research Center for Metallurgical Automation and Measurement Technology of Ministry of Education, Wuhan University of Science and Technology, Wuhan 430081, China b c
A R T IC LE I N F O
ABS TRA CT
Keywords: X-ray pulsar Delay estimation Bispectrum Wavelet X-ray pulsar-based navigation (XPNAV)
X-ray pulsar-based navigation technology provides real-time aerospace information for spacecraft through the received X-ray photons. The time-delay estimation accuracy of the pulse integrated profile is one of the determinants of navigation accuracy. To reduce the computational complexity of the time-delay estimation algorithm for X-ray pulsar integrated pulse profile, we propose a fast time-delay estimation method based on the bispectral algorithm. This method reduces computational complexity by reducing the calculated frequency bands. Wavelet is first used to decompose the pulse profiles and obtain the low-frequency component with a high signalto-noise ratio. Then, the low-frequency component is used to estimate the time delay of the integrated pulse profile via the bispectral method. Theoretical analysis and experimental results show that the bispectrum-based time-delay estimation method can ignore high-frequency components with considerable noise. Moreover, using only the low-frequency component can substantially reduce computational complexity whilst slightly improving estimation accuracy.
1. Introduction X-ray pulsar-based navigation (XPNAV) receives X-ray photons from pulsars, it relies on the X-ray detectors carried by spacecraft, records the photons’ time of arrival (TOA) and autonomously defines the speed, location, time and attitude of spacecraft [1–3]. XPNAV is a promising astronomical autonomous navigation technique that has considerable engineering application value and strategic research significance [4,5]. Pulse TOA is the most basic observation of a pulsar navigation system. In addition to the sensitivity of the X-ray detector and the time conversion model, the measurement accuracy of TOA primarily depends on the accuracy of the time-delay estimation between the standard and integrated pulse profiles. Obtaining the time delay of the integrated pulse profile quickly and accurately through short-time observations is highly significant in improving the performance of navigation systems [6]. Research on the time-delay estimation algorithm for X-ray pulsar integrated pulse profiles is relatively mature [7–17]. The maximum likelihood function has been used to estimate pulse arrival time with high computational cost [8,9]. The cross-correlation algorithm and the nonlinear least variance method used in time-delay estimation were presented in Refs. [10] and [11]. The statistical properties of the observation profile are obtained via epoch folding, and pulse phase is estimated based on fast Fourier transform (FFT) [12]. In Ref. [13], the time-delay estimation of pulsar signals was transformed into a scalar estimation method in the
⁎
Corresponding author. E-mail address: [email protected] (Z. Kang).
https://doi.org/10.1016/j.ijleo.2019.163790 Received 30 August 2019; Accepted 13 November 2019 0030-4026/ © 2019 Elsevier GmbH. All rights reserved.
Please cite this article as: Yangkun Wu, Zhiwei Kang and Jin Liu, Optik - International Journal for Light and Electron Optics, https://doi.org/10.1016/j.ijleo.2019.163790
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Y. Wu, et al.
time domain, and an artificial neural network and the particle filter algorithm were used to estimate the delay of pulsar signals in real time. In Ref. [14], the vehicle's orbital dynamics is used to estimate the phase of X-ray pulsars' pulse, and the pulse signal model is linearized to the second order around vehicle's position and velocity, which can well estimate the phase of pulsars with low-flux pulsars. Third-order cross-wavelet cumulants were used to suppress the effects of scale expansion and noise, and time delay was accurately measured via parabolic interpolation; the computational complexity of the algorithm was considerably reduced [15]. In Ref. [16], two eigenvectors of the observation equation were obtained by constructing the equation using the cross-correlation results of the standard and integrated profiles, and time delay was estimated by using its orthogonality. A time-delay estimation algorithm based on the optimal frequency band was proposed in Ref. [17], by analysing the frequency-domain features of the pulse profile, it is concluded that the high-frequency part with large noise interference can be completely discarded and we can only use the lowfrequency part with high signal-to-noise ratios (SNRs) when estimating time delay, and the optimal frequency band was determined by the relationship between the error of time-delay estimation and the frequency band used at different SNRs. Classical FFT was used to measure the phase difference of the pulse profile in the frequency domain [18]. The measurement accuracy of the time-delay estimation method in the frequency domain is independent of the temporal resolution of the integrated profile, which is only related to the SNR of the integrated pulse profile, thereby making the Taylor FFT algorithm the most popular time-delay estimation algorithm for radio astronomy. However, the Taylor FFT algorithm exhibits poor performance at low SNRs. In Ref. [19], Xie et al. proposed a bispectrum-based pulsar time-delay estimation algorithm. Given that the bispectrum can effectively suppress Gaussian noise, an estimation accuracy that is better than that using the Taylor FFT algorithm can be obtained at low SNRs. However, the bispectrumbased algorithm requires a large amount of computation due to the 2D Fourier transform in the bispectrum and the multiple iterations in the delay calculation. Excessive computation results in the long execution time of the bispectrum algorithm. A spacecraft’s computing power is limited, and navigation requires high real-time performance. Therefore, reducing the amount of calculation whilst maintaining accuracy is the key to improving the algorithm [20]. To reduce the computational complexity of the bispectrum-based pulsar time-delay estimation algorithm, this work analyses the bispectrum features of the standard and integrated pulse profiles in the XPNAV system and the spectral features of the integrated pulse profile and proposes a fast pulsar time-delay estimation algorithm based on the bispectrum with less computation. This algorithm uses only low-frequency components that contain most of the information of the pulse profile, compares the self-bispectrum and cross-bispectrum and then calculates the extreme points to estimate time delay. The experiments based on the real data obtained by the Rossi X-ray Timing Explorer (RXTE) satellite and the simulation data show that the algorithm can effectively suppress the influence of Gaussian noise and considerably reduced the calculation amount with high estimation accuracy. 2. X-ray pulsar signal profile model The principle of obtaining the XPNAV signal profile lies in the number of arriving X-ray photons within a given period following the non-homogeneous Poisson distribution. The X-ray detector on a spacecraft captures the photons emitted by the pulsar and records the time that the photons reach the detector [10,21]. Pulse epoch folding is used to obtain the pulse profile [22–24]; that is, using the time model based on the centroid of the sun as reference, the pulsar period is divided into several equal-length phase intervals, and each photon in the multiple observation periods is placed in the corresponding phase interval according to its arrival. 2.1. X-ray pulsars’ pulse profile The standard pulse profile of X-ray pulsars is obtained through long-time astronomical observations with a very high SNR. Using pulsar B0531+21 (Crab pulsar) as an example, we create its standard pulse profile using real data obtained by the RXTE satellite shown in Fig. 1(a) (the observation time is 13,580 s and the number of phase bins is 1024). The integrated pulse profile is obtained by
Fig. 1. Standard pulse profile (a) and integrated pulse profile (b) of pulsar B0531+21. 2
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Y. Wu, et al.
Fig. 2. Amplitude spectrum of the integrated pulse profile.
the photon received by the detector within a short time. Using the data from the same period and the same periodic superposition method as that described in the preceding paragraph, the integrated pulse profile shown in Fig. 1(b) (the observation time is 100 s and the number of phase bins is 1024) is obtained. 2.2. Frequency-domain features of pulse profile The frequency spectrum of the X-ray pulsar’s (B0531+21) integrated pulse profile (100 s) is shown in Fig. 2. Since the normalized amplitude when frequency greater than 100 Hz is too small and it affects viewing, the frequency in Fig. 2 is only up to 100 Hz. The convergence of the pulse profile spectrum is observed and it is mostly concentrated in the low-frequency segment, thereby indicating that the energy of the signal is concentrated in the low-frequency part; the high-frequency part has extremely low amplitude and the spectrum line is flat, and with considerable noise interference. Thus it is affected much more by noise than low-frequency part [17]. 3. Fast time-delay estimation based on bispectrum 3.1. Pulse profile’s features on the bispectral domain The standard pulse profile bispectrum map of pulsar B0531+21 is also obtained as shown in Fig. 3, where (a) is the slice diagram of the bispectrum of the standard pulse profile (ω1is fixed). The bispectral value in the low-frequency part is considerably greater than that in the high-frequency part. That is, the bispectrum of the pulse profile is mostly concentrated in the low-frequency part. Meanwhile, (b) shows the bispectral contour value map of the standard pulse profile, and the bispectrum of the pulse profile is also concentrated in the low-frequency part. By combining the amplitude spectrum of the pulse profile shown in Fig. 2, the bispectral information of the pulse profile is
Fig. 3. Bispectral map of the standard pulse profile of pulsar B0531+21: (a) bispectrum slice diagram and (b) bispectral contour map. 3
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Y. Wu, et al.
Fig. 4. Energy spectrum of the third-order wavelet of the pulse profile.
concentrated in the low-frequency component. The amplitude of the high-frequency component is small, and the high-frequency components contain noise that affects the accuracy of time-delay estimation. Thus, it is considerably affected by noise, the phase changes in the high-frequency band bispectrum cannot be accurately estimated. Consequently, this study uses wavelet to decompose the pulse profile, obtains the low-frequency components containing the bispectral information of the signal and discards the highfrequency components. This method can substantially reduce the calculation of the algorithm and improve its speed whilst ensuring estimation accuracy. From the energy perspective [25], the energy spectrum of the third-order wavelet of the X-ray pulse profile is shown in Fig. 4. The energy is concentrated in the first frequency band. The signal energy after wavelet decomposition obtained through calculations is 99.9 % of the original signal, with nearly no loss. Therefore, the low-frequency part of the pulse profile obtained via wavelet decomposition can be used to estimate time delay. Since the integrated pulsar signal is a non-stationary signal, the wavelet decomposition can express the time-frequency local properties better than Fourier transform when obtain the lowfrequency part. And wavelet transform has better frequency resolution in low frequency band. Subsequent experiments have also proved that the algorithm using wavelet to decompose pulse profile has better performance. 3.2. Fast time-delay estimation based on bispectrum The bispectrum-based time-delay estimation algorithm uses the third-order statistical characteristic. In theory, the bispectrum of the Gaussian process is zero, and thus, the bispectrum algorithm can effectively suppress Gaussian noise [26,27]. Therefore, high estimation accuracy can be achieved when the observation time of the pulse profile is short [19]. Let s(n) = {s(0),s(1),…,s(N−1)} and p(n) = {p(0),p(1),…,p(N−1)} be the standard and integrated pulse profiles, respectively. The results of the wavelet decomposition of s(n) and p(n) are:
SAj = < s, ϕj > ; SDj = < s, ψj >
(1)
PAj
(2)
= < p, ϕj > ; j
PDj
= < p, ψj >
j
where SA and PA represent the low-frequency part of the j-th order decomposition of the standard and integrated profiles, respectively; SDj and PDj represent the high-frequency part of the j-th order decomposition of the pulse profile; ϕ is the scale function and ψ is the wavelet function. We select the low-frequency part for the time-delay estimation and let XAj (ω) and YAj (ω) be the discrete Fourier transform of SAj and PAj, respectively. The self-bispectrum of the standard pulse profiles can be obtained as follows: *
Bsss (ω1, ω2 ) Aj = XAj (ω1)*XAj (ω2)*XAj (ω1 + ω2)
(3)
Meanwhile, the cross-bispectrum of the standard and integrated pulse profiles is as follows: *
Bsps (ω1, ω2 ) Aj = XAj (ω1)*YAj (ω2)*XAj (ω1 + ω2)
(4)
where ω1 = 0,1,…,N/(2^(j + 1)) and ω2 = 0,1,…,N/(2^(j + 1)). In accordance with the Taylor FFT algorithm [18], the squared cumulative error between the self-bispectrum and the crossbispectrum is defined as [19]
χ 2 (τ ) =
∑ ∑ ||Bsps (ω1, ω2 ) Aj| − |Bsss (ω1, ω2 ) Aj | *exp(iθ (ω1, ω2) − iφ (ω1, ω2) − i2π ω1*τ /N )|2 ω1
(5)
ω2
where θ (ω1, ω2) and φ (ω1, ω2) are the phases of the cross-bispectrum 4
Bsps (ω1, ω2 ) Aj
and self-bispectrum
Bsss (ω1, ω2 ) Aj ,
respectively,
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Y. Wu, et al.
Fig. 5. Flowchart of fast bispectrum-based algorithm.
and τ is the assumed delay of the integrated pulse profiles. The value of χ 2 (τ ) is smaller and the value of ∂χ 2 (τ )/ ∂τ is closer to zero when τ is closer to the real delay D:
∂χ 2 (τ )/ ∂τ = −4π / N *∑ ω1∑ (|Bsps (ω1, ω2 ) Aj|*|Bsss (ω1, ω2 ) Aj|*sin(θ (ω1, ω2) − φ (ω1, ω2) − 2π ω1τ / N )) ω1
(6)
ω2
Thus, Eq. (6) can be used to estimate time delay. The bispectrum algorithm is unaffected by the Gaussian noise in the pulsar signal. Thus, the algorithm’s accuracy is better than that of the classic Taylor FFT algorithm. However, the algorithm requires excessive computation because the bispectrum is a 2D function. Thus, the amount of calculation using Eq. (6) is proportional to the square of the data’s length. The use of wavelet decomposition to obtain the low-frequency part of the pulse contour can reduce the number of segments in one period, thereby reducing the computational complexity of the algorithm. The effect of noise on estimation accuracy in the high-frequency part is also eliminated. If we further reduce the number of segments in one period, then the calculation of the algorithm will decrease, but will also lead to a decrease in estimation accuracy. Thus, the higher the estimation accuracy, the higher the number of segments in one period will be. However, the amount of calculation will increase with a quadric function. Therefore, experiments should be conducted and the results should be analysed to consider estimation accuracy and reasonable running time and memory. Fig. 5 presents a flowchart of the method. 3.3. Computational complexity analysis of this method For the traditional bispectrum algorithm, a large amount of calculation is required because two iterations are necessary to calculate the square cumulative error in accordance with Eqs. (5) and (6). Therefore, if one period of the pulsar pulse profile is divided into N small bins, N2/2 additions, N2/2+N/2 multiplications and N2/4 trigonometric operations are required to calculate the square cumulative error after we obtain the bispectrum. Obtaining the pulse profile by n-order wavelet decomposition requires 2*N*(1-2-n) multiplications and 2*(N-1)*(1-2-n) additions [28]. When the number of order is 1, a low-frequency component with N/2 bins is obtained, then only N2/8 additions, N2/8 + N/4 multiplications and N2/16 trigonometric additions are required to calculate the square cumulative error. The computation is approximately 1/4 that in the traditional bispectrum algorithm, thereby considerably reducing the running time of the algorithm. After calculating the bispectrum of the profile, the traditional algorithm will obtain an N × N bispectrum matrix, and the size of the matrix obtained after the decomposition is only N/2 × N/2. Thus, the storage and processing of the bispectrum can be reduced. If we improve the wavelet order, then the low-frequency part with fewer segments will be obtained, and the amount of calculation and running time of the algorithm will be reduced further. 4. Experiments and analysis To verify the performance of this method, three pulsars, namely, B0531+21, B1821−24 and B1937+21, were selected as the navigation star. The experimental data of pulsar B0531+21 were obtained from the RXTE satellite, taken from the American High Energy Astrophysics Science Archive Research Centre [29], package number is 40805-01-05-00; the experimental data of pulsars B1821−24 and B1937+21 were obtained from the European Pulsar Network Database [30]. The standard pulse profile of simulation data is obtained by normalizing the data from EPN database directly and the integrated pulse profile is generated by Poisson distribution of the standard pulse profile. Follow-up experiments were conducted on a notebook computer with i5 dual-core and a main frequency of 2.6 GHz. The simulation environment was MATLAB R2015b. 4.1. Influence of different wavelet basis functions on accuracy Different wavelet basis functions may affect estimation accuracy. The commonly used wavelet basis functions are the Haar, Symlet, Coiflet, Biorthogonal and Meyer wavelets [31,32]. These five wavelet basis functions are selected to perform the preceding experiments on the three pulsars. The results are presented in Figs. 6 and 7. The experimental results show that using different wavelet basis functions exerts minimal effect on the accuracy of time-delay estimation. From the running speed perspective, the running time of each wavelet basis function is similar. The running time of Haar wavelet is slightly shorter than that of other wavelet basis functions, and Haar wavelet has the advantages of simplicity, orthogonality and symmetry. So we use the Haar wavelet in the previous experiments. 5
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Y. Wu, et al.
Fig. 6. Accuracy of different wavelet basis functions: (a) B0531+21, (b) B1821−24 and (c) B1937+21.
Fig. 7. Running time of different wavelet basis functions: (a) B0531+21, (b) B1821−24.
4.2. Influence of order of wavelet decomposition on accuracy and the selection of order The order of wavelet decomposition may affect estimation accuracy. The higher the order, the less the computation of the algorithm will be. However, it also causes the loss of bispectral information, which reduces accuracy. The experiments indicate the number of orders that wavelet decomposition can obtain for the best estimation accuracy. The three pulsars above are tested, and the number of pulse profile bins is set to 1024. The result is presented in Fig. 8. The experimental results show that when the order does not exceed two, the estimation error is small and precision is high; when the order exceeds two, the estimation error increases; when the order exceeds three, the estimation error increases sharply, and thus, is not shown in the figure. This finding is attributed to the bispectral information of the pulse profile being partially lost, which decreases the accuracy of the time-delay estimation. Fig.9 shows that when the wavelet order increases, the estimation error decreases first and then increases, and the running time decreases all the time. When the order is 2, the error is smallest and the running time is much shorter than that of the classical bispectrum algorithm. Obtaining the optimal order balanced good estimation accuracy and short running time by experiments and recorded in the database. And later experiments can directly obtain the optimal order (showed in Fig. 5).
Fig. 8. Accuracy of different wavelet orders: (a) B0531+21, (b) B1821−24 and (c) B1937+21. 6
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Y. Wu, et al.
Fig. 9. Influence of different wavelet orders on the performance of the algorithm: (a) B0531+21, (b) B1821−24 and (c) B1937+21.
4.3. Computational complexity analysis of the proposed method The following section presents a comparison of the estimation accuracy of our method and that of the traditional bispectrum algorithm and Taylor FFT algorithm. The period of the B0531+21 pulse is 33.5 ms and the number of phase bins is 1024. From Section 2.1, the standard pulse profile is obtained by integrating the data of the entire packet (integrated time: 13,580 s). The integrated pulse profile is obtained by adding time delay to the real data of different integrated times. The results are presented in Fig. 10 after 100 Monte-Carlo experiments. Fig. 10(a) and (b) show the estimation error of the algorithm when the integrated time is 1–10 s and 10–100 s, respectively (real data), and Fig. 10(c) shows the estimation error of the algorithm when the integrated time is 0.1–1 s (simulation data, the area of detector is 100 cm2). Comparing the proposed method with traditional bispectrum algorithm, Taylor FFT algorithm, generalized cross-correlation method and fast transverse recursive least squares method(FTRLS) [33], the picture shows that the accuracy of this method is better than that of other algorithm and similar to that of the traditional bispectrum method. The accuracy of the estimation using the second-order decomposition of the pulse profile by wavelet is slightly better than those of the traditional bispectrum algorithm and the first-order decomposition of the profile by wavelet when the integrated time is longer than 10 s. The simulation also exhibits similar results, the accuracy of the proposed method is better than the traditional bispectrum and Taylor FFT algorithm, and is similar to that of the method based on optimal frequency band proposed in Ref. [17], thereby confirming that the proposed method is also feasible for the simulation data. To verify the generality of the method, the same experiments are performed on low-flow pulsars B1821−24 and B1937+21. The results are shown in Fig. 11. The flow density of B1821−24 is 1.93 × 10−4 (ph/cm2/s), and the period is 3.05 ms. The flow density of B1937+21 is 4.99 × 10−5 (ph/cm2/s), and the period is 1.56 ms. The number of bins is set to 1024. The estimation accuracy of the second-order decomposed pulse profile is also slightly better than those of the traditional bispectrum algorithm, the first-order decomposition of the profile and Taylor FFT algorithm, which indicates that using the low-frequency part can achieve high estimation accuracy.
4.4. Calculation of different algorithms Running time is used as a measure for the amount of computation and running speed of the algorithm. We performed 100 MonteCarlo experiments to compare the running time of the three pulsars using our method, the traditional bispectrum algorithm and Taylor FFT algorithm. The average running time of each experiment is shown in Fig. 12. The experimental results indicate that the computational complexity of the traditional bispectrum algorithm, Taylor FFT algorithm and our method is independent of the pulse integrated time. Our method can effectively reduce the computational complexity and improve the speed of time-delay estimation. The running speed of the pulse profile after the second-order decomposition is further improved. The specific running times of these algorithm and the generalized cross-correlation method are shown in Table 1. Although the generalized cross-correlation has shorter running time, its error of time-delay estimation is much bigger than our
Fig. 10. Estimated error at different integrated times of B0531+21: (a) 1–10 s, (b) 10–100 s and (c) 0.1–1 s. 7
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Y. Wu, et al.
Fig. 11. Estimated error at different integrated times: (a) B1821−24 and (b) B1937+21.
Fig. 12. Comparison of running time using different methods: (a) B0531+21(1–10 s), (b) B0531+21 (10–100 s), (c) B1821−24, (d) B1937+21 and (e) B0531+21 (0.1–1 s, simulation data).
Table 1 Running times of the proposed method and other algorithm. Pulsar
Bispectrum algorithm/s
Taylor FFT algorithm/s
1st-order decomposed/s
2nd-order decomposed/s
Generalized cross-correlation method/ s
B0531+21 B1821−24 B1937+21
5.9737 5.7775 6.3808
0.0972 0.0954 0.0945
0.6358 0.5098 0.5666
0.0712 0.0614 0.0652
0.00085 0.00089 0.00088
algorithm (see Fig. 10).
4.5. Comparison of the low-frequency part of the bispectrum with the algorithm in this study The proposed algorithm is compared with the method of directly estimating time delay using the low-frequency part of the bispectrum. The error and running time results are shown in Fig. 13. The error in time-delay estimation using the low-frequency part of the bispectrum is slightly smaller than that using the classical bispectrum algorithm, but still larger than the error using the second8
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Y. Wu, et al.
Fig. 13. Comparison of the accuracy and time of the low-frequency bispectrum and the algorithm proposed in this study: (a), (b) B0531+21 and (c), (d) B1821−24.
order wavelet decomposition. The running times of both methods are considerably shorter than that of the classical bispectrum algorithm, and the running time of the second-order wavelet decomposition is shorter. 4.6. Effect of noise on estimation accuracy Background noise will also affect the accuracy of X-ray pulsar time-delay estimation. The comparison of the estimation errors
Fig. 14. Accuracy of time delay-estimation under different noises: (a) B1821−24 and (b) B1937+21. 9
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Y. Wu, et al.
Fig. 15. Influence of different decomposition methods on the accuracy of estimation: (a) real data and (b) simulation data.
under different background noises is shown in Fig. 14. Experiments are performed using pulsar B1821−24 and B1937+21 pulse profiles via second-order wavelet decomposition. The value of background noise exerts a considerable influence on the accuracy of the algorithm. The larger the background noise is, the greater the estimation error will be, thereby reducing estimation accuracy. 4.7. Comparison of Fourier transform and wavelet transform Using Fourier transform to obtain the low-frequency part of the signal is also a common method. Fig.15 shows the influence of the method and wavelet transform on the accuracy of the estimation algorithm. The real data and the simulated data of the pulsar B0531+21 are used in the experiments, the cut-off frequency of Fourier transform is 512 Hz. The experimental results show that the effect of Fourier transform is similar to that of the original algorithm, and the effect of wavelet transform is slightly better than that of Fourier transform. And the number of segments in one period doesn’t reduce after Fourier transform, so there is no reduction in the computational complexity compared with the original algorithm. So the wavelet transform is superior to the Fourier transform in improving the performance of the algorithm. 5. Conclusion In this study, we found that high-frequency components with minimal bispectrum values and considerable noise interference can be disregarded when estimating time delay based on bispectrum by analysing the bispectral characteristics of the X-ray pulsar pulse profile. In addition, we can only use the low-frequency components with high bispectrum values to conduct the experiment. Wavelet transform is used to decompose the pulse profile, and the low-frequency components of the profile are obtained, which not only maintain the bispectrum of the pulse profile, but also effectively reduce computational complexity. The experimental results show that this method can considerably reduce computational complexity and slightly increase the accuracy of pulse profile time-delay estimation. Funding This work is supported by the National Natural Science Foundation of China under Grant nos. 61772187 and 61873196. References [1] S.I. Sheikh, The Use of Variable Celestial X-Ray Sources for Spacecraft Navigation Dissertation, University of Maryland, Maryland, 2005. [2] J. Liu, X.L. Ning, X. Ma, J.C. Fang, Geometry error analysis in solar doppler difference navigation for the capture phase, IEEE Trans. Aeros. Electron. Syst. 55 (5) (2019) 2556–2567. [3] X.L. Ning, M.Z. Gui, J.C. Fang, G. Liu, Differential X-ray pulsar aided celestial navigation for mars exploration, Aerosp. Sci. Technol. 62 (1) (2017) 36–45. [4] J. Liu, J.C. Fang, G. Liu, Fractional differentiation-based observability analysis method for nonlinear X-ray pulsar navigation system, Proc. Inst. Mech. Eng. Part G: J. Aerosp. Eng. 232 (8) (2018) 1467–1478. [5] J. Liu, J.C. Fang, Z.W. Kang, J. Wu, X.L. Ning, Novel algorithm for X-ray pulsar navigation against doppler effects, IEEE Trans. Aeros. Electron. Syst. 51 (1) (2015) 228–241. [6] Z.W. Kang, H.C. He, J. Liu, Adaptive pulsar time delay estimation using wavelet-based RLS, Optik 171 (2018) 266–276. [7] Y.D. Wang, Z. Wei, S.M. Sun, L. Li, X-ray pulsar-based navigation using time-differenced measurement, Aerosp. Sci. Technol. 36 (2014) 27–35. [8] A.A. Emadzadeh, J.L. Speyer, Asymptotically efficient estimation of pulse time delay for X-ray pulsar based relative navigation, Aiaa Guidance, Navigation, & Control Conference (2009). [9] J.X. Li, X.Z. Ke, Maximum-likelihood TOA estimation of X-ray pulsar signals on the basis of Poison model, Chin. Astron. Astrophys. 35 (1) (2011) 19–28. [10] A.A. Emadzadeh, J.L. Speyer, On modeling and pulse phase estimation of X-ray pulsarsr, IEEE Trans. Signal Process. 58 (9) (2010) 4484–4495. [11] A.A. Emadzadeh, A.R. Golshan, J.L. Speyer, Consistent estimation of pulse delay for X-ray pulsar based relative navigation, Joint 48th IEEE Conference on Decision and Control (2009) 1488–1493.
10
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Y. Wu, et al.
[12] M.F. Xue, X.P. Li, H.F. Sun, H.Y. Fang, A fast pulse phase estimation method for X-ray pulsar signals based on epoch folding, Chin. J. Aeronaut. 29 (3) (2016) 746–753. [13] P.F. Li, G.D. Xu, L.M. Dong, T.R. Hou, A real time estimation method of time-delay for X-ray pulsar signal, Acta Aeronaut. Astronaut. Sin. 35 (7) (2014) 1966–1976. [14] Y.D. Wang, Z. Wei, Pulse phase estimation of X-ray pulsar with the aid of vehicle orbital dynamics, J. Navig. 69 (2) (2016) 414–432. [15] Z. Su, L.P. Xu, T. Wang, H. Zhang, A new time delay measurement algorithm for pulsar accumulated pulse profile, J. Astronaut. 32 (6) (2011) 1256–1261 [Chinese]. [16] Q.Q. Lin, P. Shuai, L.W. Huang, X.Y. Zhang, X.M. Bei, A high-precision TOA estimation method for X-ray pulsar-based navigation system, J. Astronaut. 37 (6) (2016) 704–710 [Chinese]. [17] H.Y. Fang, B. Liu, X.P. Li, et al., Time delay estimation method of X-ray pulsar observed profile based on the optimal frequency band, Acta Phys. Sin. 65 (11) (2016) 119701. [18] J.H. Taylor, Pulsar timing and relativistic gravity, Philos. Trans. Phys. Sci. Eng. 341 (1660) (1992) 117–134. [19] Z.H. Xie, L.P. Xu, G.R. Ni, Time offset measurement algorithm based on bispectrum for pulsar integrated pulse profile, Acta Phys. Sin. 57 (10) (2008) 6683–6688. [20] Z.W. Kang, X. He, J. Liu, X-ray pulsars time delay estimation using GSO-based bispectral feature points, Opt. Int. J. Light Electron Opt. 127 (2016) 5050–5054. [21] Y.D. Wang, W. Zheng, D.P. Zhang, X-ray pulsar-based navigation method using one sensor and modified time-differenced measurement, Proc. Inst. Mech. Eng. Part G: J. Aerosp. Eng. 233 (1) (2017) 299–309. [22] H.F. Sun, X. Sun, H.Y. Fang, et al., Building X-ray pulsar timing model without the use of radio parameters, Acta Astronaut. 143 (2018) 155–162. [23] A.A. Emadzadeh, J.L. Speyer, X-ray pulsar-based relative navigation using epoch folding, IEEE Trans. Aerosp. Electron. Syst. 47 (4) (2011) 2317–2328. [24] Y. Wang, W. Zhang, Pulsar phase and Doppler frequency estimation for XNAV using on-orbit epoch folding, IEEE Trans. Aerosp. Electron. Syst. 52 (5) (2017) 2210–2219. [25] J. Liu, Z.H. Yang, Z.W. Kang, X. Chen, Fast CS-based pulsar period estimation method without tentative epoch folding and its CRLB, Acta Astronaut. 160 (2019) 90–100. [26] Z.H. Xie, L.P. Xu, G.R. Ni, Pulsar signal recognition based on one-dimension selected line spectra, J. Infrared Millim. Waves 26 (3) (2007) 187–195 [Chinese]. [27] X.D. Zhang, Modern Signal Processing, 2nd ed., Tsinghua University Press, Beijing, 2003, pp. 328–335 [Chinese]. [28] O. Rioul, P. Duhamel, Fast algorithms for discrete and continuous wavelet transforms, IEEE Trans. Inf. Theory 38 (2) (1992) 569–586. [29] NASA’s HEASARC Data Access. [Cited 28 December 2017]. Available from: https://heasarc.gsfc.nasa.gov/docs/archive.html. [30] EPN (European Pulsar Network). [Cited 28 December 2017]. Available from: http://www.jb.man.ac.uk/research/pulsar/Resources/epn/browser.html. [31] A.S. Yaseen, O.N. Pavlova, A.N. Pavlov, Noisy signal filtration using complex wavelet basis sets, Tech. Phys. Lett. 43 (7) (2017) 641–644. [32] X.P. Liu, W. Yuan, L.L. Han, et al., X-ray pulsar signal de-noising for impulse noise using wavelet packet, Aerosp. Sci. Technol. 64 (2017) 147–153. [33] L. Ljung, M. Morf, D. Falconer, Fast calculation of gain matrices for recursive estimation schemes, Int. J. Control 27 (1) (1978) 1–19.
11
Contents lists available at ScienceDirect
Optik journal homepage: www.elsevier.com/locate/ijleo
Original research article
A fast pulse time-delay estimation method for X-ray pulsars based on wavelet-bispectrum Yangkun Wua, Zhiwei Kanga,*, Jin Liub,c a
College of Computer Science and Electronic Engineering, Hunan University, Changsha 410082, China College of Information Science and Engineering, Wuhan University of Science and Technology, Wuhan 430081, China Engineering Research Center for Metallurgical Automation and Measurement Technology of Ministry of Education, Wuhan University of Science and Technology, Wuhan 430081, China b c
A R T IC LE I N F O
ABS TRA CT
Keywords: X-ray pulsar Delay estimation Bispectrum Wavelet X-ray pulsar-based navigation (XPNAV)
X-ray pulsar-based navigation technology provides real-time aerospace information for spacecraft through the received X-ray photons. The time-delay estimation accuracy of the pulse integrated profile is one of the determinants of navigation accuracy. To reduce the computational complexity of the time-delay estimation algorithm for X-ray pulsar integrated pulse profile, we propose a fast time-delay estimation method based on the bispectral algorithm. This method reduces computational complexity by reducing the calculated frequency bands. Wavelet is first used to decompose the pulse profiles and obtain the low-frequency component with a high signalto-noise ratio. Then, the low-frequency component is used to estimate the time delay of the integrated pulse profile via the bispectral method. Theoretical analysis and experimental results show that the bispectrum-based time-delay estimation method can ignore high-frequency components with considerable noise. Moreover, using only the low-frequency component can substantially reduce computational complexity whilst slightly improving estimation accuracy.
1. Introduction X-ray pulsar-based navigation (XPNAV) receives X-ray photons from pulsars, it relies on the X-ray detectors carried by spacecraft, records the photons’ time of arrival (TOA) and autonomously defines the speed, location, time and attitude of spacecraft [1–3]. XPNAV is a promising astronomical autonomous navigation technique that has considerable engineering application value and strategic research significance [4,5]. Pulse TOA is the most basic observation of a pulsar navigation system. In addition to the sensitivity of the X-ray detector and the time conversion model, the measurement accuracy of TOA primarily depends on the accuracy of the time-delay estimation between the standard and integrated pulse profiles. Obtaining the time delay of the integrated pulse profile quickly and accurately through short-time observations is highly significant in improving the performance of navigation systems [6]. Research on the time-delay estimation algorithm for X-ray pulsar integrated pulse profiles is relatively mature [7–17]. The maximum likelihood function has been used to estimate pulse arrival time with high computational cost [8,9]. The cross-correlation algorithm and the nonlinear least variance method used in time-delay estimation were presented in Refs. [10] and [11]. The statistical properties of the observation profile are obtained via epoch folding, and pulse phase is estimated based on fast Fourier transform (FFT) [12]. In Ref. [13], the time-delay estimation of pulsar signals was transformed into a scalar estimation method in the
⁎
Corresponding author. E-mail address: [email protected] (Z. Kang).
https://doi.org/10.1016/j.ijleo.2019.163790 Received 30 August 2019; Accepted 13 November 2019 0030-4026/ © 2019 Elsevier GmbH. All rights reserved.
Please cite this article as: Yangkun Wu, Zhiwei Kang and Jin Liu, Optik - International Journal for Light and Electron Optics, https://doi.org/10.1016/j.ijleo.2019.163790
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Y. Wu, et al.
time domain, and an artificial neural network and the particle filter algorithm were used to estimate the delay of pulsar signals in real time. In Ref. [14], the vehicle's orbital dynamics is used to estimate the phase of X-ray pulsars' pulse, and the pulse signal model is linearized to the second order around vehicle's position and velocity, which can well estimate the phase of pulsars with low-flux pulsars. Third-order cross-wavelet cumulants were used to suppress the effects of scale expansion and noise, and time delay was accurately measured via parabolic interpolation; the computational complexity of the algorithm was considerably reduced [15]. In Ref. [16], two eigenvectors of the observation equation were obtained by constructing the equation using the cross-correlation results of the standard and integrated profiles, and time delay was estimated by using its orthogonality. A time-delay estimation algorithm based on the optimal frequency band was proposed in Ref. [17], by analysing the frequency-domain features of the pulse profile, it is concluded that the high-frequency part with large noise interference can be completely discarded and we can only use the lowfrequency part with high signal-to-noise ratios (SNRs) when estimating time delay, and the optimal frequency band was determined by the relationship between the error of time-delay estimation and the frequency band used at different SNRs. Classical FFT was used to measure the phase difference of the pulse profile in the frequency domain [18]. The measurement accuracy of the time-delay estimation method in the frequency domain is independent of the temporal resolution of the integrated profile, which is only related to the SNR of the integrated pulse profile, thereby making the Taylor FFT algorithm the most popular time-delay estimation algorithm for radio astronomy. However, the Taylor FFT algorithm exhibits poor performance at low SNRs. In Ref. [19], Xie et al. proposed a bispectrum-based pulsar time-delay estimation algorithm. Given that the bispectrum can effectively suppress Gaussian noise, an estimation accuracy that is better than that using the Taylor FFT algorithm can be obtained at low SNRs. However, the bispectrumbased algorithm requires a large amount of computation due to the 2D Fourier transform in the bispectrum and the multiple iterations in the delay calculation. Excessive computation results in the long execution time of the bispectrum algorithm. A spacecraft’s computing power is limited, and navigation requires high real-time performance. Therefore, reducing the amount of calculation whilst maintaining accuracy is the key to improving the algorithm [20]. To reduce the computational complexity of the bispectrum-based pulsar time-delay estimation algorithm, this work analyses the bispectrum features of the standard and integrated pulse profiles in the XPNAV system and the spectral features of the integrated pulse profile and proposes a fast pulsar time-delay estimation algorithm based on the bispectrum with less computation. This algorithm uses only low-frequency components that contain most of the information of the pulse profile, compares the self-bispectrum and cross-bispectrum and then calculates the extreme points to estimate time delay. The experiments based on the real data obtained by the Rossi X-ray Timing Explorer (RXTE) satellite and the simulation data show that the algorithm can effectively suppress the influence of Gaussian noise and considerably reduced the calculation amount with high estimation accuracy. 2. X-ray pulsar signal profile model The principle of obtaining the XPNAV signal profile lies in the number of arriving X-ray photons within a given period following the non-homogeneous Poisson distribution. The X-ray detector on a spacecraft captures the photons emitted by the pulsar and records the time that the photons reach the detector [10,21]. Pulse epoch folding is used to obtain the pulse profile [22–24]; that is, using the time model based on the centroid of the sun as reference, the pulsar period is divided into several equal-length phase intervals, and each photon in the multiple observation periods is placed in the corresponding phase interval according to its arrival. 2.1. X-ray pulsars’ pulse profile The standard pulse profile of X-ray pulsars is obtained through long-time astronomical observations with a very high SNR. Using pulsar B0531+21 (Crab pulsar) as an example, we create its standard pulse profile using real data obtained by the RXTE satellite shown in Fig. 1(a) (the observation time is 13,580 s and the number of phase bins is 1024). The integrated pulse profile is obtained by
Fig. 1. Standard pulse profile (a) and integrated pulse profile (b) of pulsar B0531+21. 2
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Y. Wu, et al.
Fig. 2. Amplitude spectrum of the integrated pulse profile.
the photon received by the detector within a short time. Using the data from the same period and the same periodic superposition method as that described in the preceding paragraph, the integrated pulse profile shown in Fig. 1(b) (the observation time is 100 s and the number of phase bins is 1024) is obtained. 2.2. Frequency-domain features of pulse profile The frequency spectrum of the X-ray pulsar’s (B0531+21) integrated pulse profile (100 s) is shown in Fig. 2. Since the normalized amplitude when frequency greater than 100 Hz is too small and it affects viewing, the frequency in Fig. 2 is only up to 100 Hz. The convergence of the pulse profile spectrum is observed and it is mostly concentrated in the low-frequency segment, thereby indicating that the energy of the signal is concentrated in the low-frequency part; the high-frequency part has extremely low amplitude and the spectrum line is flat, and with considerable noise interference. Thus it is affected much more by noise than low-frequency part [17]. 3. Fast time-delay estimation based on bispectrum 3.1. Pulse profile’s features on the bispectral domain The standard pulse profile bispectrum map of pulsar B0531+21 is also obtained as shown in Fig. 3, where (a) is the slice diagram of the bispectrum of the standard pulse profile (ω1is fixed). The bispectral value in the low-frequency part is considerably greater than that in the high-frequency part. That is, the bispectrum of the pulse profile is mostly concentrated in the low-frequency part. Meanwhile, (b) shows the bispectral contour value map of the standard pulse profile, and the bispectrum of the pulse profile is also concentrated in the low-frequency part. By combining the amplitude spectrum of the pulse profile shown in Fig. 2, the bispectral information of the pulse profile is
Fig. 3. Bispectral map of the standard pulse profile of pulsar B0531+21: (a) bispectrum slice diagram and (b) bispectral contour map. 3
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Y. Wu, et al.
Fig. 4. Energy spectrum of the third-order wavelet of the pulse profile.
concentrated in the low-frequency component. The amplitude of the high-frequency component is small, and the high-frequency components contain noise that affects the accuracy of time-delay estimation. Thus, it is considerably affected by noise, the phase changes in the high-frequency band bispectrum cannot be accurately estimated. Consequently, this study uses wavelet to decompose the pulse profile, obtains the low-frequency components containing the bispectral information of the signal and discards the highfrequency components. This method can substantially reduce the calculation of the algorithm and improve its speed whilst ensuring estimation accuracy. From the energy perspective [25], the energy spectrum of the third-order wavelet of the X-ray pulse profile is shown in Fig. 4. The energy is concentrated in the first frequency band. The signal energy after wavelet decomposition obtained through calculations is 99.9 % of the original signal, with nearly no loss. Therefore, the low-frequency part of the pulse profile obtained via wavelet decomposition can be used to estimate time delay. Since the integrated pulsar signal is a non-stationary signal, the wavelet decomposition can express the time-frequency local properties better than Fourier transform when obtain the lowfrequency part. And wavelet transform has better frequency resolution in low frequency band. Subsequent experiments have also proved that the algorithm using wavelet to decompose pulse profile has better performance. 3.2. Fast time-delay estimation based on bispectrum The bispectrum-based time-delay estimation algorithm uses the third-order statistical characteristic. In theory, the bispectrum of the Gaussian process is zero, and thus, the bispectrum algorithm can effectively suppress Gaussian noise [26,27]. Therefore, high estimation accuracy can be achieved when the observation time of the pulse profile is short [19]. Let s(n) = {s(0),s(1),…,s(N−1)} and p(n) = {p(0),p(1),…,p(N−1)} be the standard and integrated pulse profiles, respectively. The results of the wavelet decomposition of s(n) and p(n) are:
SAj = < s, ϕj > ; SDj = < s, ψj >
(1)
PAj
(2)
= < p, ϕj > ; j
PDj
= < p, ψj >
j
where SA and PA represent the low-frequency part of the j-th order decomposition of the standard and integrated profiles, respectively; SDj and PDj represent the high-frequency part of the j-th order decomposition of the pulse profile; ϕ is the scale function and ψ is the wavelet function. We select the low-frequency part for the time-delay estimation and let XAj (ω) and YAj (ω) be the discrete Fourier transform of SAj and PAj, respectively. The self-bispectrum of the standard pulse profiles can be obtained as follows: *
Bsss (ω1, ω2 ) Aj = XAj (ω1)*XAj (ω2)*XAj (ω1 + ω2)
(3)
Meanwhile, the cross-bispectrum of the standard and integrated pulse profiles is as follows: *
Bsps (ω1, ω2 ) Aj = XAj (ω1)*YAj (ω2)*XAj (ω1 + ω2)
(4)
where ω1 = 0,1,…,N/(2^(j + 1)) and ω2 = 0,1,…,N/(2^(j + 1)). In accordance with the Taylor FFT algorithm [18], the squared cumulative error between the self-bispectrum and the crossbispectrum is defined as [19]
χ 2 (τ ) =
∑ ∑ ||Bsps (ω1, ω2 ) Aj| − |Bsss (ω1, ω2 ) Aj | *exp(iθ (ω1, ω2) − iφ (ω1, ω2) − i2π ω1*τ /N )|2 ω1
(5)
ω2
where θ (ω1, ω2) and φ (ω1, ω2) are the phases of the cross-bispectrum 4
Bsps (ω1, ω2 ) Aj
and self-bispectrum
Bsss (ω1, ω2 ) Aj ,
respectively,
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Y. Wu, et al.
Fig. 5. Flowchart of fast bispectrum-based algorithm.
and τ is the assumed delay of the integrated pulse profiles. The value of χ 2 (τ ) is smaller and the value of ∂χ 2 (τ )/ ∂τ is closer to zero when τ is closer to the real delay D:
∂χ 2 (τ )/ ∂τ = −4π / N *∑ ω1∑ (|Bsps (ω1, ω2 ) Aj|*|Bsss (ω1, ω2 ) Aj|*sin(θ (ω1, ω2) − φ (ω1, ω2) − 2π ω1τ / N )) ω1
(6)
ω2
Thus, Eq. (6) can be used to estimate time delay. The bispectrum algorithm is unaffected by the Gaussian noise in the pulsar signal. Thus, the algorithm’s accuracy is better than that of the classic Taylor FFT algorithm. However, the algorithm requires excessive computation because the bispectrum is a 2D function. Thus, the amount of calculation using Eq. (6) is proportional to the square of the data’s length. The use of wavelet decomposition to obtain the low-frequency part of the pulse contour can reduce the number of segments in one period, thereby reducing the computational complexity of the algorithm. The effect of noise on estimation accuracy in the high-frequency part is also eliminated. If we further reduce the number of segments in one period, then the calculation of the algorithm will decrease, but will also lead to a decrease in estimation accuracy. Thus, the higher the estimation accuracy, the higher the number of segments in one period will be. However, the amount of calculation will increase with a quadric function. Therefore, experiments should be conducted and the results should be analysed to consider estimation accuracy and reasonable running time and memory. Fig. 5 presents a flowchart of the method. 3.3. Computational complexity analysis of this method For the traditional bispectrum algorithm, a large amount of calculation is required because two iterations are necessary to calculate the square cumulative error in accordance with Eqs. (5) and (6). Therefore, if one period of the pulsar pulse profile is divided into N small bins, N2/2 additions, N2/2+N/2 multiplications and N2/4 trigonometric operations are required to calculate the square cumulative error after we obtain the bispectrum. Obtaining the pulse profile by n-order wavelet decomposition requires 2*N*(1-2-n) multiplications and 2*(N-1)*(1-2-n) additions [28]. When the number of order is 1, a low-frequency component with N/2 bins is obtained, then only N2/8 additions, N2/8 + N/4 multiplications and N2/16 trigonometric additions are required to calculate the square cumulative error. The computation is approximately 1/4 that in the traditional bispectrum algorithm, thereby considerably reducing the running time of the algorithm. After calculating the bispectrum of the profile, the traditional algorithm will obtain an N × N bispectrum matrix, and the size of the matrix obtained after the decomposition is only N/2 × N/2. Thus, the storage and processing of the bispectrum can be reduced. If we improve the wavelet order, then the low-frequency part with fewer segments will be obtained, and the amount of calculation and running time of the algorithm will be reduced further. 4. Experiments and analysis To verify the performance of this method, three pulsars, namely, B0531+21, B1821−24 and B1937+21, were selected as the navigation star. The experimental data of pulsar B0531+21 were obtained from the RXTE satellite, taken from the American High Energy Astrophysics Science Archive Research Centre [29], package number is 40805-01-05-00; the experimental data of pulsars B1821−24 and B1937+21 were obtained from the European Pulsar Network Database [30]. The standard pulse profile of simulation data is obtained by normalizing the data from EPN database directly and the integrated pulse profile is generated by Poisson distribution of the standard pulse profile. Follow-up experiments were conducted on a notebook computer with i5 dual-core and a main frequency of 2.6 GHz. The simulation environment was MATLAB R2015b. 4.1. Influence of different wavelet basis functions on accuracy Different wavelet basis functions may affect estimation accuracy. The commonly used wavelet basis functions are the Haar, Symlet, Coiflet, Biorthogonal and Meyer wavelets [31,32]. These five wavelet basis functions are selected to perform the preceding experiments on the three pulsars. The results are presented in Figs. 6 and 7. The experimental results show that using different wavelet basis functions exerts minimal effect on the accuracy of time-delay estimation. From the running speed perspective, the running time of each wavelet basis function is similar. The running time of Haar wavelet is slightly shorter than that of other wavelet basis functions, and Haar wavelet has the advantages of simplicity, orthogonality and symmetry. So we use the Haar wavelet in the previous experiments. 5
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Y. Wu, et al.
Fig. 6. Accuracy of different wavelet basis functions: (a) B0531+21, (b) B1821−24 and (c) B1937+21.
Fig. 7. Running time of different wavelet basis functions: (a) B0531+21, (b) B1821−24.
4.2. Influence of order of wavelet decomposition on accuracy and the selection of order The order of wavelet decomposition may affect estimation accuracy. The higher the order, the less the computation of the algorithm will be. However, it also causes the loss of bispectral information, which reduces accuracy. The experiments indicate the number of orders that wavelet decomposition can obtain for the best estimation accuracy. The three pulsars above are tested, and the number of pulse profile bins is set to 1024. The result is presented in Fig. 8. The experimental results show that when the order does not exceed two, the estimation error is small and precision is high; when the order exceeds two, the estimation error increases; when the order exceeds three, the estimation error increases sharply, and thus, is not shown in the figure. This finding is attributed to the bispectral information of the pulse profile being partially lost, which decreases the accuracy of the time-delay estimation. Fig.9 shows that when the wavelet order increases, the estimation error decreases first and then increases, and the running time decreases all the time. When the order is 2, the error is smallest and the running time is much shorter than that of the classical bispectrum algorithm. Obtaining the optimal order balanced good estimation accuracy and short running time by experiments and recorded in the database. And later experiments can directly obtain the optimal order (showed in Fig. 5).
Fig. 8. Accuracy of different wavelet orders: (a) B0531+21, (b) B1821−24 and (c) B1937+21. 6
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Y. Wu, et al.
Fig. 9. Influence of different wavelet orders on the performance of the algorithm: (a) B0531+21, (b) B1821−24 and (c) B1937+21.
4.3. Computational complexity analysis of the proposed method The following section presents a comparison of the estimation accuracy of our method and that of the traditional bispectrum algorithm and Taylor FFT algorithm. The period of the B0531+21 pulse is 33.5 ms and the number of phase bins is 1024. From Section 2.1, the standard pulse profile is obtained by integrating the data of the entire packet (integrated time: 13,580 s). The integrated pulse profile is obtained by adding time delay to the real data of different integrated times. The results are presented in Fig. 10 after 100 Monte-Carlo experiments. Fig. 10(a) and (b) show the estimation error of the algorithm when the integrated time is 1–10 s and 10–100 s, respectively (real data), and Fig. 10(c) shows the estimation error of the algorithm when the integrated time is 0.1–1 s (simulation data, the area of detector is 100 cm2). Comparing the proposed method with traditional bispectrum algorithm, Taylor FFT algorithm, generalized cross-correlation method and fast transverse recursive least squares method(FTRLS) [33], the picture shows that the accuracy of this method is better than that of other algorithm and similar to that of the traditional bispectrum method. The accuracy of the estimation using the second-order decomposition of the pulse profile by wavelet is slightly better than those of the traditional bispectrum algorithm and the first-order decomposition of the profile by wavelet when the integrated time is longer than 10 s. The simulation also exhibits similar results, the accuracy of the proposed method is better than the traditional bispectrum and Taylor FFT algorithm, and is similar to that of the method based on optimal frequency band proposed in Ref. [17], thereby confirming that the proposed method is also feasible for the simulation data. To verify the generality of the method, the same experiments are performed on low-flow pulsars B1821−24 and B1937+21. The results are shown in Fig. 11. The flow density of B1821−24 is 1.93 × 10−4 (ph/cm2/s), and the period is 3.05 ms. The flow density of B1937+21 is 4.99 × 10−5 (ph/cm2/s), and the period is 1.56 ms. The number of bins is set to 1024. The estimation accuracy of the second-order decomposed pulse profile is also slightly better than those of the traditional bispectrum algorithm, the first-order decomposition of the profile and Taylor FFT algorithm, which indicates that using the low-frequency part can achieve high estimation accuracy.
4.4. Calculation of different algorithms Running time is used as a measure for the amount of computation and running speed of the algorithm. We performed 100 MonteCarlo experiments to compare the running time of the three pulsars using our method, the traditional bispectrum algorithm and Taylor FFT algorithm. The average running time of each experiment is shown in Fig. 12. The experimental results indicate that the computational complexity of the traditional bispectrum algorithm, Taylor FFT algorithm and our method is independent of the pulse integrated time. Our method can effectively reduce the computational complexity and improve the speed of time-delay estimation. The running speed of the pulse profile after the second-order decomposition is further improved. The specific running times of these algorithm and the generalized cross-correlation method are shown in Table 1. Although the generalized cross-correlation has shorter running time, its error of time-delay estimation is much bigger than our
Fig. 10. Estimated error at different integrated times of B0531+21: (a) 1–10 s, (b) 10–100 s and (c) 0.1–1 s. 7
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Y. Wu, et al.
Fig. 11. Estimated error at different integrated times: (a) B1821−24 and (b) B1937+21.
Fig. 12. Comparison of running time using different methods: (a) B0531+21(1–10 s), (b) B0531+21 (10–100 s), (c) B1821−24, (d) B1937+21 and (e) B0531+21 (0.1–1 s, simulation data).
Table 1 Running times of the proposed method and other algorithm. Pulsar
Bispectrum algorithm/s
Taylor FFT algorithm/s
1st-order decomposed/s
2nd-order decomposed/s
Generalized cross-correlation method/ s
B0531+21 B1821−24 B1937+21
5.9737 5.7775 6.3808
0.0972 0.0954 0.0945
0.6358 0.5098 0.5666
0.0712 0.0614 0.0652
0.00085 0.00089 0.00088
algorithm (see Fig. 10).
4.5. Comparison of the low-frequency part of the bispectrum with the algorithm in this study The proposed algorithm is compared with the method of directly estimating time delay using the low-frequency part of the bispectrum. The error and running time results are shown in Fig. 13. The error in time-delay estimation using the low-frequency part of the bispectrum is slightly smaller than that using the classical bispectrum algorithm, but still larger than the error using the second8
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Y. Wu, et al.
Fig. 13. Comparison of the accuracy and time of the low-frequency bispectrum and the algorithm proposed in this study: (a), (b) B0531+21 and (c), (d) B1821−24.
order wavelet decomposition. The running times of both methods are considerably shorter than that of the classical bispectrum algorithm, and the running time of the second-order wavelet decomposition is shorter. 4.6. Effect of noise on estimation accuracy Background noise will also affect the accuracy of X-ray pulsar time-delay estimation. The comparison of the estimation errors
Fig. 14. Accuracy of time delay-estimation under different noises: (a) B1821−24 and (b) B1937+21. 9
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Y. Wu, et al.
Fig. 15. Influence of different decomposition methods on the accuracy of estimation: (a) real data and (b) simulation data.
under different background noises is shown in Fig. 14. Experiments are performed using pulsar B1821−24 and B1937+21 pulse profiles via second-order wavelet decomposition. The value of background noise exerts a considerable influence on the accuracy of the algorithm. The larger the background noise is, the greater the estimation error will be, thereby reducing estimation accuracy. 4.7. Comparison of Fourier transform and wavelet transform Using Fourier transform to obtain the low-frequency part of the signal is also a common method. Fig.15 shows the influence of the method and wavelet transform on the accuracy of the estimation algorithm. The real data and the simulated data of the pulsar B0531+21 are used in the experiments, the cut-off frequency of Fourier transform is 512 Hz. The experimental results show that the effect of Fourier transform is similar to that of the original algorithm, and the effect of wavelet transform is slightly better than that of Fourier transform. And the number of segments in one period doesn’t reduce after Fourier transform, so there is no reduction in the computational complexity compared with the original algorithm. So the wavelet transform is superior to the Fourier transform in improving the performance of the algorithm. 5. Conclusion In this study, we found that high-frequency components with minimal bispectrum values and considerable noise interference can be disregarded when estimating time delay based on bispectrum by analysing the bispectral characteristics of the X-ray pulsar pulse profile. In addition, we can only use the low-frequency components with high bispectrum values to conduct the experiment. Wavelet transform is used to decompose the pulse profile, and the low-frequency components of the profile are obtained, which not only maintain the bispectrum of the pulse profile, but also effectively reduce computational complexity. The experimental results show that this method can considerably reduce computational complexity and slightly increase the accuracy of pulse profile time-delay estimation. Funding This work is supported by the National Natural Science Foundation of China under Grant nos. 61772187 and 61873196. References [1] S.I. Sheikh, The Use of Variable Celestial X-Ray Sources for Spacecraft Navigation Dissertation, University of Maryland, Maryland, 2005. [2] J. Liu, X.L. Ning, X. Ma, J.C. Fang, Geometry error analysis in solar doppler difference navigation for the capture phase, IEEE Trans. Aeros. Electron. Syst. 55 (5) (2019) 2556–2567. [3] X.L. Ning, M.Z. Gui, J.C. Fang, G. Liu, Differential X-ray pulsar aided celestial navigation for mars exploration, Aerosp. Sci. Technol. 62 (1) (2017) 36–45. [4] J. Liu, J.C. Fang, G. Liu, Fractional differentiation-based observability analysis method for nonlinear X-ray pulsar navigation system, Proc. Inst. Mech. Eng. Part G: J. Aerosp. Eng. 232 (8) (2018) 1467–1478. [5] J. Liu, J.C. Fang, Z.W. Kang, J. Wu, X.L. Ning, Novel algorithm for X-ray pulsar navigation against doppler effects, IEEE Trans. Aeros. Electron. Syst. 51 (1) (2015) 228–241. [6] Z.W. Kang, H.C. He, J. Liu, Adaptive pulsar time delay estimation using wavelet-based RLS, Optik 171 (2018) 266–276. [7] Y.D. Wang, Z. Wei, S.M. Sun, L. Li, X-ray pulsar-based navigation using time-differenced measurement, Aerosp. Sci. Technol. 36 (2014) 27–35. [8] A.A. Emadzadeh, J.L. Speyer, Asymptotically efficient estimation of pulse time delay for X-ray pulsar based relative navigation, Aiaa Guidance, Navigation, & Control Conference (2009). [9] J.X. Li, X.Z. Ke, Maximum-likelihood TOA estimation of X-ray pulsar signals on the basis of Poison model, Chin. Astron. Astrophys. 35 (1) (2011) 19–28. [10] A.A. Emadzadeh, J.L. Speyer, On modeling and pulse phase estimation of X-ray pulsarsr, IEEE Trans. Signal Process. 58 (9) (2010) 4484–4495. [11] A.A. Emadzadeh, A.R. Golshan, J.L. Speyer, Consistent estimation of pulse delay for X-ray pulsar based relative navigation, Joint 48th IEEE Conference on Decision and Control (2009) 1488–1493.
10
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Y. Wu, et al.
[12] M.F. Xue, X.P. Li, H.F. Sun, H.Y. Fang, A fast pulse phase estimation method for X-ray pulsar signals based on epoch folding, Chin. J. Aeronaut. 29 (3) (2016) 746–753. [13] P.F. Li, G.D. Xu, L.M. Dong, T.R. Hou, A real time estimation method of time-delay for X-ray pulsar signal, Acta Aeronaut. Astronaut. Sin. 35 (7) (2014) 1966–1976. [14] Y.D. Wang, Z. Wei, Pulse phase estimation of X-ray pulsar with the aid of vehicle orbital dynamics, J. Navig. 69 (2) (2016) 414–432. [15] Z. Su, L.P. Xu, T. Wang, H. Zhang, A new time delay measurement algorithm for pulsar accumulated pulse profile, J. Astronaut. 32 (6) (2011) 1256–1261 [Chinese]. [16] Q.Q. Lin, P. Shuai, L.W. Huang, X.Y. Zhang, X.M. Bei, A high-precision TOA estimation method for X-ray pulsar-based navigation system, J. Astronaut. 37 (6) (2016) 704–710 [Chinese]. [17] H.Y. Fang, B. Liu, X.P. Li, et al., Time delay estimation method of X-ray pulsar observed profile based on the optimal frequency band, Acta Phys. Sin. 65 (11) (2016) 119701. [18] J.H. Taylor, Pulsar timing and relativistic gravity, Philos. Trans. Phys. Sci. Eng. 341 (1660) (1992) 117–134. [19] Z.H. Xie, L.P. Xu, G.R. Ni, Time offset measurement algorithm based on bispectrum for pulsar integrated pulse profile, Acta Phys. Sin. 57 (10) (2008) 6683–6688. [20] Z.W. Kang, X. He, J. Liu, X-ray pulsars time delay estimation using GSO-based bispectral feature points, Opt. Int. J. Light Electron Opt. 127 (2016) 5050–5054. [21] Y.D. Wang, W. Zheng, D.P. Zhang, X-ray pulsar-based navigation method using one sensor and modified time-differenced measurement, Proc. Inst. Mech. Eng. Part G: J. Aerosp. Eng. 233 (1) (2017) 299–309. [22] H.F. Sun, X. Sun, H.Y. Fang, et al., Building X-ray pulsar timing model without the use of radio parameters, Acta Astronaut. 143 (2018) 155–162. [23] A.A. Emadzadeh, J.L. Speyer, X-ray pulsar-based relative navigation using epoch folding, IEEE Trans. Aerosp. Electron. Syst. 47 (4) (2011) 2317–2328. [24] Y. Wang, W. Zhang, Pulsar phase and Doppler frequency estimation for XNAV using on-orbit epoch folding, IEEE Trans. Aerosp. Electron. Syst. 52 (5) (2017) 2210–2219. [25] J. Liu, Z.H. Yang, Z.W. Kang, X. Chen, Fast CS-based pulsar period estimation method without tentative epoch folding and its CRLB, Acta Astronaut. 160 (2019) 90–100. [26] Z.H. Xie, L.P. Xu, G.R. Ni, Pulsar signal recognition based on one-dimension selected line spectra, J. Infrared Millim. Waves 26 (3) (2007) 187–195 [Chinese]. [27] X.D. Zhang, Modern Signal Processing, 2nd ed., Tsinghua University Press, Beijing, 2003, pp. 328–335 [Chinese]. [28] O. Rioul, P. Duhamel, Fast algorithms for discrete and continuous wavelet transforms, IEEE Trans. Inf. Theory 38 (2) (1992) 569–586. [29] NASA’s HEASARC Data Access. [Cited 28 December 2017]. Available from: https://heasarc.gsfc.nasa.gov/docs/archive.html. [30] EPN (European Pulsar Network). [Cited 28 December 2017]. Available from: http://www.jb.man.ac.uk/research/pulsar/Resources/epn/browser.html. [31] A.S. Yaseen, O.N. Pavlova, A.N. Pavlov, Noisy signal filtration using complex wavelet basis sets, Tech. Phys. Lett. 43 (7) (2017) 641–644. [32] X.P. Liu, W. Yuan, L.L. Han, et al., X-ray pulsar signal de-noising for impulse noise using wavelet packet, Aerosp. Sci. Technol. 64 (2017) 147–153. [33] L. Ljung, M. Morf, D. Falconer, Fast calculation of gain matrices for recursive estimation schemes, Int. J. Control 27 (1) (1978) 1–19.
11