An Algorithm for Mitigating Transient RFI in Pulsar Observation

Since pulsar periods are stable, an accurate ephemeris model can be obtained for a pulsar through long-term measurements. The TOAs predicted from the model are different from the observed TOAs resulting in residuals. In the case of ignoring system and model errors, the model is an approximation to the actual observation (Wu et al. 2018; Lyne & Graham-Smith 2012). Based on this, the ephemeris model and TOA formula can be used to accurately determine the pulsar phase.

From Figure 1, t₀ is the time of the reference phase in the ephemeris model; t is the start time of the current observation, and t_mod is the TOA of the current observation. The pulsar phase φ(t) is given by

$$\begin{eqnarray}\begin{array}{rcl}\varphi (t) & = & {\varphi }_{0}+60\cdot \hat{t}\cdot {f}_{0}+K(1)+\hat{t}\cdot K(2)\\ & & +{\hat{t}}^{2}\cdot K(3)+\bigtriangleup ,\end{array}\end{eqnarray} \tag{ 1 }$$

where φ₀ is the reference phase; $\hat{t}=t-{t}_{0};$ f₀ signifying the reference pulsar rotation frequency; △ is the infinitesimal of higher order, which is ignored in the following calculations. The polynomial coefficient K(n) can be obtained from the pulsar ephemeris using TEMPO2 software (Edwards et al. 2006; Hobbs et al. 2006). The φ(t) obtained from Equation (1) is the pulsar phase relative to t₀. Performing 360° modulus of φ(t), we obtain the $\hat{\varphi }(t)$ for the current observation:

$$\begin{eqnarray}&&\hat{\varphi }(t)=\varphi (t) \% 360,\end{eqnarray} \tag{ 2 }$$

where % is the modulo operation. Similarly, t_mod is obtained from the unit conversion of φ(t).

Zoom In Zoom Out Reset image size

After obtaining the specific position for the pulsar reference phase, the dispersion delay is corrected for the TOAs at each sub-band using

$$\begin{eqnarray}&&\bigtriangleup {t}_{i}^{\mathrm{DM}}=D\times \left(\displaystyle \frac{1}{{v}_{i}^{2}}-\displaystyle \frac{1}{{v}_{0}^{2}}\right)\times \mathrm{DM}.\end{eqnarray} \tag{ 3 }$$

Here, DM is the dispersion measure; $\bigtriangleup {t}_{i}^{\mathrm{DM}}$ is the dispersion delay time between v_i and v₀; D is the dispersion constant; v_i is the ith sub-band observation frequency; and v₀ is its corresponding reference frequency (Lorimer & Kramer 2012). Hence, the specific position of the phase at different frequencies can be obtained using Equation (3).

Since the pulse phase varies across the frequency channels, after obtaining the precise phase at each frequency channel, the pulse width needs to be solved separately to ensure accurate extraction of the pulse data. In our calculation, we select the standard profile phase at the center of the observing station. Then, a baseline data is extracted across a width that is two times greater than the width measured at 10% of the peak intensity (W₁₀) of the integrated profile (Lyne & Graham-Smith 2012). For a stable observing system, the baseline data can be considered as linearly correlated with the system temperature. Hence, linear fitting is used to subtract the background data. Next, the positions on the left and on the right outside of the baseline data where negative values first appear are determined and marked as the left starting position and right ending position of the pulse phase, respectively. The pulse profile is then extracted based on these two positions.

After extracting the pulsar profile, the off-pulse data are white background noise and possible RFI. We use the method of threshold setting and loop iterations to identify and label the RFI in the data.

First, we calculate the mean value, μ, and the standard deviation, σ, of the off-pulse data. The threshold range is then set as (μ − k σ, μ + k σ), where $k\in {{\mathbb{R}}}^{+}$ is the coefficient of the threshold range.

In practice, we find that the standard 3σ rule does not work well in our case. Accordingly, we reduce the coefficient k for the threshold to adapt to the RFI situation. By analyzing the observation data, we arrive at the reasonable range for k ∈ [1.8, 3] (Lyon et al. 2016).

In addition, we label the data above the threshold as RFI. In order to determine whether all RFIs have been recognized, PPSD recalculates μ and σ, with the new values represented by $\mu ^{\prime}$ and $\sigma ^{\prime}$, of the unlabeled data and obtains a new threshold. If all unlabeled data are below the new threshold, the RFI labeling is complete. Otherwise, the new threshold is used to identify and label the remaining RFI until all unlabeled data are within the latest threshold.

Next, the latest $\mu ^{\prime}$ and $\sigma ^{\prime}$ are used to generate random data points that obey Gaussian white noise (WGN) $N({\mu }_{i}^{\prime} ,{\sigma }_{i}^{\prime} )$ to substitute for the RFI. After that, all off-pulse data obeys the WGN distribution. The on-pulse data is then put back to the original position in the data, thereby completing the process for mitigating the transient RFI.

The above processing is repeated for other frequency channels until all frequencies have been completed for mitigation. Finally, the processing results after RFI mitigation are output and saved. The specific algorithm workflow diagram is shown in Figure 2.

Zoom In Zoom Out Reset image size

With the rapid development and application of electronic and communication technology, radio astronomical observatories generally face the challenge of RFI. The NanShan 25 m Radio Telescope (NSRT) is located far away from the city at the north of the Tianshan Mountains in China. In the initial stage of the construction, the radio astronomy observation environment at the site was away from man-made RFI (Liu et al. 2019). Since then, the RFI at the NSRT site has gradually worsened in partial time and in specific observing directions, and the existing mitigation method for transient RFI has become unsatisfactory for scientific purposes.

For testing of PPSD, we selected some pulsar timing observations from NSRT at L band that contained obvious transient RFI. The receiver backend for the data collection was Digital Filter Bank 3 (DFB3), which was developed by CSIRO, Australia. The kernel data parameters included the frequency range, which covered the range of 1300–1812 MHz, and the sub-integration time of 30 s. The operating system for the experimental environment was Ubuntu 20.4, and the programming language was python3.7.

We evaluate the effectiveness of the algorithm based on the signal-to-noise ratio (S/N) of the pulsar signals. The S/N is expressed as

$$\begin{eqnarray}&&{\rm{S}}/{\rm{N}}=\displaystyle \frac{{S}_{{\rm{p}}}}{{N}_{{\rm{rms}}}},\end{eqnarray} \tag{ 4 }$$

where S_p is the peak value of the original signal, N_RMS refers to the rms of the noise excluding the signal: $\sqrt{\tfrac{{\sum }_{i=1}^{n}{x}_{i}^{2}}{n}}$ (Gallager 2008; Lorimer & Kramer 2012).

We obtained a set of folded data of PSR J0332+5434 as a test example to identify and mitigate the transient RFI. We compare the results before and after the data processing, which are shown in Figure 3. The plots in the upper row are obtained from the original data, whereas the plots in the lower row are results after mitigating the transient RFI. Several obvious transient RFIs in the original data, including their shapes and intensities, are indicated by the red rectangles in the upper row. The pulsar profile is shown in the green rectangle in each plot, and the real pulsar profile after de-dispersion is displayed in the two plots in the right column.

Zoom In Zoom Out Reset image size

The original data, without mitigation, is shown in the middle plot in the upper row, which contains many strong RFIs, including transient RFI, resulting in a low S/N. However, the de-dispersion smooths the transient RFI to a certain extent as shown in the right plot in the upper row. Nevertheless, the average pulse profile has improved and clearly visible, and the S/N has also increased. However, there are still obvious transient RFIs remaining in the data. The PPSD uses TOA model to calculate the pulse phase and the width of the pulse window in every frequency subchannels, which are used to mitigate RFI by threshold (see Section 2 for the specific details on the process). From the lower row, the transient RFIs are obviously mitigated, when compare with that in the upper row, and without changing the pulsar profile. As shown in the two plots in the middle column, the S/N has increased from 2.7–26.1, which is more than eight times. From the two plots in the right column, the S/N has increased by nearly nine times from 6.9–67.2.

In the above example, the specific coefficient for the threshold is 2. The performance of RFI mitigation is related to the threshold coefficient. To further mitigate the spikes in the off-pulse baseline in the lower row in Figure 3, a smaller threshold coefficient is needed. When RFI coincides with the pulse window in the sub-band, PPSD has certain limitations. The probability of RFI being present in the pulse window for a sub-band may increase under the following conditions: (i) long transient RFI duration, (ii) large pulse window, and (iii) the source is a millisecond pulsar, etc. However, these conditions are not common.

In general, the off-pulse transient RFI is effectively mitigated after using PPSD. The threshold coefficient can change dynamically according to the number, distribution, and strength of the transient RFI. We find that the best range for the threshold coefficient is between [1.8, 3] (Lyon et al. 2016).

PSRCHIVE provides the PAZ and PAZI sub-routines to mitigate RFI in observation data. PAZ automatically zaps the RFI in frequency based on median smoothed and profile lawn mowing methods. However, it is limited to transient RFI. On the other hand, PAZI presents non-de-dispersion plot, which is helpful for manually mitigating the RFI. The user manually zaps the RFI data, and PAZI will re-assign the data. For comparison, we use PAZI to mitigate transient RFI in the same data set of PSR J0332+5434 as described in Section 3.2, and the results are shown in Figure 4. As indicated by the red dashed rectangles in Figure 4, the values used to replace the mitigated transient RFI are less than the background noise, and the transient RFI appears over zapped. Although this may not always cause obvious impact on the final result, the efficiency of manual zapping is limited. More importantly, PAZI relies heavily on user's experience, meaning that different people performing the RFI zapping may lead to different results. In addition, PAZI mitigates all data indiscriminately, even if the data is inside the pulse window.

Zoom In Zoom Out Reset image size

By comparison, PPSD algorithm remedies the shortcomings of PSRCHIVE for mitigation of transient RFI, and provides the ability to automatic identification and mitigation of the transient RFI in pulsar observation.

We applied this method to 10 other pulsars and achieved an S/N improvement ranging from 40%–800%. The stronger the RFI is in an observation, the more significant the S/N improvement is after PPSD is used. To compare the PPSD results, the average profile from three typical observations are shown in Figure 5. In the upper panel, the data is not de-dispersed, which shows the real shapes of the transient RFI. On the contrary, the shapes of the transient RFI are affected by the de-dispersion as shown in the middle panel. In the lower panel, the plots show the pulse profile after transient RFI mitigation. The overall results show that PPSD can effectively mitigate transient RFI and enhance the S/N of pulsar observations.

Zoom In Zoom Out Reset image size

In order to test the feature of PPSD algorithm on recovering the original pulsar signal, we also verified the validity of the data before and after the transient RFI mitigation on each frequency channel. We found that the original data of the pulse profile was recovered effectively.