# A Fast Detection Algorithm for the X-Ray Pulsar Signal

Authors: Hao Liang; Yafeng Zhan
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

The detection of the X-ray pulsar signal is important for the autonomous navigation system using X-ray pulsars. In the condition of short observation time and limited number of photons for detection, the noise does not obey the Gaussian distribution. This fact has been little considered extant. In this paper, the model of the X-ray pulsar signal is rebuilt as the nonhomogeneous Poisson distribution and, in the condition of a fixed false alarm rate, a fast detection algorithm based on maximizing the detection probability is proposed. Simulation results show the effectiveness of the proposed detection algorithm.


---

## Body

## 1. Introduction

X-ray pulsars are rapidly rotating neutron stars which could emit X-ray signals periodically, stably, and uniquely [1]. Since Downs first proposed interplanetary navigation using pulsars in 1974 [2], X-ray Pulsar Navigation (XPNAV) system has been considered in lots of scenarios for autonomous navigation [3–7]. XPNAV gets the spacecraft’s position vector from measuring Time Difference of Arrival (TDOA) of the X-ray pulsar signals between the spacecraft and the Solar System Barycenter (SSB). Because of the large distance (thousands of light years) between the X-ray pulsars and the solar system, the received X-ray pulsar signal is always weak and easy to be disturbed by kinds of noise [8]. It is difficult to recover the X-ray pulsar profile [9] and get high precision TDOA, and each observation time for the X-ray pulsars may be as long as hundreds of seconds [3, 10]. However, before using the X-ray pulsar signals for navigation, we should firstly determine whether we have received the desired X-ray pulsar signals, which would help the detector’s altitude controller to confirm the right direction of the X-ray pulsar. Unlike measuring TDOA, the signal detection process should be both accurate and fast. Therefore, it is a challenge to detect the weak X-ray pulsar signal with short time and limited photons.The first widely used X-ray pulsar signal detection algorithm is based on the Fast Fourier Transform (FFT), which has been proved to be invalid when the signal’s profile is nonsinusoidal or the noise is nongaussian [11]. The Gregory-Loredo method parameterizes the number of bins for epoch folding and calculates the average ratio of the probability for the periodic model to the constant model [12]. The Gregory-Loredo method is effective on discovering new pulsars when there is no prior knowledge about the signal’s shape; however, it does not use the prior information and is less efficient. The other detection algorithms, like S-transform [13, 14] and Match Filter (MF) [15], consider the noise obeying the Gaussian distribution after epoch folding. However, that happens when the number of photons in the bins of epoch folding is larger than 20 [8]. For the detection problem, the observation time is not long enough to support approximating the Poisson noise by the Gaussian [16], and the algorithms above have this deficiency. Zhang et al. built the Poisson distribution model of X-ray pulsar and proposed a detection method by using the Bayesian theory [16]. However, it simply built the signal and noise by homogeneous Poisson distribution, which is imprecision. Besides that, it did not consider the influences of the unknown signal phase and the Doppler shift, which make it hard to be used in the practical scenario.For XPNAV, the accurate prior knowledge about the X-ray pulsar signal is essential. The prior knowledge includes the X-ray background noise’s rate, the signal’s flux, period, and profile. By using the prior information, we may build the Probability Density Function (PDF) [10] and get a fast detection algorithm based on Neyman-Pearson (NP) criterion [17]. Even in the condition of few photons for detection, the proposed algorithm’s performance could meet detection requirements of XPNAV.The rest of this paper is organized as follows. Section2 introduces the X-ray pulsar signal model. Section 3 describes the fast detection algorithm in detail. Computer simulations are shown in Section 4 and some conclusions are given in the final part.


## 2. X-Ray Pulsar Signal Model

Unlike the radio signal, the X-ray pulsar signal shows the particle property, which means the output of X-ray detector on the spacecraft is not an analog signal but a discrete arriving time sequence of individual photons. Assuming the observation time isTobs, the number of received photons is NTobs, and the PDF of NTobs can be described by nonhomogeneous Poisson process [10]:(1)PNTobs=k=∫0Tobsλtdtkexp⁡-∫0Tobsλtdtk!,where λt is the arriving rate of the photons as(2)λt=λb+λshφt.Here λb is the X-ray background noise’s rate which is time-invariant in the cosmos environment, λs is the X-ray pulsar signal’s flux, hφt is the X-ray pulsar signal’s normalization standard profile template, φt=fr×t+ϕ0 is the X-ray pulsar signal’s phase, fr is the received frequency, and ϕ0 is the initial phase. Since the X-ray pulsar signal is periodic, hφt+T=hφt, where T is the period.


## 3. The Fast Detection Algorithm

### 3.1. The Optimal NP Algorithm

We divide the observation timeline intoN equal bins, which is shown in Figure 1. The duration time of each of the bins is Δt. When Δt is small enough, the arriving rate in the nth bin λnθ could be treated as a determined variable:(3)λnθ≈1Δt∫t0+n-1Δtt0+nΔtλt;θdt,where θ=frϕ0T; the supper index T represents the transpose of a vector. The number of photons in nth bin xn is a random variable with the PDF:(4)pxn=k;θ=λnθΔtkk!exp⁡-λnθΔt.Figure 1
The observation timeline.Since the random variables in the different bins are independent, the Joint Probability Density Function (JPDF) could be expressed as(5)pX;θ=∏n=1Npxn;θ=∏n=1NλnθΔtxnxn!exp⁡-λnθΔt,where X=x1,x2,…,xN.If the X-ray pulsar signal exists in the received photons,λnθ=λb+λshθ, then the JPDF of X could be expressed as(6)H1:p1X;θ=∏n=1Nλb+λshθΔtxnxn!exp⁡-λb+λshθΔt.If the X-ray pulsar signal does not exist in the received photons,λnθ=λb, then the JPDF of X could be expressed as(7)H0:p0X;θ=∏n=1NλbΔtxnxn!exp⁡-λbΔt.In the condition of a fixed probability of false alarmpFA, NP criterion maximizes the detection probability [17]. Based on NP criterion, we could do the Likelihood Ratio Test (LRT):(8)LX;θ=p1X;θp0X;θ=γ.Ifγ is larger than the threshold γth, we can determine that we have received the X-ray pulsar signal. Otherwise, the X-ray pulsar signal does not exist in the received photons. Equation (8) could be simplified as follows:(9)LX;θ=p1X;θp0X;θ=∏n=1Nλnθλbxnexp⁡λb-λnθΔt,(10)ln⁡LX;θ=∑n=1Nln⁡λnθλbxn+∑n=1Nλb-λnθΔt,(11)ln⁡LX;θ+∑n=1NλsθΔt=∑n=1Nxnln⁡λnθ-ln⁡λb.Since∑n=1N(λs(θ))Δt is independent of the measurement variable xn, the left side of (11) could be added to the threshold, and then(12)∑n=1Nxnln⁡λnθ-ln⁡λb=γ′.IfΔt→0, there is only one or none photon in each bin, and then (12) could be converted as(13)Gφ=∑i=1Mln⁡λφi-ln⁡λb=γ′,where M is the total number of the photons, φ=φ1φ2⋯φMT, φi is the phase of ith photons, which could be calculated by(14)φi=fr·ti+ϕ0.From (13), the optimal detection variable is Gφ. Thus, the steps of the optimal NP detection algorithm for the X-ray pulsar signal are as below:(I)
The X-ray detector on-board receives the photons and notes each photons’ arriving timeti.(II)
Convert the transferred arriving timeti into φi using (14).(III)
Search the arriving rateλφi of each arriving time from the XPNAV database (since λφi is periodic, the XPNAV database is a table which stores one period of the function λφi), and calculate the logarithm lnλφi.(IV)
Add up allln⁡λφi, and compare the summation Gφ with γth.The prior knowledge before the detection process includes the XPNAV database about the arriving rate functionsλt and γth. The XPNAV database could be built by observing the X-ray pulsars for months and some explorers have tried to do that, such as the Rossi X-Ray Timing Explorer (RXTE) from the National Aeronautics and Space Administration (NASA) and the Röntgen Satellite (ROSAT) from the German Aerospace Center (DLR). It is difficult to get an analytical expression for γth. Monte Carlo simulations would be used to calculate γth. The calculation process is as below:(I)
Set a fixed probability of false alarmpFA.(II)
Generate the arriving time of the photons based on the Poisson distribution with the arriving rateλb.(III)
CalculateGφ based on the above steps of the optimal NP detection algorithm.(IV)
Repeat step (II) and step (III) forK times; K is the times of the Monte Carlo simulations.(V)
Rank each simulation result in descending order, and set a thresholdγth to make the proportion of Gφ greater than γth equal pFA. γth is the threshold we pursue.


### 3.2. The Suboptimal NP Algorithm

Besides the prior knowledge, how to convert the arriving timeti into the phase φi is another important issue. ϕ0 is different due to the spacecraft’s position, and fr could be affected by the Doppler shift due to the spacecraft’s velocity:(15)φi=fr×ti+ϕ0=f0+fd×ti+ϕ0,(16)fd=v·rcf0,where f0 is the fixed frequency of the X-ray pulsar, fd is the Doppler frequency, v is the spacecraft’s velocity vector, r is the X-ray pulsar’s direction vector, and c is the light velocity. Supposing there are no prior distributions of ϕ0 and fd, we could use the two-dimension search based on the Maximum Likelihood (ML) estimation [10] to deal with the issue.Assumingϕ0’s searching range is 0,1, dividing the range into p segments, the jth searching phase is ϕ0j. Dividing the fd’s searching range into q segments, the kth searching frequency is fdk. Then, take ϕ0j and fdk into step (II) of the optimal NP detection algorithm and calculate Gj,kφ;ϕ0j,fdk. ϕ0 and fd in maxGj,kφ;ϕ0j,fdk have the maximum probability equaling the real ones, where max· is the maximum function. Therefore, we could use maxGj,kφ;ϕ0j,fdk as the detection variable to be compared with γth. We may call the algorithm resolving the unknown ϕ0 and fd as the suboptimal NP detection algorithm.For the practical navigation application,ϕ0 and fd are always unknown, so the suboptimal NP detection algorithm is more useful than the optimal one.


## 3.1. The Optimal NP Algorithm

We divide the observation timeline intoN equal bins, which is shown in Figure 1. The duration time of each of the bins is Δt. When Δt is small enough, the arriving rate in the nth bin λnθ could be treated as a determined variable:(3)λnθ≈1Δt∫t0+n-1Δtt0+nΔtλt;θdt,where θ=frϕ0T; the supper index T represents the transpose of a vector. The number of photons in nth bin xn is a random variable with the PDF:(4)pxn=k;θ=λnθΔtkk!exp⁡-λnθΔt.Figure 1
The observation timeline.Since the random variables in the different bins are independent, the Joint Probability Density Function (JPDF) could be expressed as(5)pX;θ=∏n=1Npxn;θ=∏n=1NλnθΔtxnxn!exp⁡-λnθΔt,where X=x1,x2,…,xN.If the X-ray pulsar signal exists in the received photons,λnθ=λb+λshθ, then the JPDF of X could be expressed as(6)H1:p1X;θ=∏n=1Nλb+λshθΔtxnxn!exp⁡-λb+λshθΔt.If the X-ray pulsar signal does not exist in the received photons,λnθ=λb, then the JPDF of X could be expressed as(7)H0:p0X;θ=∏n=1NλbΔtxnxn!exp⁡-λbΔt.In the condition of a fixed probability of false alarmpFA, NP criterion maximizes the detection probability [17]. Based on NP criterion, we could do the Likelihood Ratio Test (LRT):(8)LX;θ=p1X;θp0X;θ=γ.Ifγ is larger than the threshold γth, we can determine that we have received the X-ray pulsar signal. Otherwise, the X-ray pulsar signal does not exist in the received photons. Equation (8) could be simplified as follows:(9)LX;θ=p1X;θp0X;θ=∏n=1Nλnθλbxnexp⁡λb-λnθΔt,(10)ln⁡LX;θ=∑n=1Nln⁡λnθλbxn+∑n=1Nλb-λnθΔt,(11)ln⁡LX;θ+∑n=1NλsθΔt=∑n=1Nxnln⁡λnθ-ln⁡λb.Since∑n=1N(λs(θ))Δt is independent of the measurement variable xn, the left side of (11) could be added to the threshold, and then(12)∑n=1Nxnln⁡λnθ-ln⁡λb=γ′.IfΔt→0, there is only one or none photon in each bin, and then (12) could be converted as(13)Gφ=∑i=1Mln⁡λφi-ln⁡λb=γ′,where M is the total number of the photons, φ=φ1φ2⋯φMT, φi is the phase of ith photons, which could be calculated by(14)φi=fr·ti+ϕ0.From (13), the optimal detection variable is Gφ. Thus, the steps of the optimal NP detection algorithm for the X-ray pulsar signal are as below:(I)
The X-ray detector on-board receives the photons and notes each photons’ arriving timeti.(II)
Convert the transferred arriving timeti into φi using (14).(III)
Search the arriving rateλφi of each arriving time from the XPNAV database (since λφi is periodic, the XPNAV database is a table which stores one period of the function λφi), and calculate the logarithm lnλφi.(IV)
Add up allln⁡λφi, and compare the summation Gφ with γth.The prior knowledge before the detection process includes the XPNAV database about the arriving rate functionsλt and γth. The XPNAV database could be built by observing the X-ray pulsars for months and some explorers have tried to do that, such as the Rossi X-Ray Timing Explorer (RXTE) from the National Aeronautics and Space Administration (NASA) and the Röntgen Satellite (ROSAT) from the German Aerospace Center (DLR). It is difficult to get an analytical expression for γth. Monte Carlo simulations would be used to calculate γth. The calculation process is as below:(I)
Set a fixed probability of false alarmpFA.(II)
Generate the arriving time of the photons based on the Poisson distribution with the arriving rateλb.(III)
CalculateGφ based on the above steps of the optimal NP detection algorithm.(IV)
Repeat step (II) and step (III) forK times; K is the times of the Monte Carlo simulations.(V)
Rank each simulation result in descending order, and set a thresholdγth to make the proportion of Gφ greater than γth equal pFA. γth is the threshold we pursue.


## 3.2. The Suboptimal NP Algorithm

Besides the prior knowledge, how to convert the arriving timeti into the phase φi is another important issue. ϕ0 is different due to the spacecraft’s position, and fr could be affected by the Doppler shift due to the spacecraft’s velocity:(15)φi=fr×ti+ϕ0=f0+fd×ti+ϕ0,(16)fd=v·rcf0,where f0 is the fixed frequency of the X-ray pulsar, fd is the Doppler frequency, v is the spacecraft’s velocity vector, r is the X-ray pulsar’s direction vector, and c is the light velocity. Supposing there are no prior distributions of ϕ0 and fd, we could use the two-dimension search based on the Maximum Likelihood (ML) estimation [10] to deal with the issue.Assumingϕ0’s searching range is 0,1, dividing the range into p segments, the jth searching phase is ϕ0j. Dividing the fd’s searching range into q segments, the kth searching frequency is fdk. Then, take ϕ0j and fdk into step (II) of the optimal NP detection algorithm and calculate Gj,kφ;ϕ0j,fdk. ϕ0 and fd in maxGj,kφ;ϕ0j,fdk have the maximum probability equaling the real ones, where max· is the maximum function. Therefore, we could use maxGj,kφ;ϕ0j,fdk as the detection variable to be compared with γth. We may call the algorithm resolving the unknown ϕ0 and fd as the suboptimal NP detection algorithm.For the practical navigation application,ϕ0 and fd are always unknown, so the suboptimal NP detection algorithm is more useful than the optimal one.


## 4. Computer Simulations

Three weak X-ray pulsars (B0540-69, B1744-24A, and B1823-13) are selected to simulate the proposed detection algorithm’s performance. The parameters of the X-ray pulsars are shown in Table1, and the normalization standard profile templates [18] are shown in Figure 2. The X-ray background noise’s rate λb is 0.005 photons/s/cm2. The X-ray pulsar signal’s flux λs could be obtained from Table 1 [1]. The X-ray detector’s area is 1 m2, and the fixed probability of false alarm pFA is set as 0.001.Table 1
The parameters of X-ray pulsar sources.
ID Name RA(J2000) (hh:mm:ss) DE(J2000) (dd:mm:ss) Flux (photons/cm2/s) Period (s) 1 B0540-69 05:39:39 −69:44:36 5.15 e × 10 - 3 0.0505 2 B1744-24A 17:48:02 −24:46:38 1.09 × 10 - 3 0.0116 3 B1823-13 18:26:13 −13:34:47 2.63 × 10 - 3 0.1015Figure 2
The normalization standard profile templates of three X-ray pulsars.The initial phaseϕ0 obeys the uniform distribution with range 0,1. The unknown velocity v obeys the uniform distribution with range -7.9,7.9 km/s, in which 7.9 km/s is the first cosmic velocity. The Doppler frequency fd could be calculated by (16). The method to generate the photons refers to reference [10].The MF algorithm calculates the correlation peak for the X-ray pulsar signal’s profile with the standard signal template [15]. Since the MF algorithm is optimal when the noise is obeying the Gaussian distribution, we choose the MF algorithm for comparison.We divide the simulation algorithm into three groups: (I) the optimal NP detection algorithm; (II) the suboptimal NP detection algorithm, and the number of segments for searching bothϕ0 and fd is set as 10; (III) the MF algorithm with unknown information about ϕ0 and fd, and the number of bins for epoch folding is 64. We perform the Monte Carlo simulations for 1000 times to calculate the detection probability, and the results are shown from Figures 3–5.Figure 3
The detection probability for B0540-69.Figure 4
The detection probability for B1744-24A.Figure 5
The detection probability for B1823-13.From Figures3–5, the detection probability increases with the increasing observation time. Meanwhile, in the same simulation conditions, the suboptimal NP detection algorithm’s detection probability is better than the MF’s. The MF algorithm is effective when the noise is obeying the Gaussian distribution. However, in the condition of short observation time and limited photons, the noise obeys Poisson distribution, which makes the MF algorithm less effective than the proposed detection algorithm.Since the limited number of photons, the estimation forϕ0 and fd based on ML may be imprecise, and the detection probability of the suboptimal NP detection algorithm is lower than the optimal NP detection algorithm. Therefore, the optimal NP detection algorithm is a performance bound.


## 5. Conclusion

The detection for the X-ray pulsar signal is important for XPNAV and some other navigation scenarios using X-ray pulsars. In this paper, we propose a fast detection algorithm for the X-ray pulsar signal based on NP criteria. Computer simulations show that, in the condition of short observation time and limited number of photons, the detection probability of the proposed algorithm is still high. Besides that, the fast detection algorithm may be extended to other detection issues when the signals obey the Poisson distribution.



---

*Source: 9607821-2017-01-04.xml*