Correct Estimation of Higher-Order Spectra: From Theoretical Challenges to Practical Multi-Channel Implementation in SignalSnap

TL;DR

This paper introduces unbiased, consistent estimators for higher-order spectra using multivariate k-statistics, implemented in an efficient GPU-accelerated library, enabling practical analysis of large-scale, multi-channel signal data.

Contribution

It reformulates higher-order spectral estimation with modern unbiased estimators, and provides an open-source GPU library for scalable, accurate spectral analysis.

Findings

Eliminates artifacts in higher-order spectral estimates.

Enables analysis of datasets over hundreds of gigabytes within minutes.

Facilitates applications in quantum measurements and non-Gaussian data analysis.

Abstract

Higher-order spectra (Brillinger's polyspectra) offer powerful methods for solving critical problems in signal processing and data analysis. Despite their significant potential, their practical use has remained limited due to unresolved mathematical issues in spectral estimation, including the absence of unbiased and consistent estimators and the high computational cost associated with evaluating multidimensional spectra. Consequently, existing tools frequently produce artifacts, no existing software library correctly implements Brillinger's cumulant-based trispectrum, or fail to scale effectively to real-world data volumes, leaving crucial applications like multi-detector spectral analysis largely unexplored. In this paper, we revisit higher-order spectra from a modern perspective, addressing the root causes of their historical underuse. We reformulate higher-order spectral estimation using recently derived multivariate k-statistics, yielding unbiased and consistent estimators that eliminate spurious artifacts and precisely align with Brillinger's theoretical definitions. Our methodology covers single- and multi-channel spectral analysis up to the bispectrum (third order) and trispectrum (fourth order), enabling robust investigations of inter-frequency coupling, non-Gaussian behavior, and time-reversal symmetry breaking. Additionally, we introduce quasi-polyspectra to uncover non-stationary, time-dependent higher-order features. We implement these new estimators in SignalSnap, an open-source GPU-accelerated library capable of efficiently analyzing datasets exceeding hundreds of gigabytes within minutes. In applications such as continuous quantum measurements, SignalSnap's rigorous estimators enable precise quantitative matching between experimental data and theoretical models.