Millisecond Cadence Radio Frequency Interference Filters

To find outliers, we need both a central location and a measure of scale. The central location tells us the expected value of our data and is a useful replacement value for bad data. The measure of scale tells us the dispersion of the distribution, standard deviation being the most common measure of scale. RFI can be orders of magnitude brighter than the signal of interest, so we must use estimators that are robust to large outliers. There are many possible choices, see Fridman (2008) for an overview. We choose Median Absolute Deviation (MAD) because of its stability. For a set of samples, $X$, we define

Rousseeuw & Croux (1993) calculate the coefficient $\hat{k}=1.4826$ that relates MAD to standard deviation

From Eq (1) we see MAD will be stable as long as the bottom half of the distribution is unaffected by RFI. This stability gives MAD an advantage over other robust measures of scale that use data higher in the distribution. We considered other robust measures (Interquartile Range, Rousseeuw & Croux (1993) estimators $S_{n}$ and $Q_{n}$, and Biweight Midvariance (Beers et al., 1990)), but they either did not provide as stable results as MAD, had high computation cost, or proved difficult to implement.

The time domain MAD filter works to remove bright points in the dynamic spectra while preserving the astrophysical signal. Implementations of this filter have been used successfully in several experiments, for examples see, CHIME/FRB Collaboration et al. (2018); Buch et al. (2019); Rajwade et al. (2020). These implementations are bespoke to a particular instrument. Our aim is a more generic implementation that will work a high time resolution data from a variety of telescopes and sources.

Fig. 1 shows a section of dynamic spectra that includes three bursts from pulsar B0329+54 recorded using the GREENBURST backend (Agarwal et al., 2020; Surnis et al., 2019). GREENBURST is an L-Band transient search backend, the filtering performance here is typical for other telescope and backends. We want to remove the bright terrestrial sources while preserving the large pulses. We need to choose sections of dynamic spectra which are large compared to the RFI but smaller than the pulse. Such a block of data will provide robust statistics to remove and replace the RFI, while the pulse affects the statics and is preserved. We will operate two dimensional sub-blocks of the dynamic spectra.

The analogue receiver chain can add substantial frequency structure to the dynamic spectra. An example of this structure is show in the bottom right plot of Fig. 1. To perform outlier tests on two dimensional chunks of dynamic spectra, we must remove this structure, while preserving the the outlining RFI. If done correctly, this separates the RFI from the shape of the frequency response across the observing band. The RFI will then be removed as a non-Gaussian outlier.

We considered several fitting methodologies, some of which are shown in Fig. 2. The non-linear nature of RFI precludes filters such as high-pass or Savitzky–Golay (Savitzky & Golay, 1964). Iterative Chebyshev polynomial (Chebyschev, 1853) fitting, where an initial fit to all the points is used to exclude large outliers for a subsequent fit. B-splines (de Boor, 1978) fitted with a Huber loss function, where nearby points are fitted with a squared error loss and outlier points are fitted with a linear loss (Huber, 1964). We also the investigated asymmetrically re-weighted penalized least squares (arPLS) algorithm (Baek et al., 2015; Zeng et al., 2020). arPLS gives very good results, however its iterative matrix algebra is much slower than other filters. For this work, we settled on a rolling median filter. The kernel size is a adjustable parameter, which we provide a reasonable default. The median filter has the advantage of being very fast as well as being stable, were the polynomial and arPLS will sometimes diverge, leading to good data being flagged. Iterating the median filter to produce a cascade filter similar to Galassi et al. (2009) did not produce meaningful improvement.

Once the dynamic spectra has been detrended in frequency and time, we subdivide the dynamic spectra into subbands of $n$ channels. This allows for changes in dynamic spectra noise level across the bandpass, from a bandstop filter for example We then calculate the MAD and the median along the time axis of each of the subbands. We then use the rolling median filter to find the median of median and the median of the MADS along the time axis. This procedure is a recursive approximation to finding median and MAD of two dimensional blocks of spectra, while relaxing the requirement that the dynamic spectra dimensions be an integer number of subblock dimensions. We perform two iterations. The first iteration detrends using the median filter on the channel and time medians, preserving any RFI spikes while removing the intrinsic band and time structure. The MAD filter is run and pixels that are flags are replaced with Not a Numbers (NaNs). The NaNs are not considered during the subsequent detrends, this allows us to better model the good parts of the spectrum. We detrend the data again, now using only the medians along the two axes without the median filter. This iteration allows the filter to be more sensitive to smaller scale structure. The MAD filter is run again and the flagged values from both iterations are replaced by NaNs. A final detrend is followed by replacing NaNs with the desired value, which is now the median along both axes of the dynamic spectra. If the pulse has a low dispersion measure we can re-add the time trend, to preserve pulse power.

Comparing Fig. 1 and Fig. 3 shows the effects of the MAD filter. The frequency axis has been flattened, removing the effects of the signal chain. The time domain had likewise been flattened, however we re-added the total powers, to preserve the brightness of the the low-DM pulses. Almost all of the bright salt-and-pepper noise visible in Fig. 1 has been removed. The narrow band RADAR around channel 1300 has been greatly reduced, but is still present. This calls for frequency-domain filtering which we’ll discuss in the next section. Likewise, the bright narrow time RFI present around channel 3500 has been greatly reduced, but is still present. The communication band around channel 1900 has been completely removed. Importantly, much of the bright pulsar pulses remains. This indicates that the filter is working as expected.

For a given pulse morphology, larger DMs will cause the pulse to be spread out over a larger section of the dynamic spectra. This causes the MAD sections to be larger compared to the pulse. This can cause a DM dependence on the recovered pulse. However, as shown in Section 3 that for ($\mathcal{O}(100)$ S/N) candidates this effect is not measurable. As discussed in Amiri et al. (2021), FRB DM and pulse brightness are related. High DM FRBs will likely be less bright, and thus less likely to have signal excises by the MAD filter.

As seen in Figs. 1 and 3, we often encounter persistent periodic signals which can sometimes survive the time domain MAD filter. We would like to remove these anthropomorphic signals while preserving the periodic pulsar signal. This cleaning is best done in the frequency domain. We use the fast Fourier transform (FFT) along the time axis to make a representation of the dynamic spectra in Fourier space. The FFT give a real and imaginary components, or equivalently a magnitude and phase. We are interested in excising powerful periodic signals in each of the channels, but are uninterested when these signals start relative to the start of the dynamic spectra. Thus we can ignore the phase and flag the magnitudes. To flag the magnitudes, we roughly follow the same procedure as described in Section 2.2.

We take the magnitudes of the Fourier Transformed dynamic spectra, split these powers into frequency subbands. This split allows for changing noise and pulse structure across frequencies. We robustly fit in time and frequency to remove large trends show structure across the bandpass. The origin of these trends include changes in instrument sensitivity and stranding waves in the analogue chain. Red noise will leave a large spike in the Fourier domain time space. Once we have the detrended the magnitudes, we compute MAD across channels for each of the subbands. As discussed by Nita et al. (2007), for Gaussian noise the real and imaginary parts of the Fourier transform will be in turn be Gaussian distributed. The magnitudes will be $\chi^{2}$ distributed. The filter will be slightly more aggressive than pure Gaussian nose of a given sigma.

Broadband signals (such as FRBs and Pulsars) will be present across many channels. When only considering the magnitudes, broadband signals will be aligned across the band, since their arrival time is controlled by the phase. These astronomical signals will increase the median and the noise level uniformly and will be protected from the filter while narrow band RFI, and will only show in a few channels and will be excised. This filtering methodology contrasts with that of Maan et al. (2021), who also filter in the Fourier domain. Maan et al. (2021) use statistics on per-channel blocks in Fourier space flagging blocks that have outlying differences in mean or standard deviation. Our method is more forgiving to broadband signals. This is advantageous for bright pulsars, at the cost of excising broadband periodic noise (such as power lines). We will filter these broadband signals with our final filter.

Frequency channels that are used used for communication, but don’t contain periodic signals can still benefit. Digital communication consists of steps between on and off. These steps are abrupt and causes considerable ringing in the Fourier domain. The Fourier MAD filter removed the largest peaks of the ringing leading to a smoother channel in the time domain. This smoothing can be seen between Channels #1500-1700 in Figs. 1 & 3. After the filtering, these channels are much smoother, lowering the noise level of the dynamic spectra.

Once magnitudes above a defined threshold are flagged, we must make a choice about replacement. The easiest replacement option is zero. This can be problematic because sharp edges in the Fourier domain can lead to ringing in the time domain since we are computing finite FFTs. We did not find these effects noticeable, so we found zero replacement satisfactory. Alternatively we could rescale the powers or use Least-Squares Spectral Analysis (LSSA) to handle missing samples (Chapter 13 Press et al., 2007).

We see the effects of this filter on B0329+54 in Fig. 4. The bright radar pulses have been removed, additionally the communication RFI has been greatly reduced. Bright salt and pepper noise remains. The three pulses are largely the same, losing only a little power in the sub-pulse that precedes and follows the main pulse. We also see that the bandpass shape is preserved. For the time domain MAD filter, we removed this because it provided no useful information and makes replacement values difficult to determine. We did not have to do this flattening step with this filter because the total is determined by the first component for the FFT. We preserve this component and thus do not require the use of a replacement value. In this sense, the Fourier MAD filter can be seen as a smoothing filter, spreading out RFI power over the channel. This is acceptable because we are interested in change in power level, not the absolute level.

The first highpass filter for pulsar searches was described by Eatough et al. (2009), in which the dynamic spectra is collapsed over frequency to form a zero-DM time series. This time series is then subtracted off the dynamic spectra, thereby removing signals that are common across all channels at a given time sample. This is very effective at removing structureless RFI sources, such as spark plugs and unshielded electrical equipment.

This highpass filter worked well for the Parkes Multibeam Pulsar Survey (Manchester et al., 2001), which has one bit dynamic range and 288 MHz bandwidth made for uniform telescope response. If RFI is strong enough to flip one bit, the RFI was likely to show up across all the channels. More modern receivers can more than a GHz of bandwidth and often have more than eight bits of dynamic range. These combined give rise to more complex zero-dm RFI than was seen in Eatough et al. (2009), as shown in Figure 5. Bandstop filters and band roll off coupled with changing RFI intensity leads the the RFI being under-subtracted in some frequencies and over subtracted in others. To account for this we extend methodology of Eatough et al. (2009) into a more general highpass filter. We FFT along the frequency axis, replace the lowest $n$ Fourier components with zeros and FFT back, this filtering is equivalent to Eatough et al. (2009) for $n=1$. As both then remove the DC component across frequencies. As discussed by Eatough et al. (2009), this filtering will reduce the pulse energy as a function of dispersion measure. Low DM pulses occupy a large portion of the dynamic spectra. These pulses can then lose a significant portion of their energy to the filter. As discussed in Eatough et al. (2009), this energy loss can be partially recovered by using convolution kernel that takes into account the how the filter changes the shape of the time series. However, we don’t know of any search pipelines that have implemented these kernels.

This ideal high pass filter’s jump discontinuity will introduce Gibbs ringing. We tried tapered transition windows, but the improvement was minimal and this made the filter less intuitive to use. As discussed in the next section, this filter is used almost exclusively with the MAD filter, which detrends frequency structure, reducing the size of the low Fourier components before this filter is run. As shown in the right most panel of Fig. 6, some ringing is visible, however it is near the noise level, and thus not of great concern.

Each of the above filters operates in a different domain, so it is often useful to combine the three of them. We show such a case in Fig. 6 where we have multiple bright RFI sources. We begin by running the time domain MAD filter. This removes salt-and-pepper RFI and detrends the data in time and frequency. As shown in the second panel of Fig. 6, periodic RFI is still bright and present after the MAD filter. We then apply the FFT-MAD filter to remove this persistent narrow-band RFI. The third box of Fig. 6 shows the now periodic RFI has been removed. The long sweeping pulses are zero-DM RFI. They show a sweep because these blocks are dispersed. As shown in Figure 5, mean subtraction does not fully remove the effects of this signal. We then apply the high-pass filter, removed the first six modes to remove the broad-band RFI. Shown by Figure 6, this process greatly reduced the structure in the dynamic spectra, revealing candidates that would be missed if only one filter was used.

This filter order was chosen because the MAD filter detrends the dynamic spectra and remove large spikes, linearizing the dynamic spectra which make the operations in Fourier space more stable. We noticed zero-DM subtraction alone would often increase the standard deviation of de-dispersed time series. This is possibly due to very bright narrow-band RFI, which as we see in Fig. 6, can be significantly brighter than broadband RFI, being smeared out in frequency by the zero-DM subtraction. Hence we put the highpass filter last, after all the narrow-band RFI is removed.

We used will³³3https://github.com/josephwkania/will (Kania & Bandura, 2022) to create synthetic pulses and inject them into Gaussian noise that has the same MAD and median as a GREENBURST dynamic spectra. The pulse creation will be described in detail in the next section. For this section and the rest of the paper, we used the pulse detection in will. We do a pulse detection of a chunk of dynamic spectra by dispersing the dynamic spectra around the pulse, perform a boxcar convolution and detrend and measure the peak in the time series as well as the noise level. This later measurement is an additional reason why we used will, as common single pulse search pulse packages don’t report noise levels. This noise measurement is needed because we are interested in measuring the interaction between noise level and recovered pulse characteristics An additional reason for using our own search implementation is for a standard single pulse search would be extended over many DM and width trials. Candidates are then clustered in DM–time–pulse width space, with one value reported for the cluster. RFI can affect this clustering. See Aggarwal et al. (2021) for a discussion of clustering methodologies. Since we know the pulse and width of the injected pulse, so we know the optimal width and DM to search, this allows us to skip the clustering step and simplify our analysis.

We cleaned the dynamic spectra with our four filters. The effects of filters are shown in Fig. 7. We see that over DM and pulse pulse intensity, the filters behave as expected. The time domain MAD filter removes parts of very bright pulses, however some part of the pulse persists and is recoverable. The MAD filter drives the response of the composite filter bright pulse regime. The FFT MAD filter had no dependency on pulse power, the pulse magnitudes line up and the pulse protects itself. The highpass filter also behaves as expected. The signal removed by this filter is passes on how the fraction of the total bandwidth the pulse occupies and this of course is not a function of pulse power.

The response over DM similarly behaves as expected. Low DM pulses occupy a higher percentage of the band and more of the pulse gets removed by the highpass filter. At higher DM, the pulse occupies less of frequencies and less of the pulse is removed. This DM-response filter shows us that filtering introduced some variability, on order of 10%, to the recovered S/N. This variability was expected, but we are surprised by size of the effect. This variability comes from two sources. First, the background and the pulse are randomly generated. These realizations will have different random bright points both on and off the pulse that will get flagged. The flagging will cause random variation around the mean injected S/N. The second effect, for the MAD and composite filters, is as the pulse changes in DM, the pulse is spread out over different numbers of MAD sub-arrays. If the pulse only occupies a small portion of the sub-array, more of the pulse could be flagged. If the pulse occupies a large portion of the sub-array, then the pulse statistics will drive the sub-array statistics and none of the pulse will be flagged. This filter-pulse interaction is determined by the characteristics of the pulse, the pulse phase relative to the filter sub-arrays, and DM. We partially marginalize over these last two qualities in Fig. 7 by using multiple pulses randomly phase shifted relative the filter and multiple DM averaged into a bin. We use the standard deviation of S/N in the DM bins to produce one sigma errors.

We can use these two plots to determine what filtering scheme would be useful. If we are interested in low DM pulses, we should use a less aggressive high-pass filter, or forgo this filtering step. If we expect bright pulses, we would forgo MAD filtering which has a poor response to pulse where S/N is above a few hundred. The former decision can be made by even during blind searched by consulting estimated dispersion measure modules, such as Cordes (2004) and Yao et al. (2017) for a given pointing direction and estimating a minimum reasonable DM.

We used three hours of GREENBURST data from a single pointing that had been search but and did not have any bright astrophysical pulses. This observation was done at night which usually has less RFI but during a period where the telescope elevation was low which can allow more RFI into the telescope. GREENBURST records filterbank files (i.e. a discretely sampled array of frequency and time data from the telescope) that that are 597 s long. We took twenty two of these filterbank files and used them as the base to inject pulses. For each file we samples the parameter space described in Table 1, all the pulse properties except for the individual pulse powers are the some for a given file. See Agarwal et al. (2020) for the physical insight behind the parameter ranges. Each file had 288 pulses interjected. We repeated this processes 115 times for a total of 33120 pulses. This base files where used multiple times, but our analysis is focused on the effects on burst morphology and not the effects of sky position.

We create pulses that are Guassian along the time and the frequency axis. We convolve the time profile with an exponential scatter function $e^{-t/\tau}/\tau$ where $\tau$ is the scatter time in Table 1. For the frequency profile we multiply the Gaussian by a spectral index $(F/F_{\text{center}})^{\alpha}$, where $\alpha$ is the spectral index, $F$ is the channel frequency and $F_{\text{center}}$ is the reference frequency, we choose the centre channel. The radio light travels through ionized gas on its path from the source to the telescope. This plasma in plasma has random structure. The structure creates the structure, causes amplitude variations across frequency (Burke & Graham-Smith, 2009). We simulate scintillation by multiplying the frequency profile with the factor

where $N_{s}$ is half of the number of bright patches. We injected all the scintles at the same phase, since we do not anticipate this having a significant effect. Finally we take into account the receiver roll off and bandstop filter. These are frequencies where the analogue chain has low sensitivity where we are unlikely to see an astrophysical signal. We module the rolloff by calculating the medians along the band-pass, we then smooth these medians with a median filter. We calculate the standard deviation of this smoothed band-pass. We flag channels that are within 2/3 of a standard deviation of zero power. The transition between on and off regions is made using a Gaussian filter.

The previous paragraphs describe the continuous functions used to make the pulse profile. We need discrete values to inject into the the dynamic spectra. To transform from continuous to discrete, we sample continuous function for a given number of samples. The number of samples will be proportional to the pulse brightness and is reported in Table 1. We reference this number and not S/N because we do not have a way to know the noise level without RFI because we cannot perfectly remove all RFI. Based on Mickaliger et al. (2018), we used a log-normal distribution of pulse powers, but did not include an exponential tail. Since we precisely know the injected pulse powers, we can see how the recovered pulse energy distribution is affected by both RFI and RFI filters.

First we sample the space described by Table 1, which we use to create a pulse profile. From this profile we calculate the ideal boxcar width. We sample this profile and inject this into the file.

We sequentially clean this file with each one of the filters, and perform the same pulse detection as described above. We record these values and provide a summary in Table 2. The first column is the average pulse amplitude, $\overline{\text{Amplitude}}$, that is the mean amplitude of the boxcar convolved pulse. We see that each filter has some loss and composite filter inherits these losses. The next column, $\overline{\text{Noise}}$ is the average of the convolved noise levels at each pulse window. We see that both outlier filters provide a substantial reduction in the noise level. The average S/N over all the pulses is $\overline{\text{S/N}}$.

To test the distribution of injected against that of the recovered pulses we used two statistical tests, the Kolmogorov-Smirnov (KS) and the Anderson-Darling (AD) two sample tests. The goal of both these tests is to quantify the difference between the injected and recovered pulse distributions. In the ideal case, the RFI filtering should bring these distributions closer together by removing spurious signals that distort our view of the recovered pulse without degrading the pulses. Since we are injecting pulse samples we need to find a way to map these values to S/Ns. We take inspiration from Lilliefors Test (Lilliefors, 1967). Lilliefors test is a normality test for when the mean and variance is unknown. The test divides by the variance and subtracts the mean. Likewise we can make the powers and the recovered S/N unit variance and zero mean, and then apply both two sample tests. This will change the degrees of freedom, so we cannot report significance values without numerical simulation of these tests.

We see that we have mixed results. For some RFI filters, particularly the FFT-MAD, we see this values go down indicating that the two distributions are closer to each other. We see that this is not the case for the MAD filter, the Anderson-Darling value increases, while the Komogrorov-Smirnov value decreases. The Anderson-Darling statistic, $A^{2}$ for two continuous distributions $F$ and $F_{n}$ is defined by Anderson & Darling (1952) as

This equation shows that the Anderson-Darling test is weighted toward the tails of the distribution. Our test these tails will be high and low S/N pulses. This gives us an explanation why the KS and AD statistics diverge when applying the MAD filter. The MAD filter stabilized the noise level and separates low S/N candidates from the noise, this improves the KS statistic, which is equally weighted among all the pulses. This filter will clip bright pulses, reducing the returned fidelity of the S/N high tail and increasing the AD static. Which filter is most useful partially depends on the characteristics of the pulses of interest. We see that the outlier filters do a good job of increasing the occupancy, that is the percent of pulse windows that have a S/N above six. However in in this test, the composite filter is too aggressive. The % flagged tells us the amount of data that is flagged in each domain, the composite % flagged is the sum of the percent flagged of reach filter. The high cadence filters can make significant improvements to the recovered pulse characteristics while flagging only a small amount of of the data.

The final test makes use of real pulses in GREENBURST data. Pulsar will have more complex burst structures than previously modeled. We’ve chosen four pulsars that cover a range of dispersion measures and brightness as shown in Table 4. These pulsars are bright, and will have many visible single pulses.

As in the previous subsection we report the $\overline{\text{Amplitude}}$, $\overline{\text{Noise}}$, $\overline{\text{S/N}}$, and stacked S/N. We again see that the filters can significantly improve both the individual pulse S/N and the stacked S/N. However there is no one-size-fits-all solution. For bright pulsar J1935+1616, the FFT-MAD filter does the best, but the time domain MAD filter removes more of the pulse than is ideal. This result is as we expect from Figure 7. A less expected result is J0358+5413, this pulsar is about half as bright as J1935+1616 but we see the same effect with the latter. This is because much of J0358+5413’s energy is in very bright (above 200 S/N) pulses. This is discussed further in Appendix A. These bright pulses get clipped. The less bright pulses are not removed, and their S/N increases as expected. This gives mixed result for the MAD filter, the number of windows with a pulse above six sigma increases while the average pulse S/N decreases.

PSR J1705–1906 has a low DM of 23 cm^-3 pc, this causes the pulse take up a significant fraction of the bandwidth. The highpass filter subsequently removes 10% of the pulse power. This result again shows the need to consider the source while choosing the filtering. Instead of three modes, we should have excised one or zero Fourier modes.

Unlike the case with synthetic pulses, we do not know the intrinsic power distribution, so we cannot perform the KS or AD tests. Instead we will use the Spearman’s rank correlation coefficient to measure the relationship between a pulses S/N and the noise level. We choose the Spearman rank coefficient because it describes how well the relationship between two variables can be described by a monotonic function. We expect for a given signal brightness that S/N and noise level will be monotonically decreasing. In the ideal case, S/N and noise level will be uncorrelated, the noise will be Normal and any changes in pulse S/N will reflect intrinsic changes to the source brightness. If these variables are correlated, then S/N is partially measuring the evolving radio frequency environment, and the pulse measurements as reported by S/N will be less meaningful.

We see that this raw value is negative, as expected. With filtering, we can bring this metric closer to zero. This trend to zero is an indication that the S/Ns better reflect the intrinsic pulse energy and is less affected by the varying radio frequency environment.

Even though these are bright pulses, the filtering noticeably increases the percent of windows with pulses. This means that search pipelines should be able to find more pulses. These pulses can then be used to better understand the source.

The percentage flagged is higher than for the simulated pulses. This reflects that these observations took place at times throughout the day, leading to a worse RFI environment. The composite filter removes less than ten percent of the data. This small amount of flagging is an advantage of filtering on the shortest timescale. Flagging methods that remove blocks of data will often flag higher percent of the data. More flagging means more signal lost to the filter.