A FPGA-based Fast Converging Digital Adaptive Filter for Real-time RFI Mitigation on Ground Based Radio Telescopes

Radio Frequency Interference (RFI) is a growing concern in the radio astronomy community. Wireless communication devices have never been more ubiquitous, and sensing devices such as automotive radars for automation (cruise control) and safety (pre-collision airbag deployment) are using high frequency bands not used for commercial applications before. Nowadays we can find consumer applications using frequencies from a few KHz to 81 GHz and beyond, most of them on mobile devices (Hasch 2015). The radio astronomy community is well aware of the problem and has developed a number of different methods to cope with various types of RFI (refer to (Series 2013) for a comprehensive list). Interferometers have additional capabilities to mitigate RFI sources by synthesizing a beam null over the source, A technique called null-steering (Fridman 2005), but it requires knowledge of the source position in time.

Single-dish telescopes are particularly susceptible to RFI coupling through antenna sidelobes. For single-dish telescopes several methods have been described in literature based on flagging and excision, real-time blanking, and frequency flagging (Niamsuwan et al. 2005; Nita et al. 2007; Winkel et al. 2007; Offringa et al. 2010). All these methods produce some degree of data loss or require assumptions to be made on the astronomical signal. Currently, ad-hoc methods are used for each type of RFI, and no general solution has been adopted to cope with the problem.

A more ambitious approach is to identify, characterize, and subtract the RFI from the astronomical signal in real-time by using adaptive filters (Barnbaum & Bradley 1998; Briggs et al. 2000; Kesteven et al. 2005). The main pro of this approach is that it could deal with RFI with burst and persistent characteristics, broad or narrow band, and tackle fast moving sources or interference with unknown direction of arrival, provided that the adaptive filter converges in timescales shorter than the changes of the RFI. Nevertheless no reports on the speed of convergence of such adaptive filters are found in literature, and only modest bandwidths have been reported.

Here we report the development of a real-time, frequency domain digital adaptive filter implemented on a state of the art Field Programmable Gate Array (FPGA) capable of processing 4096 spectral channels in a 1 GHz of instantaneous bandwidth. The filter is able to cancel a broad range of RFI types and quickly adapt to changes on the RFI source characteristics, minimizing the data loss without any assumption on the astronomical or interfering signal properties.

Zoom In Zoom Out Reset image size

The filter architecture is based on the work of (Kesteven et al. 2010). The design assumes that the RFI signal can be independently sampled without overlap of the astronomical signal. Specifically, it is assumed that the astronomical signal and the sampled RFI are uncorrelated. Since RFI is normally many orders of magnitude stronger than the astronomical signals, and the reference antenna can be made to have high gain over the horizon with almost no gain on the sky, the above assumption is easy to meet for most astronomical sources and RFI environments (see Figure 1). Nevertheless, moving, high elevation RFI sources might require a tracking directional reference antenna.

To sample and subtract the RFI, both the primary and reference signals have to be digitized by two Analog to Digital Converters (ADCs) sharing the same clock source, i.e., with a couple of synchronized ADCs. After digitalization, two pipelined Fast Fourier Transforms (FFTs) of 4096 spectral bins are calculated. Synchronized ADCs and FFTs have been implemented for digital signal processing in radio astronomy, particularly applied to digital sideband separation (Finger et al. 2013, 2015). In our design the ADCs run at 2 GSPS, to sample a DC to 1 GHz baseband.

We use the following convention to name the different signals within the system:

• A_n: astronomical signal frequency sample n.
• I_n: interference signal frequency sample n.
• R_n: interference signal measured in the reference antenna.
• G_i: update output at cycle i (after N samples accumulations).
• F_i,n: filter output sample n after i update cycles.

Figure 2 shows the block diagram of the filter.

Zoom In Zoom Out Reset image size

After digitalization two pipelined FFTs of 4096 spectral bins are calculated. The reference channel is then conjugated and multiplied by the output signal. The filter itself has two sequential cycles: accumulation and update. During the accumulation cycle, N samples are added while the update register keeps the value calculated in the last accumulation cycle. The accumulation length N is a user selectable parameter. After each accumulation cycle, G_i is updated to be the sum of its last value plus the accumulator output, scaled by the constant $\varepsilon /N$. The parameter ε is also user selectable. The update of G_i, in turn, updates the filter's output. Finally the filter output is integrated for 536 ms to mimic the operation of an astronomical back-end.

For the mathematical analysis, it is assumed that the primary and the reference signals are purely sinusoidal, of constant amplitude and phase, which in the frequency domain can be represented by complex values.

With the above definitions we can write the following differences equations:

$$\begin{eqnarray}&&{G}_{0}=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{G}_{i+1}=\varepsilon \displaystyle \frac{1}{N}\displaystyle \sum _{j={Ni}}^{N(i+1)-1}{F}_{i,j}{R}_{j}^{* }+{G}_{i},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{F}_{i,n}={A}_{n}+{I}_{n}-{G}_{i}{R}_{n},\end{eqnarray} \tag{ 3 }$$

where (1) is the initial condition of the filter, (2) and (3) are the update and the filter output equations. To solve this system of equations, we first replace (3) in (2).

$$\begin{eqnarray}&&{G}_{i+1}=\varepsilon \displaystyle \frac{1}{N}\displaystyle \sum _{j={Ni}}^{N(i+1)-1}({A}_{j}+{I}_{j}-{G}_{i}{R}_{j}){R}_{j}^{* }+{G}_{i},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&=\,\varepsilon \displaystyle \frac{1}{N}\displaystyle \sum _{j={Ni}}^{N(i+1)-1}{A}_{j}{R}_{j}^{* }+{I}_{j}{R}_{j}^{* }-{G}_{i}{R}_{j}{R}_{j}^{* }+{G}_{i},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&=\,\varepsilon \displaystyle \frac{1}{N}\mathop{\underbrace{{\displaystyle \sum }_{j={Ni}}^{N(i+1)-1}{A}_{j}{R}_{j}^{* }}}\limits_{={C}_{i}}+\varepsilon \displaystyle \frac{1}{N}\displaystyle \sum _{j={Ni}}^{N(i+1)-1}\mathop{\underbrace{{I}_{j}{R}_{j}^{* }}}\limits_{K}\\ &&\quad -\,\varepsilon \displaystyle \frac{1}{N}{G}_{i}\displaystyle \sum _{j={Ni}}^{N(i+1)-1}\mathop{\underbrace{{R}_{j}{R}_{j}^{* }}}\limits_{r}+{G}_{i},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&=\,\varepsilon \displaystyle \frac{1}{N}{C}_{i}+\varepsilon \displaystyle \frac{1}{N}\displaystyle \sum _{j={Ni}}^{N(i+1)-1}K-\varepsilon \displaystyle \frac{1}{N}{G}_{i}\displaystyle \sum _{j={Ni}}^{N(i+1)-1}r+{G}_{i},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&=\,{G}_{i}-\varepsilon {{rG}}_{i}+\displaystyle \frac{\varepsilon }{N}{C}_{i}+\varepsilon K,\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&=\,\mathop{\underbrace{(1-\varepsilon r)}}\limits_{:= b}{G}_{i}+\displaystyle \frac{\varepsilon }{N}{C}_{i}+\varepsilon K,\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{G}_{i+1}={{bG}}_{i}+\displaystyle \frac{\varepsilon }{N}{C}_{i}+\varepsilon K.\end{eqnarray} \tag{ 10 }$$

where in (6) we use the fact that I_n and R_n are correlated by K. In (9) we assume that r changes slowly compared to the filter update cycle i, so we can define $b:= 1-\varepsilon r$. This yields a first order difference equation in (10); its solution is given by

$$\begin{eqnarray}&&{G}_{i}=\displaystyle \frac{\varepsilon }{N}\displaystyle \sum _{k=0}^{i-1}{b}^{k}{C}_{i-k-1}+\varepsilon K\displaystyle \frac{1-{b}^{i}}{1-b},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&=\,\displaystyle \frac{\varepsilon }{N}\displaystyle \sum _{k=0}^{i-1}{b}^{k}{C}_{i-k-1}+\displaystyle \frac{K}{r}(1-{b}^{i}).\end{eqnarray} \tag{ 12 }$$

Choosing ε such that $| b| \lt 1$ ensures the filter converges. Replacing the obtained expression for G_i in (3) we get

$$\begin{eqnarray}&&{F}_{i,n}={A}_{n}+{I}_{n}-{R}_{n}\left(\displaystyle \frac{\varepsilon }{N}\displaystyle \sum _{k=0}^{i-1}{b}^{k}{C}_{i-k-1}+\displaystyle \frac{K}{r}(1-{b}^{i})\right),\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&=\,{A}_{n}+{I}_{n}-{R}_{n}\displaystyle \frac{K}{r}(1-{b}^{i})-\displaystyle \frac{\varepsilon }{N}{R}_{n}\displaystyle \sum _{k=0}^{i-1}{b}^{k}{C}_{i-k-1},\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&=\,{A}_{n}+{I}_{n}-\displaystyle \frac{{R}_{n}{R}_{n}^{* }}{r}(1-{b}^{i}){I}_{n}-\displaystyle \frac{\varepsilon }{N}{R}_{n}\displaystyle \sum _{k=0}^{i-1}{b}^{k}{C}_{i-k-1},\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&=\,{A}_{n}+{b}^{i}{I}_{n}-\displaystyle \frac{\varepsilon }{N}{R}_{n}\displaystyle \sum _{k=0}^{i-1}{b}^{k}{C}_{i-k-1}.\end{eqnarray} \tag{ 16 }$$

From (16) it can be concluded that for $i\to \infty$ the interference I_n gets removed, as long as the filter converges faster than the change of the RFI. The speed of the filter is controlled by ε, where $b\approx 0$ (i.e., $\varepsilon \approx 1/r$), yields the fastest convergence. The third term represent the filter added noise. The noise term scales with $\varepsilon /N$, so added noise can be reduced at expenses of convergence speed by increasing N.

Notice that no assumption was made about the signals frequency so this monochromatic analysis can be extended to any number of spectral bins. Implementing one filter per each spectral channel allows processing any type of broad or narrow band signal.

The described adaptive filter was implemented in a ROACH2⁵ (Reconfigurable Open Architecture Computing Hardware) development board designed by the CASPER group (Collaboration for Astronomy Signal Processing and Electronics Research) led by Berkeley University. The most important components of the board are the Xilinx Virtex-6⁶ XC6VSX475T FPGA, and a couple of 5GSPS ADC boards (Jiang et al. 2014) manufactured by ASIAA (Academia Sinica Institute of Astronomy and Astrophysics), which use the e2v EV8AQ160⁷ ADC chips. The FPGA is programmed via the MATLAB Simulink software, using specialized libraries provided by Xilinx and CASPER group (Parsons et al. 2006).

A test setup was assembled to characterize the adaptive filter performance with different types of RFI. A band-limited noise signal was used as astronomical source, and an Arbitrary Waveform Generator (AWG) with an up-conversion stage was used as RFI source. Figure 3 shows the test setup block diagram.

Zoom In Zoom Out Reset image size

Two types of RFI were used to test the filter: A broadband noise-like signal and a stronger and complex multiple-carrier signal. The broadband signal was generated via software while the multiple-carrier signal was produced from RFI measured at the Five Hundred Meter Aperture Spherical Telescope (FAST, (Nan et al. 2011)) site. The RFI at FAST was recorded with a reference antenna mounted on top of one of the towers that holds the feed horn cabin. Figure 4 shows the RFI test signals while Figure 5 depicts the astronomical test signal.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 6 and 7 show the astronomical signal contaminated with RFI and the filtered output. No sign of the RFI can be seen in the filtered output for both broadband and multi-carrier RFI. For the filter parameters it was chosen N = 1 (no accumulation in the filter) in order to achieve maximum convergence speed. For ε it was given a value that roughly satisfy $\varepsilon {r}_{\max }\approx 1$, where r_max is the maximum power among all channels of the reference signal, which is different between each type of interference: 2¹⁵ for the Multiple-carrier RFI and 2¹¹ the broadband RFI.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 8 shows the difference between the astronomical signal and the filter output. The difference is in general within +/−1 dB, except for specific channels where differences above 2dB can be seen. These channels are correlated with ADC spurious signals. The filter output total power is 4% less than the astronomical signal for the broadband RFI and 9% less for the FAST RFI.

Zoom In Zoom Out Reset image size

Figure 9 shows the speed of convergence of the filter. T = 0 is defined as the moment the update loop is closed with a ${G}_{0}=0$. Convergence, measured for the average power (green) line, as the time needed to fall to a 1/e (36.8%) of the T = 0 level, is reached on 45.2 μs for the broadband RFI and 27.2 μs for the multiple-carrier RFI. Convergence to a 1% level is reached in 208.1 μs for the broadband RFI and 125.5 μs for the multiple-carrier RFI.

Zoom In Zoom Out Reset image size

• 1.  
  A fast-converging adaptive filter was implemented on a FPGA achieving the real time processing of 4096 spectral channels on a 1 GHz bandwidth.
• 2.  
  Using no assumptions on the astronomical signal properties, it was analytically demonstrated that the filter converges to an output clean of RFI.
• 3.  
  Experiments were done in laboratory using broadband and narrowband RFI test signals. In both cases, the output equals the clean astronomical signal with negligible residual contamination.
• 4.  
  The time of convergence (for a decrease to a 1%) was 208.1 μs for a broadband noise-like RFI signal and 125.5 μs for a multiple-carrier RFI signal recorded at the FAST radio telescope, demonstrating that the filter can adapt to very fast-changing RFI, like the one produced by burst transmissions from fast moving sources.
• 5.  
  Embedded in an astronomical back end, the reported filter should be able to cope with broad and narrow band RFI with negligible data loss.
• 6.  
  The adaptive filter may also find applications in communications on highly RF-polluted environments, and in defense, to counteract GPS or radio communication jammers.

This work was supported by the Chilean National Commission for Scientific and Technological Research (CONICYT) through its grants CATA-Basal PFB06, ALMA 31150012 and FONDECYT 11140428, and by the Chinese Academy of Science South America Center for Astronomy (CASSACA). D.L. and D.R. also acknowledge the support from the International Partnership Program of Chinese Academy of Sciences, Grant No.114A11KYSB20160008. We thank Xilinx Inc. for the donation of FPGA chips and software licenses.