Draft:Polyspectra

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). ^[4]

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.^[5]

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. ^[1] Fourier coefficients of ${\displaystyle z(t)}$ $${\displaystyle a_{k}={\frac {T}{N}}\sum _{j=0}^{N-1}g_{j}z_{j}e^{2\pi ikj/N}\,,}$$ are calculated from the process ${\displaystyle z_{j}=z(jT/N-t_{0})}$ which is divided into windows of length ${\displaystyle T}$ with ${\displaystyle N}$ points per window. A window function ${\displaystyle g_{j}=g(jT/N)}$ reduces spectral leakage.

Approximate polyspectra are obtained from cumulants of Fourier coefficients: $${\displaystyle {\begin{aligned}&{\text{Average: }}&S_{z}^{(1)}&\approx {\frac {NC_{1}(a_{0})}{T\sum _{j=0}^{N-1}g_{j}}}\\\\[2pt]&{\text{Power spectrum: }}&S_{z}^{(2)}(\omega _{k})&\approx {\frac {NC_{2}(a_{k},a_{k}^{*})}{T\sum _{j=0}^{N-1}g_{j}g_{j}^{*}}}\\\\[2pt]&{\text{Bispectrum: }}&S_{z}^{(3)}(\omega _{k},\omega _{l})&\approx {\frac {NC_{3}(a_{k},a_{l},a_{k+l}^{*})}{T\sum _{j=0}^{N-1}g_{j}^{2}g_{j}^{*}}}\\\\[2pt]&{\text{Trispectrum: }}&S_{z}^{(4)}(\omega _{k},\omega _{l},\omega _{p})&\approx {\frac {NC_{4}(a_{k},a_{l},a_{p},a_{k+l+p}^{*})}{T\sum _{j=0}^{N-1}g_{j}^{3}g_{j}^{*}}}.\\\end{aligned}}}$$ The distances of the spectral positions ${\displaystyle \omega _{k}=2\pi k/T}$ decrease with increasing window length ${\displaystyle T}$ resulting in a higher spectral resolution. The cumulants ${\displaystyle C_{n}}$ can be estimated by the k-statistic which provides unbiased estimators ${\displaystyle c_{n}}$ (estimator).^{[6] [7] [8]} The estimators are given by^[9] $${\displaystyle {\begin{aligned}&c_{2}(x,y)={\frac {m}{m-1}}({\overline {xy}}-{\overline {x}}\,{\overline {y}})\\\\[2pt]&c_{3}(x,y,z)={\frac {m^{2}}{(m-1)(m-2)}}({\overline {xyz}}-{\overline {xy}}\,{\overline {z}}-{\overline {xz}}\,{\overline {y}}-{\overline {yz}}\,{\overline {x}}+2{\overline {x}}\,{\overline {y}}\,{\overline {z}})\\\\[2pt]&c_{4}(x,y,z,w)={\frac {m^{2}}{(m-1)(m-2)(m-3)}}{\big [}(m+1){\overline {xyzw}}\\&\qquad \qquad \qquad \qquad \qquad -(m+1)({\overline {xyz}}\,{\overline {w}}+{\overline {xyw}}\,{\overline {z}}+{\overline {xzw}}\,{\overline {y}}+{\overline {yzw}}\,{\overline {x}})\\&\qquad \qquad \qquad \qquad \qquad -(m-1)({\overline {xy}}\,{\overline {zw}}+{\overline {xz}}\,\,{\overline {yw}}+{\overline {xw}}\,{\overline {yz}})\\&\qquad \qquad \qquad \qquad \qquad +2m({\overline {xy}}\,{\overline {z}}\,{\overline {w}}+{\overline {xz}}\,{\overline {y}}\,{\overline {w}}+{\overline {xw}}\,{\overline {y}}\,{\overline {z}}+{\overline {yz}}\,{\overline {x}}\,{\overline {w}}+{\overline {yw}}\,{\overline {x}}\,{\overline {z}}+{\overline {zw}}\,{\overline {x}}\,{\overline {y}})\\&\qquad \qquad \qquad \qquad \qquad -6m{\overline {x}}\,{\overline {y}}\,{\overline {z}}\,{\overline {w}}{\big ]},\end{aligned}}}$$ where ${\displaystyle {\overline {(..)}}}$ denotes the average of ${\displaystyle m}$ samples.