A demonstration of the effect of fringe-rate filtering in the hydrogen epoch of reionization array delay power spectrum pipeline

ABSTRACT

1 INTRODUCTION

Fluctuations in the brightness temperature of redshifted 21 cm line emission from neutral hydrogen provide a unique probe of cosmic structure in the early Universe (Scott & Rees 1990; Subramanian & Padmanabhan 1993; Furlanetto, Oh & Briggs 2006; Morales & Wyithe 2010; Pritchard & Loeb 2012; Mellema et al. 2013; Pober et al. 2014), from the cosmic ‘Dark Ages’ following last-scattering (⁠|$z \approx 30 - 1000$|⁠), through Cosmic Dawn (⁠|$z \lesssim 30$|⁠), and the Epoch of Reionization (⁠|$z \approx 6 - 15$|⁠). Other probes, such as the optical, infrared, and ultraviolet emission from early stars and galaxies (Doré et al. 2016; Robertson 2022), and molecular line emission (Kovetz et al. 2017), suffer from obscuration by dust (Fudamoto et al. 2021), interlopers/confusion (Cheng et al. 2016; Lidz & Taylor 2016; Gong, Chen & Cooray 2020; Chen & Pullen 2022), and more restrictive limits to their detectable redshift range. Developing practical approaches to measuring the 21cm fluctuations is therefore of great utility to the cosmology and galaxy formation communities. Despite over a decade of concerted effort using multiple generations of radio telescopes, a robust detection of this signal at high redshift (⁠|$z \gtrsim 6$|⁠) has remained elusive (Mirocha et al. 2019; Liu et al. 2022). The principal reason is the large dynamic range between the 21cm signal itself, which is rather faint (typically of order mK), and the radio emission from other (foreground) sources, measuring from tens of K at |$\sim$|GHz frequencies to thousands of K at |$\lesssim 100$| MHz (Santos, Cooray & Knox 2005).

In principle, the 21cm fluctuations and the foreground sources have highly distinctive spectral and angular signatures, making it possible to separate them by applying appropriate cuts or filters in these domains (Santos et al. 2005). In practice, however, the imperfect calibration of any given radio telescope gives rise to spurious couplings and modulation of the signals, blurring the distinction between them (Morales et al. 2012). Since the foregrounds are so much brighter, even very small calibration errors can leak spurious foreground emission into otherwise 21cm-dominated regions of the data, so that recovery of the 21cm signal is severely hampered.

A wide variety of signal processing methods have been used to try and make headway with this problem. A very common approach is to transform the frequency-dependent data into a spectral basis, such as a Fourier basis or principal components of the data’s frequency–frequency covariance matrix (e.g. Liu et al. 2009; Parsons et al. 2012; Wolz et al. 2014; Alonso et al. 2015). This exploits the differing spectral behaviours (and different brightness/signal-to-noise) of the 21cm signal and foregrounds to provide some degree of separation between them, albeit imperfectly due to coupling/modulation effects mentioned above. For individual interferometer baselines, this can be thought of as applying a filter in the delay (⁠|$\tau$|⁠) domain, where delay is the Fourier pair of frequency, |$\nu$|⁠, with Fourier basis functions |$e^{2\pi i \nu \tau }$|⁠. Filters to remove foreground contamination are typically most aggressive at low delay, |$\tau \sim 0$|⁠, but may also affect higher delays. This leads to some partial suppression of the 21cm signal too (otherwise known as ‘signal loss’).

The ‘foreground avoidance’ approach uses filters that separate out the most contaminated regions of the spectrum, including the 21cm signal in those regions, permitting only the nominally less-contaminated modes to be used in analysis. Conversely, the ‘foreground removal’ approach attempts to model and remove the contamination while leaving the 21cm signal as intact as possible. Both approaches have serious drawbacks; foreground avoidance is lossy, and can still allow some degree of contamination or statistical biases if couplings between foreground-dominated and 21cm signal-dominated modes are sufficiently large (Barry et al. 2016; Jensen et al. 2016; Ewall-Wice et al. 2017), whereas foreground removal can introduce additional artifacts into the data unless the subtracted model is very accurate (Pober et al. 2016; Offringa et al. 2019).

The spectral dimension of the data is not the only one that can usefully differentiate between the different signal components however. Baseline-dependent filtering of visibilities (or angular filtering of reconstructed images) can help separate contamination from our Galaxy (mostly at large angular scales) from extragalactic emission (on all scales) (Charles et al. 2023). Particular weighting and filtering schemes can also be used at the calibration step to suppress spurious features caused by antennas or baselines that are more susceptible to them.

In this paper, we focus on the design and effects of filters in fringe-rate (⁠|$f$|⁠) space (Roshi & Perley 2003; Parsons & Backer 2009; Parsons et al. 2016). This is the Fourier pair of local sidereal time (LST, denoted by |$t$| in this paper), which for a drift-scan telescope is a proxy for Right Ascension (terrestrial drift-scan telescopes observe a stripe of the sky at fixed Declination, with increasing RA as the Earth rotates). The Fourier transform of time sequences of visibilities was first used in Peckham (1973) to distinguish between fringe-rates as part of a technique for phase correcting interferometers. In 2003, filtering of fringe-rates was used by Watson et al. (2003) to remove the Sun and Moon from observations of Jupiter using the very small array (Rusholme 2005). Fringe-rate filtering was applied to observations from the precision array for probing the epoch of reionization, used for the generation of 21cm power spectrum and upper limits on the Epoch of Reionization (EoR) (Parsons et al. 2014; Ali et al. 2015). It is intended to be a used on observations from the hydrogen epoch of reionization array (HERA; DeBoer et al. 2017) as we will discuss. Other telescopes are adopting the technique, for example the Very Large Array has used fringe-rate filters for radio frequency interference (RFI) excision (Helmboldt et al. 2019).

The rate at which sky sources pass through the fringes of an interferometer baseline depends on the baseline length, and orientation of the baseline relative to the source path over time. Fringe-rate filters are typically applied to each baseline separately. Sources rotating with the sky pass through the fringe pattern of a baseline more rapidly for longer baseline lengths, which have fringes closer together; the visibilities generated by a source will oscillate (in phase) faster in time for those baselines, and will appear at high-fringe-rates. Sources further from the celestial poles have faster apparent motion on the sky, and so will have higher fringe-rates than sources closer to the poles. The fringe-rate can also be positive or negative depending on which side of the pole a source is at (with respect to the centre of the primary beam), and tends to decrease as sources get closer to the local horizon due to projection effects. Importantly, this means that fringe-rate can be used as a crude proxy for source position, making it possible to filter or down-weight emission in the sidelobes of the primary beam, which contaminate Fourier modes needed for detecting the 21cm signal (Thyagarajan et al. 2015). Removal of sky locations from the primary beam is in effect ‘beam sculpting’, demonstrated in Parsons et al. (2016). The ‘mainlobe’ filter to be described below is an example of this.

Another application of fringe-rate filtering is to separate sources of emission that are locked to the sky’s rotation from those that are not (e.g. Kern et al. 2019, 2020). Stationary terrestrial sources will not fringe for example, so should be found at around |$f \approx 0$|⁠, whereas aircraft, satellites, and other transient phenomena might be expected to give rise to signals that pass through the fringe pattern much faster than sky rotation would allow. The ‘notch’ filter, described below, removes fringe-rates around |$f \approx 0$|⁠, and the mainlobe filter can be used to deal with non-sky high-fringe-rates. These filters will also help to excise unwanted systematics that produce spurious fringe-rates, particularly if the systematic has a differentiated spectral behaviour that can also be isolated in delay space. Finally, the thermal noise contribution to the data changes randomly from one time sample to the next. For white noise, this gives rise to a signal that has similar power at all fringe-rates. By filtering out fringe-rates outside the range corresponding to sky-locked signals like the 21cm emission, some of this noise can be suppressed.

In forthcoming analyses of data from HERA, fringe-rate filtering will be adopted as a standard part of the observation analysis pipeline. This was driven by the presence of mutual coupling systematics that are well-localized in delay-fringe-rate space, and demonstrations showing that fringe-rate filters are useful in mitigating these (Kern et al. 2020). This can make calibration easier; for example, Charles et al. (2024) demonstrated its use in suppressing calibration errors arising from poorly modelled diffuse emission, showing roughly an order of magnitude improvement in calibration fidelity. More general beam-sculpting, i.e. the mainlobe filter, along with the notch filter, present interesting options for enhancing the use of fringe-rate filters in HERA data processing.

With the benefits of fringe-rate filtering also come a number of drawbacks. Applying any filter can lead to some degree of signal loss (of the 21cm signal), and introduce correlations over time, which must be rigorously characterized to avoid invalidating upper limits on the 21cm power spectrum. Applying a filter in fringe-rate space corresponds to convolution of the (transformed) filter in LST space, so signals will be smeared in time, and noise will become correlated in time. These correlations can subvert common approximations that are used in data analysis pipelines, for instance the way in which the signal-to-noise ratio is expected to grow as progressively more time samples are averaged together. Some of these issues are addressed in Cheng et al. (2018) and Kolopanis et al. (2019), and analysis needs to be updated based on the filtering described in this paper.

The temporal resolution and duration of the observations also limits the resolution of the filter, so that naive Fourier filtering gives rise to unwanted aliasing or ringing effects. We will describe below how this influenced the implementation of filtering in HERA, based on discrete prolate spheroidal sequences. Finally, flagging due to RFI and other data quality issues complicates the transformation between the LST and fringe-rate spaces; poor handling of the discontinuities that are introduced by flagging can also lead to ringing artifacts. The use of discrete prolate spheroidal sequences can tolerate some missing data (Ewall-Wice et al. 2021), but we do not examine that in this paper. Inpainting methods such as Gaussian constrained realizations are being considered to deal with flagged data (Kennedy et al. 2023).

There are other methods that, like fringe-rate filtering, use the modes present in time sequences of visibilities to manipulate those visibilities for particular aims. |$m$|-mode analysis uses time-based spherical harmonic m-modes to decompose the sky and enable foreground subtraction (Shaw et al. 2014, 2015). Correlations in visibilities over time can be used to statistically separate the EoR signal from foregrounds (Paul et al. 2014; Patwa, Sethi & Dwarakanath 2021). These and other methods provide alternatives to fringe-rate filtering.

In this paper, we apply two different types of fringe-rate filter to simulated observations. A simulation with a sky model consisting of a single point source and an East–West baseline are used to build intuition about fringe-rate filtering; these build on the examples in Parsons et al. (2016) but show the use of discrete prolate spheroidal sequences. Following on from these, full sky simulations using 10 HERA dishes are passed through certain key steps of a ‘standard’ HERA analysis, of the kind used in HERA Collaboration (2022). We examine the outputs of the steps of the analysis, including the final cylindrically averaged delay spectra, to observe the effect of the filtering.

This paper does not provide a rigorous statistical analysis of the effects of fringe-rate filtering and in particular does not address signal loss; these subjects will be dealt with elsewhere (e.g Pascua et al. 2024).

The paper is organized as follows: Section 2 provides a brief pedagogical introduction to fringe-rate filtering, using simple examples based on point source sky models to build intuition for how the filters work. Section 3 summarizes the visibility simulations we use in the rest of the paper, while Section 4 defines the mock data analysis pipeline that we apply to the simulations to act as a simplified model of the HERA pipeline. Section 5 presents our main results, which show how the products of each stage of the simplified data analysis pipeline are affected by fringe-rate filtering. Finally, we conclude in Section 6.

2 FRINGE-RATE FILTERING

Parsons et al. (2016) contains a detailed description of fringe-rate filtering theory and methods, which we will summarize here in the context of our experiments and the implementation of fringe-rate filtering in the HERA pipeline.

A single interferometer baseline correlates the voltage signal received by two antennas, producing a complex visibility for each time sample and frequency channel. We write the visibility for a baseline |$\boldsymbol{b}$|⁠, observing frequency |$\nu$|⁠, and time |$t$|⁠, as

where |$I$| and |$A$| are the sky intensity and primary (power) beam in direction |$\boldsymbol{n}$| on the celestial sphere, the integral is over all directions on the sky that are above the local horizon, and we have suppressed the frequency dependence of the terms for brevity. For a drift-scan interferometer such as HERA, sources on the sky rotate overhead as the antennas maintain a fixed pointing at the zenith. The received intensity of a source depends on the antennas’ directional sensitivity (or ‘primary beam’), and the phase of the visibility depends on the source direction relative to the baseline. It is the phase that is most important for fringe-rate filtering.

A simplified example of how the phase changes as a source drifts through the primary beam is shown in the top-left panel of Fig. 1. A single source at a declination of 0|$^\circ$| is passing over an East–West oriented baseline situated at the equator. The source rises at the East horizon at an LST of 6 h, passes through the zenith at 12 h, and sets on the West horizon at 18 h. In this example, we set the antennas to have uniform beams, so that the observed amplitude of the source is constant over time. The phase, however, changes over time, and is indicated by showing the real part of the visibility; the imaginary part changes in a similar manner. The oscillations of the phase are the fringes, and the speed of oscillation is the fringe-rate. We see from this illustration that the fringe-rate is not constant – it is lowest when the source is at the horizon, and fastest when the source is at zenith.

Figure 1.

The effect of the mainlobe filter on the visibility for a single source that transits the zenith, as observed by a purely East–West baseline at the equator, with a uniform primary beam. The left column shows the amplitude and real component of the visibility, with fringes clearly indicated by the oscillation of the real component. The right column shows the amplitude of the fringe-rates of the visibilities (the Fourier transform along the time axis). The top row shows the unfiltered visibilities, the second row shows them after filtering with the mainlobe filter, and the third row is the difference between the two. The mainlobe filter specifies that only fringe-rates within the bounds indicated by the grey vertical lines should be retained; the shaded region indicates those fringe-rates that should be excluded. The filter has been implemented using discrete prolate spheroidal sequences, leading to some leakage outside the bounds, as discussed in the text.

Open in new tabDownload slide

If we had performed a Fourier transform along the time axis in a short-time window, to a good approximation we would have recovered a single tone, i.e. a spike at a particular fringe-rate, for this source. The specific fringe-rate that is observed depends on the baseline length and orientation, but also the time |$t$| around which we have placed our short window. Fourier transforming over the whole time range would instead produce a collection of spikes that have been drawn out over a range of fringe-rates due to the variation in fringe-rate between horizon and zenith (Fig. 1, top-right panel).

We can now use the varying fringe-rate to construct a filter that depends on distance from the zenith. From our example, we can see that the fringe-rates with the highest absolute value, at |$f \simeq -1$| mHz, correspond to when the source was near zenith. By applying a tophat filter in Fourier space that select only modes around |$f \simeq -1$| mHz (indicated by the grey vertical lines), we can suppress the observed intensity of the source away from zenith, as shown in the middle panels of Fig. 1. The procedure can be generalized to any source declination, array latitude, baseline length, and orientation etc., as well as to include the effect of the primary beam, which also modulates the amplitude of sources as they rotate through the beam. In this example, the failure to remove some modes outside the filter bounds (i.e |$f \simeq -0.5$| mHz), is due to the use of discrete prolate spheriodal sequences, and will be discussed subsequently.

Following Parsons et al. (2016), we can make some approximations that permit a simplified interpretation of filtering in fringe-rate space. By relating the change in baseline orientation with time to the rotation of the Earth, and making an approximation that the fringe pattern varies faster with time than the primary beam term, we can replace the primary beam term in equation (1) with an effective ‘filtered’ beam (Parsons et al. 2016),

where

Here, |$w(f)$| is a weighting in fringe-rate space corresponding to the chosen fringe-rate filter, |$\tilde{\gamma }$| is the inverse Fourier transform of a tapering function that determines how long the time window used for the fringe-rate transform is, and |$\omega _{\oplus }$| is the angular velocity vector corresponding to the Earth’s rotation. If the time window |$\gamma$| is broad, |$\tilde{\gamma }$| is highly peaked in the fringe-rate domain, and so the convolution is localized around |$f \approx (\nu /c)(\boldsymbol{b} \times \omega _{\oplus }) \cdot \boldsymbol{n}$|⁠. For a given baseline vector |$\boldsymbol{b}$|⁠, this determines a fringe-rate for each direction on the sky, |$\boldsymbol{n}$|⁠. Fringe-rates are not unique to each direction however; instead, there are regions of constant fringe-rate which form rings on the sky (Parsons et al. 2016). These can be mapped more or less crudely to different regions of the primary beam, hence the interpretation of fringe-rate filtering as beam sculpting.

Note that some baseline orientations (i.e. mostly North–South aligned ones) have |$\boldsymbol{b}(t) \times \omega _{\oplus } \approx 0$|⁠. For these, most regions on the sky are collapsed into a narrow range around |$f = 0$| mHz, making fringe-rate filtering impractical for these baselines.

2.1 Fringe-rate filters for HERA

In this work, we will focus on two types of fringe-rate filter that will be used in forthcoming HERA analyses. The first is the ‘mainlobe’ filter, which implements a tophat window in fringe-rate space that rejects components of the observed visibilities with fringe-rates that fall outside the mainlobe of the primary beam. This down-weights the primary beam sidelobes and horizon, which suffer most from chromatic effects that are an important source of systematic contamination, particularly when bright sources such as Fornax A pass through them. By removing these fringe-rates, we can hope to reduce the amount of data that has to be flagged due to systematic contamination. This filter can be thought of as a very basic form of beam sculpting. We do not attempt to re-weight the primary beam by applying a shaped filter however; it is purely a cut on the range of fringe-rates to retain only ones that are within the beam mainlobe.

The second is a ‘notch’ filter, which is intended to remove a narrow range of fringe-rates around |$f \approx 0$| mHz. These correspond to components of the signal that are not fringing, i.e. those that are stationary with the array. While some of the sky signal can fall within this range (e.g. near the celestial poles), these fringe-rates are also contaminated by systematic effects like mutual coupling (Kern et al. 2019).

The filters are specified as rectangular tophats in fringe-rate, with centre |$f_\mathrm{ c}$| and half-width |$\Delta f$|⁠. For the mainlobe filter this specification defines the fringe-rates that are to be kept, and the filters vary with baseline and frequency. For the notch filter, this defines the fringe-rates to be rejected, typically a small range around 0 mHz, and the filter specification is the same for all baselines.

In the ideal case, the two tophat filters could be implemented by Fourier transforming the visibilities along the time dimension, masking in/out the appropriate range of fringe-rates, then and transforming back. Several complications arise for real data which means that a more sophisticated implementation must be used however.

Of most concern is missing data. For any given set of observations, a substantial fraction of times will be flagged due to RFI, poor calibration solutions, or system problems. The missing data appear as sharp features to a Fourier transform, and so generate ringing that scatters parts of the signal across fringe-rate space. A related phenomenon is due to the non-periodicity of the data, which is seen by the Fourier transform as another discontinuity at the band edges.

To reduce ringing artifacts, a tapering function such as a Blackman–Harris window (Harris 1978) can be applied to the visibility data to enforce periodicity or apodise the edges of the flagged regions. Alternatively, more sophisticated spectral analysis methods can be used that avoid direct Fourier transforms and instead use basis functions with properties that reduce the risk of ringing. HERA uses the latter kind of filtering method, based on discrete prolate spheroidal sequences (DPSS; Mathews, Breakall & Karawas 1985; Moore & Cada 2004; Ewall-Wice et al. 2021). These are a set of orthonormal basis functions that are concentrated in Fourier space for a given tophat window. In other words, they are constructed to minimize the total amount of power that is able to leak outside the tophat region in the fringe-rate domain.

However, DPSS are not all created equal, in terms of how well they concentrate power. For this reason it is possible to create a set of DPSS that are a complete basis for a sequence of visibilities; the sequence can then be exactly decomposed into DPSS, much like a Fourier transform for Fourier modes. For filtering, we generally only want to retain a few DPSS that are most concentrated within the filter window. The more DPSS that are retained, the better the approximation to a tophat filter in fringe-rate with sharp edges, but the worse the leakage outside the tophat window region. The HERA filtering pipeline has several options available to choose the number of DPSS to retain in an automated fashion, e.g. based on a minimum permitted concentration ratio (Ewall-Wice et al. 2021).

Once we have the reduced set of DPSS, we transform each of them back into the LST domain and use a linear least-squares fitting algorithm to find their best-fitting coefficients. This approach naturally handles missing data and non-periodicity at the band edges. The filter is only an approximation to a tophat, and will have some degree of leakage outside the tophat region and slight suppression of power inside the region. For clarity, note that a taper is not applied in the time dimension when performing the DPSS transform, and so some edge effects can still be present.

Because the filter is implemented through this fitting procedure, there is no simple, unique, filter profile that can be plotted. We will, however, present examples of the filter in action, which illustrate a typical filter shape when comparing unfiltered and filtered data.

Finally, we note that the shape of the filter is also limited by the resolution in the fringe-rate domain, which is determined by the length of the time/observing window, |$\Delta t$|⁠. Longer time windows permit higher resolution, i.e. |$\delta f = 1/\Delta t \approx 0.05$| mHz for |$\Delta t = 6$| h.

The advantages and disadvantages of DPSS for filtering 21cm observations is discussed more fully in (Ewall-Wice et al. 2021). A simple example of fitting DPSS, for the mainlobe filter, is given in Appendix  A. We use the eigenvalue cutoff method for selecting highly concentrated DPSS; its meaning, the value we used, and its consequences, are explored in Appendix  B.

2.2 Examples of filters applied to idealized situations

At the beginning of Section 2 (and shown in Fig. 1), we discussed the effect of a mainlobe filter on the highly simplified situation of a single point source transiting through the zenith as observed by an East–West baseline at the equator with a uniform primary beam. Returning to this example for a more detailed examination, we see that most of the signal lies between fringe-rates |$-2.2 \lesssim f \lesssim 0.5$| mHz. The predominantly negative rate is due to the direction of Earth’s rotation. The few positive fringe-rates in this example are an artifact of the Fourier transform to fringe-rate domain, as we have not used a taper to enforce periodicity, and the uniform primary beam also fails to provide any tapering of its own. In general, some positive fringe-rates are expected to be observed from sky emission, however corresponding to sources on the far side of the celestial pole relative to the beam centre.

The determination of the fringe-rate limits is based on beam power, which is low towards the horizon and highest at zenith. These also correspond to locations where the fringe-rates are low and high, respectively. The fringe-rates are plotted against the beam power, and cut based on the 5 per cent and 95 per cent percentile of the beam power. The low and high fringe-rates become the filter limits. The high limit of 95 per cent is used simply based on the evidence of Fig. 1, where it cuts out the low fringe-rates after the peak; the behaviour of the filter on the baselines of an actual telescope should be examined to see if this is reasonable. Other schemes are also possible, such as basing the low fringe-rate limit on beam nulls, but these have not yet been implemented. The same scheme and limits are used for all the mainlobe filtering experiments in this paper.

To complement this example, we also plot the effect of the notch filter on the same scenario in Fig. 2. The mainlobe filter retains fringe-rates near the maximum rate, and the notch filter rejects fringe-rates near |$f \simeq 0$| mHz, but the filtered data has features in common, for our simple examples. For the notch filter, the amplitude of the visibilities has been reduced closer to the horizon (LSTs of 6 and 18 h) where the fringe-rate tends to 0 mHz. The amplitude suppression kicks in much closer to the horizon than for the mainlobe filter, however, which is to be expected as the notch filter is much narrower and therefore has an effect that is more localized on the sky. The notch filtered signal in the fringe-rate domain (second row, Fig. 2) shows that a clear notch has been cut out around 0 mHz, but the notch feature does not have a sharp cutoff due to the use of the DPSS fitting method described above. This is analogous to the leakage outside the filter bounds in the mainlobe filter.

Figure 2.

Similar to Fig. 1 but now for the notch filter. The notch filter specifies that fringe-rates around |$f = 0$| mHz should be excluded, indicated by the shaded region. The filter was implemented using discrete prolate spheroidal sequences, as discussed in the text.

Open in new tabDownload slide

Despite using the DPSS fits, additional oscillatory structure has still been imposed on the fringe amplitudes after filtering, which can be seen clearly in the amplitude difference plots (bottom row). The bottom-right panel of Fig. 2, for example, shows the source of some of this structure – there is a low-level ripple at negative fringe-rates far outside the notch, and a tail at positive fringe-rates. These tails are an artifact of the DPSS fits, which trade-off filter sharpness and accuracy of signal recovery (i.e. lack of signal loss) inside the filter region against the level of leakage outside the filter region. This trade-off is controlled by the number of DPSS that are retained in the filter, which can be tuned depending on what one is trying to achieve with the filter.

We do not show the effect of the filters on a North–South baseline because they can be described simply. For a North–South baseline at the equator, with a source passing over the equator, their is no fringing, i.e. the fringe-rate is 0 mHz. The notch filter removes fringe-rate 0 mHz and so removes all the visibilities, setting the baseline to 0 over all time. Sky locations cannot be identified by fringe-rate for this baseline, so the mainlobe filter has filter bounds centred at 0 mHz and a half-width set by the fringe-rate resolution (or some other method); it passes all the visibilities without alteration.

Having gained some intuition as to the use and effect of filters in simple situations, we now turn to an examination of the filters that are being used in the HERA pipeline, and how these affect the 2D power spectra of more realistic visibility simulations.

3 VISIBILITY SIMULATIONS

In this section, we describe the simulations that will be used in the remainder of the paper. These are intended to be reasonably realistic as far as the visibilities are concerned; we include a full sky model containing point sources and diffuse emission, a realistic HERA beam model, and simulate a set of baselines for an HERA-like array using 30 dishes in the same location as the HERA telescope. The simulations cover a broad range and many samples in time and frequency to permit high resolution in the spectral (delay/fringe-rate) domains. For simplicity, we do not include flagging, bandpass or gain calibration errors, or non-redundancy of the antennas or baselines, nor do we attempt to model the systematic effects that are one of the principal reasons for implementing fringe-rate filtering. The parameters of the simulations are summarized in Table 1.

Simulated visibilities are generated using the hera_sim and matvis packages (Kittiwisit et al. 2023). We use a catalogue of 20 000 point sources that mimic the distribution of sources in the GLEAM catalogue (Hurley-Walker et al. 2017), but fill-in areas of the sky that GLEAM does not cover (for more details see Choudhuri, Bull & Garsden 2021). The point source flux densities for this catalogue range from |$0.2-82.4$| Jy, with spectral indices from |$-1.01$| to |$-0.63$|⁠. A diffuse emission model is obtained from the pygsm package (Price 2016), with the emission represented as a HEALPix map per frequency, with NSIDE = 128. This results in |$\sim 27$| arcmin pixels, which are sufficiently small compared to the characteristic angular scale of the longest baseline (⁠|$\lambda /D \approx 96$| arcmin for the 77m baseline at 140 MHz) that we expect discretization effects to be subdominant in our results. As an approximation, the visibility simulators treat the centre of each pixel effectively as a point source. We do not include a 21cm component in the simulations as it is much fainter than the foreground components, and we do not study the effects of the filters on it in this work.

The subset of 30 HERA dishes with its close-packed hexagonal layout shown in Fig. 3 (Top). The dish separation is 14.6 m (the shortest baseline length), with a longest baseline of 73 m. We label baselines using their constituent antennas, e.g. baseline (3, 6) runs East–West and has length 48.3 m. The (E-field) primary beam for each antenna is a frequency-dependent dipole-like model derived from electromagnetic simulations described in Fagnoni et al. (2021) (the ‘Vivaldi’ beam). The beam for |$X$| polarization, at 100 MHz, is shown in Fig. 3 (Middle). A slice from East to West gives the profile in Fig. 3 (Bottom). We generate only the |$X$| (East–West aligned) polarization channel in the simulations.

Figure 3.

(Top): The antenna layout of the HERA-like array of 30 dishes used for the simulations. Individual antennas are numbered for reference. The minimum spacing is 14.6 m. (Middle): The model HERA E-field beam (see text) for polarization X, 100 MHz, with East (E) and North (N) marked. The beam amplitude is expressed in decibels relative to a peak of 0 dB. (Bottom): A slice through the model HERA beam from East to West as a function of altitude.

Open in new tabDownload slide

To permit high-resolution studies of the filters, we simulate large ranges in frequency and time. The time axis comprises 6 h of HERA drift-scan observations, with 10 s of integration time per time sample, starting from Julian date 2458 116.48, for a total of 2160 time samples. This spans a range in LST of 7.44–13.46 h. The frequency range is |$100-140$| MHz, comprised of 408 channels of width 97 kHz. Using these frequencies, the 21cm line would be observed over a redshift range of |$ z \approx 9.1-13.2$|⁠; real-world analyses tend to use smaller bands for power spectrum estimation however, to reduce the effects of cosmic evolution.

After simulating visibilities for the sky model, uncorrelated white noise is added, using the radiometer equation with |$T_{\rm sys} \approx T_{\rm sky}$| to set the noise level (see Choudhuri et al. 2021).

4 DATA PROCESSING PIPELINE

In this section, we describe the analysis pipeline that we apply to the simulated visibilities (we do not perform imaging). This is an idealized, cut-down version of the HERA pipeline that was used in HERA Collaboration (2022) and Aguirre et al. (2022). A flowchart summarizing the steps of the pipeline is shown in Fig. 4.

Figure 4.

Schematic of the data processing pipeline used for this work, showing the steps taken to produce a 2D power spectra, and the placement of the fringe-rate filtering step.

Open in new tabDownload slide

The purpose of the simulated pipeline is to apply a given fringe-rate filter to the input visibilities and then form per-baseline delay power spectra, followed by a sequence of averaging and binning operations to increase the signal-to-noise. The end result is a 2D power spectrum indexed by baseline length and delay which can be converted to cosmological wave modes |$k_\perp$| and |$k_\parallel$|⁠, although we retain baseline length (⁠|$|b|$|⁠) and delay (⁠|$\tau$|⁠) for plots in this paper. The 2D power spectrum can be further averaged into a spherical power spectrum if desired.

Simulated visibilities are fringe-rate filtered before being coherently averaged in time, in chunks of 300 s. The visibilities are rephased to zenith before averaging the chunks, in order to prevent time smearing, and then rephased back to a drift-scan observation. This coherent averaging step is designed to average down the noise while causing minimal signal loss, as described in Aguirre et al. (2022). The optimal chunk size for coherent averaging is baseline- and frequency-dependent, and should generally be significantly shorter than the beam-crossing time-scale (the time taken for a patch of sky to drift through the mainlobe of the primary beam). Chunk lengths of around 420 s have been proposed for HERA (Aguirre et al. 2022; HERA Collaboration 2022), but we made a more conservative choice of 300 s for this paper.

The pipeline is run twice, once with the mainlobe filter and once with the notch filter. In the finalized HERA pipeline, it is unlikely that the filters will be separated in this way, they will both appear in the pipeline at different stages before coherent averaging (for example separated by an LST binning step when multiple days are used). However, so that we can examine them separately, we run them separately.

Following coherent time averaging, pairs of baselines that are redundant with each other (having the same length and orientation) are Fourier (delay) transformed and multiplied to form delay power spectra. This step is implemented using an optimal quadratic estimator (see Liu, Parsons & Trott 2014; HERA Collaboration 2022), but with a uniform weighting, which is not optimal. A Blackman–Harris taper (Harris 1978) is used for the Fourier transform. Pairs consisting of the same baseline (‘auto-baseline pairs’) are excluded, as these would be affected by a noise bias. After coherent averaging, delay spectra that are redundant with each other are averaged together. Due to the large size of the data and computing constraints, we limit the sets of redundant delay spectra to be averaged to a maximum of 48 delay spectra. Then follows time-averaging followed by binning by baseline length.

Note that we run the pipeline on the simulated data both with (‘filtered’) and without (‘unfiltered’) the fringe-rate filtering step, so we can study the effect of the different filters by comparing unfiltered versus filtered. As we are also interested in the effect of the filters on the noise, and how this affects the final signal-to-noise of the averaged data, we also run the pipeline on simulated visibilities without noise added (‘noiseless’), so the results can be compared with the ‘noisy’ data.

5 RESULTS

In this section, we present the results of applying the mainlobe and notch filters to the simulated data, examining the data at different points in the mock analysis pipeline. We begin by examining the effect on the visibilities themselves (Section 5.1), then examine different stages in the delay spectrum averaging process (Section 5.2), followed by the final 2D power spectrum results (Section 5.3).

5.1 Effect of filters on visibilities

In this section, we compare the filtered and unfiltered visibilities obtained from the output of the coherent time averaging step, which uses a 300 s averaging window and therefore reduces the number of time samples to 72 from the original 2160.

Fig. 5 shows the visibilities (left column) and noise (right column) for baseline (3, 6) (in Fig. 3, top), which is East–West orientated and has length 43.8 m. The top row shows the visibilities over time at a fixed frequency of 100 MHz before and after the mainlobe filter has been applied. The second row shows the same data in fringe-rate space (i.e. after applying a Fourier transform in the time dimension, without a taper function). The third row shows the visibilities as a function of frequency at a fixed LST of 10.45 h, and the fourth row shows the same data in delay space, that is, after Fourier transforming the visibilities by frequency (using a Blackman–Harris taper), at fixed LST. We have restricted the delays to focus on the horizon (grey vertical lines); the full delay range for our simulations is |$-5154$| to 5129 ns.

Figure 5.

The visibilities and fringe-rates for baseline (3,6) before and after application of the mainlobe filter. The top row left column shows the amplitude of visibilities over the time range of the simulation, both unfiltered and fringe-rate (FR) filtered using the mainlobe filter, at a frequency of 100 MHz. The top row right column shows the noise in the visibilities. The second row shows the fringe-rate transform of those visibilities, on the left, and the noise in the fringe-rate transform on the right. The shaded region indicates the fringe-rates that should be excluded by the mainlobe filter. The third row left column shows the amplitude of visibilities over the frequency range of the simulation, at the midpoint of the simulation time range (an LST of 10.45 h), again showing visibilities on the left and the noise in the visibilities on the right. The fourth row shows the delay transform obtained from the third-row visibilities, and the noise in the delay transform. The grey lines show the horizon delay for this baseline.

Open in new tabDownload slide

To obtain the noise plots, we processed the noisy and noiseless visibilities identically and took the difference.

We see that the visibilities are smoother in time than in frequency due to the coherent time averaging step. The mainlobe filter causes substantial additional smoothing in time, which has also reduced the noise level, as anticipated. The visibility as a function of frequency (at fixed time) shows a significantly different structure following filtering, which now lacks some of the significant structures on scales of a few MHz which are seen in the unfiltered data, while the larger-scale and smaller-scale structures are less affected. The delay-space plot shows more clearly how the delay structure of the visibility (at fixed time) has changed as a result of the mainlobe filter, with a notable suppression of power around the horizon delays. This is the expected behaviour for a mainlobe filter; the sidelobes and horizon, which are also more chromatic and therefore appear at higher delay, have been suppressed.

The fringe-rate plots show that the fringe-rates inside the tophat are recovered with hardly any loss of signal, while some fringe-rates outside the filter bounds have survived, but at a significantly reduced level. In fact, the filtered signal at fringe-rates |$f \gtrsim -0.6$| mHz is mostly caused by leakage of the signal within the tophat bounds outside the filter region due to the DPSS fitting procedure, rather than being a residual of the original unfiltered signal at these fringe-rates. It can be seen that the structure of this leakage signal is quite different as a function of |$f$|⁠, compared to the unfiltered signal. The large peak around |$f \approx 0$| mHz has been removed however, and the smaller one at |$f \approx -0.6$| mHz has also been significantly reduced. At |$f \gtrsim 0.5$| mHz, the filtered signal level is sometimes slightly above the unfiltered signal level however, again due to the DPSS leakage. As discussed above, the level of leakage can be tuned by changing the number of DPSS that are retained in the fitting procedure, at the expense of making the tophat filter boundary less sharp, and potentially increasing the level of signal loss inside the tophat region. The noise in the fringe-rate domain has also been strongly suppressed outside the tophat region, but left more or less intact inside the tophat and immediately around its boundary.

Fig. 6 shows the same quantities as Fig. 5, but now for the notch filter. This time, the main effect of the filter has been to remove the large spike around |$f \approx 0$| mHz. There is some minor signal loss outside of the filter range, particularly at large |$|f|$|⁠, and as expected for the DPSS fit, the shape of the filter is not perfectly sharp. The noise as a function of time and frequency is broadly similar before and after filtering, reflecting the narrowness of the notch filter, which can only remove a small fraction of the noise. Both the time and frequency structure of the visibilities is quite similar before and after the filtering, with only small shifts in the spectral features and a small reduction in signal amplitude observed.

Figure 6.

As for Fig. 5, but the notch filter is used.

Open in new tabDownload slide

Table 2 gives means and standard deviations for the real component of the noise in Figs 5 and 6. These show that the mainlobe consistently reduces the noise by almost one half, whereas the notch filter has little effect, as can be seen from the plots.

5.2 Effect of filters on averaged delay spectra

HERA uses a delay power spectrum approach (Parsons et al. 2012). Visibilities are Fourier transformed to delay (⁠|$\tau$|⁠) space, where we expect the foreground signal to be confined to a region within the ‘horizon delay’. This is the maximum delay in the time of arrival of a wavefront between the two antennas of the baseline, which occurs for a source at the horizon. Outside of the horizon delay, we would ideally expect only the 21cm signal and noise to be present, and so can ‘avoid’ foregrounds by analysing only these modes.

To estimate the delay spectra, cross-power spectra are formed by multiplying the delay-transformed visibilities from redundant baselines (same length and orientation), under the idealized assumption that redundant baselines observe the same sky visibilities, plus differences only due to noise. The cross-power spectra from all pairs of redundant baselines are then averaged at each observing time, then averaged over time, before finally averaging over baselines of the same length (but different orientation). This process follows a similar estimation procedure used in HERA Collaboration (2022). In this section, we will show the results of this process using unfiltered and filtered visibilities. Regarding Fig. 4, we begin at the delay spectrum estimation step.

At this step, we apply a delay transform to each baseline at each time, using a suitable tapering function (in this case a Blackman–Harris window), and then multiply with a suitable weighting defined by the optimal quadratic estimator formalism to form the cross-spectra. We follow a similar approach as in recent HERA analyses and use a uniform weighting, which does not produce optimal estimates, but which avoids signal loss issues; see Pascua et al. (2024) for a more detailed explanation. Note that baselines are not paired with themselves in this analysis, i.e. ‘auto-baseline pairs’ are excluded.

Following the initial formation of the delay spectra, we average together the cross-spectra of redundant baselines, at each time. This is the ‘incoherent redundant average’ step. Following this, the per-redundant-group delay spectra are averaged in time. Both of these averages can be performed in a way that respect the noise weights of each data sample, but for our simulations the weights are all uniform. We do not attempt to model the weighting/covariance between time samples introduced by the fringe-rate filtering procedure here, but note that any filter will induce correlations between previously independent time samples. This is important for errorbar estimation of the time-averaged delay spectrum. Errorbar estimation for the filtered data is a complex topic; a cursory examination is presented in Wilensky, Brown & Hazelton (2023) which needs to be expanded to cover the full HERA array and use realistic fringe-rate parameters. The final step of the pipeline we discuss here incoherently averages delay spectra into bins of the same baseline length.

We will track the pair of baselines (3, 6) and (8, 11) through the process, these are both 43.8 m East–West baselines in the same redundant group.

Fig. 7 shows the delay spectra for our chosen baseline-pair/redundant group for data that have had the mainlobe (left column) and notch (right column) filters applied. For these delay spectra the full delay range is |$-5154$|–5129 ns, but we show the range |$-2000$|–2000 ns as only flat noise extends beyond this range. The first row is the raw delay spectrum for the baseline pair without any averaging (apart from the coherent time average that was applied to the visibilities), and at a time half-way through the observation period; the second row is after averaging over the available baseline pairs within the redundant group at the same time; the third row shows the subsequent unweighted average of the previous step over all times; and the fourth row shows the average of the previous step and all other redundant/time-averaged cross-spectra of the same baseline length. We show both filtered and unfiltered visibilities. Delay spectra are complex valued and averaged as such, but it is an HERA convention to plot the absolute value of the real component, as is done here, for data and noise.

Figure 7.

Shows the delay spectra for baseline-pair (3,6) and (8, 11) as it proceeds through the pipeline, being combined with delay spectra from redundant baseline pairs (Group average), averaged over time (Time average), and combined with other delay spectra of the same baseline length (||$\vec{b}$|| average). The left column shows delay spectra generated using the mainlobe filter, the right column shows delay spectra generated using the notch filter. The vertical grey lines indicate the horizon delay. The stages of the pipeline are listed on the left.

Open in new tabDownload slide

For both filters we see that there is some alteration to the power within the horizon, with the peak shifting to 0 ns in the filtered data and reducing the power just outside the peak. The drop in power outside the horizon to the noise level occurs closer to the horizon for the filtered data, but likely due to other overall lower values for the filtered data. This drop is located well outside the horizon for both filtered/unfiltered cases, indicating a buffer region that should be ignored when examining power beyond the horizon.

After the power drops outside the horizon, the level of noise is fairly flat to the limit of the delay range, indicating we can ignore extreme delays and focus on the wedge for further analysis. The noise level for the filtered delay spectra is always below that of the unfiltered data for the mainlobe filter, but for the notch they are similar. For the mainlobe filter the noise level is below that of the unfiltered data at the beginning of the process (One baseline), but they tend to the same level towards the end of the process (⁠|$|\vec{b}|$| average). This is not seen for the notch filter, where the noise levels are about the same all the way through. Overall, the noise drops from one step to the next as averaging proceeds, but the area in and around the horizon, where foregrounds are present, increases.

5.3 Effect of filters on 2D power spectra

The delay spectrum from the previous section represents a single baseline length, and there will be several others corresponding to different lengths. Since different baseline lengths approximately measure different Fourier modes of the sky, when they are plotted vertically and stacked horizontally by baseline length, we obtain a 2D power spectrum, the last stage in the pipeline. The horizon delay increases by baseline length, and when plotted on the 2D power spectrum delineates a ‘wedge region’ in which the foregrounds are, ideally, contained, and outside of which exists, ideally, the 21cm signal and noise.

Fig. 8 shows the power spectra generated from the unfiltered, mainlobe-filtered, and notch-filtered visibilities. The first column contains the power spectra generated from noisy visibilities passed through the pipeline, the second column contains power spectra generated from noise only. The noise is extracted from the visibilities and passed through the rest of the pipeline. The horizon delay is indicated by the white lines. Table 3 shows the range of power in each power spectra.

Figure 8.

2D power spectra that are the end result of the data processing pipeline. The top row contains the spectra generated from unfiltered visibilities, the middle row contains the spectra generated from mainlobe-filtered visibilities, and the bottom contains the spectra from the notch filter. The left column contains the 2D power spectra, the right column contains the 2D noise power spectra, obtained from passing the noise in the visibilities through the pipeline.

Open in new tabDownload slide

The figure and table indicate that the power spectra are quite similar. Fig. 7 (last row) has hinted at the shape we expect to see in the baseline bins, with a fairly flat noise level and a roll off outside the horizons.

To more closely examine the difference between the unfiltered and filtered cases we take the ratio of the power spectra, shown in Fig. 9. The left ‘ratio spectrum’ is the result of dividing the mainlobe 2D power spectra by the unfiltered 2D power spectra. The right ratio spectra is the same, but for the notch filter case. Power ratios below 1 mean that the power is less in the filtered case than the unfiltered. The 2D power ratio spectra are similar although the notch filter shows a significant drop in power just outside the horizon, around fringe-rate 0 mHz, which the mainlobe filter does not. This may seem puzzling, since the mainlobe filter in our point source example for an East-West baseline (Fig. 1) would remove fringe-rate 0 MHz. However, the mainlobe filter bounds change depending on baseline orientation, and for baselines with a significant North-South orientation the mainlobe filter can retain fringe-rate 0 mHz.

Figure 9.

Each of these power spectra is the result of dividing power spectra that are the result of fringe-rate filtered visibilities, by power spectra produced from unfiltered visibilities. The left spectrum is the mainlobe 2D power spectra (Fig. 8, middle left), divided by the unfiltered 2D power spectra (Fig. 8, top left). The right spectrum is the same, but for notch filtered visibilities.

Open in new tabDownload slide

The ratio spectra also show that there are locations where the ratio is greater than 1, indicating that filtering has produced an increase in power. These locations are all outside the wedge and depend on the random noise added to the visibilities in the simulation; they may not appear at all, or change location.

To examine the ratios in more detail, Fig. 10 contains vertical cuts of the ratio spectra. Each panel represents a different baseline length. Within a panel, cuts for both the mainlobe and notch cases are indicated by the legend. The wedge horizon for the baseline length is shown by vertical grey lines. These show that at the horizon there is often a drop in power in the filtered spectra compared to the unfiltered; for example for baseline length 63.6 m there are two sharp drops at the horizon for the notch filter. We expect a drop at the horizon because the fringe-rates are around 0 mHz there, and the notch filter will remove these; the mainlobe filter will sometimes remove these, depending on the filter bounds, which depends on baseline orientation. We also see some instances where the ratio is above 1, for example for baseline length 63.6 m there are sharp increases just outside both horizons.

Figure 10.

Vertical cuts of the 2D power spectrum ratios in Fig. 9, at different baseline lengths. The vertical grey lines indicate the horizon delay.

Open in new tabDownload slide

The similarity between the power spectra is worthy of more investigation, and it appears that coherent averaging is at least partly the cause. Fig. 11 is generated from the same simulations and data processing pipeline as Fig. 9 except the coherent averaging step has been skipped. Similarly, Fig. 12 is the same as Fig. 10 with the coherent averaging step skipped. The mainlobe case looks similar in both figures (in that the ratios are variable), but the ratio values are at a lower level outside the wedge when coherent averaging is skipped. The notch case has a much more uniform ratio of 1 outside the wedge when coherent averaging is skipped, although there are still instances outside the wedge where the ratio is greater than 1. The uniformity outside the wedge indicates that the unfiltered and filtered power spectra are similar there. It is not surprising that the different filters have different effects on power spectra, but the interaction between the filters and time averaging requires further investigation.

Figure 11.

As for Fig. 9, generated from the same simulations and data processing pipeline, except that the coherent averaging step has been skipped.

Open in new tabDownload slide

Figure 12.

As for Fig. 10, generated from the same simulations and data processing pipeline, except that the coherent averaging step has been skipped.

Open in new tabDownload slide

5.4 Using more realistic HERA-like observations

The above results and analysis are based on simulations that use a frequency range of 100–140 MHz, |$ z \approx 9.1-13.2$|⁠, during which time the Universe will be evolving. To capture power spectra at a particular instant during the Universe’s evolution, we should use a smaller redshift range, while at the same time providing enough redshift bins for resolution of the spectra. When processing real observations, HERA will take power spectra from several frequency bands of width approximately 10 MHz, in this section we will use only the frequency range of HERA Band 4, 109.79–121.92 MHz, |$ z \approx 10.6-11.9$|⁠.

We also change the noise to simulate LST binning of visibilities over 100 d; LST binning will be used in the HERA pipeline and was part of the pipeline from which we extracted Fig. 4. At a particular LST we assume that the noiseless visibility will be the same every day, but random noise will be different. If the noisy visibility is averaged over those days, the noiseless visibility remains the same, but the variance of the noise, on average, reduces by a factor of 100. Therefore, we reduce the variance of the noise added to the visibilities by a factor of 100.

The power spectra that are produced from the reduced frequency band look almost the same as those in Fig. 8 except that the wedge is wider (due to the different band). For that reason we only show the unfiltered power spectra of the reduced band as an example, in Fig. 13. This shows that there is considerable leakage outside the wedge when the narrower band is used. Of more interest are the ratio plots (Fig. 14), and vertical cuts through the ratio plots (Fig. 15). The ratios of the filtered/unfiltered power spectra, Fig. 14, still show that the notch filter produces a significant dip in power outside the wedge that the mainlobe filter does not. Fig. 15, containing vertical cuts through the ratio power spectra, shows the dips produced by the notch filter more clearly.

Figure 13.

Shows the unfiltered power spectra that is produced from a reduced frequency band of 109.79–121.92 MHz. It is visually similar to the power spectra in Fig. 8, top left, but with the wedge widened due to the change in band.

Open in new tabDownload slide

Figure 14.

As for Fig. 9 but using the 2D power spectra from a reduced frequency band.

Open in new tabDownload slide

Figure 15.

As for Fig. 10 but the cuts are taken from the power spectra ratio plots in Fig. 14, which used a reduced frequency band.

Open in new tabDownload slide

6 CONCLUSIONS

The HERA is now the most sensitive radio interferometer targeting the 21 cm brightness temperature fluctuation signal at high redshift (HERA Collaboration 2023), and will likely remain so until the Square Kilometre Array-Low instrument comes online in several years from now. Forthcoming seasons of HERA data will unavoidably require Fourier filtering in both the frequency (delay) and time (fringe-rate) dimensions, due to the presence of foreground contamination and a variety of systematic effects such as mutual coupling that are relatively localized in these domains. The fringe-rate filters in particular can seem unintuitive, as they remove power in a way that is spread across the time axis, causing correlations between time samples that were not previously there, as well as introducing artifacts such as ringing and potentially some degree of signal loss. This interferes with traditional data analysis approaches that tend to split observations into smaller time windows (‘fields’) and average over time samples to build up sensitivity. The loss of intuition for how these averaging operations should work – that we are no longer able to simply combine independent data samples – risks making the analyses harder to interpret, or even leading to analysis choices that introduce errors, for instance in estimating the error bars or determining robust upper limits.

In this paper, we have used simple examples and a simplified mock data analysis pipeline applied to semi-realistic visibility simulations to build intuition into what the fringe-rate filters are actually doing to the data at each step of the pipeline. We consider two types of filter, both of them tophat-like – a mainlobe filter that retains parts of the signal with fringe-rates that fall only within the mainlobe of the primary beam, and a notch filter that removes fringe-rates around |$f \approx 0$| mHz corresponding to signals that are locked to the ground rather than rotating with the sky, which is symptomatic of certain classes of systematic effects.

We demonstrated the use of the filters on a simple setup of an East–West baseline on the equator, showing that the filters do as intended. The mainlobe filter alters fringe-rates to remove sky locations outside the primary beam, However, Fig. 1 shows that the profile of the resulting visibilities over time only approximately matches the beam profile in Fig. 3 (bottom), which it ideally should. Whether this is due to fringe-rate resolution, the manner of calculating the tophat width and location, the use of a tophat filter, the use of discrete prolate spheroidal sequences, or other factors, requires investigation, if the ideal is to be reached. It may be that the simple situation of a point source generating slowing changing perfect fringes with no amplitude deviations, is unsuitable for the intended application. In a real-world scenario, the fringe-rates of a baseline will have a more complex structure, generated from variable intensity sources all over the sky.

The fringe-rate filters applied to a real-world simulation demonstrate this complex structure (e.g. Fig. 5). The visibilities in this case contain realistic noise, which is smoothed over time and reduced by the filter. In this case, we have multiple frequencies which allow us to make delay transforms and observe the effect of the filter in delay space for a single time. While the fringe-rate filters operate over time, the visibilities are changed over frequency, altering the power within the horizon delays, adding power to some delays and reducing others. Of interest is that the peak of the delays within the horizon is shifted to |$\tau = 0$| ns if it is off 0 ns in the unfiltered data.

The use of discrete prolate spheroidal sequences shows that they do not function as a true tophat filter with a sharp edge, instead leaking outside the bounds of the tophat. Research is suggesting methods to deal with this, such as lessening the width of the tophat so that the leakage is within the filter width originally intended (Pascua et al. 2024). There are also parameters that can be applied to the use of these sequences, which we have not experimented with, that may mitigate the leakage.

The delay spectra through various averaging stages are surprisingly similar, indicating a peak at 0 ns delay, a drop-off outside the horizon, and a fairly flat noise floor beyond the horizon. What does change is the level of the noise floor, being reduced by the use of fringe-rate filtering. The noise, as expected, also reduces as more delay spectra are averaged at the various stages.

The 2D power spectra are also similar for both filters. To compare 2D power spectra from unfiltered and filtered visibilities, we ratio them, showing that the power outside the wedge drops more compared to the power inside the wedge, when filtering is used. The ratio also shows that outside the wedge, the power can sometimes increase compared to the filtered data, this is not everywhere, but at isolated points within the power spectrum.

In looking for the effects of the filters, we found that there can be a sharp drop in power at the horizon, more than anywhere else in delay space. This is an effect that we expect, particularly from the notch filter, which removes fringe-rates corresponding to emission from the horizon.

The similarities between the power spectra generated by the filters may be partly (or wholly) due to coherent averaging of visibilities after fringe-rate filtering. If coherent averaging is not used, then, at least comparing ratios of power spectra, the power outside the wedge differs depending on the filter. It therefore shows that the interaction of the filters with coherent averaging, and indeed any subsequent processing steps, should studied further.

In all the data products analysed, from visibilities to power spectra, it is clear that there is significant signal loss, simply due to removal of fringe-rates. Signal loss is not the subject of this paper and will be reported separately (Pascua et al. 2024); developing analytic methods to generate expected signal loss that can be combined with the analysis of errors on 2D power spectra (Tan et al. 2021), for example.

ACKNOWLEDGEMENTS

This result is part of a project that has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement no. 948764; PB and MJW). HG and PB acknowledge support from Science and Technology Facilities Council (STFC) Grant ST/T000341/1.

DATA AVAILABILITY

The total size of data generated for this paper is approximately 24 TB, if one includes all intermediate data products and the full output of the pipeline for the all the alternative versions – with coherent averaging or without coherent averaging, different eigenvalue cutoffs, and different frequency ranges. While it is not practical to release these data in full online, we can can supply components of it on request. The software used to generate the simulated data is public, and the configuration files and analysis scripts are also available on request.

References

APPENDIX A: EXAMPLES OF DPSS USED IN THE MAINLOBE FILTER

We will demonstrate filtering of the visibilities that are shown in the top row of Fig. 5 except that they have not been coherently averaged.

Given a tophat filter with a width and a centre in fringe-rate space, a set of DPSS are generated that have power concentrated within the tophat filter bounds. The mechanism for generating a set involves specifying parameters such as length, resolution, and filter requirements, then selecting DPSS that most concentrate power within the filter bounds. The selection method is based on eigenvalue cutoff, which is described in Appendix  B. Eighteen DPSS are selected to use for fitting the example visibilities; the five DPSS with greatest amplitude after fitting are shown in the left panel of Fig. A1. These are real valued, with power in fringe-rate space that is centred at fringe-rate 0 mHz, but before fitting are shifted to the centre of the tophat filter by converting to complex values and assigning appropriate phases. The right panel of Fig. A1 shows the amplitudes of the fringe-rates for the five DPSS after rephasing (coloured, corresponding to those in the left panel), and the black line shows the amplitude of their (complex) sum. The vertical grey lines indicate the bounds of the tophat filter.

Figure A1.

Examples of DPSS used for mainlobe filtering. Five DPSS are shown out of a set of 18 that are used for the complete mainlobe filter. The DPSS are shown on the left, and their corresponding fringe-rates on the right. The black line shows the amplitude of the five summed DPSS in fringe-rate space.

Open in new tabDownload slide

Fig. A2 places these five DPSS within the context of the mainlobe filtering of the chosen visibilities over time, shown in grey in the left panel. The filtered visibilities using all 18 DPSS are shown in red and the filtered visibilities using only the five DPSS are shown in blue. The right panel shows the fringe-rates of the visibilities in the left panel.

Figure A2.

The example visibility data and its fringe-rates, before and after mainlobe filtering. In the left panel, the input visibilities are shown in grey, the mainlobe-filtered visibilities are shown in red, and the filtered visibilities using only five DPSS are shown in blue, as described in the text. The right panel shows the corresponding fringe-rates.

Open in new tabDownload slide

Figure A3.

The results of an experiment to determine the consequences of different eigenvalue cutoffs when selecting DPSS for a bandpass (mainlobe) filter. The left panel shows the fringe-rates after filtering the visibilities for baseline (3, 6) at 100 MHz with different cutoffs as in the legend. The right panel shows the difference between the filtered and input fringe-rates, and the leakage outside the filter bounds, by cutoff. How these are calculated is defined in the text.

Open in new tabDownload slide

APPENDIX B: DPSS EIGENVALUE CUTOFF

DPSS arise from a system of equations whose solutions have a Fourier spectrum concentrated within certain bounds (Gruenbacher & Hummels 1994). The system of equations can be expressed as a matrix equation where the DPSS are eigenvectors. Since they are eigenvectors, they have an associated eigenvalue, with larger eigenvalues indicating a more concentrated spectrum within the filter bounds. By selecting DPSS with eigenvalues above a certain value, the ‘eigenvalue cutoff’, a limited set of DPSS can be obtained and used to bandpass filter a signal.

One then has to select an eigenvalue cutoff, and there is a trade-off between how well the spectrum is concentrated, and how much spectral leakage there is outside the bounds. Fig. A3 demonstrates this trade-off. We take a time sequence of simulated visibilities for baseline (3, 6) at 100 MHz, and filter it using DPSS, with five different eigenvalue cutoffs from |$10^{-2}$| to |$10^{-12}$|⁠. The left panel shows the fringe-rates after filtering, with the legend indicating the eigenvalue cutoff and the number of DPSS that were selected by that cutoff. The right panel shows how well the filter has performed, using two measures: (1) the RMS difference between the fringe-rates of the input visibilities and the fringe-rates of the filtered visibilities, within the filter bounds, and (2) the RMS difference between the fringe-rates of the filtered visibilities and 0, outside the bounds. In other words, (1) indicates how well the filtered spectrum matches the input spectrum within the bounds, and (2) indicates how much spectral leakage there is. In the legend, ‘Diff’ denotes measure 1, and 1‘Leakage’ denotes measure 2. We want both measures to be small. The filter bounds are indicated by the vertical grey lines.

As the eigenvalue cutoff value is reduced, more DPSS modes are selected, and the fit of the filtered spectrum within the filter bound improves (‘Diff’ goes to 0), but the leakage gets worse (‘Leakage’ increases). However, it is clear that some leakage will always be present, whereas some ‘Diff’ will not always be present, since it is possible to reduce ‘Diff’ to near 0. In the case of the mainlobe filter, this means that the fringe-rates we want to retain are preserved very well, as is seen in the second row, left panel, of Fig. 5. The notch filter is a different matter, where the match is not as good, but the notch filter suffers from issues of resolution of fringe-rates where the notch is so small, and this needs to be further investigated.

In this work we choose to minimize ‘Diff’, and use an eigenvalue cutoff of |$10^{-12}$|⁠, accepting that there will always be leakage outside the filter bounds. The actual form of this leakage can be examined in further work, in particular how the RMS difference is spread across the fringe-rates outside the bounds, and the impact of that. Within the filter bounds, at least for the mainlobe filter, we retain a good match to the signal. In the next section, we will show the effect of different cutoffs on power spectra, and that the effect is limited.

Note that we have examined the consequences of different cutoffs for only one baseline. The same analysis could be made for all baselines, with cutoffs and the number of DPSS modes varying by baseline. The effect of the cutoff, and the filters in general, should also be examined for real data, but the cutoff of |$10^{-12}$| is a staring point that allows inspection of the filter effect for the mainlobe filter within the filter bounds, and from that the cutoff can be adjusted to observe its impact.

B1 The cutoff’s effect on wedge power spectra

To compare power spectra filtered with different eigenvalue cutoffs, we use the ratio plots in Fig. B1. These ratio plots have a different meaning to the ratio plots elsewhere in this paper. Fig. 8 contains power spectra from visibilities passed through the mainlobe and notch filters using an eigenvalue cutoff of |$10^{-12}$|⁠. The same power spectra were regenerated using an eigenvalue cutoff of |$10^{-6}$|⁠. Then, the latter were divided by the former.

Figure B1.

Each of these power spectra is the result of dividing power spectra generated from visibilities filtered using different eigenvalue cutoffs. See text.

Open in new tabDownload slide

The plots show that a cutoff of |$10^{-6}$| can increase the power by a factor of |$\sim 10$| in some locations of the power spectra. It can also decrease the power by a factor of |$\sim 10^{-3}$| in other locations, although this happens less often. Fig. B2 shows that most of the locations where there are substantial increases/decreases lie outside the wedge; inside the wedge there is little change in power. The changes are also positive and negative, rather than an obvious increase or decrease over the entire spectrum.

Figure B2.

Vertical cuts through the power spectra in Fig. B1.

Open in new tabDownload slide

While this does not indicate which cutoff value is ‘best’, for example for detecting the EoR, it does show that a |$10^{-6}$| change in magnitude of the cutoff does not radically alter a power spectrum. Given that there are so many other data manipulations involved – coherent and incoherent time averaging, redundant baseline averaging – it is possible that the cutoff value is of minor concern in comparison.

Author notes