Schuster Periodogram Essentials

• The Schuster periodogram is a nonparametric estimator that computes spectral power by correlating time series data with sine and cosine functions.
• It bridges time-domain autocorrelation with Fourier-based spectral density estimation, enabling the assessment of periodic signals and frequency leakage.
• Adaptations for irregular sampling and clustered events enhance its reliability in applications across astronomy, seismology, and other scientific fields.

The Schuster periodogram, synonymous with the "classical periodogram" (CP), is a foundational nonparametric estimator for the frequency content in time series data. It quantifies the squared correlation between observed data and complex exponentials or sinusoids at each trial frequency, delivering a spectrum of power as a function of frequency. Its original formulation dates to 1898 and forms the conceptual basis for a broad range of spectral analysis methodologies, including extensions for handling uneven sampling, statistical significance evaluation, and practical adaptations across scientific domains such as astronomy and seismology (Vio et al., 2018, &&&1&&&).

1. Mathematical Formulation

Given real-valued observations $\{x(t_0), x(t_1), \ldots, x(t_{M-1})\}$ at time points $\{t_j\}$ (which may be uniformly or irregularly spaced), the Schuster periodogram at trial frequency $\nu$ is defined by

$$a(\nu) = \frac{1}{M} \sum_{j=0}^{M-1} x(t_j) \cos(2\pi \nu t_j), \quad b(\nu) = \frac{1}{M} \sum_{j=0}^{M-1} x(t_j) \sin(2\pi \nu t_j)$$

$P_x(\nu) = a(\nu)^2 + b(\nu)^2$

For regularly spaced $t_j = j\Delta t$ and discrete Fourier frequencies $\nu_k = k/(M\Delta t)$, it reduces to the squared modulus of the Discrete Fourier Transform (DFT):

$$P_x(k) = \frac{1}{M} \left| \sum_{j=0}^{M-1} x_j \exp(-2\pi i k j/M) \right|^2$$

In the context of event data, such as earthquake catalogues, the Schuster distance and periodogram are given as

$$S(\omega) = \sum_{i=1}^N e^{i\omega t_i}, \quad D^2(\omega) = |S(\omega)|^2$$

(Park et al., 2021).

2. Connection to Sample Autocorrelation and Theoretical Spectrum

For a zero-mean stationary process $x(t)$ with autocorrelation $R_x(\tau) = E[x(t)x(t+\tau)]$, the (theoretical) spectral density $S_x(\nu)$ is the Fourier transform of $R_x$:

$$S_x(\nu) = \int_{-\infty}^{\infty} R_x(\tau) e^{-2\pi i \nu \tau} d\tau$$

In discrete time, the finite-lag spectral estimate using the sample autocovariance $C(k)$ is

$$C(k) = \frac{1}{M-|k|} \sum_{j=0}^{M-|k|-1} x(t_j) x(t_{j+|k|}),\quad \hat S_x(\nu) = \sum_{k=-(M-1)}^{M-1} C(k) e^{-2\pi i \nu k \Delta t}$$

Under regular sampling, this formulation is algebraically equivalent to the periodogram definition above. In the irregularly sampled case, the periodogram can be viewed as a discrete analog to the direct correlation of data with sinusoids, providing a measure of signal strength at each $\nu$ (Vio et al., 2018).

3. Statistical Properties and Null Distributions

Regular Sampling (White Noise Null)

Under the null hypothesis that $x(t_j)$ are independent Gaussian noise with variance $\sigma_n^2$, the coefficients $a(\nu_k)$, $b(\nu_k)$ at DFT frequencies $\nu_k$ are independent $\mathcal{N}(0, \sigma_n^2/2)$ random variables. Thus,

$$P_x(\nu_k)/\sigma_n^2 \sim \mathrm{Exp}(1),\quad E[P_x(\nu_k)] = \sigma_n^2,\quad \mathrm{Var}[P_x(\nu_k)] = \sigma_n^4$$

Probabilities for family-wise errors (maximum over $N^*=M/2$ positive frequencies) are estimated as

$$\mathrm{Prob} \left[ \max_k P_x(\nu_k) > z \right] = 1 - [1-\exp(-z/\sigma_n^2)]^{N^*}$$

Irregular Sampling

For arbitrary $\{t_j\}$, the cosine and sine projections are not orthogonal, and their variances and covariance become frequency- and sampling-dependent:

$$\mathrm{Var}[a(\nu)] = \sigma_n^2 A(\nu),\quad \mathrm{Var}[b(\nu)] = \sigma_n^2 B(\nu),\quad \mathrm{Cov}[a(\nu), b(\nu)] = \sigma_n^2 C(\nu)$$

with $A(\nu), B(\nu), C(\nu)$ given by quadratic sums over trigonometric functions. The PDF of $P_x(\nu)$ is then a generalized form involving the modified Bessel function $I_0$:

$$f_P(z) = \frac{1}{2\sigma_+\sigma_-\sqrt{1-\rho^2}} \exp(-\alpha z) I_0(\beta z),\quad z\geq 0$$

where parameters ($\sigma_\pm, \rho, \alpha, \beta$) are defined in terms of these variances and covariances (Vio et al., 2018).

Event Data—Poisson and Clustering Effects

For event times (e.g., earthquake origin times), under the homogeneous Poisson null hypothesis, the Schuster spectrum $D^2(\omega)/N$ approaches an exponential distribution, so

$$p(\omega) = \exp\left[-D^2(\omega)/N\right] \approx \textrm{Uniform}[0, 1]$$

When events form clusters (e.g., aftershocks), the expected mean of $D^2(\omega)$ becomes frequency dependent and determined by the aftershock kernel, requiring explicit correction (see Section 5) (Park et al., 2021).

4. Handling Irregular Sampling and Spectral Window

Irregular sampling destroys the mutual orthogonality of trigonometric components and alters the statistical structure of the periodogram:

• The variances $A(\nu)$, $B(\nu)$ may differ, and covariance $C(\nu)$ may be nonzero, especially for low frequencies or gappy data.
• The periodogram ordinates at different frequencies become statistically dependent.
• FFT-based computation is generally not applicable; direct evaluation scales as $O(MN)$, though Non-Uniform FFT (NUFFT) methods provide an $O((M+N)\log(M+N))$ alternative.
• The spectral window $W(\nu) = \sum_j \exp(-2\pi i \nu t_j)$ characterizes frequency leakage and interprets sidelobes and frequency resolution.

Significance thresholds using the exponential approximation can be anti-conservative unless $C(\nu)$ is small across the search grid. For precise inference, the exact null distribution or empirical effective degrees of freedom (rank of the covariance matrix of periodogram bins) must be established (Vio et al., 2018).

5. Comparison with the Lomb–Scargle Periodogram

The Lomb–Scargle periodogram (LSP) modifies the classical approach by normalizing both projections to have equal variance and decorrelating their covariance at each frequency via an adaptive time shift $\tau$, so that all bins share an $\mathrm{Exp}(1)$ null-distribution. Its definition is

$$\tilde a(\nu) = \frac{\sum_j x_j \cos 2\pi\nu(t_j-\tau)}{\sqrt{ \sum_j \cos^2 2\pi\nu(t_j-\tau) } },\quad \tilde b(\nu) = \frac{\sum_j x_j \sin 2\pi\nu(t_j-\tau)}{\sqrt{ \sum_j \sin^2 2\pi\nu(t_j-\tau) } },\ P_{LS}(\nu) = \tilde a(\nu)^2 + \tilde b(\nu)^2$$

Key differences and consequences:

• LSP ensures identically distributed Exp(1) nulls at all frequencies, simplifying single-bin false-alarm calculations.
• This normalization distorts peak heights and inter-bin relationships when variances and covariances deviate from ideal values, potentially suppressing true signals in bins with strong phase correlations (e.g., frequency bins overlapping large, repeating gaps).
• The classical periodogram, while requiring more involved statistical characterization under irregular sampling, remains directly interpretable as the data–sinusoid correlation and reliably represents spectral content except in rare pathological cases (Vio et al., 2018).

6. Adaptations for Clustered Events: The Schuster Spectrum Test in Seismology

Standard Schuster spectrum testing (SST) for periodicity assumes a homogeneous Poisson process for event times. In earthquake catalogs, aftershock clustering violates this assumption, leading to inflated false positives. To correct this, the null expectation of the periodogram must incorporate the aftershock model. The expected value becomes

$$E[D^2(\omega)] = \nu_0 \nu + \nu_0|1 + \nu \chi(\omega)|^2$$

where $\chi(\omega)$ is the characteristic function of the aftershock waiting-time density $\lambda(t)/\nu$, and $\nu_0$, $\nu$ are primary and aftershock rates respectively (Park et al., 2021).

For good statistical control, a quantile regression model (QRM) is used to estimate the frequency-varying mean of $D^2(\omega)$ under the null (aftershock clusters modeled), producing corrected p-values as $$p(\omega) = \exp[-D^2(\omega)/\hat{E}[D^2(\omega)]]$$. This yields uniform null p-values and consistent detection performance, even in the presence of clustering, without the need for explicit declustering of the catalog (Park et al., 2021).

7. Implementation, Limitations, and Applications

Efficient computation depends on the sampling scheme. For regular grids, FFT delivers $O(M\log M)$ complexity. For irregular or event data, either direct evaluation or NUFFT is necessary. Smoothing or multitaper methods are often applied to reduce variance and bias caused by the spectral window.

Key limitations:

• The Schuster periodogram is a biased estimator: its mean is a convolution of the true spectrum with the squared modulus of the spectral window, $|W(\nu)|^2$.
• Frequency leakage from sidelobes may obscure weak signals near strong peaks.
• Extensive missing data or heavy clustering can violate null approximations and degrade significance control.

Applications span astrophysics and seismology, with the method extended to robustly test periodicity in clustered catalogues, as exemplified in synthetic earthquake fields and real catalogs from the New Madrid Seismic Zone and the Nepal mid-crustal cluster (Vio et al., 2018, Park et al., 2021).

Summary Table: Schuster Periodogram—Core Features and Contrasts

Aspect	Schuster/CP Definition	Implications
Mathematical Form	$a(\nu), b(\nu)$ projections, $P_x(\nu)=a^2+b^2$	Direct measure of power at $\nu$
Null Distribution	Exponential (regular), Bessel-form (irregular)	Controls false-positives with caveats
FFT Applicability	Yes (regular only); otherwise direct sum/NUFFT	O(M log M), or O(MN)/O((M+N) log(M+N))
Connection	Direct with autocorrelation and spectral window	Interpretable spectral leakage structure
Extension	QRM/mean-correction for cluster events (seismology)	Robust to aftershock-driven clustering

The Schuster periodogram remains the fundamental estimator for spectral analysis in nonparametric settings, with well-understood statistical properties under white noise and ongoing adaptations for irregular sampling and clustered event processes (Vio et al., 2018, Park et al., 2021).