* [23] Hoheisel T, Kanzow C, Mordukhovich BS, Phan H (2012) Generalized Newton’s method for nonsmooth equations based on graphical derivatives. Nonlinear Anal. 75:1324–1340.
  * [24] Izmailov AF, Solodov MV (2014) Newton-Type Methods for Optimization and Variational Problems (Springer, New York).
  * [25] Jiang H, Qi L, Chen X, Sun D (1996) Semismoothness and superlinear convergence in nonsmooth optimization and nonsmooth equations. De Pillo G, Giannessi F, eds. Nonlinear Optimization and Applications (Springer, New York), 197–212.
  * [26] Josephy NH (1979) Newton’s method for generalized equations. Technical Summary Report No. 1965. Mathematics Research Center, University of Wisconsin (Madison, WI).
  * [27] Kenderov P (1975) Semi-continuity of set-valued monotone mappings. Fundamenta Mathematicae 88:61–69.
  * [28] Klatte D, Kummer B (2002) Nonsmooth Equations in Optimization. Regularity, Calculus, Methods and Applications (Kluwer Academic Publishers, Dordrecht, The Netherlands).
  * [29] Kummer B (1988) Newton’s method for non-differentiable functions. Guddat J et al., eds. Advances in Mathematical Optimization (Akademi-Verlag, Berlin), 114–124.
  * [30] Meng F, Sun D, Zhao Z (2005) Semismoothness of solutions to generalized equations and the Moreau-Yosida regularization. Math. Program. 104:561–581.
  * [31] Mifflin R (1977) Semismooth and semiconvex functions in constrained optimization. SIAM J. Control Optim. 15:957–972.
  * [32] Mohammadi A, Mordukhovich BS, Sarabi ME (2020) Variational analysis of composite models with applications to continuous optimization. Math. Oper. Res. (to appear). arXiv:1905.08837.
  * [33] Mohammadi A, Mordukhovich BS, Sarabi ME (2020) Parabolic regularity in geometric variational analysys. Trans. Amer. Math. Soc. (to appear). arxiv:1909.00241.
  * [34] Mohammadi A, Sarabi ME (2020) Twice epi-differentiability of extended-real-valued functions with applications in composite optimization. SIAM J. Optim. (to appear). arXiv:1911.05236.
  * [35] Mordukhovich BS (1992) Sensitivity analysis in nonsmooth optimization. Field DA, Komkov V, eds. Theoretical Aspects of Industrial Design, SIAM Proc. Appl. Math. 58:32–46 (Philadelphia, PA).
  * [36] Mordukhovich BS (1993) Complete characterizations of openness, metric regularity, and Lipschitzian properties of multifunctions. Trans. Amer. Math. Soc. 340:1–35.
  * [37] Mordukhovich BS (2006) Variational Analysis and Generalized Differentiation, I: Basic Theory, II: Applications (Springer, Berlin).
  * [38] Mordukhovich BS (2018) Variational Analysis and Applications (Springer, Cham, Switzerland).
  * [39] Mordukhovich BS, Nghia TTA (2015) Second-order characterizations of tilt stability with applications to nonlinear programming. Math. Program. 149:83–104.
  * [40] Mordukhovich BS, Nghia TTA (2016) Local monotonicity and full stability of parametric variational systems. SIAM J. Optim. 26:1032–1059.
  * [41] Mordukhovich BS, Outrata JV (2001) On second-order subdifferentials and their applications, SIAM J. Optim. 12:139–169.
  * [42] Mordukhovich BS, Rockafellar RT (2012) Second-order subdifferential calculus with applications to tilt stability in optimization. SIAM J. Optim. 22:953–986.
  * [43] Mordukhovich BS, Sarabi ME (2020) Generalized Newton algorithms for tilt-stable minimizers in nonsmooth optimization. arXiv:2004.02345.
  * [44] Outrata JV, Sun D (2008) On the coderivative of the projection operator onto the second-order cone. Set-Valued Anal. 16:999–1014.
  * [45] Pang J-S (19900 Newton’s method for B-differentiable equations. Math. Oper. Res. 15:311–341.
  * [46] Poliquin RA, Rockafellar RT (1996) Prox-regular functions in variational analysis. Trans. Amer. Math. Soc. 348:1805–1838.
  * [47] Poliquin RA, Rockafellar RT (1998) Tilt stability of a local minimum. SIAM J. Optim. 8:287–299.
  * [48] Qi L, Sun J (1993) A nonsmooth version of Newton’s method. Math. Program. 58: 353–367.
  * [49] Robinson SM (1979) Generalized equations and their solutions, I: basic theory. Math. Program. Study 10:128–141.
  * [50] Rockafellar RT, Wets RJ-B (1998) Variational Analysis (Springer, Berlin).
  * [51] Rockafellar RT, Zagrodny D (1997) A derivative-coderivative inclusion in second-order nonsmooth analysis. Set-Valued Anal. 5:1–17.
  * [52] Tibshirani R (1996) Regression shrinkage and selection via the Lasso. J. R. Stat. Soc. 58:267–288.
  * [53] Ulbrich M (2011) Semismooth Newton Methods for Variational Inequalities and Constrained Optimization Problems in Function Spaces (SIAM, Philadelphia, PA).
  * [54] Wang X (2004) Subdifferentiability of real functions. Real Anal. Exchange 30:137–172.
  * [55] Yao J-C, Yen ND (2009) Coderivative calculation related to a parametric affine variational inequality. Part 1: Basic calculation. Acta Math. Vietnam. 34:157–172.


#  A Comparative Analysis To Deal With Missing Spectral Information Caused By
RFI In Cosmological Hi 21CM Observations

Arnab Chakraborty Discipline of Astronomy, Astrophysics and Space Engineering,
Indian Institute of Technology Indore, Indore 453552, India  Department of
Physics, McGill University, 3600 rue University, Montreal, QC H3A 2T8, Canada
McGill Space Institute, McGill University, 3550 rue University, Montreal, QC
H3A 2A7, Canada  Abhirup Datta  Discipline of Astronomy, Astrophysics and
Space Engineering, Indian Institute of Technology Indore, Indore 453552, India
Aishrila Mazumder  Discipline of Astronomy, Astrophysics and Space
Engineering, Indian Institute of Technology Indore, Indore 453552, India

###### Abstract

We investigate the effect of radio-frequency interference (RFI) excision in
estimating the cosmological Hi 21 cm power spectrum. Flagging of RFI-
contaminated channels results in a non-uniform sampling of the instrumental
bandpass response. Hence, the Fourier transformation (FT) of visibilities from
frequency to delay domain contaminates the higher foreground-free delay modes,
and separating the spectrally fluctuating Hi signal from spectrally smooth
foregrounds becomes challenging. We have done a comparative analysis between
two algorithms, one-dimensional CLEAN and Least Square Spectral Analysis
(LSSA), which have been used widely to solve this issue in the literature. We
test these algorithms using the simulated SKA-1 low observations in the
presence of different RFI flagging scenarios. We find that in the presence of
random flagging of data, both algorithms perform well and can mitigate the
foreground leakage issue. But, CLEAN fails to restrict the foreground leakage
in the presence of periodic and periodic plus broadband RFI flagging and gives
an extra bias to the estimated power spectrum. However, LSSA can restrict the
foreground leakage for these RFI flagging scenarios and gives an unbiased
estimate of the Hi 21 cm power spectrum. We have also applied these algorithms
to the upgraded GMRT observation and found that both CLEAN and LSSA give
consistent results in the presence of realistic random flagging scenarios for
this observed data set. This comparative analysis demonstrates the
effectiveness and robustness of these two algorithms in estimating the Hi 21
cm power spectrum from the data set affected by different RFI scenarios.

methods: data analysis-methods: interferometric-radio: cosmology: dark ages,
reionization, first stars - galaxies: evolution – large-scale structure of
universe

††journal: APJ

##  1 Introduction

The redshifted 21 cm “spin-flip” transition (Field 1958) of the neutral
hydrogen (Hi ) is considered as a promising tool to understand the global
phase transition of the Universe from neutral phase to ionized phase (Epoch of
Reionization; EoR) over a redshift range ($6<\textit{z}<15$) (Field 1958;
Madau et al. 1997; Furlanetto et al. 2006; Pritchard & Loeb 2012; Loeb &
Furlanetto 2013; Dayal & Ferrara 2018). One of the promising way to study EoR
is measuring Hi 21 cm power spectrum along with tomographic imaging of the
intergalactic medium (IGM) using low-frequency large interferometric arrays.
(Bharadwaj & Sethi 2001; Bharadwaj & Ali 2005; Morales & Hewitt 2004;
Zaldarriaga et al. 2004). Several upcoming and ongoing telescopes have been
designed to study the 21 cm power spectrum, such as Donald C.Backer Precision
Array to Probe the Epoch of Reionization (PAPER, Parsons et al. 2010), the Low
Frequency Array(LOFAR, van Haarlem et al. 2013), the Murchison Wide -field
Array (MWA, Li et al. 2018), the Hydrogen Epoch of Reionization Array (HERA,
DeBoer et al. 2017) and the Square Kilometer Array (SKA1 LOW, Koopmans et al.
2015). Also a statistical detection of post-reionization (post-EoR) Hi 21 cm
signal holds the key to understand the large scale structure formation, source
clustering and will put constrain to the cosmological parameters (Bharadwaj &
Pandey 2003; Bharadwaj & Ali 2005; Bull et al. 2015; Wyithe & Loeb 2008).

However, the redshifted Hi signal from the EoR and post-EoR epoch is hidden
behind several orders of magnitude bright radio emissions from different
Galactic and extragalactic foregrounds (Bharadwaj & Ali 2005; Zaldarriaga et
al. 2004; Jelić et al. 2008; Di Matteo et al. 2004). Several methods are
available in the literature to estimate the Hi power spectrum amidst this
bright sea of foregrounds. One can model the foregrounds accurately and then
subtract them from the data set (Yatawatta et al. 2013; Procopio et al. 2017).
Previous studies have shown that foreground is spectrally smooth, whereas the
redshifted Hi 21 cm signal has spectral structure. Due to this spectral
smoothness and chromatic primary beam response, foregrounds are localized
within a ‘wedge’ shape region in two dimensional (2D) cylindrically averaged
power spectrum. A traditional way is to avoid these foregrounds dominated
Fourier modes inside the wedge and search for cosmological Hi signal in the
so-called ‘EoR-window’. This technique is known as foreground avoidance (Datta
et al. 2010; Parsons et al. 2012; Vedantham et al. 2012; Morales et al. 2012,
2019). Along with that, there are techniques to separate spectrally smooth
foregrounds from rapidly fluctuating Hi signal without any prior modeling of
foregrounds (Harker et al. 2009; Chapman et al. 2013; Mertens et al. 2018).
However, all methods rely on the fact that the spectral and spatial properties
of foregrounds are very well known. Suppose there is some oscillating feature
across the frequency band, introduced by instrumental effects or by any post-
processing methods. In that case, it will pose a hindrance in isolating the
redshifted Hi signal from the foregrounds.

One possible cause that can introduce rapidly fluctuating components is the
excision of radio frequency interference (RFI). In general, RFI is localized
in time or frequency domain. It is usual practice to measure visibilities with
high time and frequency resolution, which helps identify RFI across time and
frequency. Then those RFI contaminated samples get flagged in post-processing
of the data analysis. This results in a non-uniform sampling of the
instrumental bandpass response, which otherwise is smooth in general. Now,
each visibility is the Fourier transform of the product of the sky signal with
instrumental bandpass response. Hence, RFI flagging introduces irregularities
in visibilities. To estimate the power spectrum of the cosmological signal,
one needs to Fourier transform the measured visibilities from frequency
($\nu$) to delay ($\eta$) domain. During this Fourier transformation, the
irregularities or non-uniform sampling of the bandpass will create rapid
fluctuation across the delay axis, which is proportional to the Fourier
$k_{\parallel}$ modes. This can be thought of as a convolution of visibilities
in the delay domain with the delay space PSF, which is the Fourier transform
of the non-uniform instrumental bandpass. The delay space PSF will have large,
rapidly oscillating sidelobes, which will create a rapidly fluctuating
component when convolved with visibilities in delay space. Hence, this RFI
excision causes ripple across $k_{\parallel}$ axis and the ‘EoR-window’ will
become contaminated by leakage of foregrounds.

There are several ways existing in literature to deal with this issue. PAPER
analyses have used a wide-band iterative deconvolution algorithm (WIDA) to
remove the foreground components in delay space and then estimate the delay
space power spectrum (Parsons et al. 2014; Jacobs et al. 2015; Kerrigan et al.
2018). Kolopanis et al. (2019) used Blackman-Harris tapering window function
during Fourier transforming to delay domain to reduce foreground leakage
during FFT. Trott et al. (2016) used least-square spectral analysis (LSSA) to
perform the line-of-sight transform from frequency to delay space instead of
Fourier transform to mitigate the issues related to non-uniform bandpass shape
for MWA observation. Patil et al. (2017),Gehlot et al. (2019), Mertens et al.
(2020) have also used LSSA to the LOFAR data set to put an upper limit on the
EoR Hi 21 cm power spectrum. Offringa et al. (2019) used Gaussian weighted
interpolation of missing visibility samples due to RFI flagging before any
averaging of data. However, they used the same LSSA during the transformation
from frequency to delay space.

Table 1: Details of Simulations Simulation parameters
---
Telescope | SKA-1 low
Min baseline | 42 m
Max baseline | 2.8 km
Frequency | 138-146 MHz
Number of channels | 126
Central redshift ($z$) | 9.0
| Cosmological signal:
| simulated Hi signal cube
Skymodel | +
| Foreground: extragalactic point
| sources and diffuse foregrounds
| +
| Gaussian noise
|
---|---
Figure 1: Left: the array layout of SKA-1 low with 297 elements. Right: the
snapshot $uv$-coverage of the simulated data for SKA-1 low.

Deshpande et al. (1996) have used 1D CLEAN (Högbom 1974) algorithm to estimate
the power spectrum of pulsar timing residual from non-uniformly sampled timing
measurement. They showed that the CLEAN algorithm could help reconstruct the
power spectra from non-uniformly sampled time sequences and dramatically
improved the dynamic range of the estimated power spectra. Deshpande et al.
(2000) have used the same CLEAN algorithm to estimate the power spectrum of
the optical depth of the cold Hi gas in the Galaxy toward Cas-A and Cygnus-A
from the incomplete sampling of the optical depth over the extent of the
source. In this work, we implement the same 1D CLEAN algorithm in the delay
space to deconvolve the effects of sidelobes of the delay space PSF (point
spread function). Then we reconstruct the delay spectra of visibilities and
estimate the 2D and 3D power spectrum of the cosmological Hi signal.

LSSA has been used in digital signal processing as well as in astronomy to
deal with the non-uniformly sampled data set Barning (1963); Lomb (1976);
Stoica et al. (2009). The spectrum of non-uniformly sampled data is difficult
to interpret as the true peak in the spectrum gives several other peaks
(aliases) of various heights distributed throughout the delay spectrum. In
other words, the alias structure of the major peak creates confusion in
estimating the true underlying spectrum Lomb (1976). This is similar to
foreground contamination of large $k_{\parallel}$ modes beyond horizon limit
as described before. LSSA is a commonly used method, where sinusoidal of
different frequencies are fitted to the non-uniform data via the traditional
least-square method, resulting in a reasonably good approximation of the
spectrum of non-uniform data. The Least Square spectrum provides the best
measure of the power contributed by different frequencies to the overall
variance of the data and can be regarded as the extension of the Fourier
methods to non-uniform data (Lomb 1976). In the limit of equal spacing data
set, this method reduces to Fourier power spectrum.

In this work, we show the application of these widely used CLEAN and LSSA
algorithms on the simulated data (SKA-1 low) for different flagging scenarios
and compare our results with that of a uniformly sampled data set (i.e., no
flagging). We also apply these to the actual observation of the upgraded Giant
Metrewave Radio Telescope (uGMRT) data set in the presence of realistic
flagging. This comparative analysis shows which algorithm gives a more robust
and improved estimation of the cosmological Hi power spectrum in the presence
of RFI flagging, without adding extra bias due to the missing samples issue.

This paper is structured as follows: The simulation of data for SKA-1 low
array configuration with observation setup and sky model has been discussed in
Sec. 2. In Sec. 3, we discuss the basic formalism to estimate the power
spectrum from the simulated observation, including missing samples due to
flagging. The application of the CLEAN and LSSA for different flagging
scenarios is shown in Sec.4. The application to the uGMRT observation is
presented in Sec. 5. Finally, we conclude in Sec. 6.

|  |
---|---|---
Figure 2:  The simulated image slices of the EoR signal (left) and foreground
(middle) emissions and noise (right) at redshift $z=9$. The images cover sky
areas of $5^{\circ}\times 5^{\circ}$.

##  2 Simulations

To investigate the effect of missing frequency channels on the estimation of
the power spectrum and to demonstrate the performance of the CLEAN and LSSA,
we use the data from simulated observations using the SKA-1 low. We follow the
procedure given in Li et al. (2019) for the simulation of EoR and foreground
signal simulation and briefly mention the procedure here. Li et al. (2019)
extrapolate the low-frequency observation to make foreground sky maps that
exhibit realistic structure of the sky. Foreground sources include the
contribution from diffuse Galactic synchrotron emission, free-emission,
extragalactic point sources, and radio haloes (Wang et al. 2010; Li et al.
2019). The Galactic synchrotron emission is simulated by extrapolating the
reprocessed 408 MHz Haslam map (Remazeilles et al. 2015) using a power-law
spectrum, the index of which is taken from Giardino et al. (2002). The
Galactic free-free radiation map has been constructed from H$\alpha$ survey
map (Finkbeiner 2003), where the close relation between H$\alpha$ and free-
free emission has been employed. Galactic diffuse emission is simulated at a
central position of (R.A., Dec.) = $(0^{\circ},-27^{\circ})$, which is at high
Galactic latitude and suitable candidate field for EoR observation (Beardsley
et al. 2016). The extragalactic point sources are simulated by using the SKADS
catalog (Wilman et al. 2008). The radio halos are simulated by generating a
sample of galaxy clusters with the Press-Schechter formalism (Press &
Schechter 1974) and then using the scaling relation (cluster mass-X-ray
temperature and X-ray temperature-radio power relations) to estimate their
radio emission. For more details regarding foreground simulation, please see
(Wang et al. 2010; Li et al. 2019).

For the simulation of the EoR signal, we have used the Evolution of 21cm
Structure 111http://homepage.sns.it/mesinger/EOS.html data set (Mesinger et
al. 2016) and extract the image slices at corresponding frequencies from the
light-cone cube of the recommended ‘faint galaxy’ models. These EoR sky maps
are tiled and re-sampled to an FoV of $5^{\circ}\times 5^{\circ}$ with a pixel
resolution of $20\arcsec$. Note that this re-scaling is done to match the same
sky coverage and resolution of the foreground sky maps.

To incorporate the instrumental effects into these simulated sky maps, we use
the SKA-1 low 222https://astronomers.skatelescope.org/DOCUMENTS array
configuration and the OSKAR 333https://github.com/OxfordSKA/OSKAR package
(Mort et al. 2017) to simulate the instrumental observation. The array
configuration consists of 512 stations, each of which consists of 256 antennas
randomly distributed within a circle of 35m in diameter. The stations are
distributed in two parts: (1) 224 stations are distributed in a core region of
diameter 1 km and the remaining distributed in three spiral arms extended up-
to 35 km. The array layout of the inner core region is shown in Fig. 1. OSKAR
uses the Radio Interferometer Measurement Equation (RIME) of Hamaker et al.
(1996) to simulate the full stokes visibility data set, however, here we only
simulate the RR component of visibilities. We provide the array configuration
of SKA-1 low to the telescope model and simulate the visibilities for each sky
map at an (R.A., Dec.) = $(0^{\circ},-27^{\circ})$. We also include the
uncorrelated Gaussian noise to our simulation corresponding to 1000 hours of
observation with SKA-1 low. The system temperature ($T_{\rm sys}$) has
contribution from both sky and the instrument receiver, i.e, $T_{\rm
sys}=T_{\rm sky}+T_{\rm rcvr}$. The assumed sky temperature is $T_{\rm
sky}=60\rm K\Big{(}\frac{300\mathrm{MHz}}{\nu}\Big{)}^{2.25}$ (Thompson et al.
2017) and the receiver temperature is $T_{\rm rcvr}=40\rm K$ (Bacon et al.
2020). We choose a 8 MHz bandwidth around the central frequency 142 MHz
($z\sim$ 9) with 64 kHz channel resolution (126 channels). The $uv$-plane
corresponding to the simulated data is being shown in Fig. 1. We then use
WSCLEAN 444https://gitlab.com/aroffringa/wsclean (Offringa et al. 2014) to
make dirty image cube of size $300\times 300\times 126$ (pixels $\times$
pixels $\times$ $N_{\rm chans}$) for the EoR signal and foreground emission.
We then convert these cubes from Jy/beam to Kelvin (K) units (Patil et al.
2017) and then merge them to make a combined dirty image cube,
$I^{D}(\theta_{x},\theta_{y},\nu)$, where $\theta_{x},\theta_{y}$ are the sky
positions and $\nu$ is the frequency. We use this image cube to estimate the
Hi 21 cm power spectrum and briefly discussed below. The image slice
corresponding to the EoR Hi 21 cm signal (left panel) and foreground (middle
panel) and instrumental Gaussian noise (right panel) at $z=9$ are shown in
Fig. 2.

##  3 Estimation of Power Spectrum in presence of missing Frequency Channel

In this section, we will briefly describe the formalism used in this work to
estimate the power spectra from the dirty image cube,
$I^{D}(\theta_{x},\theta_{y},\nu)$, which also suffers from missing frequency
channels, due to RFI flagging. We follow the methodology as described in
Morales & Hewitt (2004), Datta et al. (2010), Patil et al. (2017).

We first perform the 2D spatial Fourier transform of the image cube to make
the visibility cube

$V({\mathbb{U}},\nu)=\int
I^{D}(\bm{\theta},\nu)e^{-i2\pi{\mathbb{U}}.{\bm{\theta}}}\rm
d^{2}{\bm{\theta}},$ (1)

where, ${\mathbb{U}}=(u,v)$ and $\bm{\theta}=(\theta_{x},\theta_{y})$ are the
vectors.

We introduce the RFI flagging by multiplying the visibility cube with the
frequency dependent sample weights ($S(\nu)$), given by:

$\displaystyle S(\nu)$ $\displaystyle=$ $\displaystyle 0,\forall\rm
flagged~{}frequency~{}channel$ (2) $\displaystyle=$ $\displaystyle
1,\forall\rm channels~{}without~{}flagging$

Here $S(\nu)$ essentially corresponds to the flagging mask due to RFI, which
has been applied to the visibility data cube. Then we transform the visibility
cube to the delay domain ($\eta$-domain) by Fourier transforming the data
along the frequency axis, given by (Morales & Hewitt 2004):

$V({\mathbb{U}},\eta)=\int
V({\mathbb{U}},\nu)B(\nu)S(\nu)e^{-i2\pi\nu\eta}d\nu,$ (3)

where, $B(\nu)$ is the instrumental bandpass response, taken flat in this
case.

According to the Fourier convolution theorem, visibility
$V({\mathbb{U}},\eta)$ can be expressed as convolution of
$\mathscr{F}[V({\mathbb{U}},\nu]$ with the convolving kernel
$\mathscr{F}[B(\nu)S(\nu)]$ as:

$V({\mathbb{U}},\eta)=\mathscr{F}[V({\mathbb{U}},\nu)]*\mathscr{F}[B(\nu)S(\nu)],$
(4)

where $\mathscr{F}$ denotes the Fourier transform operator and $*$ denotes the
convolution.

We refer to the Fourier transform of $[B(\nu)S(\nu)]$ to the delay domain as
the delay space point spread function (PSF) and denote it as $D_{\rm psf}$.
Since the excision of RFI contaminated channels results in a non-uniform
sampling of the instrumental bandpass, which is encoded in $S(\nu)$, the
simple Fourier transform of visibilities (Eqn. 3) along the frequency axis
will result in spectral leakage. This causes leakage of smooth foregrounds
into the EoR window. In other words, the side lobes of the delay space PSF
($D_{\rm psf}$) due to RFI flagging will be large. Convolution of $D_{\rm
psf}$ with $\mathscr{F}[V({\mathbb{U},\nu}]$ will result in leakage of
foregrounds beyond the horizon limit into the EoR window and will show an
oscillating pattern above the horizon delay.

To understand the problem, we show the delay spectrum of a single baseline of
length $162\lambda$ (arbitrarily chosen) from our simulated data cube (see
Fig. 3). We first Fourier transform the data cube to the delay domain without
using any flagging, i.e, $S(\nu)=1,\forall\rm channels$. The resulting delay
spectrum of the baseline is shown in Fig. 3 in gray. We refer to this as no
RFI case in the plot. The vertical dashed lines are the delay horizon limits
for this particular baseline. Then we randomly flagged 10% channels and
performed the delay transformation.

One can use a window function $W(\nu)$ during the Fourier transformation to
the delay domain to reduce the spectral leakage (Thyagarajan et al. 2013;
Kolopanis et al. 2019). Then Eqn. 3 and Eqn. 4 will become:

$V({\mathbb{U}},\eta)=\int
V({\mathbb{U}},\nu)B(\nu)S(\nu)W(\nu)e^{i2\pi\nu\eta}d\nu,$ (5)

and

$V({\mathbb{U}},\eta)=\mathscr{F}[V({\mathbb{U},\nu}]*\mathscr{F}[B(\nu)S(\nu)W(\nu)],$
(6)

Here, $W(\nu)$ will act as a spectral taper during delay transformation and
will reduce the spectral leakage to some extent. We show in Fig. 3 the
resulting delay spectrum of the baseline with a Blackman-Harris window as a
spectral weighting function (in blue). As seen from the figure, there is still
significant spectral leakage beyond the horizon limit even after using the
window function. This is mainly due to the convolution of visibility with the
$\mathscr{F}[B(\nu)S(\nu)W(\nu)]$ , where the missing channels due to flagging
is encoded inside of $S(\nu)$.

We employ two previously widely used methods on data, CLEAN (Parsons et al.
2014) and LSSA (Patil et al. 2017; Trott et al. 2016), to mitigate this
spectral leakage issue due to missing samples (RFI flagging). Kolopanis et al.
(2019) used a Blackman-Harris tapering function to reduce the spectral leakage
in their reanalysis of PAPER-64 data. However, as we have seen above, a simple
spectral tapering is unable to mitigate foreground spillover due to missing
samples in the presence of only 10% RFI flagging. We briefly discuss the other
two methods, CLEAN and LSSA, below and compare our results here.

Figure 3:  Delay spectrum of an arbitrarily chosen baseline of length
162$\lambda$. The top blue curve is the delay spectrum with 10% channels are
randomly flagged. The green dashed curve is the delay spectrum after using 1D
CLEAN algorithm across delay axis. The red dashed curve is the delay spectrum
of the baseline after using Least Square Spectral Analysis (LSSA) algorithm to
this flagged data set. The gray curve is the delay spectrum of the baseline,
when there is no flagging (no RFI). The vertical dashed black lines are the
horizon limit, which is $1.14\mu$s, for this particular baseline.

CLEAN: To restrict the spectral leakage beyond the horizon limit, we
deconvolve the $D_{\rm psf}$ using a one-dimensional Hogbom CLEAN algorithm
across the delay axis (Högbom 1974; Parsons & Backer 2009; Parsons et al.
2014; Kerrigan et al. 2018). The CLEAN algorithm iteratively searches for the
maximum peak in the delay spectrum of visibility and then subtracts a fraction
(0.1; gain term) of the delay space PSF from that until a relative threshold
in the residual is reached. This threshold is determined by the ratio of the
initial and final amplitude in the spectrum and provided as the tolerance
parameter in the CLEAN algorithm (Kerrigan et al. 2018; Ewall-Wice et al.
2021). In other words, the tolerance parameter sets the degree of foreground
subtraction during CLEANing and can be set as low as we want to CLEAN deeper
(Ewall-Wice et al. 2021). The choice of this tolerance parameter for various
flagging scenarios is discussed in Sec. 4. We first zero-pad the visibilities
by the same number of channels as that of the simulation and apply a Tukey
window (with $\alpha=0.2$) as a tapering function before performing the CLEAN
(Kern et al. 2019; Ewall-Wice et al. 2021). After reaching the threshold in
the residual, we add the residual to the CLEAN model components. This delay
space CLEANing effectively deconvolves the delay space PSF and reduces the
side lobes of the delay spectrum due to non-uniform sampling of the bandpass
and gives the final delay spectra. In Fig. 3, we show the delay spectrum after
using 1D CLEAN across the delay axis in green dashed line. It is clear from
the figure that, with this simple deconvolution technique, the leakage of
foreground beyond the horizon limit due to missing channels can be restricted
to minimal.

Table 2: Details of different cases studied in this work Different Cases
Studied in this Work
---
Case | Flagging Scenario
I | No flagging or No RFI
II | 20% random channels flagged
| 20% random channels flagged
III-A | $+$
| periodic flagging, with 200 KHz flag
| in every 1.28 MHz
| 20% random channels flagged
III-B | $+$
| periodic flagging, with 300 KHz flag
| in every 0.64 MHz
| 20% random channels flagged
| $+$
| periodic flagging, with 200 KHz flag
IV | in every 1.28 MHz
| $+$
| 2 MHz wide frequency chunk from the
| center of the band has been flagged

Least Square Spectral Analysis (LSSA): Another widely used approach to do the
line-of-sight Fourier transform of the data from frequency space to delay
space in the presence of missing samples is Least Square Spectral Analysis
(LSSA) (Patil et al. 2017; Trott et al. 2016). LSSA estimates the Fourier
spectrum by fitting sinusoids to the data samples using the traditional least
square method. A discrete data set is approximated by a weighted sum of
sinusoids of progressively determined frequencies, using a standard linear
regression or least-squares fit. If one knows the covariance of the data, then
that can be used to weight the data vector to get the approximated functional
form of the discrete data samples (Stoica et al. 2009; Lomb 1976). Trott et
al. (2016) first used this method in their Hi 21 cm power spectrum estimation
pipeline and applied it to the MWA data set. Later Patil et al. (2017); Gehlot
et al. (2019); Mertens et al. (2020) also used this methodology in their
analysis of the LOFAR data set to put upper limit on Hi 21 cm power spectrum.
We have implemented this method here as well to mitigate the spectral leakage
in the delay spectrum of visibilities due to missing frequency samples. A
smooth model, which is a Fourier series up to a specified order, is being
fitted to the complex visibility data cube with flags, using a linear least-
square solver (scipy.optimize.lsq_linear). We use 300 Fourier modes to fit the
visibility spectrum in this analysis. The delay spectrum after application of
LSSA for a single baseline is shown in red dashed line in Fig. 3. We find that
similar to CLEAN, the LSSA approach is also able to reduce the spectral
leakage in the presence of random RFI, and we can go down to the signal level
(without RFI) in the delay spectrum beyond the horizon limit.

|
---|---
Figure 4:  Left: The cylindrically averaged 2D power spectrum in
$\mathbb{k}_{\perp}$ vs $k_{\parallel}$ space for the Case I, i.e, when there
is no RFI (no flagging). The black dashed line denotes the horizon limit plus
a buffer due to Blackmann-Harris window. Right: The spherically averaged 3D
power spectrum using the modes above the horizon limit with $1\sigma$ error
bars (in green). We also show, in black dashed line, the input Hi 21 cm power
spectrum of 21cmFAST.

Power spectrum estimation: After getting the delay visibility cube
($V(\mathbb{U},\eta)$), we proceed to estimate the cylindrically averaged 2D
power spectrum using the formula (Morales & Hewitt 2004) :

$P(k_{\perp},k_{\parallel})=\Big{(}\frac{\lambda^{2}}{2k_{B}}\Big{)}^{2}\Big{(}\frac{X^{2}Y}{\Omega
B}\Big{)}|V(\mathbb{U},\eta)|^{2},$ (7)

where $\lambda$ is the wavelength corresponding to the band-center, $K_{\rm
B}$ is the Boltzmann constant, $\Omega$ is sky-integral of the squared antenna
primary beam response, B is the bandwidth and, X and Y are the conversion
factors from angle and frequency to transverse co-moving distance (D(z)) and
the co-moving depth along the line of sight, respectively (Morales & Hewitt
2004). The power spectrum $P(k_{\perp},k_{\parallel})$ is in units of
$\mathrm{K}^{2}(\mathrm{Mpc}/h)^{3}$. The wave numbers $k_{\perp}$ and
$k_{\parallel}$ are related to the baseline vector ($\mathbb{U}$) and delay
($\eta$) as :

$k_{\perp}=\frac{2\pi|\mathbb{U}|}{D(z)},$ (8)

$k_{\parallel}=\frac{2\pi\eta\nu_{21}H_{0}E(z)}{c(1+z)^{2}},$ (9)

where, $\nu_{21}$ is the rest-frame frequency of the 21 cm spin flip
transition of Hi , $z$ is the redshift to the observed frequency, $H_{0}$ is
the Hubble parameter and
$E(z)\equiv[\Omega_{\mathrm{M}}(1+\textit{z})^{3}+\Omega_{\Lambda}]^{1/2}$.
$\Omega_{\mathrm{M}}$ and $\Omega_{\Lambda}$ are matter and dark energy
densities, respectively (Hogg 1999). In this work, we use the best fitted
cosmological parameters of the Planck 2018 analysis (Planck Collaboration et
al. 2018) .

The 3D power spectrum can be calculated by spherically averaging
$P(\mathbb{k}_{\perp},k_{\parallel})$ in independent $k$-bins and
dimensionless 3D power spectrum can be written as (Datta et al. 2010):

$\Delta^{2}(k)=\frac{K^{3}}{2\pi^{2}}<P(\mathbb{k})>_{k},$ (10)

where $k=\sqrt{k_{\perp}^{2}+k_{\parallel}^{2}}$. Throughout this work we will
follow this formalism to estimate the 2D and 3D power spectrum. We estimate
the uncertainty in the estimated power spectrum by calculating the variance of
the measured power spectrum and dividing it by the number of independent modes
averaged together in each $k$-bins (Tegmark 1997). We quote this $1\sigma$
uncertainty as the error in the 3D power spectrum at each $k$-modes.

|
---|---
|
Figure 5:  The ratio of the cylindrically averaged 2D power spectrum between
different algorithms and no RFI scenarios, for Case II. Top left panel is the
ratio of the flagged to no RFI scenario. Top right and bottom left panels show
the ratio after applying the CLEAN and LSSA algorithms to the RFI flagged data
set with respect to the no RFI scenario, respectively. The bottom right panel
shows the 3D power spectrum after spherically averaging the modes above the
horizon limit for CLEAN (blue) and LSSA (red) algorithms. Here we also show
the 3D power spectrum for the no RFI case (green) and power spectrum of the
injected 21cmFAST signal (black dashed curve).

##  4 Application on Simulated Data

To understand the effect of RFI flagging in estimating the cosmological Hi
power spectrum, we first flagged frequency channels before performing the FT
along the frequency axis using Eqn. 5 and then estimated the cylindrically and
spherically averaged power spectrum. Below we discuss the different cases we
consider in our analysis and the corresponding results.

###  4.1 Case Studies

To understand the effectiveness of the CLEAN and LSSA algorithms to reduce the
spectral leakage for realistic flagging scenarios, we study the following
cases:

Case I: There is no flagging, i.e., no RFI scenario. The power spectrum is
estimated from this unflagged sample.

Case II: We randomly selected 20% channels and flagged them. This is the most
simplistic scenario and can be observed in data taken with any radio
telescope. Here, we conservatively choose 20% random channels to flag.
However, note that some of the low-frequency radio telescope s were built in
radio-quiet zones, where the percentage of random flagging may be less than
this. After applying the random flags to the visibility data cube, we use
CLEAN and LSSA algorithms to go to the delay domain and estimate the power
spectrum.

Case III-A: The FPGA-based digital-receiver employed at MWA digitizes the
beam-formed signal and performs a first stage coarse frequency channelization
of the data, i.e., 1.28 MHz coarse-bands across 32 MHz bandwidth of MWA data
(Prabu et al. 2015). This step applies a filter shape and aliases channels on
the edge of each coarse-band (25 coarse bands across 32 MHz bandwidth) (Prabu
et al. 2015; Beardsley et al. 2016). To avoid this known aliasing effect,
Dillon et al. (2015) flagged the 160 KHz channel (after smoothing the data by
2 channels) in every 1.28 MHz coarse-band in their analysis of the MWA data.
Beardsley et al. (2016) also flagged two 40 KHz channels on either side of the
coarse-band edges (total of four channels per coarse-band edge) to avoid the
aliasing effect. In addition, they also flagged the central 40 KHz channel of
each coarse-band of MWA data, which corresponds to coarse-band DC mode. These
channels have been observed to contain anomalous power, likely due to very
low-level rounding errors in the polyphase filter bank of the digital
receivers. Barry et al. (2019) also used the same flagging strategy in their
updated analysis of MWA data to get a more robust upper limit on the Hi power
spectrum. The flagging of channels at either side of each coarse-bands and at
the center of the coarse-band creates a periodic flagging scenario in addition
to the random flagging due to persistent RFI. To incorporate this in our
analysis, we flag 20% randomly selected channels, and also 200 KHz flags are
introduced in every 1.28 MHz across the 8 MHz bandwidth of our simulated data.
We call this flagging scenario as random plus periodic flagging.

Case III-B: This is similar to the previous Case III-A, with the change in
periodicity and increase of flagging channel width. In this case, in addition
to random flagging, we introduced 300 KHz flags in every 0.64 MHz across the 8
MHz bandwidth.

Case IV: In this case, we introduce the low-level (faint) broadband RFI
flagging corresponding to digital TV (DTV) signal as reported in MWA
observation (Wilensky et al. 2019; Barry et al. 2019). Wilensky et al. (2019)
developed a pipeline, called SSINS, which can identify those ultra-faint RFI
and remove the contaminated data (Wilensky et al. 2019). Barry et al. (2019)
found a significant improvement in the estimated power spectrum after removing
the low-level broadband RFI contaminated data using SSINS in their re-analysis
of the MWA data set. To incorporate this broadband RFI in our analysis, we
flag a 2 MHz wide frequency chunk at the center of the band, in addition to
random and periodic flagging.

The different cases are also summarized in Table 2. We compare the results
from these different cases in the following section.

|
---|---
|
Figure 6:  This figure is same as Fig. 5, but for the Case III-A. Here the
periodic flagging is being applied to the data on top of random flagging, i.e,
200 KHz flags in every 1.28 MHz across the 8 MHz bandwidth (see text).

###  4.2 Results

Case I: We show in Fig. 4 (left panel), the cylindrically averaged 2D power
spectrum in $\mathbb{k}_{\perp}$ and $k_{\parallel}$ space (as given by
Equation 7) when there is no RFI, i.e, no flagging involved. We find that the
spectrally smooth bright foregrounds coupled with instrumental response are
confined within a ‘wedge’ shape region in this 2D space. The black dashed line
is the horizon limit with an additional buffer due to the Blackman-Harris
window used during FT along the frequency direction. The region above this
line contains less power than the foreground wedge and is popularly known as
the ‘EoR window’. One can avoid the foreground-dominated modes inside the
‘wedge’ and use only the modes inside the ‘EoR window’ to search for Hi 21 cm
signal. This technique is called ‘foreground avoidance’. The spherically
averaged 3D power spectrum using the modes outside the foreground wedge has
been shown in the right panel of Fig. 4. The black dashed line shows the
injected Hi 21 cm power spectrum from 21cmFAST modulated by the instrument
beam. Here, we assume a simple Gaussian beam of the instrument. We find that
the estimated 3D power spectrum by avoiding the foreground for no RFI scenario
is consistent with the input Hi 21 cm signal. This case study also provides
the null test, where we can recover the injected Hi 21 cm power spectrum using
our power spectrum estimation methodology as described earlier.

|
---|---
|
Figure 7:  This figure is same as Fig. 6, but for the Case III-A, i.e, where
we introduced 300 KHz flag in every 0.64 MHz across the 8 MHz band in addition
to the random flagging (see text).

Case II: In Fig. 5, we show our findings for random flagging scenario. We plot
the ratio of the 2D cylindrically averaged power spectrum for different
algorithms applied to the randomly flagged data. At the top left of Fig. 5, we
show the ratio of the 2D power spectrum of flagged to no RFI scenario in log-
scale. We find a significant spillover of foregrounds into the ‘EoR window’
when there are missing samples due to RFI flagging. Note that, although we use
the BH window during FT along the frequency axis in this case, it completely
fails to restrict the bleed-in of foregrounds into the ‘EoR window’. The use
of a tapering window function during the FT helps to suppress the effects of
the sharp edge of the finite bandpass (Thyagarajan et al. 2013), but it can
not deal with the missing samples due to RFI flagging. At the top right and
bottom left panels of Fig. 5, we show the ratio after applying the CLEAN and
LSSA algorithms to the RFI flagged data set with respect to the no RFI
scenario, respectively. Here, we use the tolerance parameter in CLEAN to be
$10^{-9}$. We find that there are only a few modes outside the foreground
wedge, which have a little foreground contamination in comparison with the no
RFI case when we use the CLEAN algorithm. The contamination is even less when
applying LSSA to the flagged data set during FT. Nevertheless, these two plots
show that one can restrict the foreground spillover to a minimum after using
CLEAN and LSSA algorithms to the flagged data set in estimating the power
spectrum. After applying these algorithms, one can choose a foreground-free
‘window’ from this 2D cylindrical power spectrum and spherically average the
modes to estimate the 3D power spectrum. At the bottom right of Fig. 5, we
show the 3D power spectrum after spherically averaging all the modes above the
horizon line, i.e., the modes inside the ‘EoR window’. We compare the result
of the 3D power spectrum for CLEAN and LSSA algorithms to the no RFI scenario
(Case I) and the input 21cmFAST Hi power spectrum. We find that we can still
recover the $k-$modes of the Hi 21 cm power spectrum in the presence of the
random flagging after using these two algorithms.

Case III-A: The results for random and periodic flagging are shown in Fig. 6.
Similar to the findings of Case II, here, we also find that the presence of
periodic flag along with the random flagging severely contaminates the ‘EoR
window’ when only the BH tapering function is being applied (top left panel of
Fig. 6). At the top right panel, we show the ratio of the 2D power spectrum
after the application of CLEAN to the periodic plus randomly flagged data set.
We find that the Fourier modes inside the EoR window are still contaminated by
foreground spillover after applying the CLEAN algorithm. The extra bias in
CLEAN arises due to the un-subtracted foregrounds in the residual, which still
contains side-lobes from RFI flagging (Ewall-Wice et al. 2021). To CLEAN the
spectrum much deeper than Case II, we choose the tolerance parameter to be
$10^{-11}$. But, even after lowering the tolerance parameter, we find
significant bias due to foreground leakage into the EoR window. However, when
we apply the LSSA to the flagged data set (bottom left panel), we can restrict
the foreground contamination of Fourier modes beyond the ‘wedge’ in the
presence of periodic flagging. The 3D power spectrum at the bottom right panel
also shows that there is an extra bias to the power spectrum for the CLEAN
algorithm with respect to the no RFI scenario. But, LSSA can mitigate this
foreground leakage in the presence of periodic flagging and gives an unbiased
estimate of the 3D power spectrum.

|
---|---
|
Figure 8:  This figure is same as Fig. 6, but for the Case IV. Here we
introduced ultra-faint broadband RFI flagging in addition to the random and
periodic flagging. A 2 MHz chunk at the center of the band is being flagged,
which corresponds to broadband RFI (DTV signals).

Case III-B: The results for the increased flagging channel width (300 kHz) and
a decreased periodicity (in every 0.64 MHz) than Case III-A are shown in Fig.
7. The results are nearly similar to the Case III-A when we apply BH and LSSA
algorithms to the flagged data set. The application of the BH tapering
function shows a significant foreground leakage above the horizon. However,
applying the LSSA algorithm shows that this leakage is negligible. We also
find that, unlike Case III-A, the application of CLEAN gives higher foreground
leakage into the EoR window. Note that, here also, we use the tolerance
parameter to be $10^{-11}$ in the CLEAN algorithm to reduce the foreground
residual. The spherically averaged 3D power spectrum also shows the same
feature. We get a large bias due to the foreground leakage at almost all
scales when applying the CLEAN algorithm. However, we do not find any such
bias in the spherically averaged power spectrum when applying the LSSA
algorithm. The analysis of these two cases, Case III-A and Case III-B, show
that CLEAN is not an optimal way to restrict foreground leakage in the
presence of periodic flagging (with any periodicity) along with random
flagging, and it gives a significant bias to the estimated 3D power spectrum.
However, the foreground leakage due to the periodic RFI flagging can be
restricted with the application of LSSA, and the modes within the EoR window
become available again for detection of the cosmological Hi power spectrum.

Case IV: We show the results for the broadband RFI flagging scenario in Fig.
8. We find that when a large frequency chunk is missing due to ultra-faint
broadband RFI flagging (Wilensky et al. 2019), the CLEAN algorithm shows a
significant amount of foreground leakage in comparison with all other cases.
The application of the BH tapering function shows similar results as with
previous cases. However, applying the LSSA algorithm constrains the leakage
during the FT along the frequency axis in this case. The spherically averaged
power spectrum also shows the presence of a substantial bias in the estimated
power spectrum in the case of CLEAN. In contrast, we can still recover the
Fourier modes inside the EoR window in the case of LSSA. The LSSA can give an
unbiased and robust estimate of the underlying power spectrum in the presence
of broadband RFI flagging in addition to the random and periodic flag.

The different flagging scenarios used in this comparative analysis prove that
LSSA gives a more sensitive and robust way to estimate the cosmological Hi
power spectrum from the observed data. However, a sufficiently deep CLEANing
in the delay space can mitigate the foreground leakage issue in the presence
of random RFI flagging only.

##  5 Application to the uGMRT DATA

Here, we extend our comparative analyses to a real observation of the ELAIS-N1
field observed with the uGMRT (in GTAC cycle 32) during May 2017. The details
of the observation, a brief discussion on data analysis and the estimation of
the power spectrum from the observed visibility data set are discussed in the
following sub-sections.

Table 3: Observational details of the target field ELAIS N1 for four observing sessions Project code | 32_120
---|---
Observation date | 5, 6, 7 May 2017
| 27 June 2017
Bandwidth | 200 MHz
Frequency range | 300-500 MHz
Channels | 8192
Integration time | 2s
Correlations | RR RL LR LL
Total on-source time |  $\sim$13 h (ELAIS N1)
Working antennas | 26
Pointing centre |  $16^{h}10^{m}01^{s}$ $+54^{d}30^{m}36^{s}$ (ELAIS N1)

###  5.1 Observations and Data Analysis

The data was taken over four nights, where the total integration time,
including calibrators, is 25 hours. The total bandwidth (200 MHz) between
300-500 MHz is divided into 8192 channels, which results in a 24 kHz frequency
resolution. The raw data has a time resolution of 2 sec. The observational
details are tabulated in Table 3. The details of data analysis are presented
in (Chakraborty et al. 2019). The raw data was flagged for RFI with AOFLAGGER
(Offringa et al. 2012) and then averaged to 16 sec time resolution to reduce
the data volume. About 40% of the data was flagged after RFI flagging,
calibration, and self-calibration steps. We have not averaged the data across
frequency coordinate to get the maximum $k_{\parallel}$ modes and work with 24
kHz frequency resolution. This choice will allow us to access a larger region
of the EoR window above the foreground wedge.

|
---|---
Figure 9: The cylindrically averaged 2D power spectrum after applying the
CLEAN (left panel) and LSSA (right panel) to the one night’s data set of the
ELAIS N1 field observed with the uGMRT. The data spans a 8 MHz bandwidth
around 310 MHz ($z=3.58$). The black dashed line in both panels is the horizon
limit.

###  5.2 Power Spectrum estimation from the data

In this work, we show the result of a small 8 MHz chunk around 304 MHz ($z\sim
3.6$) of one night’s data ($\sim$ 6 hours including calibrators). This
demonstrates the effect of CLEAN and LSSA algorithms in estimating the power
spectrum with realistic flagging applied to the observed visibilities. In our
previous work Chakraborty et al. (2021), we have already analyzed the whole
data set and put the upper limit on the post-EoR Hi 21 cm power spectrum at
redshifts $z=1.96,2.19,2.62,3.58$. In Chakraborty et al. (2021), we have used
the CLEAN algorithm to mitigate the missing samples issue due to RFI flagging.
However, here we also apply the LSSA along with CLEAN on a single night data
set to compare the effectiveness of these two algorithms in a realistic
scenario.

The details of the estimation of the power spectrum from the observed uGMRT
data set are mentioned in Chakraborty et al. (2021) and we briefly mention the
methodology here. The simple squaring of visibilities as described in Eqn. 7
to estimate the power spectrum will result in positive noise bias since the
observed data also contains noise (Bharadwaj & Pandey 2003; Ali et al. 2008).
To avoid this, we grouped the visibilites within a $uv$-cell whose dimension
is governed by the inverse of the half-power beamwidth of the primary beam
($\theta_{\mathrm{HPBW}}$). We then cross-correlate the visibilities
associated with a particular $uv$-cell. This results in a correlation matrix
for each $uv$-cell. The diagonal terms of that matrix are the self-correlation
of visibilities and contain the positive noise bias. The off-diagonal terms
are the cross-correlation of visibilities that lie inside the $uv$-cell. The
noise is expected to be uncorrelated between two visibilities, and so the off-
diagonal terms do not contain the positive noise bias (Ali et al. 2008). We
take the average of the off-diagonal terms of a correlation matrix of a
particular $uv$-cell and quote it as an estimated power for that baseline. The
underlying assumption is that the cosmological Hi signal is correlated among
visibilities that fall within a $uv$-cell whose dimension is inverse of the
primary beam. This formalism is proposed by Bharadwaj & Pandey (2003) and
investigated in detail using observation by Ali et al. (2008), Ghosh et al.
(2012), Choudhuri et al. (2016), Bharadwaj et al. (2019), Chakraborty et al.
(2021), Pal et al. (2021).

We have applied both CLEAN and LSSA to transform the observed visibility data
from frequency to delay domain in the presence of realistic flagging and then
estimate the power spectrum following the method as described above. In Fig.
9, we show the cylindrically averaged 2D power spectrum after applying CLEAN
(left panel) and LSSA (right panel) algorithms to the observed visibilities.
We find that both of CLEAN and LSSA can restrict the foreground leakage above
the foreground wedge and give nearly identical results. Above the horizon
line, there is a region where the value of the power spectrum is much less in
comparison to the power inside the wedge. This region is relatively free of
foreground contamination, and one can spherically average the
$\mathbb{k}_{\perp}$ and $k_{\parallel}$ modes of that region to estimate the
3D power spectrum. In Fig. 10, we show the ratio of the 2D power spectrum
estimated using CLEAN and LSSA algorithms. The ratio is close for most of the
$k$-modes above the horizon limit. This shows that in the presence of
realistic random flagging of the uGMRT data set, both CLEAN and LSSA give
nearly identical results and are consistent with each other. Note that, here,
we do not estimate the spherically averaged power spectrum of cosmological Hi
signal with this small data set. We only demonstrate the effectiveness of
CLEAN and LSSA algorithms in case of realistic flagging due to RFI. The 3D
spherically averaged power spectrum using the modes free of foreground
contamination above the horizon line for different night’s data set as well as
for the combined data set are already presented in Chakraborty et al. (2021).

Figure 10: The ratio of the 2D power spectrum between LSSA and CLEAN
algorithms after applying those to the data set of uGMRT. This plot is
basically the ratio of the 2D power spectrum shown in Fig. 9.

##  6 Summary and Conclusion

The flagging of RFI-contaminated channels can cause excess power leaking into
the EoR window (above the horizon limit) in the 2D cylindrical power spectrum.
The non-uniform sampling of the bandpass due to RFI excision causes large
sidelobes in the delay space PSF. The convolution of this non-ideal delay
space PSF with the delay spectrum of the data causes foreground leakage into
the EoR window. One dimensional CLEAN algorithm in delay space and Least
Square Spectral Analysis (LSSA) are two widely used algorithms, which have
been applied to the observed data set to mitigate this foreground leakage
issue. In this work, we have done a comparative analysis of these two
algorithms using simulations in the presence of different realistic flagging
scenarios.

We have simulated the data for SKA-1 low array configuration with cosmological
Hi 21 cm signal and foreground signals as sky model at redshift $z\sim 9$. We
have compared three different RFI flagging scenarios. At first, we randomly
flagged 20% channels of the simulated data and then applied CLEAN and LSSA to
restrict the foreground leakage into the EoR window (Case II). We estimate the
2D cylindrically averaged power spectrum and compare the result with the no
RFI scenario, i.e., no flagging. We have found that both of these algorithms
can restrict the foreground leakage into the EoR window for random flagging of
channels. We spherically averaged the modes inside the ‘EoR window’ to
estimate the 3D power Hi spectrum. We have compared the 3D power spectrum with
the no RFI scenario as well as with the injected Hi 21 cm signal to our
simulation, as taken from 21cmFAST. The results are consistent with each
other. We conclude that in the presence of random flagging, the application of
CLEAN and LSSA can recover the $k-$modes to extract the cosmological Hi 21 cm
signal power spectrum.

Then we study the case where we have flagged 200 KHz channels in every 1.28
MHz across the whole 8 MHz band of our simulated data set in addition to 20%
random flagging (Case III-A). This is identical to observed data with the MWA
telescope. We have also changed the periodicity and flagging channel width
(Case III-B), i.e., we flag 300 KHz channels in every 0.64 MHz across 8 MHz
bandwidth. We applied CLEAN and LSSA and compared the 2D and 3D power spectrum
with no RFI scenario and with the input Hi 21 cm power spectrum. We have found
that the CLEAN algorithm fails to restrict the foreground leakage for both
cases, and there is an extra bias in the 3D power spectrum. However, the LSSA
method can limit the foreground leakage even in the presence of periodic
flagging. We can average the modes above the ‘wedge’ and reconstruct the Hi 21
cm power spectrum using the LSSA algorithm. We conclude from this analysis
that LSSA works better than CLEAN and gives an unbiased estimate of the power
spectrum in the presence of periodic flagging.

Lastly, we study the scenario where a 2MHz wide frequency chunk has been
flagged from the center of the band in addition to the random and periodic
flagging. This scenario corresponds to the ultra-faint broadband RFI due to
the dish TV signal. After applying different algorithms for this case, we
found that CLEAN gives a significant bias in the power spectrum and fails to
restrict the foreground leakage. But, LSSA can constrain the leakage and
provide an unbiased estimate of the Hi power spectrum.

We have also applied these algorithms to the observed data set of the ELAIS N1
field using uGMRT to verify these findings with realistic flagging. However,
we only choose a small 8 MHz chunk around 310 MHZ of one night’s data (6
hours) and estimate the 2D power spectrum using CLEAN and LSSA algorithms.
Here also, we have found that both of these algorithms give identical results
in restricting the foreground contamination above the ‘wedge’, since the uGMRT
data contains the random flagging only. We can limit the foreground spillover
with these algorithms and extract the cosmological Hi signal from the uGMRT
observation. The first upper limit on post-EoR Hi signal using this data set
has already been presented in Chakraborty et al. (2021), where we have used
the CLEAN algorithm to mitigate the missing sample issue. However, this
analysis with the data shows that the power spectrum estimated using the LSSA
algorithm is consistent with the CLEAN algorithm.

In summary, this paper explores the effectiveness of 1D CLEAN and LSSA
algorithms to deal with missing frequency channel issues in CD/EoR
observations. We show that both of these algorithms are effective in the
presence of random flagging. However, CLEAN fails in the presence of periodic
and broadband RFI flagging. But, the LSSA can still restrict the foreground
spillover beyond the horizon limit and give a robust estimate of the Hi 21 cm
power spectrum. This methodology gives a sensitive estimation of Hi power
spectrum even in the case where 20% random channels are flagged, as well as a
periodic flagging of channels have been applied across the band. The
comparative analysis presented here establishes the advantages and necessity
of using one algorithm over the other to estimate Hi power spectrum from the
observed data infested with RFI. This work is relevant for upcoming telescopes
with CD/EoR science goals like the SKA 555https://www.skatelescope.org/, HERA
666https://reionization.org/ and also for the ongoing projects like the LOFAR
777http://www.lofar.org/astronomy/eor-ksp/lofar-eor-project/lofar-epoch-
reionization-project.html, MWA 888https://www.mwatelescope.org/science/epoch-
of-reionization-eor, PAPER 999http://eor.berkeley.edu/, etc. However, it
should be noted that the missing channels issue is not only limited to CD/EoR
frequencies but also to post-EoR redshifted frequencies which are important
for the Hi 21cm intensity mapping experiments like the CHIME
101010https://chime-experiment.ca/en , TIANLAI
111111https://tianlai.bao.ac.cn, HIRAX 121212https://hirax.ukzn.ac.za/, the
SKA-1 mid 131313https://www.skatelescope.org/mfaa/ and OWFA
141414http://www.ncra.tifr.res.in/ncra/ort.

ACKNOWLEDGEMENTS
We thank the anonymous referee and the scientific editor for helpful comments
and suggestions that have helped to improve this work. AD would like to
acknowledge the support of EMR-II under CSIR No. 03(1461)/19. AM would like to
thank Indian Institute of Technology Indore for supporting this research with
Teaching Assistantship.

DATA AVAILABILITY

The simulated data cube used in this analysis will be shared on reasonable
request to the corresponding author. The point source catalog of ELAIS-N1
field at 400 MHz compiled from uGMRT observation can be found here
https://github.com/Arnabiitindore/Catalog-of-ELAIS-N1-field and in the link
https://doi.org/10.5281/zenodo.5822497.

$Software$ :

This work heavily relies on the Python programming language
(https://www.python.org/). The packages used here are pyuvdata
(https://github.com/RadioAstronomySoftwareGroup/pyuvdata; aipy
(https://github.com/HERA-Team/aipy/tree/main/aipy); Hazelton et al. 2017),
astropy (https://www.astropy.org/; Astropy Collaboration et al. 2013, 2018),
numpy (https://numpy.org/) van der Walt et al. 2011; Harris et al. 2020, scipy
(https://www.scipy.org/) Virtanen et al. 2020, matplotlib
(https://matplotlib.org/) Hunter 2007, OSKAR
(https://github.com/OxfordSKA/OSKAR),
WSCLEAN(https://gitlab.com/aroffringa/wsclean).

## References

  * Ali et al. (2008) Ali, S. S., Bharadwaj, S., & Chengalur, J. N. 2008, MNRAS, 385, 2166, doi: 10.1111/j.1365-2966.2008.12984.x
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Bacon et al. (2020) Bacon, D. J., Battye, R. A., Bull, P., et al. 2020, Publications of the Astronomical Society of Australia, 37, e007, doi: 10.1017/pasa.2019.51
  * Barning (1963) Barning, F. J. M. 1963, Bull. Astron. Inst. Netherlands, 17, 22
  * Barry et al. (2019) Barry, N., Wilensky, M., Trott, C. M., et al. 2019, ApJ, 884, 1, doi: 10.3847/1538-4357/ab40a8
  * Beardsley et al. (2016) Beardsley, A. P., Hazelton, B. J., Sullivan, I. S., et al. 2016, ApJ, 833, 102, doi: 10.3847/1538-4357/833/1/102
  * Bharadwaj & Ali (2005) Bharadwaj, S., & Ali, S. S. 2005, MNRAS, 356, 1519, doi: 10.1111/j.1365-2966.2004.08604.x
  * Bharadwaj et al. (2019) Bharadwaj, S., Pal, S., Choudhuri, S., & Dutta, P. 2019, MNRAS, 483, 5694, doi: 10.1093/mnras/sty3501
  * Bharadwaj & Pandey (2003) Bharadwaj, S., & Pandey, S. K. 2003, Journal of Astrophysics and Astronomy, 24, 23, doi: 10.1007/BF03012189
  * Bharadwaj & Sethi (2001) Bharadwaj, S., & Sethi, S. K. 2001, Journal of Astrophysics and Astronomy, 22, 293, doi: 10.1007/BF02702273
  * Bull et al. (2015) Bull, P., Ferreira, P. G., Patel, P., & Santos, M. G. 2015, ApJ, 803, 21, doi: 10.1088/0004-637X/803/1/21
  * Chakraborty et al. (2019) Chakraborty, A., Roy, N., Datta, A., et al. 2019, MNRAS, 490, 243, doi: 10.1093/mnras/stz2533
  * Chakraborty et al. (2021) Chakraborty, A., Datta, A., Roy, N., et al. 2021, ApJ, 907, L7, doi: 10.3847/2041-8213/abd17a
  * Chapman et al. (2013) Chapman, E., Abdalla, F. B., Bobin, J., et al. 2013, MNRAS, 429, 165, doi: 10.1093/mnras/sts333
  * Choudhuri et al. (2016) Choudhuri, S., Bharadwaj, S., Chatterjee, S., et al. 2016, MNRAS, 463, 4093, doi: 10.1093/mnras/stw2254
  * Datta et al. (2010) Datta, A., Bowman, J. D., & Carilli, C. L. 2010, ApJ, 724, 526, doi: 10.1088/0004-637X/724/1/526
  * Dayal & Ferrara (2018) Dayal, P., & Ferrara, A. 2018, Phys. Rep., 780, 1, doi: 10.1016/j.physrep.2018.10.002
  * DeBoer et al. (2017) DeBoer, D. R., Parsons, A. R., Aguirre, J. E., et al. 2017, PASP, 129, 045001, doi: 10.1088/1538-3873/129/974/045001
  * Deshpande et al. (1996) Deshpande, A. A., Alessandro, F. D., & McCulloch, P. M. 1996, Journal of Astrophysics and Astronomy, 17, 7, doi: 10.1007/BF02709341
  * Deshpande et al. (2000) Deshpande, A. A., Dwarakanath, K. S., & Goss, W. M. 2000, ApJ, 543, 227, doi: 10.1086/317104
  * Di Matteo et al. (2004) Di Matteo, T., Ciardi, B., & Miniati, F. 2004, MNRAS, 355, 1053, doi: 10.1111/j.1365-2966.2004.08443.x
  * Dillon et al. (2015) Dillon, J. S., Neben, A. R., Hewitt, J. N., et al. 2015, Phys. Rev. D, 91, 123011, doi: 10.1103/PhysRevD.91.123011
  * Ewall-Wice et al. (2021) Ewall-Wice, A., Kern, N., Dillon, J. S., et al. 2021, MNRAS, 500, 5195, doi: 10.1093/mnras/staa3293
  * Field (1958) Field, G. B. 1958, Proceedings of the IRE, 46, 240, doi: 10.1109/JRPROC.1958.286741
  * Finkbeiner (2003) Finkbeiner, D. P. 2003, ApJS, 146, 407, doi: 10.1086/374411
  * Furlanetto et al. (2006) Furlanetto, S. R., Oh, S. P., & Briggs, F. H. 2006, Phys. Rep., 433, 181, doi: 10.1016/j.physrep.2006.08.002
  * Gehlot et al. (2019) Gehlot, B. K., Mertens, F. G., Koopmans, L. V. E., et al. 2019, MNRAS, 488, 4271, doi: 10.1093/mnras/stz1937
  * Ghosh et al. (2012) Ghosh, A., Prasad, J., Bharadwaj, S., Ali, S. S., & Chengalur, J. N. 2012, MNRAS, 426, 3295, doi: 10.1111/j.1365-2966.2012.21889.x
  * Giardino et al. (2002) Giardino, G., Banday, A. J., Górski, K. M., et al. 2002, A&A, 387, 82, doi: 10.1051/0004-6361:20020285
  * Hamaker et al. (1996) Hamaker, J. P., Bregman, J. D., & Sault, R. J. 1996, A&AS, 117, 137
  * Harker et al. (2009) Harker, G., Zaroubi, S., Bernardi, G., et al. 2009, MNRAS, 397, 1138, doi: 10.1111/j.1365-2966.2009.15081.x
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hazelton et al. (2017) Hazelton, B., Beardsley, A., Pober, J., et al. 2017, Hera-Team/Pyuvdata: Version 1.2, v1.2, Zenodo, doi: 10.5281/zenodo.1044022
  * Högbom (1974) Högbom, J. A. 1974, A&AS, 15, 417
  * Hogg (1999) Hogg, D. W. 1999, arXiv e-prints, astro.  https://arxiv.org/abs/astro-ph/9905116
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Jacobs et al. (2015) Jacobs, D., Pober, J., Parsons, A., & Paper Team. 2015, in American Astronomical Society Meeting Abstracts, Vol. 225, American Astronomical Society Meeting Abstracts #225, 403.01
  * Jelić et al. (2008) Jelić, V., Zaroubi, S., Labropoulos, P., et al. 2008, MNRAS, 389, 1319, doi: 10.1111/j.1365-2966.2008.13634.x
  * Kern et al. (2019) Kern, N. S., Parsons, A. R., Dillon, J. S., et al. 2019, ApJ, 884, 105, doi: 10.3847/1538-4357/ab3e73
  * Kerrigan et al. (2018) Kerrigan, J. R., Pober, J. C., Ali, Z. S., et al. 2018, ApJ, 864, 131, doi: 10.3847/1538-4357/aad8bb
  * Kolopanis et al. (2019) Kolopanis, M., Jacobs, D. C., Cheng, C., et al. 2019, ApJ, 883, 133, doi: 10.3847/1538-4357/ab3e3a
  * Koopmans et al. (2015) Koopmans, L., Pritchard, J., Mellema, G., et al. 2015, in Advancing Astrophysics with the Square Kilometre Array (AASKA14), 1.  https://arxiv.org/abs/1505.07568
  * Li et al. (2018) Li, W., Pober, J. C., Hazelton, B. J., et al. 2018, ApJ, 863, 170, doi: 10.3847/1538-4357/aad3c3
  * Li et al. (2019) Li, W., Xu, H., Ma, Z., et al. 2019, MNRAS, 485, 2628, doi: 10.1093/mnras/stz582
  * Loeb & Furlanetto (2013) Loeb, A., & Furlanetto, S. R. 2013, The First Galaxies in the Universe
  * Lomb (1976) Lomb, N. R. 1976, Ap&SS, 39, 447, doi: 10.1007/BF00648343
  * Madau et al. (1997) Madau, P., Meiksin, A., & Rees, M. J. 1997, ApJ, 475, 429, doi: 10.1086/303549
  * Mertens et al. (2018) Mertens, F. G., Ghosh, A., & Koopmans, L. V. E. 2018, MNRAS, 478, 3640, doi: 10.1093/mnras/sty1207
  * Mertens et al. (2020) Mertens, F. G., Mevius, M., Koopmans, L. V. E., et al. 2020, MNRAS, 493, 1662, doi: 10.1093/mnras/staa327
  * Mesinger et al. (2016) Mesinger, A., Greig, B., & Sobacchi, E. 2016, MNRAS, 459, 2342, doi: 10.1093/mnras/stw831
  * Morales et al. (2019) Morales, M. F., Beardsley, A., Pober, J., et al. 2019, MNRAS, 483, 2207, doi: 10.1093/mnras/sty2844
  * Morales et al. (2012) Morales, M. F., Hazelton, B., Sullivan, I., & Beardsley, A. 2012, ApJ, 752, 137, doi: 10.1088/0004-637X/752/2/137
  * Morales & Hewitt (2004) Morales, M. F., & Hewitt, J. 2004, ApJ, 615, 7, doi: 10.1086/424437
  * Mort et al. (2017) Mort, B., Dulwich, F., Razavi-Ghods, N., de Lera Acedo, E., & Grainge, K. 2017, MNRAS, 465, 3680, doi: 10.1093/mnras/stw2814
  * Offringa et al. (2019) Offringa, A. R., Mertens, F., & Koopmans, L. V. E. 2019, MNRAS, 484, 2866, doi: 10.1093/mnras/stz175
  * Offringa et al. (2012) Offringa, A. R., van de Gronde, J. J., & Roerdink, J. B. T. M. 2012, A&A, 539, A95, doi: 10.1051/0004-6361/201118497
  * Offringa et al. (2014) Offringa, A. R., McKinley, B., Hurley-Walker, N., et al. 2014, MNRAS, 444, 606, doi: 10.1093/mnras/stu1368
  * Pal et al. (2021) Pal, S., Bharadwaj, S., Ghosh, A., & Choudhuri, S. 2021, MNRAS, 501, 3378, doi: 10.1093/mnras/staa3831
  * Parsons & Backer (2009) Parsons, A. R., & Backer, D. C. 2009, AJ, 138, 219, doi: 10.1088/0004-6256/138/1/219
  * Parsons et al. (2012) Parsons, A. R., Pober, J. C., Aguirre, J. E., et al. 2012, ApJ, 756, 165, doi: 10.1088/0004-637X/756/2/165
  * Parsons et al. (2010) Parsons, A. R., Backer, D. C., Foster, G. S., et al. 2010, AJ, 139, 1468, doi: 10.1088/0004-6256/139/4/1468
  * Parsons et al. (2014) Parsons, A. R., Liu, A., Aguirre, J. E., et al. 2014, ApJ, 788, 106, doi: 10.1088/0004-637X/788/2/106
  * Patil et al. (2017) Patil, A. H., Yatawatta, S., Koopmans, L. V. E., et al. 2017, ApJ, 838, 65, doi: 10.3847/1538-4357/aa63e7
  * Planck Collaboration et al. (2018) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2018, arXiv e-prints, arXiv:1807.06209.  https://arxiv.org/abs/1807.06209
  * Prabu et al. (2015) Prabu, T., Srivani, K. S., Roshi, D. A., et al. 2015, Experimental Astronomy, 39, 73, doi: 10.1007/s10686-015-9444-3
  * Press & Schechter (1974) Press, W. H., & Schechter, P. 1974, ApJ, 187, 425, doi: 10.1086/152650
  * Pritchard & Loeb (2012) Pritchard, J. R., & Loeb, A. 2012, Reports on Progress in Physics, 75, 086901, doi: 10.1088/0034-4885/75/8/086901
  * Procopio et al. (2017) Procopio, P., Wayth, R. B., Line, J., et al. 2017, PASA, 34, e033, doi: 10.1017/pasa.2017.26
  * Remazeilles et al. (2015) Remazeilles, M., Dickinson, C., Banday, A. J., Bigot-Sazy, M. A., & Ghosh, T. 2015, MNRAS, 451, 4311, doi: 10.1093/mnras/stv1274
  * Stoica et al. (2009) Stoica, P., Li, J., & He, H. 2009, IEEE Transactions on Signal Processing, 57, 843, doi: 10.1109/TSP.2008.2008973
  * Tegmark (1997) Tegmark, M. 1997, Phys. Rev. D, 55, 5895, doi: 10.1103/PhysRevD.55.5895
  * Thompson et al. (2017) Thompson, A. R., Moran, J. M., & Swenson, George W., J. 2017, Interferometry and Synthesis in Radio Astronomy, 3rd Edition, doi: 10.1007/978-3-319-44431-4
  * Thyagarajan et al. (2013) Thyagarajan, N., Udaya Shankar, N., Subrahmanyan, R., et al. 2013, ApJ, 776, 6, doi: 10.1088/0004-637X/776/1/6
  * Trott et al. (2016) Trott, C. M., Pindor, B., Procopio, P., et al. 2016, ApJ, 818, 139, doi: 10.3847/0004-637X/818/2/139
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * van Haarlem et al. (2013) van Haarlem, M. P., Wise, M. W., Gunst, A. W., et al. 2013, A&A, 556, A2, doi: 10.1051/0004-6361/201220873
  * Vedantham et al. (2012) Vedantham, H., Udaya Shankar, N., & Subrahmanyan, R. 2012, ApJ, 745, 176, doi: 10.1088/0004-637X/745/2/176
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Wang et al. (2010) Wang, J., Xu, H., Gu, J., et al. 2010, ApJ, 723, 620, doi: 10.1088/0004-637X/723/1/620
  * Wilensky et al. (2019) Wilensky, M. J., Morales, M. F., Hazelton, B. J., et al. 2019, PASP, 131, 114507, doi: 10.1088/1538-3873/ab3cad
  * Wilman et al. (2008) Wilman, R. J., Miller, L., Jarvis, M. J., et al. 2008, MNRAS, 388, 1335, doi: 10.1111/j.1365-2966.2008.13486.x
  * Wyithe & Loeb (2008) Wyithe, J. S. B., & Loeb, A. 2008, MNRAS, 383, 606, doi: 10.1111/j.1365-2966.2007.12568.x
  * Yatawatta et al. (2013) Yatawatta, S., de Bruyn, A. G., Brentjens, M. A., et al. 2013, A&A, 550, A136, doi: 10.1051/0004-6361/201220874
  * Zaldarriaga et al. (2004) Zaldarriaga, M., Furlanetto, S. R., & Hernquist, L. 2004, ApJ, 608, 622, doi: 10.1086/386327