What is the Fourier Transform of a Spatial Point Process?

Spatial point processes are an important form of observational data structure, for example in forest ecology [77], communications networks [57, 49], epidemiology [31], social science [72], pharmacology [35] and medicine [2] amongst many other application fields. Understanding the properties of a point process can be approached from many different perspectives [24], and the aim of this paper is to determine how to extract the frequency (or scale/wavenumber) behaviour of an observed point process, as well as connect that to its spectral representation.

The most common assumption used when analysing spatial processes is that of spatial homogeneity (or stationarity). This means that if you shift all of the observations by a fixed spatial shift, the distribution of those observations does not change from the distribution of the original sample; and if the observational window is changed, then understanding the full set of observations remains tractable. The consequence of this probabilistic invariance in distribution is the spectral representation of a stochastic process. The existence of the spectral representation of a stochastic process means the Fourier transform fully characterises the second order properties of that stochastic process. The Fourier transform also characterises the second order properties of a point process, see e.g. [20]. Yet unlike random fields and time series, spectral analysis of point processes is still in its infancy, see also [5, 6, 56], and critically, the digital processing of a point process remains fully outstanding. Recent interest in Fourier features in machine learning based approaches for point patterns such as  [38, 44] show the potential of using Fourier based information as features for estimation and detection. The work in this manuscript both establishes what Fourier features to calculate for homogeneous processes from a sampling perspective, and their large but finite sampling area properties, just like [33] determined the large but finite properties of Fourier representations for random fields. We note that machine learning techniques [73, 48] utilised Fourier features in learning algorithms for parametric models before precise and general understanding was established.

Thus despite inspirational work by the aforementioned authors, several key aspects of discrete spectral analysis lie unresolved when applied to observed point processes. In particular, 1) given the Fast Fourier Transform is unavailable (unlike the case for regularly sampled time series or random fields), how do we efficiently implement the calculation of whichever discrete Fourier transform we chose to define? 2) How does that discrete transform relate to the underlying spectral measure of the process? 3) Can the discrete transform be improved by linear operations such as tapering as is the case for other stochastic processes? 4) If yes, how do we then select a taper? 5) How do we define a radial or isotropic transformation to simplify the representation of the spatial point process? There are many, many more questions to answer before bringing the spectral analysis of spatial point processes to the sophistication of the analysis of random fields, but these currently unresolved questions are positioned as the first hurdles to overcome in the path of whomever wishes to develop more sophisticated spectral analysis techniques. Once the fundamental properties of the spectral representation of a stationary point process have been developed and understood, results could be extended to inhomogeneous point processes, but that is not the aim of this work.

To put this in the context of information theory; understanding how to compute the Fourier transform from a spatially compact signal had already sparked a lengthy debate from the first introduction of tapering [70], and selection of optimal tapers using localization operators [67, 66, 17]. While some aspects of tapering in the context of point patterns is reflected by topics in mature (multi)dimensional tapering and spectral estimation [1, 45, 46, 33], understanding how to adapt these ideas to spatial point processes remains fully outstanding, and must be answered using our understanding of irregularly sampled processes [11]. Understanding and characterising point patterns is naturally a problem of general interest, e.g. [12, 13].

Why then do we want to understand the spectral information of a point pattern? The spectral information of any stochastic process is the same as the spatial information as there is a bijection between the two, but the former is directly showing the scale of variation of the process. To illustrate, let us simulate a complex point pattern and compare its spectral content with its usual spatial representation. Log-Gaussian Cox processes [55], popularly used in applications, are marginally stationary processes whose patterns are more variable than models that do not depend on latent random variables. The pattern (a single realization), and its estimated spectrum (the subject of this paper) are shown in Figure 1. We have chosen a process which is isotropic, a choice we have made so that its usual spatial representation and spectrum is isotropic also, and thus easier to interpret. In the left hand subplot we see the pattern itself, in the middle subplot we show a common spatial summary, the pair correlation function (pcf) comparing both the estimate and truth, and in the right hand subplot we show the spectrum on a decibel scale ($10\log_{10}$), both estimated and true. While the raw pattern (left) looks mainly unremarkable, naturally its periodicity is present in the pcf (middle), but is more obvious in the spectrum (right). We can even estimate the periodicity as having frequency 1. This is immediately recoverable directly from its estimated spectrum.

To get to the point where we can estimate the spectrum of a point pattern, we need to answer the five questions posed earlier. We will therefore start by determining how to compute Fourier transforms of point processes, discussing the question of how to form direct spectral estimators [58] in this setting. This will be answered in terms of the mean and (co)-variance of the different discrete Fourier transforms that we define. The first surprising result is that the natural direct spectral estimator is biased unless it is mean corrected (this bias was noted by [23], but only corrected in an ad-hoc fashion). We show how to perfectly eliminate this bias.

Second, to understand the variance and covariance of a direct spectral estimator, we need to calculate a fourth order moment of our choice for the discrete spectral transform. This is complicated since Campbell’s formula is required to derive its form for a point process. This approach can be contrasted with using Isserliss’ theorem to determine higher order moments for a Gaussian Process [40]. Calculating covariances between direct spectral estimators gives us a way to determine what grid of wavenumbers the spectral estimator is uncorrelated at, and thus where to evaluate it.

Third, tapering [70] is required to define the direct spectral estimator and avoid leakage; but how to taper a point pattern is an open question. Most tapers are defined for regularly spaced stochastic processes; but there are some continuous tapers. We choose to use the continuous tapers of Riedel and Sidorenko [60] to construct the direct spectral estimator, and use multitapers to reduce the variance of such estimators. We will determine how to implement both.

Fourth, many spatial models are radially symmetric, or isotropic. This prompts us to describe how to construct an isotropic representation of the spectral content of a point processes. There are two possible approaches, namely using the Bessel function [23] to do the transformation; or isotropizing a general spectral representation [25]. We describe how to define the appropriate tapering in this instance, inspired by other isotropic harmonic decompositions [41].

Calculating the Fourier transform of a point process is useful beyond the point pattern itself. If we sample a stochastic process with a point process then the spectrum of a point pattern is equally important to that of the stochastic process to determining the spectrum of the observations, e.g. [50]. For this reason; the study of the frequency information of an observed process is of interest in its own right, and beyond, and we will discuss the implications that arise for point processes.

In Section 2 we define the basic concepts required by the second order representation of a point pattern. In Section 3 we discuss tapering of the DFT and the direct moments of the DFT. In Section 4 we determine how to form a spectral density estimator from our understanding of how to form the DFT and its mean and variance. The next step is to study the covariance structure of such spectral density estimators, see Section 5. We then use that understanding to define how to do linear smoothing in Section 6. For multidimensional spectral representations, the full anisotropic structure can be hard to interpret; we therefore propose 1d isotropic or isotropised summaries to characterise such processes, in Section 7. We then present a simulation study in Section 8, and demonstrate the potential of the proposed methodology in an example from forest ecology in Section 9. We conclude with a discussion in Section 10.

In this Section we shall give the basic notation necessary for the spectral analysis of a spatially homogeneous point process. We assume that we observe a spatially homogeneous point process $X$ with intensity $\lambda=\rho^{(1)}>0$ and with $\rho^{(2)}(z)$ defined as the (second order) product density of $X$. Assuming $X$ is spatially homogeneous means that $\rho^{(2)}(z)$ only has one argument ($z\in{\mathbb{R}}^{d}$, namely the spatial shift) rather than depending on two local spatial variables (say $x$ and $y$ rather than just $z=x-y$). We assume that $X$ is a simple point process on ${\mathbb{R}}^{d}$, meaning no duplicate points are allowed. We take $d=2,3,\dots$, and so exclude the case of $d=1$, which is relatively well-studied, see for example the references in [19, 43, 75]. For a (Borel) set $A$ in ${\mathbb{R}}^{d}$, we define $|A|$ for the volume of $A$ and $N(A)=|X\cap A|$ for the number of points of $X\subset{\mathbb{R}}^{d}$ in $A$. For set $A$ and $z\in\mathbb{R}^{d}$, denote the $z$ shifted set by $A_{z}=\{x+z:x\in A\}$. We denote by $B$ any subset of ${\mathbb{R}}^{d}$ where points are observed (or equivalently our observation domain). For any complex number $z$ we use the superscript $z^{\ast}$ to denote the conjugate of $z$.

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 [6]. The complete covariance function of a stationary point process is therefore

$$\gamma(z)\equiv\lambda\delta(z)+\rho^{(2)}(z)-\lambda^{2},\quad z\in{\mathbb{R}^{d}},$$

(1)

where $\delta(\cdot)$ is the Dirac delta function. The spectral density function (sdf) [6] of the point process $X$ is then defined as the Fourier transform of the complete covariance function:

$$f({k})\equiv\mathcal{F}[\gamma]({k})\ =\int_{{\mathbb{R}}^{d}}e^{-2\pi i{k}\cdot z}\gamma(z)\,dz\ =\ \lambda+\int_{{\mathbb{R}^{d}}}e^{-2\pi i{k}\cdot z}[\rho^{(2)}(z)-\lambda^{2}]dz,\qquad{k}\in{\mathbb{R}}^{d}.$$

(2)

This representation is in direct analogy with the corresponding spectral decomposition of a random field or a time series. The symbol $\mathcal{F}$ denotes the Fourier transform, with $i$ as the imaginary unit. We refer to the argument ${k}$ as the “wavenumber” rather than frequency to acknowledge its multi-dimensionality. This is a natural choice of a Fourier transform in analogy to the analysis of random fields [68]. To avoid additional constants in the inverse Fourier transforms, we parameterise the Fourier transform with wavenumber instead of the customary angular frequency ${\omega}=2\pi{k}$ used by [6] and some others.

For a time series or random field there are a number of spectral results, ranging from the spectral representation theorem (I) [58, p130] and to Bochner’s (Herglotz) theorem (II) [58]. Both these results are not exactly available in the point process setting, but for a time series $X_{t}$ we can note them:

	$\displaystyle X_{t}$	$\displaystyle\overset{(I)}{=}\mu+\int_{-\frac{1}{2}}^{\frac{1}{2}}e^{2\pi i{k}t}dZ({k}),\quad{\mathrm{\bf Var}}\{dZ({k})\}\overset{(II)}{=}f({k})d{k},$		(3)
	$\displaystyle\gamma(\tau)$	$\displaystyle\overset{(III)}{=}\int_{-\frac{1}{2}}^{\frac{1}{2}}f({k})e^{2\pi i{k}\tau}\,d{k},\quad f({k})\overset{(IV)}{=}\sum_{\tau}\gamma(\tau)e^{-2\pi i{k}\tau}.$		(4)

Equation (3) decomposes $X_{t}$ into random contributions associated with each frequency ${k}$. 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({k})$ being considerably larger relative to other $f({k}^{\prime})$ as in that scenario ${\mathrm{\bf Var}}\{dZ({k})\}$ is bigger than ${\mathrm{\bf Var}}\{dZ({k}^{\prime})\}$ and therefore $|dZ({k})|$ likely to be larger than $|dZ({k}^{\prime})|$. Spectral representations of point processes are also discussed in [20, Ch. 8]. We can yet again decompose the (complete) covariance in terms of a spectral representation as in (II), and it takes the form in terms of a $d$-dimensional Fourier transform of

$$f({k})=\int e^{-2\pi i{k}\cdot z}\gamma(z)\;dz=\lambda+\int e^{-2\pi i{k}\cdot z}\left\{\rho^{(2)}(z)-\lambda^{2}\right\}\;dz.$$

(5)

We refer to $f({k})$ as the spectrum of $X$. We note that the Fourier transform can be inverted to yield the equality of (an analogy of III):

$$\rho^{(2)}(z)=\lambda^{2}+\int_{{\mathbb{R}}^{d}}f({k})e^{2\pi i{k}\cdot z}\,d{k}.$$

(6)

We assume the spectrum is the primary object of interest in our study of point patterns, as the covariance $\gamma(z)$ can be fully determined from it, and the covariance characterises the point process. We see directly from (5) that $\rho^{(2)}(z)-\lambda^{2}$ and not the complete covariance $\gamma(z)$ is playing the role of a time series covariance. From a nonparametric perspective the spectrum characterises what wavenumbers are more notable (distinct) in the process. As (5) specifies a constant level $\lambda$ to all wavenumbers, the notable (distinct) wavenumbers are determined from the Fourier transform of $\rho^{(2)}(z)-\lambda^{2}$.

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 from a simple spatial understanding. In time or space the notion is associated with the support of $\gamma(z)$ or $\rho^{(2)}(z)$. Being supported over a wavenumber ${k}$ means variation over scales $1/\|{k}\|$ is present in the correlation function, or to represent the variability in $\gamma(z)$ we need wavenumbers $\|{k}\|$ present in the spectral representation. The function $\rho^{(2)}(z)-\lambda^{2}$ is often approached in terms of what scales it is non-zero at, but the variation in the function can also be associated with long or short scales. So for a covariance function, we now have two notions of what it means to possess scales $\|{k}\|$; either $\gamma(\mathbf{u}/\|{k}\|)$ is non-zero for unit vector $\mathbf{u}$ or $\gamma(z)$ is variable (changing) and the scale of the change is $\|{k}\|$. 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$.

This is clear from Figure 2 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’ complete covariance function is discontinuous, its Fourier transform is supported over all scales (due to the Heisenberg-Gabor uncertainty principle [14]). This can be deduced from (6). To reproduce the discontinuity we need all scales in the Fourier representation.

We see from (5) that not all wavenumbers are given equal weighting. First all wavenumbers are given an equal weighting $\lambda$ and then the Fourier transform of $\rho^{(2)}(z)-\lambda^{2}$ determines the wavenumbers that are up–weighted (added) or down–weighted (subtracted) relative to the overall level of $\lambda$. 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 $\mathbb{Z}^{d}$ the simplest random process, white noise, is constant across wavenumbers yielding a spectrum that takes the form of $\sigma^{2}$ on $[-\frac{\pi}{\Delta},\frac{\pi}{\Delta}]^{d}$, where $\Delta>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 $\lambda\delta(z)$ in space. However once we remove this term, we expect a decay of the remaining spectrum $f({k})-\lambda$ as $\|{k}\|\rightarrow\infty$. Otherwise  (5) would contain additional singularities. Furthermore at ${k}=0$ we retain

$$f(0)=\lambda+\int\left\{\rho^{(2)}(z)-\lambda^{2}\right\}\;dz.$$

(7)

The second term in this expansion is not required to be zero, and for the Thomas process for example, it is not zero, as it takes the form of [56, p 55]:

$$f({k})=\lambda+\lambda\mu e^{-4\pi^{2}\|{k}\|^{2}\sigma^{2}}$$

(8)

where $\mu$ is the per-cluster expected point count and $\sigma$ is the Gaussian dispersal kernel standard deviation. It is clear from this expression that we cannot arrive at a zero contribution due to the exponential.

As we shall be studying the moments of a point process, it is convenient to restate Campbell’s theorem [10, Sec. 4.3.3], and we shall use this result multiple times. The theorem applies to any measurable function $h:R^{nd}\mapsto R^{+}$ with $n=1,2,...$, and states that (assuming product densities of order $n$ $\rho^{(n)}$ are well defined for any given $n$)

$${\mathbf{E}}\sum_{x_{1},...,x_{n}\in X}^{\neq}h(x_{1},...,x_{n})=\int_{{\mathbb{R}}^{nd}}h(x_{1},...,x_{n})\rho^{(n)}(x_{1},...,x_{n})dx_{1}...dx_{n},$$

(9)

where the summation is over distinct point tuples. Note that if the point pattern $X$ is stationary then $\rho^{(1)}(x)=\lambda$ is constant for any $x$, and $\rho^{(n)}(x_{1},...,x_{n})=\rho^{(n)}(x_{1}-x_{n},\ldots,x_{n-1}-x_{n},0)$. With some abuse of notation, we shall use the same symbol for non-stationary as well as stationary product densities, where the latter function has $n-1$ arguments for an $n$th order product density, e.g. we write $\rho^{(n)}(z_{1},...,z_{n-1})=\rho^{(n)}(z_{1},\ldots,z_{n-1},0)$.

We define a new parameterisation to capture how the process’s $f({k})$ differs from that of a Poisson process. We define the $n$th deviation of the $n$th order product density as

$$\widetilde{\rho}^{(n)}(z_{1},\dots,z_{n-1})=\frac{{\rho}^{(n)}(z_{1},\dots,z_{n-1})-\lambda^{n}}{\lambda^{n}},\quad n=2,3,\dots,$$

(10)

where the argument of the product density $z_{l}\in{\mathbb{R}}^{d}$ for $l=1,\dots,n-1$. For a Poisson process these $n=2,3,\dots$ deviations from Poissonianity will all be identically zero. By understanding deviations from Poisson behaviour, we get greater insight into the underlying process of interest. For completeness we also define the Fourier transform of the deviations:

$$\widetilde{f}^{(n)}({k})={\cal F}\left\{\widetilde{\rho}^{(n)}(z)\right\},$$

where if $n=2$ we suppress the superscript. With this definition we find that

$$\widetilde{f}({k})=\frac{f({k})-\lambda}{\lambda^{2}},\quad w\in{\mathbb{R}}^{d},\quad f({k})=\lambda+\lambda^{2}\widetilde{f}({k}).$$

(11)

Thus in a sense, $\widetilde{f}({k})$ encapsulates the deviation of the function from a constant spectral density of $\lambda$ via the term $\lambda^{2}\widetilde{f}({k})$.

Already [20, 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 $\lambda$ has been subtracted (so the decay of $\widetilde{f}({k})$ is reasonable to assume for large magnitude wavenumbers). Having established the theoretical spectral description of point processes, we now turn to their sampling properties.

In this section, we shall revisit the possible definitions of a spectral estimator for a point pattern. We shall start by defining a direct Fourier transform and taper this definition to ameliorate edge effects, resulting in families of spectral estimators. We also discuss other choices of spectral estimators used in the literature already. To see the spectral characteristics of $X$ we start from what is known as Bartlett’s periodogram estimator [23] based on observing a point pattern in region $B\subset{\mathbb{R}}^{d}$, written as

$\displaystyle I_{0}({k})$

$\displaystyle:=\hat{\lambda}+|B|^{-1}\sum_{x,y\in X\cap B}^{\neq}e^{-2\pi i{k}\cdot(x-y)},\quad{k}\in{\mathbb{R}^{d}},$

(12)

where we have set $\hat{\lambda}=|B|^{-1}|X\cap B|$. This definition uses all the data (or points) available to us, $\{x:\;x\in X\cap B\}$, but is simply a possible choice amongst all possible direct spectral estimators. Normally a direct spectral estimator is one which is bilinear in the DFT of the data, see for example the discussion in [58, Ch 5–6]. We recall a bilinear form is a function $D({\mathbf{x}},{\mathbf{y}})$ satisfying in the first argument $D({\mathbf{x}}_{1}+{\mathbf{x}}_{2},{\mathbf{y}})=D({\mathbf{x}}_{1},{\mathbf{y}})+D({\mathbf{x}}_{2},{\mathbf{y}})$, and equivalently in the second argument. Bilinear forms have been thoroughly discussed in signal processing for the analysis of time series  [52]. 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({\mathbf{x}},{\mathbf{y}})$ simply generalises a bilinear form to the Hermitian symmetry, additionally requiring $D({\mathbf{x}},{\mathbf{y}})=D({\mathbf{y}}^{\ast},{\mathbf{x}})$ as well as $D({\mathbf{x}},{\mathbf{y}})$ satisfying $D({\mathbf{x}}_{1}+{\mathbf{x}}_{2},{\mathbf{y}})=D({\mathbf{x}}_{1},{\mathbf{y}})+D({\mathbf{x}}_{2},{\mathbf{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. Guyon’s spectral estimator [34]. We define the tapered DFT of a point process $X$ for a specified general (square) integrable function $h(x)$ (referred to as a ‘data taper’ by  [58]) with Fourier transform $H({k})$ and unit norm (i.e. $\|h\|_{2}=1$), to be

$$J_{h}({k}):=\sum_{x\in X\cap B}h(x)e^{-2\pi i{k}\cdot x},\quad{k}\in{\mathbb{R}}^{d}.$$

(13)

Spectral estimators in time series are bilinear in the observed real-valued process $X_{t}$ and sesquilinear in its Fourier transform $J_{h}({k})$. 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 to achieve. We shall also define

$$\widetilde{J}_{h}({k})=J_{h}({k})-\lambda H({k}).$$

Because $x$ appears in the argument of the complex exponential, we cannot define an estimator that is bilinear directly in the data. In practice also one has to decide what taper function to use. For 1D point processes tapering has been used [15], but also multitapering has been used on interpolated data [21]. Here, we do not interpolate, do not implement localised analysis, and do not implement analysis in 1D. It is more complex to implement interpolation in higher dimensions, as they are not orderly, unlike a time series time argument.

Based on $J_{h}({k})$ the natural spectral estimator becomes its modulus square or

$$I_{h}({k}):=\left|J_{h}({k})\right|^{2}=\sum_{x\in X\cap B}h(x)h^{*}(x)+\sum_{x,y\in X\cap B}^{\neq}h(x)h^{*}(y)e^{-2\pi i{k}\cdot(x-y)}.$$

(14)

If we take $h(z)=h_{0}(z)=|B|^{-\frac{1}{2}}\mathbf{1}(z\in B)$ then we recover Bartlett’s periodogram in (12) from (14). If we do not take the $h_{0}(z)$ taper then Bartlett’s periodogram is not exactly recovered. 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 [58]. We study the whole family of estimators $I_{h}({k})$, where a new member of the family is defined for each choice of $h$. We also note that (14) is a direct spectral estimator [58, p. 207]. Multidimensional tapers are available to us [36, 65], and can be pre-computed. Our choice of tapering corresponds to using continuous tapers evaluated at random locations. We have deduced an asymptotic distribution for $J_{h}({k})$ using [7] but this result cannot be applied to $I_{h}({k})$ even if it superficially may seem to be of the correct form.

Note that the sum in (14) cannot be evaluated in a computationally efficient manner, unlike the DFT, as the locations $\{x\}$ are not regularly spaced, which is unfortunate, but unavoidable. Finally (14) is bilinear in $J_{h}({k})$ 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 [23, Eqn. 3.3] in 2D, namely

$\displaystyle\bar{I}_{0}(\|{k}\|)$

$\displaystyle:=\hat{\lambda}+\frac{1}{|B|}\sum_{x,y\in X\cap B}^{\neq}\mathcal{J}_{0}\left(2\pi\|{k}\|\|x-y\|\right),$

(15)

where $\mathcal{J}_{0}(x)$ denotes the Bessel function of order 0, specified in Section 7. The $d$-dimensional extension is given later in equation (18). $\bar{I}_{0}$ is less clearly bilinear in the data, but if we start from the estimator $I_{0}$ of (12) and average it over orientations analytically, then we arrive at this form. We therefore with a slight abuse of our terminology also refer to it as a bilinear estimator. There is one additional correction made by [23]. As we shall see, for low wavenumbers there is a bias inherent in (15). To address this problem [23, Eqn. 3.4] suggests taking for some choice of lower bound $t_{0}>0,$

$$I_{D}(\|{k}\|)=\left\{\begin{array}[]{lcr}\bar{I}_{0}(t_{0})&{\mathrm{if}}&\|{k}\|\leq t_{0}\\ \bar{I}_{0}(\|{k}\|)&{\mathrm{if}}&\|{k}\|>t_{0}\end{array}\right..$$

(16)

The authors of [23] suggest taking $t_{0}$ so that $\bar{I}_{0}(t_{0})$ is at a minimum, and also suggest smoothing the estimated ${I}_{D}$, which we clarify in Section 6. The authors also discuss iterated methods of bias correction in their Section 5, and correcting biased estimators of the correlation.

Note also that we have modified the estimator in (15) to divide by $|B|$ rather then the observed number of points in the region, $N(B)$, as do for example [23], as it is much preferable not to have a random denominator. We shall discuss the usage of $I_{0}(\cdot)$ versus $I_{D}(\cdot)$ further on in Section 7, and the implications of when the process is truly isotropic or not.

As can be surmised from (14) the expectation of the general estimator ${\mathbf{E}}\left\{I_{h}({k})\right\}$ is a convolution between the Fourier transform of the observation window $B$ 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. Inspection of (14) also raises a second problem, namely that there is a constant contribution $|B|^{-1}|X\cap B|$ which does not give any wavenumber specific information, and is correlated between wavenumbers. 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 [50, 11].

Should we choose to taper isotropically in 2D then we instead use an isotropic taper $h_{I}(\cdot)$

$$\bar{I}_{h}(\|{k}\|):=\hat{\lambda}+\sum_{x,y\in X\cap B}^{\neq}h_{I}\left(\|x-y\|\right)\mathcal{J}_{0}\left(2\pi\|{k}\|\cdot\|x-y\|\right).$$

(17)

Note that this cannot necessarily be constructed. Say $J_{h}(\cdot)$ 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}(\cdot)$? Radial tapers can be determined as described in [65], whether in 2D continuous space or discrete space. Equation (17) can be of course be extended into higher dimensions. In general we would have in dimension $d=1,2,3,4,\dots$

$$\bar{I}_{h}(\|{k}\|):=\hat{\lambda}+\sum_{x,y\in X\cap B}^{\neq}h_{I}\left(\|x-y\|\right)\mathcal{J}_{d/2-1}\left(2\pi\|{k}\|\|x-y\|\right)\|x-y\|^{-(d/2-1)}.$$

(18)

We can also use more than one taper, and will use a sequence of orthogonal functions $\{h_{{m}}\}$ [58], 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_{{m}}||_{2}=1$ and

$$\int_{{\mathbb{R}}^{d}}h_{{m}}(z)h_{{m^{\prime}}}(z)\,dz=\delta_{{{m}}{{m^{\prime}}}}.$$

(19)

We shall study the properties of direct spectral estimators of the form of Eqn (14). This will help us to characterise the second order properties of the process $X$. Most of time series analysis is based on discrete regular sampling for instance, and so most tapers are designed for that scenario. We chose to use continuous tapers, rather than interpolating the points to a regular grid. This leaves as possible tapers to use the spheroidal wavefunctions (continuous but hard to compute), as well as the cosine tapers of [60]. These correspond to separable taper choices; non-separable choices with be discussed in Section 7.

In this Section we derive the asymptotic marginal distribution of the tapered periodogram, and its asymptotic expectation. We shall start from the simplest spectral estimator, namely the Bartlett periodogram as specified by (12) and calculate its properties. We also define the transfer function corresponding to a taper $h$ to be

$$H=\mathcal{F}[h],\quad H_{0}=\mathcal{F}[h_{0}].$$

(20)

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:

$${B}_{\square}(\bm{l})=[-{l_{1}}/{2},{l_{1}}/{2}]\times\dots\times[-{l_{d}}/{2},{l_{d}}/{2}].$$

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}\lambda^{2}T(B,{k})$ has been added. To get there, we need to understand the DFT $J_{h}({k})$ better.

We note that in general $h(x)e^{-2\pi i{k}\cdot 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 [7]. If we compare the quantity in (13) with [7], then we see that $q\in{\mathbb{N}}$ and $p=1$, in their notation. We do not have to worry about $e^{-2\pi i{k}\cdot x}$ being a function of several of the points, and this is why $q=1$. Citing [7], we can deduce from their Theorem 1 that as $|B|$ diverges $J_{h}({k})$ is asymptotically Gaussian. We recover the second moments by direct calculations:

Why do we bother with determining the first two moments of the point pattern? The spectrum characterises any stationary/homogeneous process and so is a key summary to estimate [8]. Estimation will be linear in the periodogram, but once $J_{h}({k})$ is approximately Gaussian, the distribution of the periodogram is determined from these moments, using theory developed for quadratic forms in normal variates [42].

In the above Proposition the term ${\mathrm{Rel}}\left\{\cdot\right\}$ denotes “relation”, also known as the “complimentary covariance”, see e.g. [54]. With these moments, and the result of  [7] we deduce the following corollary.

We can view $J_{h}({k})$ as a complex-valued scalar or a real-valued two-vector. $N_{C}(\mu,\Sigma,C)$ is the general complex normal and its arguments are its mean $\mu$, it covariance $\Sigma$, as well as its relation or complimentary covariance $C$. It should be contrasted with the complex proper normal $N_{C}(\mu,\Sigma),$ that has zero relation.

What can we then say about $I_{h}({k})$? Using the continuous mapping theorem [22] we can deduce from (25) that if we consider only arguments ${k}$ so that the complementary covariance is negligible then

$$|J_{h}({k})|^{2}\overset{L}{\rightarrow}\frac{f({k})}{2}\chi^{2}_{2}\left(\frac{\lambda^{2}|H({k})|^{2}}{f({k})}\right),$$

(26)

where the parameters of the non–central $\chi^{2}_{\nu}(\mu^{2})$ ($\nu$ denoting the degrees of freedom, and $\mu^{2}$ the non-centrality parameter). We give the form of $H({k})$ in (20). We can also use the uniform integrability of the variable $J_{h}({k})$, to get the asymptotic moments of $I_{0}({k})$ from these limits. In fact, the expectation of any member of that family of tapered estimates takes the form of:

$${\mathbf{E}}\left\{I_{h}({k})\right\}=\int_{{\mathbb{R}^{d}}}|H({k}^{\prime}-{k})|^{2}f({k}^{\prime})d{k}^{\prime}+\lambda^{2}|H({k})|^{2}.$$

(27)

This can be derived straightforwardly by direct computation from (14) by taking expectations with Campbell’s formula and using the convolution theorem. We see immediately the bias of this estimator, namely the $\lambda^{2}|H({k})|^{2}$ term. To remove the non-centrality bias, we define the bias-corrected periodogram from the de-biased discrete Fourier transform

$$\widetilde{I}_{h}({k})=\left|\widetilde{J}_{h}({k})\right|^{2}.$$

(28)

Note also the discussion in [20, p. 292], where the theory of signed measures is used to make the equivalent definition, if with additional mathematical sophistication. Therefore $J_{h}({k})-\lambda H({k})=\widetilde{J}_{h}({k})$ can be called the signed measure. The quantity $\widetilde{J}_{h}({k})$ is the DFT of the process $X^{0}$ dual to the mean-corrected random measure $N^{0}(dx)=N(dx)-\lambda dx$ where $N$ is the dual counting measure of the point process $X$. As we have subtracted $H({k})$ 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 ${k}=0$. For the periodogram we have $h_{0}(x)=\mathbf{1}_{B}(x)/\sqrt{|B|}$. Thus we have

$$J_{0}({k})=\frac{1}{\sqrt{|B|}}\sum_{x\in X\cap B}e^{-2\pi i{k}\cdot x},\quad J_{0}(0)=\frac{1}{\sqrt{|B|}}N(B)=\sqrt{\|B\|}\widehat{\lambda}.$$

(29)

This removes the problem of a non-zero mean of the DFT, as long as we assume that we know the intensity $\lambda$, 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|\rightarrow\infty$), by just counting the number of points and dividing by the area. For completeness we also define the signed measure DFT as

$$\widetilde{J}_{h}({k})=J_{h}({k})-\hat{\lambda}H({k})\quad.$$

(30)

We note directly from (25), that

$$\widetilde{J}_{h}({k})\rightarrow N_{C}\{0,f({k}),r({k})\}.$$

(31)

Also it follows

$$\widetilde{J}_{0}(0)=\sqrt{|B|}\widehat{\lambda}-\widehat{\lambda}H_{0}(0)=\sqrt{|B|}\widehat{\lambda}-\sqrt{|B|}\widehat{\lambda}=0,$$

(32)

trivially. Thus we cannot estimate the DFT at wavenumber 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 wavenumber zero. In a time series setting the periodogram at wavenumber zero is often not plotted.

We note that a second way to do bias correction is by means of

$$\breve{I}({k};h)=\sum_{x,y\in X}h(x)h(y)e^{-2\pi i{k}\cdot(x-y)}-\widehat{\lambda^{2}}\left|H({k})\right|^{2}.$$

(33)

The advantage of the first estimator in (28) 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, where as $\hat{\lambda}$ is an unbiased estimator for $\lambda$, $\widehat{\lambda^{2}}$ can be a biased estimator for $\lambda^{2}$.

For additional generality, we could consider a more sophisticated estimator than (33) or (28) see e.g. [14, Ch. 6], and define

$$I_{g}({k})=\sum_{x\in X\cap B}\sum_{y\in X\cap B}g(x,y)e^{-2\pi i{k}\cdot\left(x-y\right)}-\widehat{\lambda^{2}}G({k},{k}),$$

(34)

where $g(x,y)$ is a kernel which may have a specified support, that has to be selected, and $G({k},{k})=\sum_{x}\sum_{y}g(x,y)e^{-2\pi i{k}\cdot\left(x-y\right)}$. This estimator suffers from not being positive by design. For example we could define $g(x,y)=h(x)h^{\ast}(y)$ with $h(x)$ a multi-dimensional data taper  [36], 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}({k})$ may satisfy a number of desirable characteristics such as positivity, asymptotic unbiasedness and computational efficiency.

Looking at the debiased periodogram, we can now determine further properties of $I_{h}({k})$, or properties of $\widetilde{I}_{h}({k})$. To produce an estimator that is smoothed we need to study the covariance and variance of $\widetilde{I}_{h}({k}_{1})$ and $\widetilde{I}_{h}({k}_{2})$. We first look at ${\mathbf{E}}\left\{\widetilde{I}_{h}({k})\right\}$. We find that its form for large spatial regions is specified by the following theorem.

Thus as long as $H({k})$ is getting more concentrated in wavenumber (e.g. $H({k})\rightarrow\delta({k})$), this is asymptotically in $|B|$ an unbiased estimator of $f({k})$. Using the signed measure DFT of (30) to study $X$ is more convenient, as this lets us avoid the contribution that affects the low wavenumbers.

Let us write down what grid we get large sample unbiasedness 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 wavenumber elements is per unit length, such as $m^{-1}$. Referring to Lemma 4.1 we see that the zero wavenumber, or taper $h_{0}$ related bias term $T(B,{k})$, vanishes using this grid not only removing any asymptotic bias at low wavenumbers, but also giving us a sampling of the wavenumbers that approximately leads to independent periodogram ordinates, rather like in the random field case, see Corollary 5.1.

These results establish what wavenumber grid we should evaluate a standard spectral estimator of a time series at, in analogy to the the Fourier frequencies [58, p. 197–198] in time series. 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 wavenumber [58, p. 98]. This does not exist for point processes. It may seem counter intuitive that there is no upper limit to the wavenumbers 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. For random locations of the point process the pairwise distances can be any real-valued value, so the Nyquist wavenumber does not (in a sense) exist.

This theorem establishes that the tapered periodogram is a biased estimator of the spectrum of the point pattern. Previous authors made up an ad-hoc corrections [23] to remove the bias. Our result will suggest a method to arrive at a positive spectral density estimator that is not biased. This is the basic and important result; equivalent to [8, Section V]. In fact, asymptotics are well–understood both for time series and random fields. Our (previous) understanding of asymptotics for spectral estimation of point processes hearkens back to [20] or [23]. In the latter reference the authors informally refer back to Chemistry theory by Hansen and McDonald, discussing the properties of gasses.

This corollary then informs us how to do estimation, by removing the bias. This corollary shows the importance of removing the spectral bias before squaring – otherwise it gets hard to isolate and remove the effect at zero wavenumber, which is exacerbated at higher intensities (growing $\lambda$), as is clear from (36). Note that the estimator $\widetilde{I}_{h}({k})$ is debiased relative to ${\mathbf{E}}\left\{I_{h}({k})\right\}$, but is still guaranteed to be non-negative. This can be compared to removing $\widehat{\mu}$ from a time series before calculating the periodogram.

Finally, to employ smoothing for purposes of variance reduction, we will extend the tapered estimator to include multiple tapers. Define ${M}\geq 1$ estimates of the spectrum via

$$I_{{{{m}}}}^{t}({k})=\sum_{x,y\in X\cap B}h_{{{m}}}(x)e^{-i2\pi{k}\cdot x}h^{\ast}_{{{m}}}(y)e^{i2\pi{k}\cdot y},\quad{{m}}=1,\dots,{M},$$

(38)

which is not bias-corrected but where we assume

$$\int_{{\mathbb{R}^{d}}}h_{{{m}}}(x)h^{\ast}_{{m^{\prime}}}(x)dx=\delta_{{{m}}{{m^{\prime}}}},$$

where $\delta_{{{m}}{{m^{\prime}}}}=1$ if ${{m}}={{m^{\prime}}}$ and $0$ otherwise. Of course (38) still suffers from low-wavenumber bias. To define a debiased estimator we take

$$\widetilde{I}_{{{{m}}}}^{t}({k})=\left(\sum_{x\in X\cap B}h_{{m}}(x)e^{-i2\pi{k}\cdot x}-\lambda H_{{m}}({k})\right)\left(\sum_{y\in X\cap B}h^{\ast}_{{m}}(y)e^{i2\pi{k}\cdot y}-\lambda H^{\ast}_{{m}}({k})\right),\quad{{m}}=1,\dots,{M}.$$

(39)

We subsequently average these estimates over ${{m}}$ to reduce variance.

In (38) we have a product of $h_{{m}}(x)h^{\ast}_{{m}}(y)$. This is an inevitable consequence of the bilinear form of the periodogram from the Fourier transform. The most natural way of making multidimensional and separable tapers would be for $x=\begin{pmatrix}x_{1}&\dots&x_{d}\end{pmatrix}^{T}$ to define the component-wise product taper $h_{{m}}(x)=\tilde{h}_{{{m}}_{1}}(x_{1})\dots\tilde{h}_{{{m}}_{d}}(x_{d})$. We will then need to re-enumerate ${{m}}_{1},\dots,{{m}}_{d}$ into a single index. For instance in 2D if we have three tapers in 1D, then we will end up with nine tapers in 2D, and ${{m}}$ will range from one to nine. These estimates can also be arrived at starting from (13) and $J_{h}({k})$. As we shall use ${M}$ tapers linearly, and ${M}^{2}$ for the quadratic estimators, we shall use the subscript ‘0’ to refer to the untapered periodogram, and ${{m}}=1,\dots,{M}$ 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 [18, 76]. The idea of tapering is to ameliorate edge effects that leads to the asymptotically leading bias term. This is important in 1D [70, 76] but of greater importance in 2D and higher [18]. In 2D and higher the edge effects become more pronounced, and indeed asymptotically dominant [47, 53]. We also use multitapering to stabilise the variance of any estimator, as described for time series in [76]. We shall describe the necessary steps to perform linear estimation of the spectrum of a 2D and higher dimension point process in Section 6.

Let us set theory aside for a moment and consider some spectral estimates for archetypal point patterns. Figure 3 demonstrates the debiased periodogram, the debiased multitapered periodogram, with sine-tapers, and a rotationally averaged 1D summary of the peridogram (defined in Sections 6&7) for four point patterns exhibiting different structural behaviour. The aspect ratio of the figure panels are kept the same only to have an orderly figure: Due to debiasing we can estimate the periodogram on wavenumbers of our own choosing, and not just on the Fourier grid which in these examples would be different for different tapers. 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 $\|{k}\|\approx 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 2). 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 $\|{k}\|=0$ for the non-debiased version, which illustrates why the proposed debiasing step is relevant when applying the method in practice.

In this Section we determine the covariance of the periodogram, in anticipation of smoothing the periodogram. If we compare the developments of the previous sections to that of spectral analysis, it would seem that we are in a good position to estimate the spectrum, especially as we can assume the spectral deviation function $\widetilde{f}({k})$ is smooth. The smoothness of $\widetilde{f}({k})$ follows from the decay of $\widetilde{\rho}(z)$. Simple Fourier theory stipulates that the decay of $\widetilde{\rho}(z)$ yields the smoothness of $\widetilde{f}({k})$. From a mean-square error perspective on estimating the spectrum, to understand how to smooth the periodogram away from zero wavenumber, and to do so we need to further study the variance and covariance of direct spectral estimators.

Our main concern is : 1) Is the variance finite? 2) Can we find a grid of wavenumbers so that the estimated spectrum is uncorrelated at these points? To be able to answer such questions we must study what the variance and covariance of the DFT is. 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

$$\phi_{k}(x)=e^{-2\pi i{k}\cdot x},\quad{k},x\in{\mathbb{R}}^{d}.$$

(40)

With this notation we have

$$I^{t}_{{m}}({k})=\sum_{x,y\in X\cap B}h_{{{m}}}(x)h_{{{m}}}^{\ast}(y)\phi_{k}(x-y),\quad{{m}}=1,\dots M.$$

(41)

Above we note that we have the same value of ${{m}}$ across $x$ and $y$ in (41) as this corresponds to a modulus square.