{'fraction_non_alphanumeric': 0.028524222800387414, 'fraction_numerical': 0.007254496078394136, 'mean_word_length': 4.610270615810782, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Imaging the Epoch of Reionization: limitations from foreground confusion and imaging algorithms

Tomography of redshifted 21 cm transition from neutral Hydrogen using Fourier synthesis telescopes is a promising tool to study the Epoch of Reionization (EoR). Limiting the confusion from Galactic and Extragalactic foregrounds is critical to the success of these telescopes. Instrumental response or the Point Spread Function (PSF) of such telescopes is inherently 3 dimensional with frequency mapping to the Line of Sight (LOS) distance. EoR signals will necessarily have to be detected in data where continuum confusion persists; therefore, it is important that the PSF has acceptable frequency structure so that the residual foreground does not confuse the EoR signature. This paper aims to understand the 3 dimensional PSF and foreground contamination in the same framework. We develop a formalism to estimate the foreground contamination along frequency, or equivalently LOS dimension, and establish a relationship between foreground contamination in the image plane and visibility weights on the Fourier plane. We identify two dominant sources of LOS foreground contamination-'PSF contamination' and 'gridding contamination'. We show that 'PSF contamination' is localized in LOS wavenumber space, beyond which there potentially exists an 'EoR window' with negligible foreground contamination where we may focus our efforts to detect EoR. 'PSF contamination' in this window may be substantially reduced by judicious choice of a frequency window function. Gridding and imaging algorithms create additional 'gridding contamination' and we propose a new imaging algorithm using the Chirp Z Transform (CZT) that significantly reduces this contamination. Finally, we demonstrate the analytical relationships and the merit of the new imaging algorithm for the case of imaging with the Murchison Widefield Array (MWA).


Introduction
Several synthesis telescopes such as MWA 1 , LO-FAR 2 , and PAPER 3 are currently being built to study the Epoch of Reionization (EoR) by observing the redshifted 21 cm line transition from neutral hydrogen. While future instruments may map the 21 cm brightness in three dimensions (frequency corresponds to line of sight), current experiments lack the required sensitivity and instead hope to estimate the power spectrum of 21cm brightness fluctuations. The brightness of these fluctuations is about 5 orders of magnitude smaller than the Galactic and Extragalactic foregrounds. Moreover, removing the effects of the telescope point spread function (PSF) from the synthesis images may not be possible to the level of sensitivity demanded by EoR power spectrum measurements. Consequently, foreground contamination and instrumental effects present a major hurdle to 21 cm tomography.
A host of simulations have demonstrated the potential in current instruments to recover the EoR power spectrum despite thermal uncertainties, foreground contaminants, and instrumental effects. Among instrumental effects, simulations have shown the frequency dependence of PSF to be of particular concern in EoR power spectrum estimation (Bowman, Morales & Hewitt 2009). In 21 cm tomography, frequency maps to the line of sight, and any instrumental variation in frequency can potentially 'mimic' line of sight 21 cm brightness fluctuations. Most of the proposed foreground subtraction algorithms rely on the ability to distinguish smooth spectral features of Galactic and Extragalactic contaminants from the relatively rapid spectral fluctuations (expected from theory) of the 21 cm EoR signal. Spectral features in the telescope PSF may affect our ability to 'fit' and remove undesirable foregrounds. The sidelobe response of confusing foreground contaminants may vary with frequency and generate an instrumen-tal spectral structure at each image pixel-a phenomenon that has been dubbed 'mode-mixing' in literature (Bowman, Morales & Hewitt 2009). Apart from mode-mixing, calibration errors and choice of gridding algorithms can also lead to undesirable instrumental response in frequency domain.
While simulations have evaluated such modemixing effects for realistic arrays like MWA (Liu et al. 2009;Bowman, Morales & Hewitt 2009), this paper presents an approach that is analytical, and hence, more general. In particular, the relationship between the PSF in image coordinates (lmf ) and the telescope sampling function in Fourier coordinates (uvf ) is established, and the mode-mixing problem is cast in Fourier space where the telescope measurements are really made. In this formalism, implications for choice of experimental parameters and foreground subtraction algorithms are discussed. Apart from the analytical treatment, various practical data processing operations that influence the PSF are discussed using an example of snapshot imaging with MWA.
The remainder of the paper is organized as follows. Section 2 defines the problem analytically and derives an expression relating the line of sight foreground contamination to the weighting function in Fourier space. To better understand the relationships established in Section 2, Section 3 applies them to the simple case of imaging a point source with a complete visibility distribution. The analytical treatment in Section 2 does not take into account the effects of discretization of visibilities in the re-gridding and imaging algorithms. Section 4 discusses these effects and complements the results of Section 2. In particular, Section 4 describes a new imaging algorithm based on the Chirp Z Transform (CZT) that alleviates the problem of foreground contamination. Section 5 demonstrates the concepts developed in Sections 2, 3, and 4 through simulating a snapshot observation with MWA. Section 6 draws up conclusions and proposals for future work.

Analytical treatment
We first analytically evaluate the effect of a frequency variant PSF at a given image pixel due to an unsubtracted point source by evaluating its contribution to foreground contamination in LOS wavenumber space. We call this contamination 'PSF contamination'. A single point source is considered here to attain an intuitive understanding of 'PSF contamination'. Generalization for a realistic sky will be made in subsequent sections. Array element locations are assumed to be coplanar.

Notation
We define the co-ordinate system to avoid ambiguity. The baselines are defined on a local plane whose origin is at the latitude and longitude of the array center. A two dimensional vector x = [x y] (in meter units) represents the projected baseline on the local plane. The corresponding baseline u = [u v] (in wavelengths) at a frequency f is then given by u = x.(f /c), where c is the speed of light. A two-dimensional direction vector l towards any direction is represented by its direction cosines: l = [l m]. An interferometer observation gives a visibility cube with two spatial axes representing the baseline and a frequency axis. An image cube on the other hand has two axes representing direction cosines on the sky and a third frequency axis. Assuming co-planarity of array elements, within the framework of Fourier synthesis imaging, u and l are Fourier conjugates.

Weighting function
The telescope sampling function is represented as W (x, f ), which is the weight assigned to the baseline represented by the vector x, and includes any visibility taper applied prior to imaging. In general,W (x) depends on instrumental, observational and processing parameters. Among instrumental parameters is the antenna layout that determines the set of available baselines. Observational parameters, which include pointing direction, duration of rotation synthesis, and mode of observation-drift scan or track mode, define the uv coverage. Observational parameters along with the antenna layout determines the 'natural' weight at each x-number of visibilities falling into the pixel centered at x. Among processing parameters are the visibility grid size and extent, convolution kernel used to grid the visibilities, and any visibility mask and/or taper used during imaging. The weighting function may be factored into a frequency independent geometric component, a frequency dependent 'taper' component, and a frequency window function: Instrumental and observational parameters are geometric in nature and hence frequency independent. They may be sufficiently represented by W (x). Processing parameters like visibility 'taper' function and visibility masks are frequency dependent and are represented by T (x, f ). B(f ) is the window function in frequency assigned to each visibility.

Derivation
An image made with a weighted baseline distribu-tionW (x, f ) has a frequency dependent PSFÃ (l, f ), where l is the two dimensional unit vector representing direction on the sky. Tilde notation is used for image domain functions as opposed to their Fourier domain counterparts. Consider the sidelobe contribution at zenith (l = [0 0]) from a source with flux density of unity and spectral index of zero located at any l. This contributionS(l, f ) varies with frequency due to a changing PSF. We wish to compute the energy in these spectral variations on different scales in terms of the telescope weighting functionW (x, f ).
PSFÃ(l, f ) is the Fourier Transform of the uvweighting function A(u, f ). We may represent this asÃ where the denominator is included to obtainÃ(l, f ) as the normalized PSF for all frequencies, and essentially does bandpass calibration of visibilities. The inclusion ofÃ(l, f ) essentially does bandpass calibration of visibilities and results in a PSF that is normalised at all frequencies. Any additional frequency weighting may be absorbed into the frequency window function B(f ). The function A is just a scaled version of the functionW through the scaling relationship given by u = x.(f /c). We may represent this scaling relationship as The quantity uc/f has units of distance and is frequency invariant: it is simply the baseline distance in meters. We use this to make the substitution x = uc/f , y = vc/f and get where is the area un-derW . The importance of this substitution is that it expresses all baselines in units of meter as opposed to wavelength and makes them independent of frequency.
Consider a point source with a flux density of unity and a spectral index of zero located at l. The response at zenith due to a sidelobe on this source is given bỹ We are interested in the energy inS(l, f ) at different frequency scales. In other words, to evaluate 'PSF contamination', we have to compute the power spectrum S(l, η) by Fourier transformingS(l, f ) with frequency offset ∆f and time η as Fourier conjugates. We use ∆f instead of f because the spatial EoR fluctuations are evaluated about an origin r 0 corresponding to a frequency f 0 and redshift z 0 .
S(η, l) may be expressed as Using ∆f = f − f 0 , f 0 = (f L + f H )/2, where f L and f H are the lowest and highest frequencies in the instrument bandpass, and writing the conjugate variable as f instead of ∆f we get Substituting from Equation 4 forÃ(l, f ) and interchanging the order of integrations we get The inner integral, in general, may have to be evaluated numerically. However, to have an insight into the nature ofS(l, η), we evaluate Equation 8 for the case where there is no taper or mask applied in the uv domain. In this case, T (x, f ) ≡ 1, and A W is frequency independent. The inner integral is now just a scaled and shifted version of the Fourier transform of the frequency window function. Consider a rectangular frequency window defined by where B = f H − f L is the frequency bandwidth. The inner integral now evaluates tõ This gives us the relationship we are looking for: The power spectrumS(l, η) is a convolution of the weighting function W (x) and the convolving kernelB, which is a shifted version of the Fourier transform of the frequency window function B(f ).

Ideal case of complete visibility coverage
This section provides an insight into the terms and equations of the analytical derivation, by understanding them in the case of a visibility distribution with complete coverage. Complete coverage refers to the case where there exists at least one visibility measurement in every pixel in the uv grid at all frequencies within the instrument bandwidth, and over the spatial frequency range spanned by the array configuration. A key idea in the analytical treatment of Section  2 is the invariance of the PSF in lf . This concept is manifested in the invariance of Equation 4 in lf for a frequency invariant weighting. Frequency invariant weighting implies a 'taper' function T (x) that is invariant in frequency and dependent only on x. Figure  1 shows an image of the PSF in lf domain for a uniformly weighted visibility distribution with complete coverage. The plot extends over a 10% spread in both frequency and direction cosine (∆l/l o = ∆f /f o = 0.1). Due to the lf invariance, we expect the PSF to have the same value at all points defined by lf = l o f o . Consequently, the variation in an interval δl of the PSF evaluated at f o , is the same as the variation in an interval f o of the PSF evaluated at l o , so long as δl/l = −δf /f or approximately δl/l o = −δf /f o . This can be obtained by differentiating lf = l o f o and then approximating f to f o and l to l o , which hold for δl/l << 1 and δf /f << 1. Figure 1 also evaluates 'PSF contamination' in the zenith pixel in η space due to PSF sidelobes on a 1 Jy point source at l o . The localization of this contamination in light of Equations 10 and 11 is pivotal to our understanding of 'PSF contamination' and deserves a detailed explanation.

Localization of 'PSF contamination'
Consider a single point source at direction cosine l o entering the zenith pixel through the PSF sidelobe at l o . This is depicted in Figure 2. Due to the lf invariance, the frequency structure of the PSF at the zenith pixel is same as the sidelobe structure of the PSF in l evaluated at f o . The frequency structure of the PSF in a bandwidth δf is identical to the sidelobe structure of the PSF (evaluated at . This truncated version of the PSF can be extracted by multiplying the PSF with a window centered at l o and having width δf l o /f o . A Fourier transform of the product of these two functions is the convolution of their individual Fourier transforms which are the baseline weighting function W (x) and the convolution kernel from Equation 11. It may be seen from Equation 10 that the sinc-function envelope of the convolution kernel has a main lobe width equal to 2c/(l o ∆f ), which in general is far narrower than the visibility range 2x max . The rapid variation in the convolution kernel is a complex sinusoid in x with frequency l o f o /c. At any η, the convolution in Equation 11 is sensitive to a small range in W (x) that is centered at x = cη/l o with a width c/(l o ∆f ) which is the main-lobe width of the convolving kernel. For a visibility distribution with uniform weight and complete coverage like the one discussed, owing to the rapid fluctuations of the complex sinusoid, there exists significant response only around x max . This explains the delta-function like 'PSF contamination' in η centered at η = x max l o /c (see Figure 1). Sources from different directions in the sky have different values of l o and hence contribute to foreground contamination at different fairly localized values of η. The relative magnitude of these impulses is determined by the relative sidelobe levels at the respective directions. Figure 3 shows the delta-function like 'PSF contamination' in η for point sources entering the imaged pixel at zenith through different parts of the PSF (different l o ). A rectangular window has been assumed in frequency domain. The finiteness of W (x) gives a 1/l roll-off in sidelobe levels which corresponds to a 1/η roll-off in the contamination at different values of η. The magnitudes of such impulses are equal to the sidelobe level at the value of l o where they are evaluated. The solid black line is an approximate 1/η fit to the peaks of the delta-function like impulses. For η > x max l/c, the sidelobes of the convolving kernel give a response that further falls off as 1/η. The Blackman-Nuttall window substantially lowers contamination in the 'EoR window' (η > x max l o /c).

Sidelobes of the convolving kernel
We have used the phrase 'delta-function like' to describe the 'PSF contamination' for a source entering the zenith pixel from the direction l o . The exact profile is described by the sinc envelope of the convolving kernel. Owing to this, the η space may be divided into two regions. For l o ∈ [0 1], the delta-function like peaks are confined to η ∈ [0 x max /c]. This defines a region of relatively high foreground contamination, say R 1 . However, for η > x max /c, there still exists 'PSF contamination' due to the sidelobes of the convolving kernel. The contamination in this regime, say R 2 , falls off as 1/η since the sidelobes of a sincfunction envelope of the convolving kernel have a 1/x roll off. While EoR power spectrum estimation for η ∈ [0 x max /c] will suffer from stronger contamination, the level of contamination for η > x max /c may be reduced by a judicious choice of a bandpass window B(f ) that gives low sidelobe levels in R 2 , and hence, R 2 defines an 'EoR window' where we may focus our efforts to detect EoR.
In essence, a Fourier synthesis array gives a PSF that varies with frequency. If the baseline extent of the array is x max , the spectrum of frequency variations will have significant energy only in R 1 defined by η ∈ [0 , x max /c]. The instrument will sample the η-space up to η max = 1/B. In R 2 defined by η ∈ [x max /c , 1/B] the foreground contamination is due to the 1/η-form sidelobes of the contaminants in R 1 . These 1/η sidelobes are the result of a 'top hat' frequency window. A smooth window function in frequency like, for example, a Blackman-Nuttall window (Nuttall 1981) will substantially reduce the sidelobe contamination in R 2 (by more than three orders of magnitude) at the cost of reduced resolution in η space (see Figure 3). .5 x 10 −4 . The region of high foreground contamination R 1 extends up to η = x max /c = 3.33 µsec that corresponds to k ||m of about 1.2 Mpc −1 . The region between k ||m and k ||max is the 'EoR window' R 2 and is a region of low foreground contamination. The unfilled black circles define the point (k ⊥max , k ||m ) for redshifts of 3, 4, 5, 6, 7, 8, 9, and 10. Each concentric curve is the local of constant k = k 2 ⊥ + k 2 || .

Instrumental k space
The instrumental k space is an effective space to represent regions and extent of foreground contamination. The instrumental k space is 2 dimensional with the comoving line of sight wavenumber k || and comoving transverse wavenumber k ⊥ as axes. If k x , k y , and k z are comoving wavenumbers along three spatial axes with origin at some point about which the EoR fluctuations will be estimated, then  have shown that where the approximation holds for a small bandwidth (small redshift range). Here f HI = 1420 MHz is the rest frequency of 21 cm HI line emission, H 0 = 70 km sec −1 Mpc −1 is the present value of the Hubble con- Figure 4 shows a plot of the measurement space (in k ⊥ −k || domain) for a Fourier synthesis array with complete visibility coverage. A Fourier synthesis array is sensitive to measurements in a region in the k space bounded by some maximum and minimum values of k ⊥ and k || . The highest k ⊥ mode sampled by the instrument is related to the uv extent of the array and is given by while the smallest k ⊥ probed by the instrument is related to the antenna extent d, and is given by The highest value of k || sampled by the instrument is related to the frequency resolution and is given by The lowest value of k || sampled by the instrument in a snapshot observation is related to the frequency bandwidth and is given by The final comment in this section is related to the implications of lf invariance of the PSF to our understanding of the location and level of contamination in the instrumental k space. Consider the foreground contamination in the neighborhood (δl << 1) of some image pixel P due to a point source at some l o removed from P. Due to the lf invariance, the transverse structure of contamination around P is identical to the line of sight contamination at P. Since δl lo = δf fo , the energy in transverse fluctuations around P at wavenumber k ⊥• will be the same as the energy in line of sight fluctuations at P at wavenumber η o = lo fo k ⊥• . Using Equation 12, the (k ⊥ , k || ) pairs with identical contamination in them are defined by For a uniformly weighted visibility distribution with complete coverage, there is significant energy only in the transverse mode k ⊥ = x max f o /c. This gives us a corresponding η o = l o /f o k ⊥ = x max l o /c: a result we established in Section 3.1.

Effect of gridding and imaging algorithms
While Section 3 described the nature of foreground contamination along the line of sight and transverse wavenumber axes, its results are a good description of the contamination in the special case of an ideal array, and hence lay the foundations needed to appreciate issues related to EoR detection with practical array geometries and data processing algorithms. Real arrays seldom have complete visibility coverage, and the Fourier relationships between visibility and image space, and between frequency and line of sight wavenumber space are, in practice, evaluated as discrete transforms. Computationally efficient transforms like the Fast Fourier Transform (FFT) are often used in Fourier synthesis arrays and almost always involve re-gridding the visibility data on a regular uv grid. This section discusses the effects of such gridding and Fourier transform algorithms used for imaging on the nature of the 3 dimensional PSF and the resulting foreground contamination in k space. The first generation EoR experiments propose to differentiate confusing foregrounds from the EoR signal by attempting to synthesize a confusion-free 'EoR window' in the line of sight wavenumber dimension. We now describe the effect of gridding and imaging algorithms on the foreground contamination in η.
We have shown in Section 3 that contamination in η space is substantially different in the two regions (see Figures 3 and 4) R 1 and R 2 , where R 1 extends from η = 0 to η = x max /c, and R 2 (or 'EoR window') extends from η = x max /c to η = 1/∆f . We have also shown that contamination in R 2 is due to the sidelobes of the convolution kernelB (see Equation 11) and the delta-function like peaks from sources at various directions are confined to R 1 . Consequently, if region R 2 exists (1/∆f > x max /c), as it should in any well-designed EoR instrument that aims to synthesize a confusion-free 'EoR window', then PSF variation in frequency at any direction cosine l is oversampled by the instrument spectral resolution ∆f . This implies that, in the ideal case, there will not be significant channel to channel variation in flux versus frequency at any sky pixel. However, processes like gridding may induce channel to channel 'jitter' in the PSF and potentially contaminate R 2 . This stochastic component of foreground contamination arising from the gridding process may be termed 'gridding contam-ination'. We now describe such stochastic effects of gridding in terms of channel to channel migration of baseline vectors from one uv pixel to another. For simplicity, we have assumed that both depicted visibility measurements have equal weights. The upper panel shows the visibility axis for some frequency channel n, and the lower panel shows the visibility axis for the adjacent frequency channel n + 1. The depicted visibility segment is shown to have a fairly sparse u coverage for clarity. Notice how the baselines migrate on the visibility grid. While grid pixel number 1 and 8 go from being 'filled' to being a 'hole' between channel n and n + 1, visibility pixels 5 and 12 go from being a 'hole' to being 'filled'. Pixels like 2 and 10 continue to be filled but experience a 'step' change in their weighted visibility value.

Baseline migration
Fourier synthesis arrays have traditionally formed image cubes in lmf space by (i) computing the baselines (in wavelength units: u = x.f /c) generated by a given antenna distribution at every frequency channel, (ii) weighting the visibilities to get desirable PSF characteristics, (iii) performing gridding convolution on the visibilities to re-grid them on a regular uv grid at every frequency channel, and (iv) using an efficient FFT algorithm to compute a 'dirty' sky map at every frequency channel. Among these signal processing algorithms, gridding convolution has not been considered in the analytical treatment of Section 2. We now describe the influence of gridding convolution on foreground contamination in η space.
A baseline vector x (in meter units) will suffer a displacement of δu = xδf /c (in wavelength units) between adjacent frequency channels in the uv plane.
Here δf is the inter-channel frequency spacing. We call this phenomenon 'baseline migration'. This displacement results in 'step' changes in gridded visibility values of pixels within the support of the convolving kernel. Many such 'step' changes in visibilities due to the displacement of many baseline vectors results in a stochastic 'jitter' in the frequency structure of the PSF. PSF 'jitter' due to baseline migration may induce channel to channel fluctuations in flux versus frequency on the image plane, thereby significantly contaminating the 'EoR window' R 2 . This PSF 'jitter' due to 'baseline migration' is expected to have high contribution (i) from longer baseline vectors due to the relatively higher quantum of baseline displacement δu, and (ii) from regions of low uv density. In regions of higher uv density a given uv pixel, during gridding convolution, receives contributions from several neighboring baselines. The visibility contribution to a given uv pixel then depends on the relative displacement of the neighboring baselines, rather than their actual displacement. On the other hand, consider the extreme case of very low uv density as depicted in Figure 5. Here the actual inter-channel displacement of the baseline vectors induces relatively high 'steps' in the post-gridding convolution visibilities. Natural weighting results in significantly lower weight assigned to regions with (i) longer baseline vectors and (ii) low uv density and is expected to result in low PSF 'jitter' in frequency. Uniform weighting, on the other hand, increases the weighting for longer baselines and is expected to result in substantially high PSF 'jitter'.

The central uv void
The angular power spectrum of Extragalactic sources has been observed to have a power-law like form (Blake & Wall 2002) with substantially greater power on small angular scales; the universe is homogeneous and isotropic on large scales. Consequently, short baselines are not expected to have substantially greater response to Extragalactic continuum sources. However, the mean sky brightness will have a response at zero spacing (u = 0, v = 0). Energy in various spatial Fourier modes is seldom perfectly isolated and the zero spacing component will leak into short baselines. Consequently, channel to channel changes in the visibility distribution close to u = 0, v = 0 may result in 'step' changes in flux versus frequency at image pixels thereby contaminating higher values of η in R 2 . While Extragalactic sources are expected to have a smooth distribution on large angular scales, flux from the Galactic synchrotron emission is expected to follow a power law with significant spatial distribution power at large angular scales. Consequently short baselines close to u = 0, v = 0 will have significant response to our Galaxy and 'baseline migrations' in this part of the uv plane will create significant PSF 'jitter' that may contaminate R 2 .
Fourier synthesis arrays usually have a central 'uv void' as the antenna elements have a minimum spacing due to constraints imposed by their own physical size. Antenna cross-talk considerations often result in array configurations with a central 'uv void' that is several times the antenna size. In such cases, the central uv void spans several uv pixels. We mentioned in Section 4.1 that 'baseline migrations' in regions of high uv density do not contribute significantly to interchannel PSF 'jitter'. Though short spacing baselines lie in a densely populated part of the uv plane, with increasing frequency, most baseline vectors are displaced radially away from any uv grid point close to u = 0, v = 0, as opposed to other uv grid points where baseline vectors also get displaced towards them. This asymmetric neighborhood near short spacing visibilities results in relatively higher inter-channel 'step' changes in the visibilities on short baseline vectors. To alleviate this problem, the central uv void may be chosen to be larger than that imposed by antenna cross-talk constraints.

The Chirp Z transform
We have seen how channel to channel changes in the visibility weighting W (u, f ) contaminates the line of sight wavenumber dimension. Here W (u, f ) includes the effect of all visibility weighting and tapers used, and W (u, f ) determines the PSF via the imaging transform. If the visibility distribution has complete coverage, maintaining invariance of W (u, f ) with frequency is possible through application of a frequency independent visibility taper function. In such a case, W (u, f ) is same at all frequencies and there is no frequency dependence in the PSF. In practice, there are many holes in the uv grid and Labropoulos et al. (2009) have suggested that a visibility mask may be used that rejects the uv pixels that are not filled at all (or most) frequencies prior to image formation. At one end, we may exploit the complete instrument data and sensitivity at the cost of increased contamination in R 2 , and at the other end we may have a frequency independent PSF with some loss in information and sensitivity. However, there exists a 'middle path' that assures that no data is rejected while restricting the foreground contamination in R 2 -image formation using the Chirp Z Transform, or CZT (Rabiner, Schafer & Rader 1969). We first present a brief description of CZT.
The Z transform of any finite discrete sequence x n , n = 0, 1, 2....N − 1 is defined as Z transform is a function of the complex variable z. When evaluated at discrete points along the unit circle defined by z k = exp (i2πk/N ), where  /(N α)). In Case 1, the generalized CZT reduces to DFT, and for Case 2 the generalized CZT reduces to a special case of CZT that will be used in the proposed imaging algorithm of Section 4.4. k = 0, 1, 2....N − 1, Z Transform of the finite discrete sequence x n reduces to its Discrete Fourier Transform (DFT). The Chirp Z Transform or CZT is another special case of the Z transform that is evaluated at discrete points on the locus given by z k = AW −k where A and W are complex constants. Such a locus describes a spiral in the complex z plane. Notice how for A = 1 and W = exp (−i2π/N ) the CZT again reduces to a DFT.
Consider a case where |A| = 1 and W = exp (−i2π/(N α)) where α is a real constant. The locus of points on the complex z plane where the transform is evaluated now lie on a unit circle like in the case of a DFT. However, angular spacing between the points on the locus is now 2π/(nα) as opposed to 2π/N in case of a DFT. The first point on the locus is also shifted by the argument of the complex constant A. The loci for the two cases are shown in Figure 6. Though |A| = 1, W = exp (−i2π/(N α)) is a special case of the generalized Chirp Z Transform, we refer to it hereinafter as Chirp Z Transform or CZT with a chirp factor α. The only difference between DFT and CZT is the presence of the chirp factor α and an offset due to the presence of A. α simply 'stretches' the locus on which the transform is evaluated and the CZT of a sequence is essentially a DFT of a 'stretched' version of the sequence-the stretching factor being α. We now describe how visibilities can be gridded in meter units and how a chirp factor can be used to compensate for the 'stretching' of baseline vectors with frequency.

A CZT based imaging algorithm
As described in Section 4.1, the principal cause of channel to channel PSF variations in conventional gridding and imaging processes is baseline migration due to independent visibility gridding at every channel. However, the baseline vector generated by a given antenna pair moves across the visibility grid with change in frequency in a very predictable manner. The baselines in wavelength units at any frequency f are merely the baselines in meter units 'stretched' by a factor α, where α = f /c. In principle, we may grid the visibilities only once in meter units and use the CZT with chirp factor α to form the image at every frequency channel. Given an image extent, the grid size in meter units is chosen to avoid aliasing even at the highest frequency in the visibility cube. Gridding the visibilities only once in meter units eliminates the problem of non-linear baseline migrations and consequently ameliorates the problem of contamination in R 2 . In summary, we grid the baselines in meter units, and thereby do not let the baselines scale with frequency, and a given antenna pair will fall on exactly the same visibility grid point at all frequencies. We then compute the image on an lm grid that is constant over frequency and use the chirp factor α in a CZT to compensate for an unchanging baseline. Such a CZT based imaging transform is given by where i, j are uv (meter units) grid pixel indicies, m, n are image grid pixel indicies, function V 1 (x ij , f ) represents the visibilities at frequency f gridded on to the common grid, and α(f ) is the chirp factor given by α(f ) = c/f . The transform in Equation 18 is similar to a DFT and may be evaluated by the same efficient FFT algorithms with a minor modification to the exponential terms. The CZT image formation results in reduced contamination in R 2 and, moreover, at a substantially reduced computational cost of gridding convolution because the gridding convolution at all frequency channels is done to a single common grid with little additional computational cost for imaging.

The MWA example
We have so far discussed the nature of foreground contamination in η space for an ideal array with complete visibility coverage. We have also outlined some of the practical considerations applicable to realistic arrays and in the process proposed a new imaging technique using the Chirp Z Transform. In this section we demonstrate the relationships and results of Sections 2, 3, and 4 using the case of snapshot imaging with MWA. MWA consists of 512 antenna elements called 'tiles'. Each tile consists of a phased array of 16 bowtie antennas. The tiles are distributed to form baselines up to about 1 km. The instantaneous bandwidth is about 30 MHz, and we assume a frequency resolution of about 40 kHz. MWA will produce snapshot images every 8 seconds. The instrumental response for this instrument configuration in the k ⊥ − k || plane is shown in Figure 4. Details of the MWA telescope design may be found in Lonsdale et al. (2009).

Confusing sources and sidelobes
To estimate 'blending confusion', we use the source counts given by Subrahmanyan & Ekers (2002). The threshold flux density S C above which we expect to find one source per beam of full width at half maximum θ arcmin at frequency f GHz is given by S C is a confusion limit for a survey with excellent visibility coverage and all sources much stronger than S C can, in principle, be identified and subtracted. The stochastic response to sources weaker than S C contribute to sky pixel rms that is close to S C . S C for MWA evaluates to 2.3 mJy and is about 3 orders of magnitude higher than the expected thermal noise in 8-second snapshot images: MWA snapshot images will be confusion limited. For our simulations we assume successful peeling of point sources above 5S C ≈ 10 mJy. The expected residual rms confusion noise after successful subtraction of sources above 10 mJy is about 3.6 mJy. Due to the dominance of confusion over thermal noise, we will neglect the effects of thermal noise hereinafter. However, it must be noted that post foreground modeling and successful subtraction, the residual image that will be used for EoR power spectrum estimation is expected to be dominated, among other residual fluctuations, by thermal noise.
The threshold to which continuum sources may be identified and subtracted also depends on our ability to distinguish spillover flux due to sidelobes on the sky. We now roughly estimate the expected rms fluctuation at any image pixel due to the sidelobes on all the confusing sources in the sky-a quantity we call 'sidelobe confusion'. If we assume (i) the sidelobes of the PSF to be a zero mean random variable, and (ii) the confusion flux in each pixel to be a random variable that is uncorrelated with the sidelobe level, then the approximate rms sidelobe confusion can be estimated fairly easily: C S = i P i S i where P i is the sidelobe level at pixel i, and S i is the primary beam weighted true sky flux due to confusing sources at pixel i. Since P i = 0, the variance in sidelobe confusion is The far sidelobes of MWA are expected to have an rms variation of about 0.05% giving P 2 i = 0.25x10 −6 , and the RMS confusion is about 3.6 mJy giving S 2 i = 13 mJy 2 . The sidelobe level assumed here is true for sidelobes far from the phase center at which 'sidelobe confusion' is being computed. In reality sidelobe levels vary from sky pixel to pixel and a constant value has been assumed here for simplicity. Far sidelobe level has been chosen since far sidelobes result in foreground contamination at higher values of η where we expect to detect the EoR with minimal contamination, and this calculation will present the worst case. While P i is assumed to be independent of i, S i is assumed to be the confusion noise weighted by the primary beam gain, and the summation is performed using the expected primary beam of the MWA primary antenna element (Full width to first null of 60 • ) to evaluate 'sidelobe confusion'. The resulting 'sidelobe confusion' from the entire sky is C S ≈ 1.9 mJy. The similarity in the values of 'blending confusion' and 'sidelobe confusion' is an interesting consequence of the similarity between the number of resolution elements in the image (weighted by the primary beam) and inverse of the sidelobe level of the synthesized beam. This implies that after subtracting bright point sources from MWA snapshot images, the resulting residual image has statistical fluctuations with roughly equal contributions from the stochastic response to the Extragalactic sources and the sidelobe levels.
The EoR brightness fluctuations are expected to have an rms of about 25 mK on scales comparable to the MWA 5 arcmin resolution. This corresponds to an rms flux of about 36 µJy. Hence, apart from confusing foreground sources, 'sidelobe confusion' can confuse any attempt to detect the EoR signal through a simple measurement of image plane variance. That said, flux from confusing sources and PSF sidelobes have structure in the frequency domain which can be exploited to detect the EoR fluctuations despite a high rms fluctuation in the image plane. While spectral structure of confusion Extragalactic sources follow relatively simple power laws, spectral structure of instrumental sidelobe confusion is more complex as described in Sections 2, 3, and 4. We now quantitatively estimate the foreground contamination in frequency space (and eventually in η space) for the case of snapshot imaging with MWA.

A snapshot observation
This subsection describes a simulation to estimate the foreground contamination in η space for the case of snapshot imaging of Extragalactic point sources using MWA. The simulation assumes a frequency invariant sky model and assumes that all sources brighter than 10 mJy have been successfully subtracted. The sky model is then interpolated into an lm grid with grid size at least 2 times finer than the MWA instrument resolution. The grid extent is covers the interval [−0.5 0.5] in direction cosines that corresponds to the MWA primary beam full width at first null, which is 60 • at 150 MHz. Gridding convolution of visibilities is performed using a triangular kernel like the one shown in Figure 5. A simple kernel was chosen purely to illustrate the nature of contamination in η space. At every frequency channel, three separate PSFs are generated using natural weighting, uniform weighting, and uniform weighting with a Gaussian taper (hereinafter referred to simply as Gaussian weighting). The simulated sky in lm is then multiplied with the PSF at zenith and the result is summed. This accumulated value at every frequency channel gives the dirty image flux density at the zenith pixel. This zenith pixel flux versus frequency is used to evaluate the spillover of foreground contamination by Fourier transforming to η domain. The simulation was performed for three cases each involving a different gridding and imaging algorithm; we discuss these below.
Case 1 employed independent gridding convolution at every frequency channel, and the dirty image was formed using an efficient FFT algorithm. Figure 7 shows the flux density versus frequency and contamination power in η for this case. We have shown earlier that for η < x max /c in R 1 , the exact structure of spillover depends on the nature of relative weightings given to visibilities. Weighting here includes the effect of uv plane distribution of baselines, user defined weighting and any tapers applied. For the case of uniformly weighted visibilities with complete uv coverage, the η spillover from a point source is expected to be localized since the convolution in Equation 11 has a significant contribution only due to the discontinuity in the weighting function at |x| = x max . In case of MWA snapshot visibilities, the weighting function is not 'smooth' for |x| < x max and the convolution in Equation 11 has a significant contribution from all |x| < x max . Though this results in a distributed spillover, a significant part of the spillover energy is still restricted to R 1 , and R 2 will, in the absence of baseline migrations, be contaminated only due to the sidelobes of the convolving kernel from Equation 11 as discussed earlier. This is the 'PSF contamination' that is expected to persist in R 2 . Any additional channel to channel variations in the PSF will contaminate R 2 beyond the 'PSF contamination'. Such channel to channel variations are evident in the upper panel of Figure 7 in the form of 'jitter', especially for the cases of uniform weighting and Gaussian weighting. This 'jitter' was found to be due to inter-channel baseline migrations as discussed in Section 4. To quantitatively estimate the amount of 'gridding contamination' due to baseline migration, the Blackman-Nuttall window with very low sidelobes was used as a frequency window function B(f ). As described in Section 3.2, such a windowing function will significantly suppress 'PSF contamination' in R 2 (ηc > 500 meter), and the dominant residual contamination in R 2 will be due to the inter-channel 'jitter' arising from the gridding process.
It is interesting to note that baseline migrations have induced relatively higher jitter for the cases of uniform weighting and Gaussian weighting as opposed to natural weighting. As discussed in Section 4.1 this is due to two reasons. First, for an inter-channel frequency spacing of δf , a baseline vector x (in meter units) suffers an inter-channel migration on the uv plane of length xδf /c (in wavelength units). Second, and more importantly, a relatively higher amount of inter-channel 'jitter' due to baseline migrations arises from regions of lower uv density (longer baseline vectors) compared to regions of higher uv density (smaller baseline vectors). Since uniform weighting, and in this case Gaussian weighting, place higher emphasis on longer baseline vectors as compared to natural weighting, the amount of 'jitter' in Figure 7 goes in decreasing order of magnitude from uniform weighting to Gaussian weighting to natural weighting. Earth rotation synthesis and drift scan-synthesis will increase the uv density and is expected to ameliorate the 'jitter'.
Case 2 involves image formation with the CZT algorithm with a similar image resolution at every frequency. As described in Section 4.4, image formation with the CZT algorithm is expected to lower 'gridding contamination' in R 2 arising from inter-channel jitter due to 'baseline migrations'. Figure 8 shows the line of sight contamination for the case of image formation using the CZT algorithm. The 'jitter' in flux density versus frequency (upper plot) has significantly reduced due to the absence of baseline migrations inherent to independent gridding at every frequency channel-an expected outcome of gridding the visibilities in meter units. Similar image resolution at every frequency channel was achieved by scaling the Gaussian taper to have the same form in wavelength units, and restricting the longest baseline vector that was used in the gridding and imaging process to the same value (in wavelength units) at every frequency channel. This resulted in inter-channel 'steps' in the weighting function at the far end of the uv plane that resulted in a small but significant amount of 'jitter' and hence a contamination in R 2 . Nevertheless, the contamination in R 2 for this case is reduced by a factor of 1.7, 6.4, and 21.7 for natural weighting, uniform weighting, and Gaussian weighting respectively. Relative contamination for natural weighting has not improved significantly as regions in the uv plane that contribute most to 'jitter' are down-weighted by natural weighting.
Case 3 involves image formation with the CZT algorithm with a frequency dependent image resolution. In this case, the set of baselines that are used in the gridding and imaging routine are identical at all frequencies. It may be noted that despite a frequency dependent image resolution, the CZT algorithm computes the sky flux density on the same lm grid at all frequencies. The results of this simulation are shown in Figure 9. The contamination in R 2 has dropped considerably owing to the elimination of channel to channel 'jitter' in the PSF by this method. The extremely low sidelobes of the Blackman-Nuttall window and the absence of 'baseline migration' in the CZT based imaging algorithm we have proposed herein result in extremely low levels of contamination in R 2 . Quantitatively, the contamination in R 2 for this simulation as compared to Case 1 has reduced by more than three orders of magnitude for all types of weighting schemes.
It is important to note that the amount of interchannel 'jitter' is largely dependent on the visibility distribution and in particular on the filling fraction in different parts of the uv plane. Therefore, the relative merits of the three different weighting schemes discussed will depend on the array configuration and observing strategy. It is however expected that the contamination in R 2 will be orders of magnitude lower if gridding and imaging are done using the proposed CZT algorithm along with a good frequency window. The simulations and results in this section are not exhaustive, and are included as a means to appreciate the various factors that influence the frequency dependence of the PSF and illustrate the potential reduction in foreground contamination.

Conclusions and future work
Redshifted 21 cm tomography using Fourier synthesis arrays has emerged as a promising tool for EoR studies. Mitigating the confusion effects of Galactic and Extragalactic foregrounds is a major challenge to such EoR experiments. Through this work we have furthered understanding of contamination arising from Extragalactic continuum foreground and instrumental effects using a common framework. In particular, we have derived analytical expressions for contamination along frequency, or equivalently line of sight (LOS) dimension, which relates the frequency-axis contamination in the image cube to the visibility weighting in Fourier space. In doing so, we have cast the problem in Fourier space where the telescope measurements indeed lie, thus enabling us to draw conclusions that are directly applicable to instrument design and data processing.
We identified two major sources of foreground contamination-'PSF contamination' and 'gridding contamination', and described their structure in LOS wavenumber space. The 'PSF contamination' due to any point source is localized in LOS wavenumber space around η = x max l/c where l is the direction cosine of the source with respect to the imaged pixel, and x max is the maximum baseline length in meters. Consequently, when imaging with continuum confusion that is offset in direction cosine upto a maximum of l max from an imaged pixel, most of the 'PSF contamination' will be confined to a regime R 1 defined by η ∈ [0 x max l max /c]. It is useful to note that there does exist a regime R 2 defined by η ∈ [x max l max /c 1/∆f ], where ∆f is the frequency resolution, where finite 'PSF contamination' persists. This contamination in R 2 may be suppressed by three orders of magnitude by judicious choice of a window function in frequency. For this reason, R 2 is an 'EoR window' where we may focus our efforts to detect EoR.
We have shown that gridding and imaging routines result in additional 'gridding contamination' in LOS wavenumber space. In particular, we have shown that independent gridding at every frequency channel may lead to a stochastic inter-channel 'jitter' in the PSF that potentially contaminates the 'EoR window' R 2 . To ameliorate problems associated with inter-channel 'jitter' due to gridding, we proposed an alternative gridding and imaging algorithm based on the Chirp Z Transform (CZT). In this algorithm the visibilities are gridded in meter units rather than wavelength units, thereby eliminating 'gridding contamination'. The CZT compensates for the stretching of baselines (in wavelength units) with frequency by introducing a 'chirp' term-a scaling factor by which the frequency of the Fourier sinusoids are multiplied. We finally demonstrated the localization 'PSF contamination' and the elimination of 'gridding contamination' in the CZT based imaging algorithm using simulations of imaging with MWA.
Our ongoing and future work is directed towards using the concepts and ideas in this paper in an all sky simulation of MWA images that incorporates different gridding convolution kernels, Earth rotation and driftscan synthesis, effects of array non-coplanarity, and the Galactic synchrotron emission. We are currently developing specific simulation pipelines that will accomplish the all sky simulations with the CZT based imaging algorithm. Such simulations will quantitatively estimate not only the merit of the CZT based imaging algorithm for EoR detection, but also the expected contamination from foreground continuum in future MWA data products. In generating these representations of foreground contamination, we populate sources in the sky out to |l| max = 0.5. The zenith pixel flux is plotted against frequency (upper panel) and η (lower panel). Notice the jitter in the upper panel for the case of uniform weighting and Gaussian weighting as opposed to the case of natural weighting. This jitter results in significant 'gridding contamination' at higher values of η as seen in the lower panel. Also shown in the lower panel is radial uv density as a function of |x||l| max for the case of natural, uniform and Gaussian weighting. The boundary between R 1 and R 2 lies at ηc = |x| max |l| max = 500 meter.  Figure 7 has significantly reduced in the zenith pixel flux density versus frequency plot (upper panel). Consequently, foreground contamination at higher values of η (lower plot) is significantly lower as compared to that in Figure 7. A finite amount of undesirable 'jitter' still persists since a different set of baselines vectors were used in the gridding and imaging process at every frequency channel so as to keep the final image resolution constant at all frequencies.