Compressed Sensing Based RFI Mitigation and Restoration for Pulsar Signals

Suppose a digital filter bank has P equal-width frequency channels. Upon saving the data to disk, the observation program divides the time samples in one pulse period into T phase bins. Let an observed 2D pulsar profile be a digital matrix F ₀ within an observing bandwidth $[{\nu }_{1}^{1},{\nu }_{P}^{2}]$. Each frequency channel is tuned to receive a signal ${{\boldsymbol{f}}}_{p}^{\mathrm{ch}},1\leqslant p\leqslant P$ within a narrow bandwidth $[{\nu }_{p}^{1},{\nu }_{p}^{2}]$ with ν_p as the central frequency, where the vector ${{\boldsymbol{f}}}_{p}^{\mathrm{ch}}={({f}_{{pt}}^{\mathrm{ch}})}_{1\leqslant t\leqslant T}$. In this way, the broadband F ₀ is split into adjacent narrowband signals. The dispersion frequency relationship (Lorimer & Kramer 2005) is calculated as

$$\begin{eqnarray}&&{\rm{\Delta }}{t}_{p}\simeq \mathrm{DC}\times \left({\nu }_{\mathrm{ref}}^{-2}-{\nu }_{p}^{-2}\right)\times \mathrm{DM},\end{eqnarray} \tag{ 1 }$$

where DC = 4.15 × 10³ (MHz² pc⁻¹ cm³ s) is the dispersion constant, DM (parsecs per cubic centimeter) is the dispersion measure, and ${\nu }_{\mathrm{ref}}=({\nu }_{1}^{1}+{\nu }_{P}^{2})/2$ is the reference frequency and often chosen to be the central frequency of the broadband. Frequencies are all in megahertz. Let d be the distance to the pulsar, and n_e be the electron number density. DM is defined as $\mathrm{DM}={\int }_{0}^{d}{n}_{{\rm{e}}}{\rm{d}}l$, where l is the line of sight. For the profiles collected by the Nanshan 26 m Radio Telescope, we set P = 1024 and T = 512 for calculation.

A necessary three-step preprocessing is carried out. Step 1: baseline removal by a channel mean subtraction. Step 2: zero-setting for saturated channels, resulting in an F ₂. Step 3: a simple threshold-to-zero to remove the huge-amplitude RFI, resulting in an F . Let γ be the threshold of Step 3. The threshold-to-zero means that for an ${{\boldsymbol{f}}}_{p}^{\mathrm{ch}},1\leqslant p\leqslant P$, if $\max ({{\boldsymbol{f}}}_{p}^{\mathrm{ch}})-\min ({{\boldsymbol{f}}}_{p}^{\mathrm{ch}})\gt \gamma$, then ${{\boldsymbol{f}}}_{p}^{\mathrm{ch}}={\bf{0}}$. Suppose the median absolute deviation (MAD) of a matrix X is MAD( X ) = median∣ X − median( X )∣, then the standard deviation $\sigma =\tfrac{\mathrm{MAD}}{0.6745}$. For the case of Gaussian white noise, γ can be selected as 3σ (Starck et al. 2010). In Step 3, a preliminary removal for huge-amplitude RFI is carried out by setting $\gamma =6\sim 8\tfrac{\mathrm{MAD}\left({{\boldsymbol{F}}}_{2}\right)}{0.6745}$ or selecting a reasonable numerical number. Figure 1 takes PSR J1645-0317 (B1642-03), for example, and shows the process.

Zoom In Zoom Out Reset image size

Consider the standard CS problem with missing channels

$$\begin{eqnarray}&&{\boldsymbol{f}}={\boldsymbol{\Phi }}{\boldsymbol{u}}+{\boldsymbol{\epsilon }},\end{eqnarray} \tag{ 2 }$$

where ${\boldsymbol{\Phi }}\in {{\mathbb{R}}}^{L\times N}(L\leqslant N)$ is the CS measurement matrix and ${\boldsymbol{\epsilon }}\in {{\mathbb{R}}}^{L}$ is the measurement error. Signal ${\boldsymbol{u}}\in {{\mathbb{R}}}^{N}$ will be recovered via the noisy linear measurements ${\boldsymbol{f}}\in {{\mathbb{R}}}^{L}$, and f is the vector form of F and formed by concatenating every column f _t , 1 ≤ t ≤ T of F . Let = σ ω , where ω is the vector of random Gaussian white noise and RFI, and σ controls the magnitude of . Unfortunately, there is no purely statistical way to decide the statistics of ω and the detailed information about σ. For the purpose of reconstruction, is a combined factor encompassing both Gaussian noise and RFI that can be executed with or without statistical characteristics. Φ has many more rows than columns. Suppose A is the number of the zero set or missing channels, then L = (P − A)T and N = PT. Given Φ, the recovery of u from f is an underdetermined linear system leading to an ill-posed problem. To regularize it, assume that u has the property of "sparsity," i.e., u is s-sparse in ${\boldsymbol{\Theta }}={{\boldsymbol{\Psi }}}^{\dagger }\in {{\mathbb{R}}}^{M\times N}$ (M = N for a basis or M > N for a frame), such that u can be expressed by s nonzero coefficients in Θ, i.e., ∥Θ u ∥₀ = ∥ α ∥₀ = s, where ${\boldsymbol{\alpha }}\in {{\mathbb{R}}}^{M}$, M ≥ N, and † denotes the inverse transform. Then, u can be reconstructed with high accuracy from the incomplete measurements f . For most cases, signals of interest are not strictly sparse, but the assumption can be relaxed to "compressibility" if most of the energy in Θ u = α is contained in its largest s coefficients, i.e., amplitudes of the sorted scalar coefficients exhibit a polynomial decay: $| {\alpha }_{i}| \leqslant c\cdot {i}^{-\tfrac{1}{b}}$, 1 ≤ i ≤ s, b ∈ [0, 1], where c is a constant. The smaller the b is, the faster the amplitudes decay, and the more compressible a signal is. Supposing Φ and Θ are two orthogonal basis or frames, and Φ is noise-like incoherent in the Θ domain, i.e., Φ Θ satisfies a sufficient condition named restricted isometry property (Candés & Tao 2005; Candés & Wakin 2008), CS tells us that the sparse coefficients α can be accurately recovered. Thus,

$$\begin{eqnarray}&&{\boldsymbol{f}}={\boldsymbol{\Phi }}{\boldsymbol{\Psi }}{\boldsymbol{\alpha }}+{\boldsymbol{\epsilon }}.\end{eqnarray} \tag{ 3 }$$

Physically realizable random matrices that are incoherent to most transforms can be used as measurement matrices. One of the classical methods to solve Equation (3) is the generalized formulation of the unconstrained regularized LS based reconstruction, which reads

$$\begin{eqnarray}&&{\boldsymbol{u}}=\mathop{\mathrm{argmin}}\limits_{{\boldsymbol{u}}\in {{\mathbb{R}}}^{N}}\ \left\{\displaystyle \frac{1}{2}\parallel {\boldsymbol{f}}-{\boldsymbol{\Phi }}{\boldsymbol{\Psi }}{\boldsymbol{\alpha }}{\parallel }_{2}^{2}+\lambda J({\boldsymbol{\alpha }})\right\},\end{eqnarray} \tag{ 4 }$$

where λ is a nonnegative regularization parameter balancing the data fidelity against the sparsity. $J:{{\mathbb{R}}}^{N}\to {\mathbb{R}}$ is a fixed nonnegative valued cost function, and it is also a penalization term of α to preserve the sparsity property. J could be various choices. For example, to solve the NP-hard problem of the ℓ₀ pseudo-norm, one can use its optimal convex approximation, i.e., a sparsity-promoting ℓ₁ minimization instead to assist the reconstruction. For another example, the total variation (TV) regularizer assumes that the pulse portions are sparse, and the usual discretizations of TV are assumed to be convex. The reason to apply the regularized LS framework of Equation (4) is that the statistics of α and w, and the detailed information about σ are not required (Vehkaperä et al. 2016). Thus, the measurement error is considered as a combination of the RFI and random Gaussian white noise. A suboptimal synthesis problem termed least absolute shrinkage and selection operator (Candés & Tao 2006; Tibshirani 1996; Chen et al. 1998) is considered:

$$\begin{eqnarray}&&{\boldsymbol{\alpha }}=\mathop{\mathrm{argmin}}\limits_{{\boldsymbol{\alpha }}\in {{\mathbb{R}}}^{M}}\left\{\displaystyle \frac{1}{2}\parallel {\boldsymbol{f}}-{\boldsymbol{\Phi }}{\boldsymbol{\Psi }}{\boldsymbol{\alpha }}{\parallel }_{2}^{2}+\lambda \parallel {\boldsymbol{\alpha }}{\parallel }_{1}\right\},\end{eqnarray} \tag{ 5 }$$

where $\parallel {\boldsymbol{\alpha }}{\parallel }_{1}={\sum }_{i=1}^{N}| {\alpha }_{i}|$. The fast iterative soft-thresholding algorithm (FISTA; Beck & Teboulle 2009) is derived from the proximal forward–backward iterative scheme (Bruck 1977; Passty 1979) in the general framework of splitting methods (Facchinei & Pang 2003). For solving Equation (5), the FISTA will be applied. Set $\lambda =\tfrac{a\cdot \mathrm{MAD}\left({\boldsymbol{F}}\right)}{0.6745}$, where a is a parameter that can be adjusted.

The deleterious smoothing effect from RFI mitigation can cause precise astrophysical measurements (e.g., timing) to be negatively affected. It is necessary to perform signal retrieval after RFI mitigation to account for these effects. The signal retrieval will be carried out by a further MCA (Starck et al. 2005) decomposition of the measurement error . The MCA seeks the sparsest of all representations over the sparse dictionaries representing the morphological components. Suppose ${\boldsymbol{\epsilon }}\in {{\mathbb{R}}}^{L}$ contains three components of the removed signal details ${{\boldsymbol{u}}}^{{\rm{M}}}\in {{\mathbb{R}}}^{L}$, the RFI ${\boldsymbol{i}}\in {{\mathbb{R}}}^{L}$ and the Gaussian white noise ${\boldsymbol{w}}\in {{\mathbb{R}}}^{L}$, i.e., = u ^M + i + w . The MCA assumes that i can be considered as the texture-like features, which will be represented by, e.g., the local discrete cosine transform, and u ^M is regarded as the piecewise smooth component represented by the curvelet. After the default MCA decomposition, the components will also be extended to matrix forms ${{\boldsymbol{U}}}^{{\rm{M}}},{\boldsymbol{I}},{\boldsymbol{W}}\in {{\mathbb{R}}}^{P\times T}$ by supplementing zero vectors to the corresponding positions of their zero set channels. To retrieve finer details of the main peak, the upper root-mean-square envelope with a sliding window of length 128 points is applied, and U ^M is taken for example to interpret the retrieval strategy. Suppose the maximum value of the upper envelope of the signal after dedispersion of U ^M is e^u, and phases of the signal samples above the threshold e^u can be recorded. After a dispersion process, phases of samples in every channel can always be obtained by a simple phase shifting based on the recorded phases. In every channel, samples located on these phases should be recycled, and other samples will be set to zeros, resulting in a new profile U ^R with most entries being zeros. For components I and W , the thresholds are 3eⁱ and 3e^w, and the new profiles are I ^R and W ^R, respectively. Let R = U ^R + I ^R + W ^R and the matrix form of ${\boldsymbol{u}}\in {{\mathbb{R}}}^{N}$ be ${\boldsymbol{U}}\in {{\mathbb{R}}}^{P\times T}$; then a final profile P = U + R can be obtained. Figure 2 shows a flowchart for the algorithm of CS-Pulsar explaining the overall execution process from input of an observed signal profile F ₀ to output of a final signal profile P . A normalization is carried out for two reasons. (1) The advantage of restoration efficiency should be clearly clarified. (2) Shape comparisons should be carried out between the reference signals and the restored ones. Several dedispersed signals are defined to avoid unnecessary confusions.

• (1)  
  y ₀: the standard reference signal.
• (2)  
  y ₂: dedispersion result of F ₂ (after Steps 1 and 2).
• (3)  
  y : dedispersion result of F (after Steps 1, 2, and 3).
• (4)  
  y _F, y _H, y _D, y _C and y _Cu: dedispersion results of U by the Fourier, Haar and Daubechies wavelets, DCT, and curvelet operators. m _F, m _H, m _D, m _C, and m _Cu: dedispersion results of RFI mitigation by the five operators without considering the restored channels. p _F, p _H, p _D, p _C, and p _Cu: dedispersion results of P by the five operators.

Zoom In Zoom Out Reset image size

Because there is no unified calibration standard at present, the reference signals are applied only for shape comparisons. Without loss of generality, the maximum value of y ₀ will always be set as 1. Then, the maximum value of y is normalized as 1, and the other signals will be normalized according to their intensity ratios to y . The peak signal-to-noise ratio (PS/N: in dB) is used to measure the similarity of the restored signals. In order to show the restoration accuracy, a restoration ratio B_r /B is applied as an assessing criterion, where B is the area of a restoration result, and B_r is the area of restored signal contents.

Figure 3 shows PSR J1645-0317 as an example for mitigation and restoration. The number of the zero set channels is 74 with γ = 25, and a is set to 4. The maximum iteration number and the termination tolerance are set to 500 and 10⁻⁵, respectively. By visual inspection, after mitigation and channel restoration, flux intensity obtained by the curvelet operator shows a better continuity in frequency distribution, and the missing channels are sufficiently restored. Figure 4 shows comparisons between the different operators applied in the CS-Pulsar test in Figure 3 after mitigation, restoration, and signal retrieval. The standard signal of PSR J1645-0317 is the summation of all observed profiles between 2003 January 3 and 2009 December 31 with the total observing time 78,594 s. Signals y ₀ (top left) and y (top right) are normalized to 1, and y ₂ (top middle) and other signals are normalized according to their intensity ratios to y . Obviously, y ₂ widens the signal and shows false peaks due to the on-pulse RFI. In this test, y can also be considered as a result of the traditional incoherent dedispersion. From middle left to bottom middle: dedispersion results of RFI mitigation (black), mitigation with restoration (blue), and restoration and signal retrieval (red) by the five operators. The PS/N values of the pulses after mitigation are 26.8353, 26.3932, 26.8876, 28.1907, and 30.1723, respectively. Bottom right: dedispersion results of all the restored 74 missing channels by the five operators, and the restoration ratios are 0.1321, 0.1088, 0.1135, 0.1329, and 0.1375, respectively. The purpose is to show the significant effect of the function of signal restoration by visual inspection.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 5 shows the restoration and retrieval results of PSR J0332 + 5434 with the mode of three peaks. The number of the zero set channels is 97 with γ = 80, and a is set to 7. The maximum iteration number and the termination tolerance are set to 1000 and 10³, respectively. In cases of multipeaked pulse profiles, the on-pulse RFI will severely disturb the recognition of lower peaks after dedispersion (see Figure 5 left and right top). The main purpose is to show the mitigation efficiency in removing on-pulse RFI leading to false peaks. In the right bottom panel, the false peaks are successfully removed and the multipeaked signal shape is preserved in the process of mitigation and restoration. Another purpose is to demonstrate the retrieval result for the strongest peak. In the middle panel, the removed signal contents of the strongest peak can be reasonably retrieved with the aid of the upper root-mean-square envelope.

Zoom In Zoom Out Reset image size

The restoration ratios with γ = 40 using the curvelet operators for PSRs J1935 + 1616 (B1933+16), J0332 + 5434 (B0329+54), J0814 + 7429 (B0809+74), J0826 + 2637 (B0823+26), and J1136 + 1551 (B1133+16) are 0.0995, 0.1243, 0.0362, 0.0278, and 0.0286, respectively. The purpose is to further reveal the efficiency of CS-based restoration. These values also imply considerable information loss that should have been considered in astrophysical measurements. For a profile, the larger the value of γ is, the less channels are removed and the lower the restoration ratio is, and vice versa. For PSR J0814 + 7429, the corresponding restoration ratios of γ = 30, 40, and 50 are 0.1683, 0.0362, and 0.0104, respectively.

High-precision timing observations of millisecond-period pulsars have diverse applications, such as testing relativistic gravity (Kramer et al. 2006), investigating binary stellar evolution (Freire et al. 2011), and using pulsar timing arrays (Shannon et al. 2013) to detect gravitational waves emitted from supermassive black hole binaries. Astrophysical measurements require RFI mitigation as a necessary routine. If the CS-Pulsar method is applied to improve timing accuracy, it must not cause issues with timing-related analysis further down the pipeline. It must be proven to not cause systematic biases or underestimated uncertainties. To be an improvement on current methods, it should reduce the timing residuals and the estimated error.

Let v ₁ and v ₂ (in seconds) be residual vectors of the original profile and a processed one, and e ₁ and e ₂ (in seconds) be length vectors of their error bars, respectively. A residual vector reveals the goodness of fit of the rotation model of an estimated template to a standard one. The rms of a group of residuals can be applied to judge the respective timing accuracy. The standard templates corresponding to different a and operators are chosen as the summations of the observations without outliers, and the observations only contain residuals distributed like white noise. An error bar judges the estimated error or goodness of fit of the signal shape of a processed signal to a standard one, and its length is directly related to the S/R. The median error (ME) is defined as the median value of the lengths of a group of error bars to reveal the respective estimated error of TOA. In this subsection, the standard template-matching approach (Taylor 1992) is applied as the TOA measurement algorithm, and the toolboxes of PSRCHIVE and TEMPO2 are used.

The profiles of PSR J1935 + 1616 were collected from 2011 to 2014, and γ is set to $\tfrac{6\mathrm{MAD}\left({{\boldsymbol{F}}}_{2}\right)}{0.6745}$ in preprocessing. Each panel of Figure 6 shows a comparison of timing results of two groups of profiles, i.e., the profiles F ₀ and P . Parameter a of the soft threshold is selected as 1, 2, 3, 5, 7, and 9, and only the results of a = 1, 3, and 9 are shown in the left, middle, and right panels, respectively. Figure 6 reveals a reasonable trend that the residuals of P are reduced with the increase of a, and the advantage of P over F ₀ will be stabilized when a is large. Table 1 shows the values of the assessing criteria. When a is small, e.g., a < 3, the Gaussian noise and RFI still deteriorate the timing results. With the increase of a, the goodness of fit of the rotation model between the estimated and standard templates has been improved, and the processed signals have higher S/N values. In general, timing accuracy benefits from channel restoration and signal retrieval due to recovered key signal details, and the detrimental smoothing effect of mitigation can be largely avoided. Table 2 shows the timing performance of the profiles P by the five operators with a = 5.

Zoom In Zoom Out Reset image size