{'fraction_non_alphanumeric': 0.041110895001381965, 'fraction_numerical': 0.01653239427968921, 'mean_word_length': 4.175248993430811, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 129, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Spectral estimation for spatial point patterns

This article determines how to implement spatial spectral analysis of point processes (in two dimensions or more), by establishing the moments of raw spectral summaries of point processes. We establish the first moments of raw direct spectral estimates such as the discrete Fourier transform of a point pattern. These have a number of surprising features that departs from the properties of raw spectral estimates of random fields and time series. As for random fields, the special case of isotropic processes warrants special attention, which we discuss. For time series and random fields white noise plays a special role, mirrored by the Poisson processes in the case of the point process. For random fields bilinear estimators are prevalent in spectral analysis. We discuss how to smooth any bilinear spectral estimator for a point process. We also determine how to taper this bilinear spectral estimator, how to calculate the periodogram, sample the wavenumbers and discuss the correlation of the periodogram. In parts this corresponds to recommending suitable separable as well as isotropic tapers in d dimensions. This, in aggregation, establishes the foundations for spectral analysis of point processes.


Introduction
On one hand, point processes are a very important special form of stochastic processes. They have been the focus of numerous books, and are applied to model patterns in criminology, ecology and earth science, just to list a few concrete examples. On the other hand spectral/harmonic analysis of stochastic processes have led to making key insights into the nature of spatial and temporal processes. Despite the long development of spectral analysis in time series, hearkening back to the 1890's, their usage is far from standard in the spatial point process literature, in fact the lack of methodological development has virtually made application non-existent. This paper aims to develop the fundamentals of spectral analysis methodology for stationary point processes. Many of the basic choices for time series and random fields are not well developed for point processes, which we answer head on.
If one only observed a stochastic process once over a given domain, then additional assumptions are required to enable estimation, basically to enable averaging within rather than across ensembles. The standard assumption that enables averaging for estimation, is that of spatial homogeneity or stationarity. We shall also make this assumption throughout this article. Without this assumption any spectral decomposition needs additional care in interpretation, and as many questions remain outstanding even in this case, we confine our interest to stationary processes. In addition, it is possible to assume isotropy, so that dependence becomes characterised in one variable. We shall discuss isotropy, but not confine our interest to solely this special case.
The spectral decomposition of a random process on one hand gives a decomposition of variance of a stochastic process into a scale-by-scale contribution of energy; on the other hand, it gives the representation of the process as a weighted average of complex exponentials so that given behaviour can be associated with a given physical scale (the inverse of its wavenumber or magnitude wavenumber). The latter notion of scale is key to many further developments in the representation of a particular stochastic process. For point processes the corresponding results are harder to derive as our observations do not take Euclidean values at fixed regular locations but are either absent or present at that given (but random) location. This means that the random location enters non-linearly in the Fourier transform, as the argument of the complex exponential, rather than as in a linear combination.
Spectral analysis of regularly sampled processes has already been studied extensively in the 1-D case (Percival and Walden, 1993b), but even as we study results into 2-D understanding the corresponding developments become more challenging, especially as non-parametric methods can become hard to interpret. For random fields (Stein, 2012;Adler and Taylor, 2010;Rosenblatt, 2012;Deb et al., 2017;Guillaumin et al., 2019) spectral analysis does resemble the analysis of time series, but the additional challenges of embedded simplified structure (lower dimensional structure) in 2-D and higher D are harder to address, and the appropriate visualization remains an outstanding challenge, see for example (Errico, 1985;Durran et al., 2017;Osgood, 2019). Also, there is some study of random fields sampled by 2-D point processes, as for example by (Chou and Lii, 1992), in 1-D (Lii and Masry, 1994).
Despite the potential promise of spectral analysis of point processes, after the initial flurry of development (Bartlett, 1963(Bartlett, , 1964Mugglestone and Renshaw, 1996b), also summarized in (Daley and Vere-Jones, 2003, Ch. 8), some results followed for its practical utility (Mugglestone, 1990;Mugglestone andRenshaw, 1996b, 1998), but its promise was never realised into broad usage. Somehow spectral analysis was never converted into a well-used practical tool set for data analysis, and normally it is referred to only when modelling and analysis depend on it (Lavancier et al., 2015;Møller and Torrisi, 2007). We suspect the reason why is, at least partly, that a number of practical choices necessary for constructing spectral summaries have never been substantiated by the proponents of spectral analysis from a theoretical or methodological point-of-view, and therefore remain somewhat mysterious and unclear. Naturally in the limit of low wavenumbers spectral analysis of point processes must coincide what with spectral analysis of random fields to be self-consistent, but despite this fact, even at low wavenumbers, there are many unexpected features as we shall show. We shall assume that the estimation of the spectrum is the end point of the analysis of the point process.
We can simply contrast this state-of-the-art with the results available for random fields (Rosenblatt, 2012;Adler and Taylor, 2010;Stein, 2012). The usage of spectral methods for point processes is limited in parts just because the Cooley-Tukey FFT algorithm (Cooley and Tukey, 1965) cannot be used when performing analysis when by definition, the spacing of observations is not regularly spaced.
The aim of this paper is to understand the properties of any direct spectral estimator defined from a spatial point process, assuming we observe a spatially stationary process. The spectral density captures all relevant second order features of the process, and so this is a fundamental characterization of the process in question. Taking inspiration from standard practice in signal processing, we shall focus on the class of bilinear spectral estimators (Cohen, 1995) that form a quadratic either in the Direct Fourier Transform (DFT) of the data or the data directly, or a linearly modified version, bias-correcting or such as using tapering (Walden, 2000).
The well-known periodogram non-parametric spectral estimator is naturally a special (and important) case of a bilinear spectral estimator. Two important tools from time series, of de-trending (Percival and Walden, 1993b;Velasco and Robinson, 2000) and tapering (Thomson, 1982) will be transferred to point processes, and they will have a unique form in this context. It goes without saying, that the point process setting is much more involved than that of time series, and we introduce appropriate separable and radial tapers, drawing on our expertise from signal processing.
To understand the spectral analysis of point processes therefore we will address the following key points: 1. How do we define spectral summaries for a point process? We recall existing frameworks for spectral analysis of a point process in Section 2.
All these points will be clarified in order to enable spectral estimation, and determine the mean square error of a spectral estimator.
Note that the Direct Fourier Transform (DFT) of the point process is not zero-mean under the assumption of stationarity, unlike for random fields, as has been noted by (Diggle et al., 1987), and discussed more abstractly in terms of random measures by (Daley and Vere-Jones, 2003). We refer to work by (Biscio and Waagepetersen, 2019) and deduce Gaussianity of the DFT. We show that the DFT of a point process will not be zero-mean, no matter its marginal distribution (compare with discussion in Diggle et al. 1987 for small wavenumbers). To be concrete, we shall determine the first and second moments of the DFT, which, combined with the Gaussian limiting results, will completely specify the distribution (under growing-domain asymptotics, compare with Stein (1995) for fixed-domain asymptotics).
The Fourier Transform (FT) will be (modulus) squared to yield the periodogram, which will have asymptotically a non-central chi-square distribution, but after de-trending, this will be a chi-square just like the modulus square of the DFT of a random field. This shows that we cannot simply assume that the properties of DFT of a random field will be directly inherited by a FT of a point process, and shows the challenges in using Diggle's (Diggle et al., 1987) isotropic estimator, as the problem of its sampling in wavenumber is much harder to treat.
Tapering (Walden, 2000) is the second method that we shall explore more fulsomely. It is not clear how to taper in the point process context, or what taper to use, or indeed how to use them. Tapers are normally defined for regularly spaced data. We shall use more rarely applied continuous-space tapers, and show for what classes of processes, multi-tapering produces reduced variance estimators. We note that other choices of implementing tapering could have been made; and this will be discussed, such as tapers for irregularly spaced stochastic processes. We shall finish up with an extensive simulation study to illustrate the performance of our proposed estimators.
Therefore, we aim to settle many basic outstanding issues in spectral analysis of point processes. We also connect to related work on sampling a continuously indexed random field, see (Chou and Lii, 1992) for example (already developed fully for 1-D processes by Lii and Masry 1994). Our end result is the arrival at an appropriate smoothing method for the raw spectral estimator of a point process, and also suitable 1-D summaries of such a process, whether the process is isotropic or not. Section 7 discusses isotropic processes and isotropic summaries. This all combines to underpin our understanding of the spectral representation of point processes. Section 8 explores the implications of our results via simulations. We conclude the paper in Section 9 with a discussion of the derived results.

Notation and Definitions
We assume that we observe a spatially homogeneous point process X with intensity λ = ρ (1) > 0 and with ρ (2) (z) defined as the (second order) product density of X. This means that ρ (2) (z) only has one argument (z ∈ R d , namely the spatial shift) rather than depending on two local spatial variables (say x and y rather than z = x − y). We assume that X is a simple point process on R d , meaning no duplicate points are allowed. We take d = 2, 3, . . . , and so exclude the case of d = 1, which is relatively well-studied, see references in e.g. (Daley, 1971;Jowett et al., 1972;Vere-Jones, 1974). For a (Borel) set A in R d , we define |A| for the volume of A and N (A) for the number of points of X ⊂ R d in A. We denote by B any subset of R d where points are observed (or equivalently our observation domain).
In complete analogue with a random field we shall define the spectral density of X as the Fourier transform of the complete covariance function of X. This is the standard approach, and was first proposed by Bartlett (1964). The complete covariance function of a stationary point process is therefore where δ(·) is Dirac's delta function. The spectral density function (sdf) (Bartlett, 1964) of the point process X is then defined as the Fourier transform of the complete covariance function: This representation is in direct analogy with the corresponding spectral decomposition of a random field or a time series. The symbol F denotes the Fourier transform, with i as the imaginary unit. We refer to the argument w as the "wave number" rather than frequency to knowledge its multi-dimensionality. This is a natural choice of a Fourier transform in analogy to the analysis of random fields (Stein, 2012). To avoid additional constants in the inverse Fourier transforms, we parameterise the Fourier transform with wave number instead of the customary angular frequency ξ = 2πw used by Bartlett (1964) and some others. For a time series or random field there are a number of spectral results, ranging from the spectral representation theorem (I) (Percival and Walden, 1993a, pp.130) and to Bochner's (Herglotz) theorem (II) (Brockwell and Davis, 1991;Percival and Walden, 1993a). Both these results are not exactly available in the point process setting, but for a time series X t we can note them: Equation (3) decomposes X t into random contributions associated with each angular frequency ω. The covariance in Eqn (4) is also decomposed into weighted frequency contribution. Both of these observations are important for the interpretation of the Fourier representation of X t . "Important" contributions correspond to f (ω) being considerably larger relative to other f (ω) as in that scenario Var{dZ(ω)} is bigger than Var{dZ(ω )} and therefore |dZ(ω)| likelier to be larger than |dZ(ω )|. Spectral representations of point processes are also discussed in (Daley and Vere-Jones, 2003, Ch. 8).
We can yet again decompose the (complete) covariance in terms of a spectral representation like in (II), and it takes the form in terms of a d-dimensional Fourier transform of f (w) = e −2πiw·z γ(z) dz = λ + e −2πiw·z ρ (2) (z) − λ 2 dz.
We refer to f (w) as the spectrum of X. We note that the Fourier transform can be inverted to yield the equality of (an analogy of III): We assume the spectrum is the primary object of interest in our study of point patterns, as the covariance γ(z) can be fully determined from it, and the covariance specifies the point process. We see directly from (5) that ρ (2) (z) − λ 2 and not the complete covariance γ(z) is playing the role of a time series covariance. This is as when transforming γ(z) rather than ρ (2) (z) − λ 2 we otherwise get a constant boost of λ that spreads across all frequencies. This does not indicate to us which frequencies are "more" important unlike for a random field. By important we here mean contribute to the process more than any other wavenumber.
We have to exercise some caution when interpreting the Fourier transform as a bijective transform. Yes, it should contain the same information about scale, but the meaning of the word scale will be different. In time or space the notion is associated with the support of γ(z) or ρ (2) (z). In frequency being supported over a wavenumber w means variation over scales 1/ w is present in the correlation function, or to represent the variability in γ(z) we need frequencies w . The function ρ (2) (z) − λ 2 is often approached in terms of what scales it is non-zero at, but the variation in the function can be associated with long or short scales. So for a covariance function, we now have two notions of what it means to possess scales w ; either γ(1/ w u) is non-zero for unit vector u or γ(z) is variable (changing) and the scale of the change is w . Imagine observing a sinusoid with period k. The function will hit unity at a regular set of z values, and so it will be supported at those z values, but also its Fourier transform will be supported at 1/k. MaternII r2 Thomas FL Thomas MS Figure 1: Pair correlation pcf = ρ (2) /λ 2 (pcf, left) and corresponding scaled spectral density function f /λ (sdf/λ, right) for two stationary and isotropic non-Poisson point process models, two variants each, with known pcfs. For the Thomas processes the sdf is given in the text. For the regular Matern II process the sdf was numerically approximated using the Hankell transform (cf. Section 7). See details of the processes and the variants in Section 8. This is clear from Figure 1 where covariance is present at smaller or medium scales smoothly for the clustered Thomas processes, or repelled for the Matérn hard-core processes. As the Matérn process's complete covariance function is discontinuous, its Fourier transform is supported over all scales (due to the Heisenberg-Gabor uncertainty principle (Cohen, 1995)). This can be deduced from (6). To reproduce the discontinuity we need all scales.
We see from (5) that not all frequencies or wavenumbers are given equal weighting. First all wavenumbers are given an equal weighting λ and then the Fourier transform of ρ (2) (z) − λ 2 determines the wavenumbers that are up-weighted (added) or down-weighted (subtracted) relative to the overall level of λ. We then observe what wavenumbers are important to the representation of the point process, which gives more information, conveniently decomposed on a scale-by-scale manner.
For a random field or a stochastic process in d-dimensions a few cartoon characteristics are important. For a discretely sampled process in Z d the simplest random process, white noise, is constant across wavenumbers yielding a spectrum that takes the form of σ 2 on [− π ∆ , π ∆ ] d , where ∆ > 0 is the sampling period, and zero otherwise. For a point process the choice of definition of the spectral density does not imply a decay because of the inclusion of the term λδ(z) in space. However once we remove this term, we expect a decay of the remaining spectrum f (w) − λ as |w| → ∞. Otherwise (5) would contain additional singularities.
Furthermore at w = 0 we retain The second term in this expansion is not required to be zero, and for the Thomas process, it is not zero. For a Thomas process the spectrum is (Mugglestone and Renshaw, 1996b, p 55) It is clear from this expression that we cannot arrive at a constant contribution due to the exponential. As we shall be studying the moments of a point process, it is convenient to restate Campbell's theorem (Chiu et al., 2013, Sec. 4.3.3), and we shall use this result multiple times. The theorem applies to any measurable function h : R nd → R + with n = 1, 2, ..., and states that where the summation is over distinct point tuples. Note that if the point pattern X is stationary then ρ (1) (x) = λ is constant for any x, and ρ (n) (x 1 , ..., x n ) = ρ (n) (o, x 2 − x 1 , ..., x n − x 1 ) ≡ρ (n) (x 2 − x 1 , ..., x n − x 1 ). We shall use theρ (n) notation to explicitly avoid confusion when parameterizing in terms of the shift z 1 = x 2 −x 1 , . . . , z n−1 = x n − x 1 (and the stationary parameterization) rather than the basic coordinates. We define a new parameterisation to capture how the process's f (w) differs from that of a Poisson process. We define the nth deviation as For a Poisson process these n = 2, 3, . . . deviations from Poissonianity will all be identically zero. By understanding deviations from Poisson behaviour, we get greater understanding of the underlying process under study. For completeness we also define the Fourier transform of the deviations: where if n = 2 we suppress the superscript. With this definition we find that Thus in a sense, f (w) encapsulates the deviation of the function from a constant spectral density of λ via the term λ 2 f (w). Daley and Vere-Jones (2003, p. 337) discusses some differences in spectral representation of time series versus that of point processes. We would argue that decay of the spectrum is still reasonable to assume once λ has been subtracted (so the decay of f (w) is reasonable to assume). Having established the theoretical spectral description of point processes, we now turn to their sampling properties.

Direct Spectral Summaries of the Data
We shall re-examine spectral estimation of point processes, and clarify a number of points that are currently unclear. We start from Bartlett's periodogram estimator based on observing a point pattern in region B ⊂ R d which can be written as: where we have setλ = |B| −1 |X ∩ B|. This definition uses all the data available to us, {x : x ∈ X ∩ B}, but is simply a possible choice amongst all possible direct spectral estimators. In general a direct spectral estimator is one which is bilinear in the DFT of the data, see for example the discussion in (Percival and Walden, 1993a, Ch 5-6). We recall a bilinear form is a function D(x, y) satisfying D(x 1 + x 2 , y) = D(x 1 , y) + D(x 2 , y), and equivalently in the second argument. Bilinear forms have been thoroughly discussed in signal processing for the analysis of time series (Loughlin et al., 1993). As a point process consists of locations the best we can hope for in terms of bilinearity will be in terms of sesquilinearity of the Fourier transform as the point locations will appear in the argument of the complex exponential. Sesquilinearity of D(x, y) simply generalizes a bilinear form to the Hermitian symmetry, additionally requiring D(x, y) = D(y * , x) as well as D(x, y) satisfying D(x 1 + x 2 , y) = D(x 1 , y) + D(x 2 , y).
In time series the bilinear form has been chosen to ensure that spectral estimators are real-valued, and often non-negative though some bilinear estimators are not, see e.g. the Kirkwood-Rihaczek spectrum (Hindberg et al., 2006), and Guyon's spectral estimator (Guyon, 1982). Spectral estimators in time series are bilinear in the observed real-valued process X and sesquilinear in its Fourier transform J h (w). We shall still require that the form of the estimator is sesquilinear in the DFT of the point process, but the estimators will not be bilinear in X, the point pattern, as this will not be possible. We define the tapered DFT of X of a point pattern for a specified general (square) integrable function h(x) (referred to as a 'data taper' by Percival and Walden 1993b) with Fourier transform H(w), to be for any unit norm h, e.g. h 2 = 1. Because x appears in the argument of the complex exponential, we cannot define estimator bilinear in the data. In practice one has to decide what taper function to use. In 1D for point processes tapering has been used (Cohen, 2014;Taleb and Cohen, 2018) with continuous wavelets, but also multitapering has been used on interpolated data (Das and Babadi, 2020). Here, we do not interpolate, do not implement localized analysis, and do not implement analysis in 1D. Higher dimension is more complex to interpolate, as it is not orderly.
We note that in general h(x) exp (−2πiw · x) is identically distributed but not independent for many choices of point processes x, making the choice of a Central Limit Theorem (CLT) a bit more complex. For a large class of processes one such CLT is given by Biscio and Waagepetersen (2019). If we compare the quantity in (13) with (Biscio and Waagepetersen, 2019), then we see that q = 1 and p = 2, in their notation. We do not have to worry about e −2πiw·x being a function of several of the points, and this is why q = 1. Requirement H1 is therefore trivially not a problem as an assumption, H2 is about the scaling of the observation area, as we assume the sides scale the same not a problem, H3 (higher order centred powers) are limited through the usage of Campbell's theorem and as q = 1 we do not have to worry about pathological scaling between components in H4. Therefore, citing (Biscio and Waagepetersen, 2019), we can deduce from their Theorem 1 that as |B| diverges, and recover the second moments by direct calculations in Section 4: Here |B| → ∞ in equal speed in all directions of B and B is the cuboid, and we define the tapered related term H(w) in (21). We deduce the moments from Proposition 5.1 combined with Theorem 5.1. We can view J h (w) as a complex-valued scalar or a real-valued two-vector. N C (µ, Σ, C) is the general complex normal (Schreier and Scharf, 2003) and its arguments are its mean µ, it covariance Σ, as well as its complimentary covariance C. It should be contrasted with the complex propert normal N C (µ, Σ), that has zero relation. Based on J h (w) the spectral estimator becomes its modulus square or If we take h(z) = h 0 (z) = 1 √ |B| 1(z ∈ B) then we recover Bartlett's peridogram in (12) from (16). The point of using a general function h(z) is that abruptly ending the inclusion of points when we leave the region B causes ripples in wavenumber, and therefore a worse estimation of the spectrum, as is also the case of time series when not using a taper (Percival and Walden, 1993b). We have deduced an asymptotic distribution for J h (w) using (Biscio and Waagepetersen, 2019) (14) that if we consider only arguments w so that the complementary covariance is negligible then where the parameters of the non-central χ 2 k (m 2 ) (k denoting the degrees of freedom, and m 2 the non-centrality parameter). We give H(w) in (21). We can also use the uniform integrability of the variable J h (w), to get the asymptotic moments of I 0 (w) from these limits.
Note that the sum in (16) cannot be evaluated in a computationally efficient manner (compare to (Cooley and Tukey, 1965)), as the locations {x} are not regularly spaced, which is unfortunate, but unavoidable. Finally (16) is bilinear in J h (w) but it is not bilinear in X, unlike the commensurate expressions for the DFT of time series and random fields.
For completeness we also note the isotropic estimator of (Diggle et al., 1987, Eqn. 3.3) in 2D, namelȳ where J 0 (x) denotes the Bessel function of order 0, specified in Section 7. The d-dimensional extension is given later in equation (24).Ī 0 is less clearly bilinear in the data, but if we start from the estimator of (12) and average it over orientations analytically, then we arrive at this form. We therefore abuse our terminology to also refer to it as a bilinear estimator. There is one additional correction made by (Diggle et al., 1987). As we shall see, for low wave numbers there is a bias inherent in (18). To address this problem (Diggle et al., 1987, Eqn. 3.4) suggests taking for some choice of lower bound t 0 Diggle et al. (1987) suggest taking t 0 so thatĪ 0 (t 0 ) is at a minimum, and also suggest smoothing the estimated I D , which we clarify in the linear estimation section. The authors also discuss iterated methods of bias correction in their Section 5, and correcting biased estimators of the correlation. They discuss difficulties of preserving short range features in the spectrum. Note also that we have modified the estimator in (18) to divide by |B| rather then the observed number of points in the region, N (B), as do for example (Diggle et al., 1987), as it is much preferable not to have a random denominator. We shall discuss the usage of I 0 (·) versus I D (·) further on in Section 7, and the implications of when the process is isotropic or not.
As can be surmised from (16) the expectation of the general estimator E {I h (w)} is a convolution between the Fourier transform of the window and the true spectrum, and if h(z) does not go to zero nicely at the boundary of B then the periodogram becomes quite complicated in terms of its expectation. Looking at (16) also raises a second problem, namely that there is a constant contribution |B| −1 |X ∩ B| which does not give any wave number specific information, and is correlated between wave numbers. This term is new to time series and random fields, if not new to the expectation of the periodogram of randomly sampled stochastic processes, see (Lii and Masry, 1994;Chou and Lii, 1992).
We can generalize the family of estimators from Eqn. (16) to take the form of We study the family of estimators I h (w). We also note that (20) is a direct spectral estimator (Percival and Walden, 1993b, p. 207). Multidimensional tapers are available to us (Hanssen, 1997;Simons and Wang, 2011), and can be precomputed. Our choice of tapering corresponds to using continuous tapers evaluated at random locations. Should we use no data taper, or equivalently use the classical taper h 0 from (8), then we introduce the notation of writing I 0 (w) for the periodogram for clarity. We then also define The expectation of any member of that family takes the form of: This can be derived straightforwardly by direct computation from (20) by taking expectations with Campbell's formula and using the convolution theorem. Should we chose to taper isotropically in 2-D-space then instead we take with isotropic taper h I (·) Note that this cannot necessarily be constructed. Say J h () with an isotropic taper would be made with the largest spherical domain and a compact support could be constructed, but what does that imply for I h ()? Radial tapers can be determined as described in (Simons and Wang, 2011), whether in 2-D continuous space or discrete space. Equation (23) can be of course be extended into higher dimensions. In general we would have in dimension d = 1, 2, 3, 4, . . .
We can see directly that As the change of variable to polar coordinates means we will multiply by z d−1 , and this will correspond to a d/2 − 1 Hankel transform of the taper h(·) multiplied by the product density. For more details on the Hankel transform see (Grafakos and Teschl, 2013;Cressie, 2015). We can also use more than one taper, and will use a sequence of orthogonal functions {h k } (Percival and Walden, 1993b), that will be used to get a new estimator. We will assume they are all of unit energy and are all (pairwise) orthogonal. This means we will require ||h k || 2 = 1 and Tapering has been viewed with variable approbation. Criticism comes especially for throwing away data (Fougere, 1977;Yuen, 1979). We shall study the properties of direct spectral estimators of the form of Eqn (20). This will help us to characterise the second order properties of the process X. Just as a brief note: it is not obvious to know how to taper. Most of time series analysis is based on discrete regular sampling for instance. We chose to use continuous tapers, rather than interpolating the points to a regular grid. This leaves the spheroidal wavefunctions (continuous but hard to compute), as well as the cosine tapers of Riedel and Sidorenko (1995). These correspond to separable taper choices; non-separable choices with be discussed in Section 7.

Expectation of Bilinear Spectral Estimators
We shall start from the simplest spectral estimator, namely the Bartlett periodogram as specified by (12). We shall now be quite concrete and understand some common cases of spatial data, and their sampling, and choose as early special cases Cartesian product domains. We shall focus on the cuboid sampling domain: Lemma 4.1. Assume that X is a homogeneous point process with intensity λ and twice differentiable spectrum f (w) at all values of w = 0. Assume X is observed in a cuboid domain B with a centroid at 0, which is growing in every dimension. Then the expected value of the periodogram I 0 (w) satisfies the equation for This is not what we would expect as a large sample expectation of the spectrum, given our experience for time series and random fields in that the term |B| −1 λ 2 T (B, w) has been added. To get there, we need to understand the DFT J h (w) better.
Proposition 4.1. Assume that X is a stationary point process with intensity λ and with spectrum f (w), and that h is a unit energy taper supported in the domain B itself only. Then the first moments of the direct spectral estimator J h (w) are given by: Proof. See Appendix C and notice the definition of (21).
In the above Proposition the term Rel {·} denotes "relation", also known as the "complimentary covariance", see e.g. (Miller, 1974). These moments will have implications for the distributional properties of I h (w) = |J h (w)| 2 . As we consider w = 0 then the variance and the relation are the same. For w small if J h (w) is Gaussian but the mean is non-zero and the moments are different than from the moments of a zero-mean time series and likely to lead to a non-centred χ 2 distribution.
As the DFT is the sufficient statistic (asymptotically) of a stationary process, in a sense, determining its first two moments are pivotal. Perhaps most important is (27) because of the non-central nature of the distribution. Understanding such distrubutions, e.g. the non-central χ 2 of the square of a Rician distribution, is not uncommon in medical imaging (Olhede and Whitcher, 2011). Inference is often not facilitated subsequently to making this observation-instead we shall propose to remove the mean, removing the non-centrality parameter.
We thus define a bias-corrected periodogram (see Daley and Vere-Jones 2003, p. 292, where the signed measure is used to make the equivalent definition), or Thus J h (w) − λH(w) is often called the "signed measure". This is the DFT of the process X 0 dual to the mean-corrected random measure N 0 (dx) = N (dx) − λdx where N is the dual of X. As we have subtracted H(w) off the discrete Fourier transform, it is no longer strictly positive, but once we take the modulus square we are guaranteed to arrive at a real-valued positive quantity. There is one wavenumber which is problematic, namely w = 0. For the periodogram we have h 0 (x) = 1 B (x)/ |B|. Thus we have This removes the problem of a non-zero mean of the DFT, as long as we assume that we know the intensity λ, and so are in the position to remove this effect. Assuming knowledge of this quantity is not a major hurdle, as it can be estimated consistently for largish areas (|B| → ∞), by just counting the number of points and dividing by the area. For completeness we also define the signed measure DFT as trivially. Thus we cannot estimate the DFT at wave number zero. For a time series analysis when calculating DFTs we subtract off the sample mean, and then also the mean corrected DFT is zero at wave number zero. In a time series setting the periodogram at frequency zero is often not plotted.
We note that two ways of bias subtraction for the periodogram is possible, namely I(w) as described by (30) and additionally we could have defined The advantage of the first estimator is that it is naturally non-negative and removes the error before we take the modulus square. This also should have a positive impact on the variance. In addition, λ 2 can be a biased estimator for λ 2 . If we are looking at a more sophisticated estimator than the separable estimator of (34) or (30) see e.g. (Cohen, 1995, Ch. 6), then we have where g(x, y) is a kernel which may have a specified support, that has to be selected. For example we could define g(x, y) = h(x)h * (y) with h(x) a multi-dimensional data taper (Hanssen, 1997), or we could make a non-separable choice of g(x, y) = h( x − y ). Depending on the choice of g(x, y) the quantity I g (w) may satisfy a number of desirable characteristics. For example positivity, asymptotic unbiasedness as an estimator, and computational efficiency. The quantity I h (w) will be referred to as the "debiased periodogram", (and the subscript unity will be removed) and is preferable to use to the normal periodogram as it removes a number of unpalatable effects.
Lemma 4.2. Assume that X is a homogeneous point process with a spectral density f (w). Then the biascorrected tapered periodogram has a first moment given by: Proof. See Appendix D.
Thus as along as H(w) is getting more concentrated in wave number (e.g. H(w) → δ(w)), this is asymptotically in |B| an unbiased estimator of f (w). Using the signed measure DFT of (32) to study X is more convenient, as this lets us avoid the contribution that affects the low wave numbers.
We now look to determine further properties of I h (w), or properties of I h (w). To produce an estimator that is smoothed we need to study the covariance and variance of I h (w 1 ) and I h (w 2 ). We first look at E I h (w) . We find that its form for large spatial regions is specified by the following theorem. Let us write down what grid we get large sample unbiasesness for, customarily called the Fourier grid, and this is useful when the observational domain B is a box.
Note that the physical units of the wave number elements is per unit length, such as m −1 . Referring to Lemma 4.1 we see that using this grid both removes any asymptotic bias at low wave numbers, but also gives us a sampling of the wave numbers that approximately leads to independent periodogram ordinates, rather like in the random field case. Note that the centering of the box B at zero, e.g. requiring the box to take the form is required for the stated result. If the observational domain is not centred, then additional phase factors enter the expression. We can compare these results to wave number grid we chose to evaluate a standard spectral estimator of a time series, the Fourier frequencies (Percival and Walden, 1993b, p. 197-198). Their basic importance follows because we can expect the direct Fourier transform to be Gaussian and so uncorrelated implies independence.
A second feature of time series is the notion of the Nyquist wave number (Percival and Walden, 1993b, p. 98). This does not exist for point processes. It may seem counter intuitive that there is no upper limit to the wave numbers we can estimate. When analysing a process that has been sampled in space, such as a random field, we expect to see aliasing. Aliasing is when variation that is happening very rapidly is confused with slower cycles, because when we have sampled more sparsely, rapid variation cannot be resolved. Because the distances can be any real-valued value, the Nyquist wave number does not exist.
Theorem 4.1. Large-Domain Expectation of the tapered periodogram. Assume that X is a stationary point process with intensity λ and twice differentiable spectrum f (w) at all values of w = 0, and that h is a unit energy taper (e.g. h = 1) supported in the cuboid domain B only. Assume X is observed in B , which is growing in every dimension, that is min j → ∞. Then the expected value of the tapered periodogram I h (w) satisfies the equation where H(w) is the Fourier transform of h(x), as defined by Eqn. (21).
Proof. See Appendix E.
In fact, asymptotics are well-understood both for time series and random fields. Our understanding of asymptotics for spectral estimation of point processes falls back to Daley and Vere-Jones (2003) or Diggle et al. (1987). In the latter reference Diggle et al. (1987) refers back to Chemistry theory by Hansen and McDonald, discussing the properties of gasses.
Corollary 4.1. Assume that X is a stationary point process with intensity λ and twice differentiable spectrum f (w) at all values of w = 0, and that h is a unit energy taper supported in the domain B only. Assume B is growing in every dimension. Then the expected value of the tapered periodogram I h (w) satisfies the equation where H(w) is the Fourier transform of h(x).
Proof. The proof can be found in Appendix G.
This result is not unexpected and generalizes Theorem 4.1. The main difference here is that we consider a general domain B rather than a cuboid domain B . This theorem allows us to determine the expectation of I h (w), in general, as is made clear in the following corollary.
Corollary 4.2. Bias-corrected periodogram. Assume that X is a stationary point process with intensity λ and twice differentiable spectrum f (w) at all values of w = 0, and that h is a unit energy taper supported in the cuboid domain B only. Assume X is observed in B , which is growing in every dimension. Then the expected value of the tapered periodogram I h (w) satisfies the equation This corollary shows the importance of removing the spectral bias before squaring-otherwise it gets hard to isolate and remove the effect at zero wave number-which is exacerbated at higher intensities (growing λ), as is clear from (38). So the estimator I h (w) is debiased relative to E {I h (w)} (as we shall see from Theorem 4.1), but is still guaranteed to be non-negative. This can be compared to removing µ from a time series before calculating the periodogram.
This gives us insight into the act of debiasing. In a sense our understanding of the spectral nature of point processes is informed by the spectral representation of randomly sampled processes (Chou and Lii, 1992). We can continue to explore the differences and similarities between spatial point processes and random fields.
Finally, to enable smoothing, define P ≥ 1 estimates of the spectrum via where we assume where δ pl = 1 if p = l and 0 otherwise. These estimates can also be arrived at starting from (13) and J h (w). As we shall use P tapers linearly, and P 2 for the quadratic estimators, we shall use the subscript '0' to refer to the untapered periodogram, and p = 1, . . . , P to refer to the subsequent tapers. Tapering is very common for stochastic processes. Initially the idea was hotly contested, but its utility is now firmly established both for time series and random fields (Dahlhaus and Künsch, 1987;Walden, 2000). The idea of tapering is to ameliorate edge effects that leads to the asymptotically leading bias term. This is important in 1-D (Thomson, 1982;Walden, 2000) but of greater importance in 2-D and higher (Dahlhaus and Künsch, 1987). In 2-D and higher the edge effects become more pronounced, and indeed asymptotically dominant (Kent and Mardia, 1996;Mardia and Marshall, 1984). This paper also uses multitapering to stabilize the variance of any estimator, as described for time series in (Walden, 2000). We shall describe the necessary steps to perform linear estimation of the spectrum of a 2-D and higher dimension point process.

Regular
Completely Rotation averaged periodogram (centered and scaled) Figure 2: Top row : Four example data patterns in a rectangle of area 2e4 units, exhibiting regularity, complete randomness, small scale clustering, and clustering with a left-right axis directionality. Second row : Their debiased periodograms excluding w = 0. Third row : Their multitapered debiased periodograms, 3 × 3 sinetapers. Bottom row : Their rotationally averaged periodograms, including the case without debiasing. The 2D summaries use a w-grid with a step size 0.005 in both directions. The rotational averaging was estimated at |w| = 0.0050, 0.0075, ..., 0.2500 using a box kernel having a radius 1.5 times the w-grid step size 0.005.
Let us set theory aside for a moment and have a look at some spectral estimates: Figure 2 demonstrates the debiased periodogram, the debiased multitapered periodogram with sine-tapers and a rotationally averaged 1D summary of the peridogram (Sections 6&7) for four point patterns exhibiting different structural behaviour. The aspect ratio of the plots are kept the same only to have an orderly figure: Due to debiasing we can estimate the periodogram on a wavenumber grid of our own choosing, and not just on the Fourier grid (which in these finite examples would be different for different tapers, cf. Section 5). The wavenumber-grid in this illustration is a regular grid with a step size 0.005 in both dimensions.
We see that information is present at wavenumbers up to |w| ≈ 0.2. The tapering smooths the periodogram, as expected. There is a dark well near the origin for the regular pattern, and a bright hump for the clustered pattern, and both features transfer to the rotationally averaged curves (compare with Figure 1). The anisotropy of the fourth pattern is hinted by the anisotropy visible in the 2D periodograms as elongation of the "hump" in the second dimension. The elongation is perpendicular to the elongations of the clusters in the data because wavenumbers relate inversely to spatial units. The rotationally averaged periodogram was also computed from the (not shown) non-debiased periodogram for comparison. There is a very prominent bias near |w| = 0 for the non-debiased version, which hopefully illustrates why the proposed debiasing step is relevant when applying the method in practice on freely chosen wavenumber grids.

Variance and Covariance of Bilinear Spectral Estimators
Before determining how to do bilinear estimation of the spectral density let us determine the second order properties of spectral estimators. Core to our understanding of smoothing will be the variance and covariance of the tapered functions. This will be established in the following theorem. For brevity define With this notation we have Above we note that we have the same value of p across x and y in (43) as this corresponds to a modulus square.
Theorem 5.1. Let I t p (k) denote the (tapered) periodogram given by (41). The covariance between the (tapered) periodogram and itself across two wavenumbers w 1 , w 2 ∈ R d then takes the form of Proof. The proof is in Appendix H.
Note that we can in fact write: The expression derived in Theorem 5.1 is quite general, and only requires the assumptions of homogeneity in space, but it is hard to understand, especially as it is hard to order the terms in size without putting very restrictive assumptions in place. We can still consider the untapered case by setting p = l and taking h(x) = |B| −1/2 1 B (x), to arrive at an expression. Let us study this general correlation in the instance of the Poisson process where we know that and consider the case of w 1 = w 2 = w. For completeness we also define the spectral bandwidth b h by Proposition 5.1. Let I t p (w) denote the tapered periodogram given by (41) for a Poisson process with intensity λ. The covariance between the tapered periodogram with one taper and the tapered periodogram with another taper then takes the form for k > b h of The form for w < b h is much more complex, and is provided in Appendix J.
Proof. The proof can be found in Appendix J.
Remark 5.1. At this stage we have simplified our assumptions to the Poisson process which is quite disappointing. The reason why we have chosen to do so is clear from Theorem 5.1. If we assume the process is Poisson then a lot of terms cancel. The correct treatment of non-Poisson processes would correspond to transforming the terms into a wavenumber representation, doing a Taylor series because of the concentration of the tapers, and then implementing the same calculation after the Taylor series.
Proposition 5.2. Let I t p (w) denote the tapered periodogram given by (41) for a Poisson process with intensity λ. The covariance between the (tapered) periodogram and itself at different wave numbers takes the form of Proof. The proof can be found in Appendix K.
Fix 0 < < 1. This is so that we can approximately find an uncorrelated grid. Recall that h(x) is unit energy in space or wave number.
Definition 5.1. The tapered Fourier wave number grid for a point process observed on a cuboid domain B (l) corresponds to the points where τ is defined from the equation These results hint at producing general methodology; we can for the Poisson process figure out how to do non-parametric estimation. This is however not enough. We need to determine the variance of the periodogram in more general. The following proposition establishes its form.
Theorem 5.2. Assume X is a stationary process. The variance of the (tapered) periodogram then takes the form of Proof. The proof can be found in Appendix L.
With these moments we can now consider the problem of estimation. We might feel an increasing degree of unease as we have succumbed to the usual fallacy of only deriving the most general results for a Poisson process. However as the observational domain becomes larger we would expect |H p (w)| to concentrate and so the lack of constancy will become less of an issue.

Linear Smoothing
Having derived the first and second order properties of the periodogram, the next task is then estimation of f (w), as the primal quantity of interest in studying the point process X. If we fully characterise f (w) then we have fully characterised γ(z), as the Fourier transform is a bijection. Just like for a time series or random field, a key question concerns the correct method for quadratic characterisation and estimation.
For the process itself (X), then its first moment is encapsulated by λ for a stationary (or homogeneous) processes. We recall the spectrum uniquely maps to the complete covariance function γ(z). We want to start from the random variables {I h (w)} (or { I h (w)}) and then non-parametrically estimate the spectral density directly from those quantities. As long as we seek to estimate f (w) at points w where the function f (w) is continuous, averaging {I 0 (w)} (or { I h (w)}) locally seems like a sensible strategy. (Diggle et al., 1987) proposed smoothing I D (w) as defined in (19). They just note that they apply a weighted average smoother with user tuned weights. This will not determine how to sample wavenumbers, or remove correlation, and we also have the alternative form of bias removal to them before smoothing.
We shall only treat growing domain scenarios (|B| → ∞ (one can contrast growing domain observations and increasing density of observations), see e.g. Zhang and Zimmerman (2005) or Stein (1995). Of course when faced with a real point pattern we cannot tell if asymptotic methods offer a good approximation to the finite sample distribution.
For sufficiently high wave numbers the constant effect in the spectral density with expectation λ will dominate both in the mean and variance of I 0 (w). Thus for higher wave numbers in w the periodogram will eventually become perfectly correlated over wavenumbers, and will be converging to λ. For a growing domain we can show thatλ converges in probability to λ.
Our understanding of the correlation of higher wave-numbers should be contrasted with the discussion of Mugglestone and Renshaw (1996a). The latter authors remark that there is an upper limit on the number of independent periodogram ordinates, see either of Mugglestone andRenshaw (2001, 1996b), but do not justify that statement (which is based on random field analysis rather than point processes). It is not clear from these manuscripts exactly how this upper limit is arrived at, as Mugglestone and Renshaw (2001) it is based on heuristic (or one could say practical) arguments on the nominal degrees of freedom present in the data, in the random field limit.
We have already described how to de-bias J h (w) into J h (w), which is a properly zero-mean quantity. We clearly do not need to do this as long as we stick to the Fourier grid. Mean correction has also been discussed by (Daley and Vere-Jones, 2003, p.332) to achieve analogous results to the spectral representation of time series. This mean correction has the advantage of letting us determine the properties of J h (w) directly, and as we chose to evaluate Bartlett's rather than Diggle's spectral estimator, we have the complex-valued quantity available to us, and can do the mean correction while still retaining a real-valued positive spectral estimate. Mean correction for the periodogram after squaring is more precarious, as we are not guaranteed a positive estimate.
Having understood the periodogram point-wise in terms of its expectation, to understand how to smooth it, we also need to make sense of its variance/covariance structure. Smoothing requires an explicit understanding of local degrees of freedom. For this purpose we shall use tapered estimates, see also (Walden, 2000). We refer to the previous section and (41) for the basic taper definitions.
We have two possible estimators, namely the bias-corrected multi-taper estimator of (Walden, 2000) as well as a multi-dimensional kernel density estimator (Duong and Hazelton, 2005). They are both linear in the estimated spectrum. We shall use the periodogram rather than another estimator for the kernel density estimation, as it is easier to keep track of bias and correlation.
We therefore define the multitaper estimator (Walden, 2000) (and refer to (Hanssen, 1997;Liu and Van Veen, 1992)) where I hp (w) is defined in (30). Choosing to use continuous tapers in (30) is not the only option, and this definition allows for nonseparable tapers (Simons and Wang, 2011), as it does not loop over Cartesian product tapers. We could have calculated tapers for our spatial samples, see (Wingham, 1992;Bronez, 1988), or we could have interpolated the points (Das and Babadi, 2020). Neither solution is palatable, and so we have stuck with continuous space tapers evaluated at random locations. This is a deliberate choice, as otherwise our method would become opaque to analysis. If we use tapers then to avoid bias, it is necessary to debias I hp (w). If we stick to the Fourier grid, then the bias is zero and not necessary to remove. Thus starting from multi-dimensional kernel density estimators (Duong and Hazelton, 2005) we instead chose a window function W (x) and matrix bandwith Ω. We then define where I h (w) is defined by (30). The 'h' just denotes the taper, and the tilde captures that we are using a bias-corrected periodogram. A useful modification is to consider one-dimensional window W and Ω = Ω and consider the averaging over directions, thus reducing the statistic to magnitude t = w only, which we refer to as the rotationally averaged periodogram. Mugglestone and Renshaw (1996b) call a version of it the R-spectrum. Alternatively, one could average also over the magnitudes to arrive at a summary function in direction. In 2D the convenient argument would be the polar angle; Mugglestone and Renshaw (1996b) call this version the θ-spectrum.
Basically both the multi-tapering and the local averaging smooth any direct spectral estimate (or bilinear estimator). Both of them have a smoothing parameter, namely the number of tapers P or the smoothing matrix bandwidth Ω. Let us determine how we select them, and discuss how this might be done in theory as well as practice. We shall focus on 2D as this is the best explored dimensionality after 1D.

Bandwidth selection in 2D
We have already discussed the tapering in (43). Some care has to be taken over the indices we loop over, especially the spatial index and the taper index. Eqn (43) has to contain the same index over the taper as it is a modulus square. On the other hand in (50) we loop over the taper indices to combine the 1-d Cartesian tapers to make 2-d Cartesian tapers, e.g. we would use 00 then 01, 10 as well as 11. The two types of indices should not be confused.
With this, we are now ready to turn our theoretical understanding of I t p (k), the tapered periodogram, into practical methodology. Starting from the multiple estimates of the spectrum we define the multitaper estimates for Note that when looking at 1D tapers p ranges from 1 to P 1 but when going higher dimension we need to loop over each dimension in turn. The higher the dimension, the more complex are our options. In 2D we shall simply explore Cartesian products of 1D tapers. The number of tapers P 1 and P 2 in each of the two dimensions in (50) needs to be selected to balance variance and bias. For intuition we shall first understand the mean square error of a local averaging procedure as that defined in (48). We shall thus be inspired by regular d-dimensional smoothing to determine how to smooth. Starting from standard understanding available for d-dimensional density estimation (Duong and Hazelton, 2005;Wand and Jones, 1993) we use the kernel function W (k) with matrix-bandwidth Ω. Scaling the kernel function using the bandwidth matrix we obtain The bandwidth function satisfies Note that this is not the same as the continuous version of the taper integrating to unity, see also the Section discusses the tapers.
With this definition we get the kernel density estimator Note that (50) is not a special case of (53). As each quadratic estimate I t p (w) is arrived at by linear operations on the DFT, and then implementing a modulus square, not all direct spectral estimates I t p (w) can be arrived at from smoothing the periodogram I 0 (w).
Lemma 6.1. Assume f (w) for all w = 0 is twice differentiable in w, and denote the mixed second derivative by ∂ 2 f ∂wi∂wj = f ij (w). Then the absolute bias of (53) as an estimator of f (w) is given for w = 0 by The next step is the variance and the mean square error of our smoothed spectral estimator.
Lemma 6.2. Assume f (w) for all w = 0 is twice differentiable. Then the variance of (53) as an estimator of f (w) is given for w = 0 by Proof. See Appendix O.
We put these results together and arrive at the optimal bandwidth.
Theorem 6.1. We assume that X is a point process with spectrum f (w) that is twice differentiable for all w = 0. Assume we smooth the periodogram using a box-car window with width σ, and takeσ = σ h as the bandwidth measures in terms of the length . In this case the mean square error takes the form of This mean square error is minimized at Proof. See Appendix P.
Naturally for a real data set, the parameters in this expression must be estimated, and the second derivative depends on the quantity we want to estimate. Given the theorem shows that an optimal smoothing bandwidth exists, it makes sense to pick one, and we propose, as we shall demonstrate in simulation studies, that visually a suitable degree of smoothing can be arrived at in examples. Thus the main point of this section is that an optimal wavenumber can be found, and it depends on the intensity of the point process as well as the roughness of the spectrum.

Multitaper number
This is all very well to find the bandwidth for the locally averaged estimator specified by (53), but we also want it for the multi-taper estimator in (50). The two most popular sequences of tapers to use are the prolate spheroidal wave sequences (Landau and Pollak, 1962) and the sinusoidal functions (Riedel and Sidorenko, 1995). The former is used because it is maximally concentrated in wave number whilst being compact in time/space; the second because it produces the least bias formally.
Note that our choice is not trivial; some authors chose to interpolate the point process to a regularly spaced process, and then taper. We could chose to design taper to match the observed locations, e.g. Bronez (1988), but the analysis of that choice would be complicated.
We shall use P d tapers in d dimensions, and in this case we can derive the same expression for the bias as in (54) but we now need to use the properties of the discrete prolate spheroidal eigensequences-the number of well-concentrated tapers with a bandwidth σ is P = 2 σ for a cuboid data analysis window (Percival and Walden, 1993b, p 369). We then get the same bias expression for each tapered periodogram as for k > σ of We do this by mapping the numbering of each dimension into a single enumeration. Having determined the expectation of each term we find that this defining the spectral window of the taper by Using spheroidal wave functions, assuming the observational spatial domain is a Cartesian product of intervals of length , we find that (Percival and Walden, 1993b, p 342) with P = 2 σ then This follows directly from the choice of direct product domains, and the use of the spheroidal wave functions. We again expand the spectral density function to We may as before note that We may argue about the bias in a manner analogous to before. We find: We now use that P = 2 σ. For this reason our bias remains the same as before: The variance of the estimator comes from individual contributions from covariance expressions: To approximate the covariance we observe Theorem 5.2 which gives the form of the variance, where we ignore the higher order terms. We then note Because of the scaling we get and so putting these together we get that We can then write We can yet again minimize in P and so use calculus to do so Again we see how the optimal taper number scales. There is a clear problem with the "curse of dimensionality", as the range of P decays with d. This indeed gives a very similar answer to our previous smoothing, and so we can estimate how many taper we need to use using the previous parameters. Again, the choice of taper necessitates knowledge of the second derivative of what we are trying to estimate-if we truly knew the second derivative we would not need to implement that estimation. Instead what the proof here shows is the scaling in the fundamental parameters of the process, and that a balance between variance and bias can be made. Furthermore this is not trivially done at the boundary of using P d or 1 taper. In practice we do not have the unknowns in equations (56) and (62). But the advantage of determining their value is that we know that an optimum exists: we shall therefore empirically look for it using techniques applied in time series. An approach in time series is to pre-select the wave number resolution 2W , see for example (Percival and Walden, 1993a, p. 369). The authors note that if the basic resolution is W then we select K = 2N W ∆t − 1. In fact for a time series sampled at ∆t steps apart we call 2N W ∆t the Shannon number (Percival and Walden, 1993a, p. 335). If 2W is the bandwidth in wave number of the datataper concerntation operator, then k = 0, ..., K − 1 where K < 2N W ∆t will yield a set of tapers {h k (·)} that are optimally concentrated to the wave number band [−W ∆t, W ∆t]. So one method of selecting the number of data tapers is to just select the wave number resolution w, and then take the value of K that matches the Shannon number.

Plotting Spectral Estimates
For regularly spaced time series it is very clear what frequencies/wave numbers we plot. We stick to the Fourier wave numbers, and plot up to the Nyquist wave number. But for a point process no such upper limit exists, as is also the case for randomly sampled stochastic processes, see (Lii and Masry, 1994).
What we will find at larger wave numbers is that the baseline contribution of λ dominates the spectrum. Unlike time series that are regularly sampled in time, and it is natural to assume that the spectrum of the process is decaying in wave number, we can only assume the spectral deviations f (w) are decaying with w .
How do we then determine what is the highest wave number we can do spectral estimation? To determine this number we must guess at the decay of f (w). We shall assume that If we knew w then our approach would be simple. We would take We can then assign w 1 as the wave number taking the value 1 2 I( 1 1 , 1 2 ) up to that taking the value 2 λ using the Fourier grid and then solve for (C 0 , α) using least squares, this giving us (Ĉ 0 ,α). With these estimates we solveĈ and use this as the highest wave number, w h .

Isotropy versus Anisotropy
We now consider the special case of isotropic point processes for which the complete covariance function depends only on distance, γ(z) = γ I (||z||) (see also Daley and Vere-Jones, 2003, p. 310-311). Such radial functions (Grafakos and Teschl, 2013) are special as we do not require any orientational specificity in our spectral representation. The isotropy of γ(z) transfers through the Fourier transform to the sdf so f (w) = f I (||w||). We shall refer to the treatment of isotropic random fields when implementing analysis (Ponomarenko and Perun, 2006). An alternative is discussed by (Errico, 1985) and (Durran et al., 2017), for random fields. The orientation invariance leads to dimensional reduction of the multidimensional Fourier transform. Three basic concepts are useful for the isotropic analysis: Recall the Bessel function of order m > −1/2 It is connected to the d-dimensional Fourier transform via (Vembu, 1961, Sec 3) we can also note relationships between radial and Cartesian Fourier representations, see for example Grafakos and Teschl (2013), or Estrada ( Then if we assume X is a stationary and isotropic process with a radial complete covariance function γ I , the sdf at wavenumber w with t = w takes the form i.e. the sdf is related to the Hankel transform of the complete covariance function. And finally recall the radially averaged set covariance function, To gain insight into how we might estimate the radial sdf besides numerically rotation-averaging the ddimensional periodogram using Eq. (49), let us analytically rotation-average Bartlett's periodogram (Eq. 20): which in the planar case simplifies to This explains the motivation behind the isotropic estimator discussed by Bartlett (1964) and Diggle et al. (1987). Them authors prefer to normalize by N (B) rather than by |B|, corresponding to further dividing (67) byλ.
Like the classical periodogram, the isotropic "shortcut" estimator of Eq.(67) is highly biased near 0, as has been observed by Diggle et al. (1987). The biases are described in the following proposition.
Proposition 7.1. Expectation of the isotropic Estimator, I. With the definition of the isotropic estimator given by (67), Proof. Follows by applying the Campbell's theorem and adding and subtracting If we compare the formulas for the isotropic sdf in Eq.(65) and the expectation in Proposition 7.1, we recognise two sources of bias, just like with the periodogram: A convolution with a functionν B /|B|, coming from the finite observation window like in the periodogram but this time radially averaged, and a centering bias term that is a radially averaged version of the d-dimensional bias. Diggle's truncation-to-local-minimum construction, see (19), tries to correct for the biases, but among other things, it will not work well for a clustered process as the spectrum near zero will be underestimated.
To work towards a tapered estimator, consider rotation averaging the debiased, tapered estimator of (30), in expectation, The problem is that a radial function h I might not exist for which h I (||x − y||) = h(x)h(y). But if we consider data-tapering the observed differences {x − y : x, y ∈ X ∩ B} instead of in B, we can device a taper h I in B ⊕ B = x∈B {B + x}, thus we propose the following isotropic estimator: where the bias correction term consists of the rotation average of both the taper and the set covariance of B, and depends only on t = |w|. We warn the reader that the estimator of (69) is not necessarily positive, but like previous authors (Lii and Masry, 1994;Bandyopadhyay and Lahiri, 2009), we sacrifice this to remove large part of the bias. Additional issue arises from estimation of the required λ 2 , as the standard estimatorλ 2 is biased with a term depending on λ/|B| and a term depending on the unknown second order product density of the underlying process. In our examples we use the slightly less biased λ 2 = N (B)[N (B) − 1]/|B| 2 which removes the first order bias. On the positive side, do note that the debiased isotropic estimator without particular tapering (i.e. h I ≡ 1) can be estimated without any tuning parameters.
We may now ask what radial taper to use. Tapers are normally designed for processes observed in discrete time or space. We have discussed separable filters, and now note that the choice of filter needs to match the observational domain. For example it is difficult to match a square observational domain with a radial filter, unless one chooses a compact filter inscribed by the box. There are therefore three considerations 1) demanding the taper be separable in space, 2) isotropic in space, or 3) exactly compact in space. Discrete and isotropic filters have been been determined numerically by Simons and Wang (2011). This will not be an option for us as we need to evaluate the tapers at random locations.
Slepian and co-authors have studied the design of tapers in arbitrary dimensions (Slepian, 1964). Often solutions are in terms of prolate spheroidal wavefunctions. Their discrete analogue is the prolate spheroidal wave sequences, which are similar to the set of Hermite function (Xu et al., 1999). However, the isotropic estimator can be seen as an estimator for the Fourier transform of the (non-stationary) first moment of the difference process Assuming the observation window is a cuboid B = We see that if j = for all j, then h I becomes radial. This is related to circular harmonic decompositions, see (Jacovitti and Neri, 2000). The authors of (Xu et al., 1999) recommend using a scaling factor which we have adjusted for a sampling region of z j ∈ (− j , j ), but which they match to the spheroidal wave functions. We propose setting a = 25, as this numerically seems to not down-weight too much data, or have a too large jump near the border of B ⊕ B.
In a sense our above construction is a "fix" as it only works for the 0th Hermite function, and using more than one taper will even if j = for all j not lead to isotropic functions. This is not unsurprising as the sampling domain will leave an inprint, and capturing this set of information requires using more than the first taper. Also, if we have a spatial sampling that is highly elongated rectangle rather than a square then another solution needs to be chosen.
Finally we may ask ourselves, what happens if we have a point process X which is anisotropic but we still evaluate (18) in 2D for the pattern. Because it incorporates angular averaging we expect the RHS of the expression to only depend on w as it does.
Proposition 7.2. Expectation of Diggle's 2D Estimator II. The large area expectation of Diggle's estimator given in (18) when B is a cuboidal box when the spectrum is not necessarily isotropic takes the form in terms of f (w) as defined in (11) of where T (B, w) is defined in Lemma 4.1.
Proof. The proof is provided in Appendix R.

Simulation studies
To study the behaviour of the estimators we have discussed, and to verify some of the debiasing results, we conducted a simulation study.

Details for executing the simulations study
We simulated three stationary models in 2D to represent the archetypal spatial patterns of regularity, complete randomness and clustering. Throughout the simulations the intensity, i.e. expected number of points per unit area, was kept constant at λ ≡ 0.01. The processes were simulated in increasing square observation windows B n = [− n /2, n /2] 2 such that on average we observed λ|B n | = λ 2 n = n = 25, 50, 100, 200, 400 or 800 points. For example, n=100 = 100 and n=400 = 200 spatial units. We chose the regime both to explore the often encountered small data scenarios and to ascertain asymptotic behaviour. The models below are tuned so that they produce patterns with visually distinct structures at n=100 .
Complete randomness i.e. the Poisson process has no adjustable parameters after λ was fixed. Spatial regularity is represented by the Matérn type II process (Illian et al., 2008, Section 6.5.2) in two variants, one with hard-core radius 2 (variant "r2") and one with radius 5 (variant "r5") which correspond to about 40% and 90% of the maximum allowed radius for the process given λ, respectively. The clustering behaviour is represented by the Thomas process (Illian et al., 2008, Section 6.3.2) with cluster intensity κ and Gaussian dispersal kernel standard deviation σ and again in two variants, one with many small or tight clusters (κ = 0.6λ, σ = 2, variant "MS") and one with a few large clusters (κ = 0.3λ, σ = 6, variant "FL"). The per-cluster expected point count µ got then fixed via the property µκ = λ to 1.67 and 3.33, respectively. The theoretical pair correlation functions and spectral density functions of the models are depicted in Figure 1. Example patterns of each model and observation window combination are given in Supplement Figure 6.
For the estimation of the 2D spectra we fixed the wavenumber grid to w ∈ [−0.3, 0.3] 2 with 101 equidistant steps in each dimension, giving the step size 0.006 in both dimensions. This scale was chosen to cover the interesting range of non-constant values for all models (cf. Figure 1). To reduce the effect of high wave number noise only the values on the sub-grid w ∈ [−0.2, 0.2] 2 were considered when integrating over w for the quality metrics discussed below. When a rotationally averaged curve was to be computed using Eq. (49), the averaging was done over the magnitude-grid |w| = 0.003, 0.006, ..., 0.300 using a box kernel and if not otherwise stated, radius of 1.25 · 0.006. The radial Hermite taper parameter was fixed to a = 25.
We The estimator for a periodogram is given in equation (30) with the debiasing term included. Bartlett's periodogram ("periodogram") has the taper h 0 (x) = 1 Bn (x)/ |B n | for window B n . For the multitapering ("mt P "), defined in equation (41) with P 2 tapers each having a parameter p = (p 1 , p 2 ) ∈ {1, ..., P } 2 , we used the orthogonal sine-tapers h p (x) = 1 Bn (x) 2 j=1 sin[πp j (x j + l n /2)/l n ]. For the kernel smoothed estimators ("smoothed m"), we first compute the Bartlett's periodogram and then, given the estimateÎ 0 on the wavenumber gird, we convolute it with a m × m discrete template having approximately Gaussian weights.
All computations were done using the R-software (R Core Team, 2020, v.3.6.3). Simulations were done with the help the R-package spatstat (Baddeley et al., 2015, v.1.64-1). Each of the discussed estimators were programmed into an R-package ppspectral, available on request from the first author.

The effect of debiasing
As we saw in Figure 2, the debiasing term has a large effect on the quality of the estimates. Table 1 provides the fractions of bias removed by the debiasing term for the Bartlett's periodogram and P = 3 multitapered periodogram. Nearly all bias is removed, with the only notable exception being the 65% bias left in the multitapered periodogram when data is a (arguably very) small sample of Thomas FL variant. We illustrate the form of the biases in Figure 3 for two model variants. The centering bias is always positive, as expected from equation (22). Bartlett's periodogram bias (cf. Lemma 5.1) is concentrated along the axis where one of the sinc-functions is constant 1 and near the origin where the second sinc function grows to 1. The multitapered periodogram has a square-shape around the origin.  Figure 4 summarizes the estimated integrated variance, squared bias and MSE of the estimators with various parameters. Only the debiased estimators are included, and the values are given in log 10 -scale and relative to the baseline given by Bartlett's periodogram. Multitaper with P = 1 (i.e. only one taper) does not differ in iMSE from the periodogram as the benefit of averaging does not apply. With more than one taper the variance is reduced, e.g. P = 3 (9 tapers) it goes down by 90%, and the bias stays at baseline. But adding more and more tapers (P = 6, 36 tapers) starts to introduce bias especially for low sample sizes. The bias is more prominent for the clustered cases, and based on the earlier discussion this is due to oversmoothing at small wavenumbers. Both the variance reduction and the oversmoothing happens with the post-hoc smoothing as well.

Overall quality of the estimators
From the results we can confirm that some smoothing is beneficial but again a balance must be struck to avoid oversmoothing (cf. Section 6). The regular process is more resilient to oversmoothing with reasonable sample sizes.

Rotational averaging and isotropic estimation
Isotropic estimation, either by averaging radially or by the isotropic shortcuts discussed in Section 7, is an important dimension reduction step for easier visual inspection of the periodogram (particularly for 3D data). The kernel-based rotation averaging depends on a bandwidth parameter, and we studied how this affects the quality metrics. We set the bandwidth radius to BF · 0.006 and varied the bandwidth factor BF from 1 to 10. Figure 5 provides the iBias 2 and iM SE for the debiased Bartlett's periodogram when comparing the rotationally averaged estimates to the true 1D isotropic sdfs (cf. Figure 1). The Poisson process has a constant sdf so oversmoothing will not be penalised. For the clustered processes even a small bandwidth introduces bias, but this is typically more than offset by the reduction in variance. The error for the few large clusters variant is more sensitive to changes in the bandwidth than the many small clusters. This is because the effective range of wavenumbers where the sdf exhibits structure is more concentrated, and thus easier to smoothen out, than the few large cluster variant (cf. Figure 1).
The strongly regular variant "r5" is similar to the clustered cases in that an optimal bandwidth is easy clearly present. For the less regular "r2" variant the error is like that of the Poisson process, possibly due to the target sdf being a weak oscillation around a constant and this being hard to detect. If we look the integrated squared bias we however do see that with large enough data oversmoothing is detectable.
We note that since multitapered periodograms and locally averaged periodograms are already smoothed, their rotational averaged estimates exhibit smoothing bias at smaller bandwidth factors (results not worth showing here).
To compare the radial averaging to direct isotropic estimation, we estimated the debiased isotropic periodogram and the squared-exponential tapered isotropic periodogram with a = 25. Table 2 provides the relative quality scores for the isotropic and the tapered isotropic periodograms relative to radially averaged periodogram with the best bandwidth for each model and data size, cf. Figure 5. The values are medians over the sample sizes n ≥ 50.
In iMSE sense the radially averaged periodogram performs the best, chiefly because iMSE is of order iVar and the radial averaging greatly reduces the variance. However, in terms of bias the isotropic estimator provides more accuracy, except for Poisson when no oversmoothing is possible. Adding a smooth taper to the isotropic estimator further reduces the bias, but can increase the variance, likely due to some data near the edges is being filtered out. The bias occurs near |w| = 0, mostly when |w| < 1/( √ 2l n ) (illustration in Supplements Figure 7).

Discussion
The Fourier transform of a point process has been treated by numerous authors (Mugglestone and Renshaw, 1996a;Bartlett, 1964Bartlett, , 1963Diggle et al., 1987). Despite this fact we have seen limited methodological development, leaving methodological understanding lagging much behind that for time series and random fields. There are a number of barriers for using and understanding the DFT of a point process. The first is computational: q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q MaternII r5   Table 2: Scores iM SE, iBias and iV ar, relative to radially averaged debiased periodogram with optimised smoothing bandwidth.
unlike the FFT the act of calculating the DFT is not O (N log N ), where N are the number of observed points. Second, there has been an often said, that Fourier analysis of point processes are not very different from that of random fields. However, the results we have developed here, show that the spectral moments are different between random fields and point processes, e.g. compare to (Stein, 2012). Following on from this realisation we have found a way of avoiding the treatment of low wave numbers differently. An important aspect of the spectral moments is the impact of the spatial domain we observe the point process in: the geometry of this object determines what wave numbers we should evaluate the Fourier transform at. Unlike for a regularly sampled stochastic process like a random field, the point process is biased at low wave numbers. This is because, as we demonstrate, we need to de-mean the spectral measure. By understanding the moments of the Fourier transform we show that we may remove this bias completely. Our realisation in this issue is that the sampling domain determines both bias and second order structure for a point process and needs to be studied directly. This has implications for our sampling of the wave number grid, where once we have debiased estimates can be formed at any wave numbers.
Third, we discuss various forms of variance reduction of a spectral variance estimator (the periodogram). We propose either a kernel smoothing of the debiased periodogram, or multitapering combined with a debiasing step. Our studies agree with the theory that unless the second derivative of the spectrum is zero, there is an optimal smoothing width. We also discuss constructing one-dimensional summaries from the two-dimensional spectrum, and tapering in many dimensions.
We have here taken the first steps to fully developing the methodology for spectral analysis of point processes, by calculating the moments of spectral summaries, deducing Gaussianity of the DFT, and fully determining how the wave number grid is influenced by the geometry of the region we observe the point process in. We imagine that authors starting from results such as those provided by this paper will seek to fill in the space remaining in between technical understanding of point processes (Daley and Vere-Jones, 2003) to their practical application (Lavancier et al., 2015;Møller and Torrisi, 2007).

B Proof of Lemma 4.1
Proof. We start the proof by noting that (using that B is centered) by the convolution theorem, after we have defined the functions: Note that from Gradshteyn 3.741 with the change of variables x j = πw j j that we get (recalling |B| = j j ): We have defined We now go back to our definition of EI (2) 0 (w), and explore this more carefully. We also fix 0 < β < 1, and ii expand the expectation as follows (realising that care needs to be taken around the pole): We can note directly, that taking a supremum over G (1) * (w − w ) in the range of integration, and assuming the integral of G (2) (w) is finite (which is true from Parseval-Rayleigh relationships

EI
(2) which requires β > 1/2 to be o(1). The same holds for the other elements apart from EI (2) 2...2 (k). We therefore have We assume thatf (w) is twice differentiable in w at w (everywhere possibly but w = 0) so that we get for all w = 0 and for w − w < w − w this implicitly defining ∇ f (w) the gradient, andH f (w) as the Hessian matrix.
We can then substitute this back in and need to calculate moments of G (1) (w) to understand the bias. Substituting this expansion back into (77) we get with a change of variables x j = π j w j : We can then note that By direct calculation we then find iii Note that this is now consistent with (73), and achieves a value of 1. Second, implies that this function |B × B −z | is symmetric (even) in z, and real-valued, which implies that its Fourier transform also is. The symmetric integral of an even times an odd function is always zero. We now also determine Let us start by a change of variables. Let us first define the matrix Thus the change of variables is With this change of variables we get that We then find if all the j = then if 0 < β This decreases with increasing as long as β < 2 (which follows as we have assumed β < 1), and as the next error term is O( −β ). We therefore want to chose β → 1, or for some fixed > 0 take β = 1 − . Putting all of these components together we get that iv C Proof of Proposition 4.1 Proof. We start from using Eqn. (9) in Section 2 and determine Subsequently we find using Eqn. (9) in Section 2 We now calculate the variance as We now define the renormalized quantity The expression in (86) then can be simplified to where we define To simplify this expression we note that v This yields the desired expression for the variance. We return to the expression of the relation of the DFT. We aim to show that We have We then note that We now seek to simplify this and write We additionally have from (84): Putting all of this together we get that We now implement a change of variables, and set z = x − y. We then have and so with the definition of the expression is a consequence of the above expressions.

D Proof of Lemma 4.2
Proof. We start from and so take We now use Campbell's theorem to evaluate these expectations. We already have from Proposition 4.1 that (2) (x − y)h(x) exp(−2πiw · x)h * (y) exp(2πiw · y) dx dy.
vi Also we note that We also note E {J h (w)} = λH(w).
Therefore we have We additionally note This can be noted from calculations in Section 2. We note that the multiplication in (95) leads to a convolution in the wave number domain. If we rewritē then rewriting (95) we get as given.

E Proof of Theorem 4.1
Proof. The expectation of the tapered Fourier transform is (see (9)) this defining the quantity EI t,(2) (w). The expression EI t,(2) (w) mirrors what we can see in (72). We see directly the benefits of tapering-|H(w)| will decay rapidly from zero so there will be no effect from the third term, and there will be less blurring in the second term apart from very close to w = 0.
From which we see that we can re-express this understanding in the Fourier domain, again using the convolution theorem. We define vii To explore this further we again do a change of variables: We define the inner product window using that h(x) is compactly supported to get We could also divide by |B ∩ B −z |h(x)h(x − z) as long as we are inside the domain. Outside B we cannot, and so this is why we cannot remove all bias in the spectrum even if we bias-correct as suggested by Guyon (1982). With this window we can note that Ideally we would like |H(w)| 2 = δ(w), but this is not possible as |B| is finite. If we assume f (w) is smooth and possessing two derivatives then if B is growing we can obtain nearly unbiasedness. We now assume that |H(w)| 2 is concentrated to a region in wave number W and again use a Taylor expansion. There are two approaches to this, either Riedel (1996) or Thomson (1990). The systematic bias: λ 2 h(B, w)h(B, −w) = |H(w)| 2 replaces the sinc functions for the square box. As we have chosen the function h to be well-concentrated Thomson (1990); Riedel and Sidorenko (1995); Riedel (1996); Walden (2000), the effect of this is limited.
We assume that We also assume thatf (w) is upper bounded by f 0 . We then have We note that We can yet again utilize the Taylor expansion of (78) inside W. We find We note that We have that |H(w)| 2 decays in |w|. If w / ∈ W then the influence of this term becomes negligible.
viii F Proof of Corollary 4.2 Proof. We start from We calculate the expectation of this object. This takes the form of We now use Theorem 4.1 to determine We now need to determine the middle term of (103). We need to determine (using Campbell's theorem) Also note Therefore we find that G Proof of Theorem 4.1 Proof. The proof of this result is somewhat more challenging. Assume that the domain B is symmetric with respect to the axis, so that the centre of mass of the region is at the origin. Assume the regions have a positive Euler characteristic, or that there are no holes in the observation region. Determine the mean length in each dimension and overlay a box B . Assume this box has width i in dimension i. Let us first describe this set-up in 2 dimensions and then what modifications are needed in higher dimensions. Now imagine B, the two dimensional observational domain need δB l removed and δB u added ti make up B . Thus B ∪ δB l = B ∪ δB u .
Note that we assume B ∩ δB l = ∅, and and also assume that i = Θ( ), to slave all dimensions to one scale.
It follows that we can construct 3 different cuboids (boxes), one inscribed by the shape, one with the average distance to the origin the same as with the shape (almost overlapping the shape), and one which inscribes the shape. We shall write these cuboids as B 0 ( ), B 1 ( ) and B 2 ( ). We note that the expectation of the modulus square of the Fourier transform is to leading order the same between the three regions, as long as they grow according to the order d .
This allows us to treat any region, as long as it is not too skinny or asymmetrical. This is why we slave the growth in any dimension to .

H Proof of Theorem 5.1
Proof. From first principles using Campbell's theorem we may deduce that the uncentered second moment of I t p (w 1 ) with I t l (w 2 ) takes the form of This latter equation defines the three terms V 1 (w 1 , w 2 ), V 2 (w 1 , w 2 ) and V 3 (where the latter does not depend on w 1 and w 2 even if V 2 is the sum of two terms, one only depending on w 1 and the other only depending on w 2 ). We now need to further split these terms up in order to learn about the cases when we can have repeats of locations or not. We start with V 1 . First, we note: x Using the simplification implied by (9) we obtain the expression of We will be able to use the orthogonality between the data tapers {h l (x)} but this will be easier in the Fourier domain where we can assume local smoothness of the Fourier transform of ρ (n) (. . . ), and carry out the usual Taylor series. Now we start by implementing the calculations for the Poisson case, to see the mechanics. The next term in the expansion is then The final term is in turn using (9) We can put these terms together to determine the covariance between the periodogram at two different wave xi numbers and using two different tapers as follows: Figuring out the size of this contribution is a little bit more complex. The integral will do a Fourier transform in y of h p (y) in w 1 and w 2 , and also the conjugate will be taken. The Fourier transform in y will be supported when |w 1 | < b w and |w 2 | < b w , but will not be supported when either magnitude gets too large. The second Fourier transform in x will be somewhat more spread-this is because we are Fourier transforming h 3 p (x), which will be a third order convolution once Fourier transformed. But as long as we can use the joint concentration in |w 1 | < b w and |k 2 | < b k this will still be negligible.
The next term is given by To see why this is true, note that if h p (x) is constant then we can easily calculate the higher order norms. We find by using the Cauchy-Schwarz inequality that if p = q that We also note that We now need to determine an upper bound. We then look at Then we have We can for any taper h(x) calculate B h 4 (x) dx explicitly. For the constant function we get 1 |B| . For tapers well-concentrated we would expect a similar decrease, but for any choice of taper we can calculate the value of the 4th norm explicitly.
Thus we understand a bit more about this term. Moving on to the next aspect of the computation, finally, This shows each individual contribution as long as |w 1 − w 2 |, |w 1 + w 2 | > b w . This completes the proof.

L Proof of Theorem 5.2
We start from Theorem 5.1. We use that we have set φ w (x) = e −2πiw·x . We then have: h(x)h(y)φ w (x − y)ρ (2) (x, y)dxdy .
We wish to simplify this expression to get a sense of the major contribution to the variance. We recall that we are analyzing a stationary process and so with some abuse of notation ρ (n) (x 2 − x 1 , . . . , x n − x 1 ) = ρ (n) (x 1 , x 2 , . . . , x n ).
These are higher order generalizations of the common spectrum. Apart from possibly at w = 0 we assume that these functions are twice differentiable, just like we did for the centred spectrum. We can now exchange the representation of Var{I h (w)} using f (n) (w), rather like our treatment of E{I h (w)} in Theorem 4.1; and then as a last step correct with the effect of a Poisson process. Referring back to (120) we see that all terms are multiplications of ρ (n) (·) with products of h(·) in space. These in the Fourier domain become convolutions between the spectrum and higher order generalizations and the window functions (products of h(x)). This follows from the convolution theorem.
To arrive at a simplified expression, we first transform betweenρ(τ ) and ρ(τ ). This adds exactly the variance term of a Poisson process, in all places but the final line of (120) which is the only place ρ(τ ) does not appear linearly. Once the representation is in terms of ρ(τ ), then we transfer the product to a convolution in the Fourier domain.
We now need to use some of the large sample characteristics. (125) We are left with terms T 4 and T 6 as well as T 0 and T 8 .
We see that the variance increases strictly with the rate λ > 1.

M 4th Order Moment of Taper
We start by using the Cauchy-Schwarz inequality to obtain We also note that We now need to determine an upper bound. We then look at Then we have We can for any taper h(x) calculate B h 4 (x) dx explicitly. For the constant function we get 1 |B| . For tapers well-concentrated we would expect a similar decrease, but for any choice of taper we can calculate the value of the 4th norm explicitly.

N Proof of Lemma 6.1
The estimate in (53) smooths in a local region of w as determined in shape by W , and scale by Ω. We assume either that either W (w) = W (w 1 ) . . . W (w d ) and that W (w j ) = W (−w j ) for 1 ≤ j ≤ d (e.g. that the function is separable and even) or that W (w) = W (|w|) (e.g. that the function is isotropic). We can then determine the bias and variance of this estimator using Theorems 4.1 and 5.2, which will be the next step in our study.
As a first step we calculate First we need to settle the normalization of W (w), and make that consistent when going between continuous and discrete operations. We note that if w = 0 then we may assume that f (w) is twice differentiable and so we may appeal to Theorem 4.1 to quantify E I h (w; Ω) . We need to be a little bit more careful with the order term because it will matter.
Referring back to Corollary 4.2 we may note E I h (w) = f (w) + e(w).
Here e(w) is defined as the error, but we can deduce its order from the proof of Theorem 4.1.
We recall the constraints of producing an unbiased estimator, namely referring back to (52) w W (w) = 1.
O Proof of Lemma 6.2 For the periodogram we approximate the second order properties by We derive these results from Theorem 6.2 combined with Proposition 6.2. Starting from these results, we shall now calculate the variance of Again using a separable and square window we have w ∈K(B (l)) |H| −2 W 2 H −1 (w − w) = ( hσ) −1 .
Differentiating we then get We have therefore checked this is a minimum. We now need to estimate these parameters in order to get a suggested bandwidth. Some of the parameters are known, e.g. d. The curvature is a little bit more complex. We shall motivate our curvature estimate by a single length scale in the spectrum.

Q Proof of Proposition 7.1
Note that by direct computation

R Proof of Proposition 7.2
An interesting question is naturally; what if we use Diggle's estimator, even if the point process is not isotropic. We recall that the estimator takes the form We can still compute the estimator for any observed point-process X, even if X was not generated by an isotropic process. The estimator I D (|w|) has expectation EI D (|w|) = λ + |B| −1 x =y EJ 0 (2π|w| · x − y ) xxii the latter defining the form of these three contributions. Just like before we shall explicitly demonstrate the effects of this convolution. We have that the Fourier transform if the Bessel function is F {J 0 (2π|w| · |z|)} (u) = J 0 (2π|w| · |z|) e −2πiz T u dz = ∞ 0 2π 0 J 0 (2π|w|r) e −2πir|u| cos(φ−φu) 2π|w|d|w|dφ = ∞ 0 J 0 (2π|w|r) 2π|w|d|w| 2π 0 e −2πi|w||u| cos(φ−φu) dφ.
We now use the convolution theorem to deduce that: Thus the low frequency bias is determined by this term. Thus the bias is worse for this version of the wave number domain representation. We then see how this "blurs out" the sinc kernel. This implies that the second part of the integral becomes

EI
(2) We therefore need to compute With these pieces we have that

EI
(2) We now see that further blurring is present in (147). In the standard non-isotropic case this was just a convolution of f (w) with |B| −1 T (B, w). To get a feeling for its behaviour we note that as |B| −1 T (B, w) → δ(w)/(2π). In this case we get for d = 2

EI
(2) This shows that asymptotically we would recover the isotropic spectral density from this component. Note that the first term, or EI (2) D (|w|) is just λ and already determined. By combining the three terms, we are left with the expectation of the Diggle isotropic estimator. By combining the terms with arrive at the result stated in the proposition.