Polyspectra

Polyspectra, also known as higher-order spectra, are frequency-domain statistical quantities that generalize the spectral density (power spectrum) to higher orders. Polyspectra of a given signal can reveal non-Gaussian behavior, time-inversion asymmetries, and phase correlations between different frequency contributions to the signal. While a standard power spectrum quantifies the distribution of intensity across frequencies, polyspectra characterize higher-order phase and amplitude correlations between frequency components. Higher-order spectra include the third-order bispectrum, the fourth-order trispectrum, and their multivariate generalizations.

Definition

In his original publications, Brillinger considers a stationary process (the signal) $z(t)$ and its multi-time cumulant. $$f({\tau }_1,\cdots ,{\tau }_{n-1})=C_n(z(t),z(t+{\tau }_1),\cdots ,z(t+{\tau }_{n-1})).$$ The cumulant does not depend on $t$ since $z(t)$ is stationary. This allows for a definition of $n$-th order polyspectra $S_z^{(n)}$ of a stationary signal $z(t)$ by $$C_n(z({\omega }_1),\cdots ,z({\omega }_n))=2\pi \delta ({\omega }_1+\cdots +{\omega }_n)S_z^{(n)}({\omega }_1,\cdots ,{\omega }_{n-1})$$ with the Fourier transform $z(\omega )={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{e}^{i\omega \tau }z(\tau )d\tau ,$ and the Dirac delta function $\delta$ which appears due to the independence of $f$ on $t$.

The second order spectrum can be rewritten as $$S_z^{(2)}(\omega )={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{e}^{i\omega \tau }{C}_2(z(t+\tau ),z(t))d\tau ,$$ where $C_2(x,y)=⟨xy⟩-⟨x⟩⟨y⟩$ is the covariance. The representation of $S_z^{(2)}$ is equivalent to the definition of the spectral density of $z$ (power spectrum) as Fourier transform of the autocorrelation function of $z(t)$.

History of polyspectra

Higher-order spectra were formalized in the works of Brillinger starting in the 1950s. Mendel and Nikias contributed the first major reviews on the topic in the early 1990s, resulting in the "Higher Order Spectral Analysis Toolbox" (HOSA), a software library that included the bispectrum and some fourth order statistics. However, Birkelund stated in 2003: "In theory, polyspectra can be applied to solve many important problems in signal processing and data analysis. In practice, however, one has been discouraged by the poor statistical properties of most polyspectral estimators." In 2026, unbiased cumulant-based estimators for higher-order spectra were proposed and analyzed. They are implemented in GPU-accelerated software libraries based on ArrayFire and alternatively based on PyTorch.

Numerics of polyspectra

The definition of polyspectra assumes a stationary process $z(t)$ that is known with an infinite temporal resolution in an infinite temporal interval. Real world data $z(t)$ allow therefore only for the calculation of estimates of polyspectra. These estimates have finite spectral resolution and finite accuracy.

Classical estimators

Direct Fourier methods usually start with the Fourier transform of the entire signal resulting in a high frequency resolution with large variance (error). Subsequently, a smoothing kernel can be applied to trade resolution for variance. Alternatively, spectra of segments of the signal are averaged (Welch-type method).

The lag window method relies on estimating cumulants in the time-domain and subsequent application of Fourier transforms.

Multitaper methods consider the entire signal tapered (windowed) with orthogonal tapers. Spectra for separate tapers are then averaged.

The classical methods rely on cumulant estimators that are only asymptotically unbiased.

Unbiased cumulant-based estimators

The following equations provide an unbiased estimate (estimator) of polyspectra up to fourth order. Fourier coefficients of $z(t)$ $$a_k=\frac{T}{N}\sum _{j=0}^{N-1}{g}_j{z}_j{e}^{2\pi ikj/N},$$ are calculated from the process $z_j=z(jT/N-t_0)$ which is divided into windows of length $T$ with $N$ points per window. A window function $g_j=g(jT/N)$ reduces spectral leakage.

Approximate polyspectra are obtained from cumulants of Fourier coefficients: $$\begin{array}{cccc}& \text{Average: }& S_z^{(1)}& \approx \frac{N{C}_1(a_0)}{T\sum _{j=0}^{N-1}{g}_j}\\ \\ & \text{Power spectrum: }& S_z^{(2)}({\omega }_k)& \approx \frac{N{C}_2(a_k,a_k^{\ast })}{T\sum _{j=0}^{N-1}{g}_j{g}_j^{\ast }}\\ \\ & \text{Bispectrum: }& S_z^{(3)}({\omega }_k,{\omega }_l)& \approx \frac{N{C}_3(a_k,a_l,a_{k+l}^{\ast })}{T\sum _{j=0}^{N-1}{g}_j^2{g}_j^{\ast }}\\ \\ & \text{Trispectrum: }& S_z^{(4)}({\omega }_k,{\omega }_l,{\omega }_p)& \approx \frac{N{C}_4(a_k,a_l,a_p,a_{k+l+p}^{\ast })}{T\sum _{j=0}^{N-1}{g}_j^3{g}_j^{\ast }}.\end{array}$$ The distances of the spectral positions ${\omega }_k=2\pi k/T$ decrease with increasing window length $T$ resulting in a higher spectral resolution. The cumulants $C_n$ can be estimated by the k-statistic which provides unbiased estimators $c_n$ (estimator). The estimators are given by $$\begin{array}{cc}& c_2(x,y)=\frac{m}{m-1}(\overline{xy}-\bar{x}\bar{y})\\ \\ & c_3(x,y,z)=\frac{m^2}{(m-1)(m-2)}(\overline{xyz}-\overline{xy}\bar{z}-\overline{xz}\bar{y}-\overline{yz}\bar{x}+2\bar{x}\bar{y}\bar{z})\\ \\ & c_4(x,y,z,w)=\frac{m^2}{(m-1)(m-2)(m-3)}[(m+1)\overline{xyzw}\\ & -(m+1)(\overline{xyz}\bar{w}+\overline{xyw}\bar{z}+\overline{xzw}\bar{y}+\overline{yzw}\bar{x})\\ & -(m-1)(\overline{xy}\overline{zw}+\overline{xz}\overline{yw}+\overline{xw}\overline{yz})\\ & +2m(\overline{xy}\bar{z}\bar{w}+\overline{xz}\bar{y}\bar{w}+\overline{xw}\bar{y}\bar{z}+\overline{yz}\bar{x}\bar{w}+\overline{yw}\bar{x}\bar{z}+\overline{zw}\bar{x}\bar{y})\\ & -6m\bar{x}\bar{y}\bar{z}\bar{w}],\end{array}$$ where $\overline{(\mathrm{..})}$ denotes the average of $m$ samples.

Applications

Bispectra find medical applications in the analysis of electroencephalograms and are used in engineering for fault detection in rotating machinery components. Polyspectra have also become a tool in quantum physics for analysing records of quantum measurements. The applications of quantum polyspectra range from quantum transport over single photon detection to spin noise measurements.