The Bispectrum and Its Relationship to Phase-Amplitude Coupling

TL;DR

This paper explores the relationship between the bispectrum and phase-amplitude coupling (PAC) in biological signals, revealing their connection and limitations, and proposing a bispectral framework for analyzing cross-frequency interactions.

Contribution

It demonstrates the fundamental link between bispectrum and PAC measures, highlights their limitations, and introduces a bispectral framework for better interpretation of cross-frequency coupling.

Findings

Bispectrum and PAC are closely related, with PAC being a smoothed bispectral estimate.

Conventional PAC measures have bias, symmetry issues, and spectral constraints.

Bispectral analysis can distinguish nested oscillations from transient signal features.

Abstract

Most biological signals are non-Gaussian, reflecting their origins in highly nonlinear physiological systems. A versatile set of techniques for studying non-Gaussian signals relies on the spectral representations of higher moments, known as polyspectra, which describe forms of cross-frequency dependence that do not arise in time-invariant Gaussian signals. The most commonly used of these employ the bispectrum. Recently, other measures of cross-frequency dependence have drawn interest in EEG literature, in particular those which address phase-amplitude coupling (PAC). Here we demonstrate a close relationship between the bispectrum and popular measures of PAC, which we relate to smoothings of the signal bispectrum, making them fundamentally bispectral estimators. Viewed this way, however, conventional PAC measures exhibit some unfavorable qualities, including poor bias properties, lack of correct symmetry and artificial constraints on the spectral range and resolution of the estimate. Moreover, information obscured by smoothing in measures of PAC, but preserved in standard bispectral estimators, may be critical for distinguishing nested oscillations from transient signal features and other non-oscillatory causes of "spurious" PAC. We propose guidelines for gauging the nature and origin of cross-frequency coupling with bispectral statistics. Beyond clarifying the relationship between PAC and the bispectrum, the present work lays out a general framework for the interpretation of the bispectrum, which extends to other higher-order spectra. In particular, this framework holds promise for the detailed identification of signal features related to both nested oscillations and transient phenomena. We conclude with a discussion of some broader theoretical implications of this framework and highlight promising directions for future development.

Figures

Equations

$$\mathrm{E}\left[X(t_1)X(t_2)...X(t_k)\right]$$

$${\hat{\mu }}_2(\tau )=\frac{1}{T}{\int }_0^{T}X(t)X(t+\tau )\mathrm{dt}$$

$$X(t)$$

$$\mathcal{F}\left\{X(t+\Delta t)\right\}=\tilde{X}(\omega )e^{i\omega \Delta t}$$

$$\begin{array}{rl}{\mu }_2(\tau )& =\mathrm{E}\left[X(t+\Delta t)X(t+\Delta t+\tau )\right]\\ & =\frac{1}{2\pi }\int \mathrm{E}\left[\tilde{X}({\omega }_1)\tilde{X}({\omega }_2)\right]e^{i({\omega }_1+{\omega }_2)t+i{\omega }_2\tau }\mathrm{E}\left[e^{i({\omega }_1+{\omega }_2)\Delta t}\right]d{\omega }_{1}d{\omega }_2\\ & =\frac{1}{2\pi }\int \mathrm{E}\left[\tilde{X}({\omega }_1)\tilde{X}({\omega }_2)\right]e^{i({\omega }_1+{\omega }_2)t+i{\omega }_2\tau }\delta ({\omega }_1+{\omega }_2)d{\omega }_{1}d{\omega }_2\\ & =\frac{1}{2\pi }\int \mathrm{E}\left[{\left|\tilde{X}({\omega }_2)\right|}^2\right]e^{i{\omega }_2\tau }d{\omega }_2\end{array}$$

$$\begin{array}{rl}{\tilde{\mu }}_K& =\mathrm{E}\left[\tilde{X}({\omega }_1)\tilde{X}({\omega }_2)\dots \tilde{X}({\omega }_K)\right]\end{array}$$

$$\sum _{i=1}^K{\omega }_i=0$$

$$\begin{array}{rl}W(\omega ,t)& =\int X\left(t+\frac{\tau }{2}\right)X\left(t-\frac{\tau }{2}\right)e^{-i\omega \tau }\mathrm{d\tau }\end{array}$$

$$B({\omega }_1,{\omega }_2)=\mathrm{E}\left[\tilde{X}({\omega }_1)\tilde{X}({\omega }_2){\tilde{X}}^{\ast }({\omega }_1+{\omega }_2)\right]$$

$$\begin{array}{r}B=\mathrm{E}\left[\iiint X(r)X(s)X(t)e^{-i{\omega }_1(r-t)-i{\omega }_2(s-t)}\mathrm{dr}\mathrm{ds}\mathrm{dt}\right]\\ =\iiint \mathrm{E}\left[X({\tau }_1+t)X({\tau }_2+t)X(t)\right]e^{-i{\omega }_1{\tau }_1-i{\omega }_2{\tau }_2}d{\tau }_{1}d{\tau }_2\mathrm{dt}\end{array}$$

$$\begin{array}{rl}{W}_3& \left({\omega }_1,{\omega }_2,t\right)=\\ & \iint X\left(\frac{2}{3}{\tau }_1-\frac{1}{3}{\tau }_2+t\right)X\left(\frac{2}{3}{\tau }_2-\frac{1}{3}{\tau }_1+t\right)X\left(t-\frac{1}{3}\left({\tau }_1+{\tau }_2\right)\right)e^{-i{\omega }_1{\tau }_1-i{\omega }_2{\tau }_2}d{\tau }_{1}d{\tau }_2\end{array}$$

$$\mathrm{E}\left[X({\tau }_1+t)X({\tau }_2+t)X(t)\right]={\mu }_3({\tau }_1,{\tau }_2)$$

$$\begin{array}{rl}B& =\iint {\mu }_3({\tau }_1,{\tau }_2)e^{-i{\omega }_1{\tau }_1-i{\omega }_2{\tau }_2}d{\tau }_{1}d{\tau }_2\end{array}$$

$$\tilde{X}(\omega )=\delta (\omega -\theta )+\delta (\omega -\gamma )+\delta (\omega -\gamma -\theta )/2+\delta (\omega -\gamma +\theta )/2$$

$${\tilde{f}}_{\mathrm{FO}}(\omega ){\tilde{f}}_{\mathrm{FO}}^{\ast }(\omega +\xi )\approx {\left|{\tilde{f}}_{\mathrm{FO}}(\omega )\right|}^2{e}^{-i\xi \tau }$$

$$\tilde{f}({\omega }_1)\tilde{f}({\omega }_2){\tilde{f}}^{\ast }({\omega }_1+{\omega }_2)\approx {\tilde{f}}_{\mathrm{SO}}({\omega }_1){\left|{\tilde{f}}_{\mathrm{FO}}({\omega }_2)\right|}^2{e}^{i{\omega }_1\tau }$$

$$\begin{array}{rl}\mathrm{E}\left[\tilde{X}({\omega }_1)\tilde{X}({\omega }_2){\tilde{X}}^{\ast }({\omega }_1+{\omega }_2)\right]& =\\ & ⟨{\tilde{f}}_{\mathrm{SO}}({\omega }_1){\tilde{f}}_{\mathrm{FO}}({\omega }_2){\tilde{f}}_{\mathrm{FO}}^{\ast }({\omega }_1+{\omega }_2)⟩\\ & +⟨{\tilde{f}}_{\mathrm{SO}}({\omega }_1)⟩⟨{\tilde{f}}_{\mathrm{FO}}({\omega }_2){\tilde{f}}_{\mathrm{FO}}^{\ast }({\omega }_1+{\omega }_2)⟩{\tilde{\mu }}_2({\omega }_1)\\ & +⟨{\tilde{f}}_{\mathrm{FO}}({\omega }_2)⟩⟨{\tilde{f}}_{\mathrm{SO}}({\omega }_1){\tilde{f}}_{\mathrm{FO}}^{\ast }({\omega }_1+{\omega }_2)⟩{\tilde{\mu }}_2({\omega }_2)\\ & +⟨{\tilde{f}}_{\mathrm{SO}}({\omega }_1){\tilde{f}}_{\mathrm{FO}}({\omega }_2)⟩⟨{\tilde{f}}_{\mathrm{FO}}^{\ast }({\omega }_1+{\omega }_2)⟩{\tilde{\mu }}_2({\omega }_1+{\omega }_2)\\ & +⟨{\tilde{f}}_{\mathrm{SO}}({\omega }_1)⟩⟨{\tilde{f}}_{\mathrm{FO}}({\omega }_2)⟩⟨{\tilde{f}}_{\mathrm{FO}}^{\ast }({\omega }_1+{\omega }_2)⟩{\tilde{\mu }}_3({\omega }_1,{\omega }_2)\end{array}$$

$$\begin{array}{rl}\mathrm{E}& \left[\tilde{X}({\omega }_1)\tilde{X}({\omega }_2){\tilde{X}}^{\ast }({\omega }_1+{\omega }_2)\right]\\ & =⟨{\tilde{f}}_{\mathrm{SO}}({\omega }_1){\tilde{f}}_{\mathrm{FO}}({\omega }_2){\tilde{f}}_{\mathrm{FO}}^{\ast }({\omega }_1+{\omega }_2)⟩+⟨{\tilde{f}}_{\mathrm{SO}}({\omega }_1)⟩⟨{\tilde{f}}_{\mathrm{FO}}({\omega }_2){\tilde{f}}_{\mathrm{FO}}^{\ast }({\omega }_1+{\omega }_2)⟩{\tilde{\mu }}_2({\omega }_1)\\ & \approx ⟨{\tilde{f}}_{\mathrm{SO}}({\omega }_1){\left|{\tilde{f}}_{\mathrm{FO}}({\omega }_2)\right|}^2{e}^{-i{\omega }_1\Delta \tau }⟩+⟨{\tilde{f}}_{\mathrm{SO}}({\omega }_1)⟩⟨{\left|{\tilde{f}}_{\mathrm{FO}}({\omega }_2)\right|}^2{e}^{-i{\omega }_1\Delta \tau }⟩{\tilde{\mu }}_2({\omega }_1)\end{array}$$

$$f_{\mathrm{FO}}=h(t)\cos (\beta (t+\tau ))$$

$$\mathrm{E}\left[{\tilde{f}}_{\mathrm{FO}}({\omega }_1){\tilde{f}}_{\mathrm{FO}}({\omega }_2){\tilde{f}}_{\mathrm{FO}}^{\ast }({\omega }_1+{\omega }_2)\right]=\tilde{h}({\omega }_1-\beta )\tilde{h}({\omega }_2-\beta ){\tilde{h}}^{\ast }({\omega }_1+{\omega }_2-\beta )\mathrm{E}\left[e^{i\beta \tau }\right]$$

$$\begin{array}{rl}I(\tau ,{\omega }_2)& ={\int }_{SO}B({\omega }_1,{\omega }_2)e^{i{\omega }_1\tau }d{\omega }_1\\ & \approx f_{\mathrm{SO}}(\tau +\Delta \tau ){\left|{\tilde{f}}_{\mathrm{FO}}({\omega }_2)\right|}^2+{\tilde{f}}_{\mathrm{SO}}({\omega }_2)e^{{\omega }_2\Delta \tau }{f}_{\mathrm{FO}}^{(2)}(\tau )\end{array}$$

$$\mathrm{E}\left[{\tilde{X}}^2({\omega }_1)\tilde{X}({\omega }_2){\tilde{X}}^{\ast }(2{\omega }_1+{\omega }_2)\right]$$

$$\mathrm{E}\left[{\tilde{X}}^k({\omega }_1)\tilde{X}({\omega }_2){\tilde{X}}^{\ast }(k{\omega }_1+{\omega }_2)\right]$$

$$c_1=\cos \left(\theta /2t\right)\ast e^{i(\gamma -\theta /2)t}\mathrm{and}{c}_2=\cos (\theta /2t)\ast e^{i(\gamma +\theta /2)t}$$

$$f_{\mathrm{AM}}(t)=c_1+c_2=2{\cos }^2(\theta /2t)e^{i\gamma t}$$

$${\tilde{f}}_{\mathrm{FM}}(\omega )=\delta (\omega -\gamma )-\delta (\omega -\gamma -\theta )/2+\delta (\omega -\gamma +\theta )/2$$

$$s_1=i\sin \left(\theta /2t\right)\ast e^{i(\gamma -\theta /2)t}\mathrm{and}{c}_2=\cos (\theta /2t)\ast e^{i(\gamma +\theta /2)t}$$

$$\left|X_{\mathrm{FM}}(t)\right|=\left|s_1+c_2\right|=\sqrt{1+{\sin }^2(\theta t)}$$

$$S_q[k,m]=\sum _{n=1}^N{h}_{qk}[n]x[n-m]\exp \left[-i\omega [k]\frac{n}{N}\right]$$

$$B[j,k]=\sum _{m=0}^{M-1}{S}_1[j,m]S_2[k,m]S_3^{\ast }[l,m]$$

Full text

We propose guidelines for gauging the nature and origin of cross-frequency coupling with bispectral statistics. Beyond clarifying the relationship between PAC and the bispectrum, the present work lays out a general framework for the interpretation of the bispectrum, which extends to other higher-order spectra. In particular, this framework holds promise for the detailed identification of signal features related to both nested oscillations and transient phenomena. We conclude with a discussion of some broader theoretical implications of this framework and highlight promising directions for future development.

1 Introduction

Many interesting properties of signals in nature relate to nonlinear, non-Gaussian and non-stationary dynamics, which are poorly indexed by second-order measures such as power and cross spectra. Because the spectrum of a stationary Gaussian process lacks any statistical dependence across frequencies, measures of frequency-domain dependence can be extremely useful for gauging the presence and nature of higher order dynamics [51]. For stationary signals, higher order spectra, “polyspectra,” which capture such dependence, are the frequency-domain representations of higher moments [43], in which they mirror the relationship between the power spectrum of a signal and its autocorrelation. The bispectrum is the third-order polyspectrum, making it an obvious place to start in the approach towards higher-order dynamics, at least from a statistical standpoint. Bispectral analysis has proved its practical worth in applications to EEG [18, 53], most notably in gauging the depth of anesthesia [5, 21, 33, 42]. Nevertheless, a significant drawback for the non-statistician remains its lack of an obvious connection to any simple physical interpretation [19]. This point has been the source of some confusion in applied literature; for example, as recently reviewed by Hyafil [28], various authors have suggested incorrectly that bicoherence relates to phase entrainment across pairs of frequencies. While it is true that bispectral measures cannot easily be reduced to any single interpretation, a goal of the present work is to show that specific forms of dependence leave easily recognized signatures in the bispectrum, making bispectral measures invaluable for distinguishing between a variety of phenomena related to cross-frequency coupling.

A separate body of work has recently emerged from EEG literature, which examines the role of another form of cross-frequency coupling, phase-amplitude coupling (PAC) [12, 29]. PAC refers to dependence between analytic amplitude at one frequency and analytic phase at another. Rather than any statistical first principle, interest in PAC has been motivated by the empirical discovery of PAC in signals recorded from the brain, alongside emerging computational and physiological models of how it arises within populations of interacting neurons [61, 32, 2, 29]. In this literature, the measurement of PAC is usually approached with second-order statistics, such as coherence, applied towards comparing analytic phase in one band with analytic amplitude extracted from another, treating the two as separate signals [49, 12, 56]. The situation with PAC is therefore reversed from bispectral measures: its physical meaning is more evident than its relationship to any general body of statistical theory.

The present work aims to bridge this gap and resolve ambiguities of meaning in both directions. The central result is that second-order measures of PAC may be fundamentally understood as estimators of the bispectrum. In the same way that windowed stationary estimates of the power spectrum can be equated to smoothings of the true signal spectrum, windowed bispectral estimators, which include those underlying measures of PAC, amount to different ways of smoothing the true signal bispectrum [22, 44, 55, 25]. In both cases, differences between estimators relate to properties of the respective smoothing kernels [14]. This observation demonstrates conclusively the meaning of PAC measures as they relate to the bispectrum and vice versa and establishes that second-order measures of PAC provide no unique information beyond what can be obtained from the bispectrum.

While PAC measures are fundamentally measures of the bispectrum, the reverse is not true; it is not correct to conclude that the bispectrum is principally a reflection of phase-amplitude coupling. Forms of “spurious” PAC, related, for example, to spectrally broad signal features, may be traced to the bispectral nature of PAC measures. One practical implication is that the superior resolution and lower bias of standard bispectral measures, in comparison to PAC measures, allows them to retain information critical for distinguishing between nested oscillations and other sources of apparent phase-amplitude coupling [37]. Following a brief review of bispectral and PAC estimation, we will observe how different regions of the bispectrum may be taken to reflect either phase-amplitude coupling or consistency of phase across a range of frequencies. It is shown that the bispectrum can be highly useful for ascertaining the presence and origin of phase-amplitude coupling. Properly applied and interpreted, bispectral statistics may overcome a number of recently highlighted limitations and ambiguities of existing PAC measures[4, 28, 37, 41, 50, 57]. For example, with conventional measures of PAC, the observable range of phase-providing frequencies is restricted by the bandwidth of the amplitude-providing band [4], yet no such limitation applies to bispectral estimates. In light of the relationship to the bispectrum, it becomes clear that this constraint is an artifact of the estimator rather than anything inherent in the quantity measured.

Many of the questions that arise in considering the relationship between PAC and the bispectrum prove to be of much more general relevance for understanding a range of signal properties that are neglected by traditional spectral measures, a topic whose importance is becoming increasingly clear [15]. In particular, we develop a model that uses the bispectrum to capture spectrally complex signal features, of which PAC is only one example. The concluding sections are devoted to a prospective review of the application of the bispectrum towards understanding the nature and functional significance of non-oscillatory and transient sources of cross-frequency coupling, beyond nested oscillations.

1.1 Organization

The overall aim of the current work is threefold; first, a general introductory background is provided to motivate applications of higher-order spectra in signal analysis, accompanied in the appendices by a more focused and technical review of the bispectrum and its estimation. The second aim is to describe a formal equivalence between bispectral estimators and measures of phase-amplitude coupling, details of which are presented in C. Finally, building on this formal relationship, a framework is developed to guide the interpretation of the bispectrum. In most places, more technical development has been left to the appendices, the results of which are summarized in the main text alongside some background explanation.

Section 2, in particular, focuses on introducing readers who have had little exposure to the ideas and applications of higher-order spectra to some general motivating principles, thus it deals with the subject very broadly. Section 3 gives a short summary of the main relevant properties of the bispectrum described in more detail in A and B. Section 4 similarly introduces phase-amplitude coupling, which is given a more detailed treatment in C. Section 5 develops a framework for understanding the relationship between PAC and the bispectrum, building on the proof provided in C. Section 6.2 briefly considers an extension to phase-frequency coupling. The concluding sections set the framework within a broader theoretical context; section 8.1 describes a model of signal encoding and transmission that exploits properties of the bispectrum, while section 8.2 anticipates further development of the framework, especially in the direction of feature extraction and identification.

2 A primer on higher-order spectra in signal processing

Abstractly, signal processing is about making sense of sequential or otherwise meaningfully ordered quantities; it considers the question, what kinds of measures best capture the essential nature of any dependence among the quantities as it relates to their position in the ordering space (henceforth assumed to be time)? An endless variety of different measures can be conceived, the usefulness and interpretation of any of which naturally depends on very concrete particulars of the application. But to a theoretically-minded statistician, armed with the hammer of distribution theory, a signal is a nail drawn from a distributional pile, and signal processing is merely the application of a general form of statistical inference to the special case of ordered data. Our slightly naive theoretician might be forgiven for turning to measures related to moments, because, at least in principle, moments completely characterize a distribution, meaning that they collectively harbor all the information there is to harbor about the process generating the signal, irrespective of its source or nature.

Moments are most familiar in the form of the mean, variance, skewness and kurtosis of univariate distributions, related to the expected values of the first four powers of the random variable. The idea generalizes to the expected value of any power; the $k^{\mathrm{th}}$-order moment is the expectation of the $k^{\mathrm{th}}$ power. The extension to multivariate distributions adds cross-terms to the picture, which are most commonly encountered in the off-diagonal values of a covariance matrix. At higher orders, cross-terms within the $k^{\mathrm{th}}$ moment may take the form of the product of $k$ distinct variates or any combination of powers of the variates that sum to $k$. For example, the $k^{\mathrm{th}}$ moment of a distribution over $(x,y,z)$ includes all of the combinations of powers that appear in the expansion of $(x+y+z)^{k}$.

In applications to signal processing, the sampling space is neither univariate nor in general multivariate, rather it is the space of possible signals, which, for theoretical purposes, occupies a continuously infinite number of dimensions. This fact causes moments to take on a decidedly intimidating air in signal-processing contexts, but the truth is that they behave more-or-less exactly like their tamely finite cousins. The ordering space of the signal (in our case, time) takes the place of the integer-valued index of multivariate distributions. In analogy to the cross-terms of multivariate distribution, moments of a continuous signal involve the expected values of the products of the signal with itself at specified points in time. For example, if $X(t)$ is a bounded signal, meaning here the relevant moments are defined, its first moment at time $t_{1}$ is the expectation of the signal itself, $\mathrm{E}\left[X(t_{1})\right]$, the second moment at $(t_{1},t_{2})$ is $\mathrm{E}\left[X(t_{1})X(t_{2})\right]$, while the expectation of the product

[TABLE]

is the $k^{\mathrm{th}}$ moment at $\left(t_{1},..,t_{k}\right)$.

In spite of their considerable theoretical appeal, the application of moment-related statistics to real-world problems in signal processing invariably runs into two obstacles. First, their estimation puts severe demands on data because moments may assume a set of distinct values over a parameter space that grows exponentially with order, starting from one that is already at least as large as the observed signal. Second, the question of how to interpret higher moments also seems to become exponentially less clear with increasing order. It is therefore impossible to estimate moment-related statistics without some judicious assumptions that rein in the parameter space, but the exercise is still bound to prove hollow if one lacks any insight into real-world meaning.

The first point is addressed easily enough by making assumptions about how the dependence relates to ordering; for example, windowing is commonly used to enforce the assumption that any dependence grows negligible outside some local neighborhood according to the window. Combined with the assumption of time invariance (reviewed in the following section) windowing and similar procedures make estimation of higher moments a tractable problem. Even so, the practical application of moment-related statistics in signal-processing contexts rarely involves orders above 4, while applications involving 3rd order dependence, of the type that will be given special attention here, already border on exotic. Moment-related statistics of order 1 and 2 are perfectly familiar as mean and covariance or autocorrelation; the leap from 2 to 3 introduces qualitatively new properties with potentially far-reaching implications, a few of which we will encounter in the following discussion.

2.1 Time invariance

A common and sometimes necessary simplifying assumption is that the statistics governing the signal remain constant over time. A very useful byproduct of time invariance, and motivation for giving special consideration to invariant measures, is that estimation may resort to averaging over time; for example, the estimate of the fixed first moment becomes the average $\hat{\mu}_{1}=\frac{1}{T}\int_{0}^{T}{X(t)\mathop{dt}}$, which ought to converge with the true expectation, $\mu_{1}$, as $T$ increases (provided again the signal is suitably bounded). Similarly, time invariance allows the second moment to vary only with respect to the delay between two points in time, $\mu_{2}(\tau)=E\left[X(t)X(t+\tau)\right]$. The second moment may therefore be estimated by averaging over time in the same way:

[TABLE]

giving the familiar autocorrelation. Likewise, for higher moments, time-invariant statistics are those obtained by considering the lags between the indexed times, collapsing over the absolute time.

2.2 Spectra and time-invariant moments

A central result in signal-processing theory concerns the relationship between time-invariant statistics and the spectral representation of the signal, given by its Fourier transform: $\tilde{X}(\omega)=\mathscr{F}\left\{X(t)\right\}$. From the inverse Fourier transform:

[TABLE]

one can readily appreciate that any shift of time, $\Delta t$, separates out into an exponential term in the spectrum of $X$:

[TABLE]

Consider the effect of a random shift, $\Delta t$, with uniform distribution (and otherwise independent of the unshifted signal) on the second moment:

[TABLE]

In other words, the time-invariant second moment, $\mu_{2}(\tau)$, is the time-domain representation of the power spectrum, a result known as the Wiener-Khintchine theorem (technically speaking, what is given here is a special case valid for integrable functions).

More generally, paralleling the situation in the time domain, moments are represented in the spectral domain as the expectation of products of the signal spectrum across different frequencies. One can conclude from a glance at Eq. (4) that any time-invariant products must arise from those combinations of frequencies for which the exponential terms cancel. For the first moment, the exponential term vanishes only at $\omega=0$, hence time invariance applies only at the origin, $\tilde{X}(0)$, related to the fixed constant mean. For the second moment, the exponent vanishes when $\omega_{1}=-\omega_{2}$, as then $\tilde{\mu}_{2}=\tilde{X}(\omega)e^{i\omega\Delta t}\tilde{X}(-\omega)e^{-i\omega\Delta t}=\left|\tilde{X}(\omega)\right|^{2}$. The idea extends naturally to higher moments as well; in general, given the product

[TABLE]

time invariant statistics relate to those terms for which

[TABLE]

Such products go by the names higher order spectra or polyspectra and it is not very difficult to show that, as with the second-order case, they are related to the time-domain moment of the corresponding order by a multi-dimensional Fourier transform. The 3rd order polyspectrum is known as the bispectrum while the 4th order is the trispectrum.

2.3 Time-Varying Spectra and the Wigner-Ville Distribution

In the previous section we equated time-invariant spectra to the spectral representation of a time-varying moment integrated over time, which we may likewise interpret as an estimate of the time-invariant spectrum obtained by averaging over time. But there are plenty of applications in which it must be assumed thatspectra vary over time. If the time-invariant estimate averages over time, shouldn’t it be possible to observe the time-varying spectrum simply by not averaging? An important class of time-frequency representations does just that, retaining the Fourier transform but skipping the integration over time; under a suitable change of variables, the outcome represents the original time-varying moment as a time-frequency distribution. The best-known example of such an object is the Wigner-Ville distribution (WVD):

[TABLE]

The WVD applies a Fourier transform along the lag dimension of the time-varying second moment, combined with a trivial change of variables that centers the lag, $\tau$, on absolute time, $t$. Though sometimes used as a time-frequency representation in its own right, it suffers from the drawback that “spurious” energy may arise as a result of interference between spectrally isolated components in regions where neither component separately contains energy.

The WVD is, arguably, more useful as a conceptual starting point in developing a general theory of power-spectral estimation: any of the common stationary and time-varying Fourier-derived power spectra can be gotten by smoothing or integrating over the distribution, which has the effect of smoothing away energy related to interference [27, 14]. Distinctive properties of any estimator are a consequence of the particular window it applies to the distribution. This idea generalizes naturally to higher moments, whose estimators may be likewise equated to smoothings of higher-order Wigner-Ville distributions, with distinguishing properties set by the choice of smoothing window.

2.4 Polyspectra and Gaussian Processes

One of the important qualitatively new properties of third-order and higher polyspectra arises from the fact that Gaussian signals are fully characterized by second-order statistics, meaning that any linear time-invariant (LTI) Gaussian process is completely determined by its mean and power spectrum. All forms of cross-frequency dependence relate to polyspectra of order 3 or higher (those with two or more frequency dimensions), therefore LTI Gaussian processes exhibit no dependence across frequencies. One of the main selling-points of higher-order spectra is that they can reveal the presence of non-Gaussian behavior as might result, for example, from non-linear signal dynamics. At the conclusion of this work we will see that it is possible to isolate such non-Gaussian signals from a background of Gaussian noise on the basis of higher-order spectra. This is an especially compelling application of higher-order spectra because the most interesting signal components are also often the least Gaussian.

2.4.1 Moments and Cumulants

Cumulants are closely related to moments, but for Gaussian processes they are identically zero at all orders above two. They may therefore be constructed by subtracting the moment of a Gaussian process with a matching spectrum from the corresponding signal moment. For higher-order spectra, this generally involves subtracting a term within those sub-domains wherein the polyspectrum reduces to a product of power spectra and/or any non-zero mean. For example, the cumulant trispectrum subtracts $\mathrm{E}\left[\left|\tilde{X}(\omega_{1})\right|^{2}\right]\mathrm{E}\left[\left|\tilde{X}(\omega_{2})\right|^{2}\right]$ from the moment trispectrum within the subspace where $\omega_{3}=-\omega_{1}$, while outside this region the moment and cumulant trispectra are equal. The cumulant bispectrum and moment bispectrum diverge only along $\omega_{1}=0$, $\omega_{2}=0$ and $\omega_{1}=-\omega_{2}$ for a non-zero mean process but are identical everywhere for a zero-mean process, a property that extends to all odd-ordered cumulants.

2.5 Cross Polyspectra

In the polyspectral moment products described by Eq (6), each term is drawn from the same signal, but there is nothing to stop us from mixing terms across multiple signals. Such cross polyspectra, which are the higher-order generalizations of the cross spectrum, describe various forms of dependence between multiple signals and are cross terms in the moments of the corresponding multivariate signal distributions. The practical application of these measures to EEG has gained some recent attention [52, 13].

Third-order cross polyspectra also arise implicitly when considering the dependence between second-order statistics, and independent measures of interest. For example, it is often of interest to relate ordinary power spectra or cross spectra to experimental variables that change over time. Standard regression models which treat power as the dependent measure implicitly model a subset of cross-bispectral interactions with the independent variable. The possible value of a more explicit treatment of third-order interactions in such analyses is a question worth pondering.

2.6 A Philosophical Aside on Time Invariance

Much of signal-processing theory is preoccupied with time-invariant measures, but the matter is usually presented as an assumption about the signal, and a rather restrictive, too-often unrealistic one, at that. It can be motivated in a different way: in many settings the origin of the time scale has no intrinsic bearing on the signal and might be chosen arbitrarily. If the observer lacks any other information by which to meaningfully situate the signal in time, such as the time of some relevant external event, then from the observer’s perspective, time is subject to some large random shift, reflecting the arbitrariness of the origin. The question becomes, what properties of the signal can such an observer measure? As illustrated in Eq. (5), a lack of timing information makes it possible to observe only time-invariant moments. Even if essential properties of the signal are not strictly constant over time, in the perspective of such an observer, the moment expectations become those of a stationary signal, whose bounded statistical properties do remain constant. Time-shift invariant properties are particularly relevant because uncertainty of timing is, more often than not, inherent in the problem of signal identification.

One may arrive at a similar idea by noticing that the time-dependent moment of order $k$ considered at the outset in Eq. (1) can be trivially equated to a time-invariant moment of order $k+1$ which includes a cross-term dependence on the impulse function, $\delta(t)$; this is so because the spectrum of the latter is constant and unit valued, $\tilde{\delta}(\omega)=1$, allowing it to partake implicitly in any product. “Situating a signal in time” therefore amounts to observing the dependence between the signal and some time-anchoring event, encoded as an impulse. But nothing gives such a function special importance next to any other function of doubtful relevance to the signal. One may conclude that time-invariant moments give a completely general description of the signal-generating process when the full structure of dependence is taken into account, whether or not the full structure includes signals encoding such “time-anchoring” event(s).

In this view, one sense of what is conventionally meant by “stationary” becomes a statement about dependence between the signal and particular time-anchoring events. There is another commonly used sense, which refers to whether low-order statistics, especially 2nd, vary over time. Such fluctuations may however be described with higher-order time-invariant statistics and so do not technically require signal statistics to change over time. For example, the tendency for signal power to fluctuate over time at a particular modulation frequency may be captured by 4th-order statistics, which encompass the power spectrum of bands within the WVD. What about the case when a glance at the spectrogram makes it plainly obvious that the signal is not stationary? In such cases, one’s eye has situated the signal in time according to some clearly identifiable series of events within the signal itself. But this is also, in a sense, what time-shift invariant statistics do: they take the signal itself as its own “time-anchoring” input. For example, in the case of PAC, phase at one frequency defines the time window over which power is observed to vary at another. In this way, time-shift invariant 3rd-order statistics may account for a particular form of 2nd-order non-stationarity.

One important qualification should be mentioned. It has been noted more than once now that moment estimation assumes the distribution to be bounded at the corresponding order, but it is possible to construct distributions whose moments are unbounded at or above some order (a property of Student’s t-distributions, for example); thus “stationary” might be applied more generally to signals arising from distributions with bounded moments, while signals with unbounded moments are properly non-stationary at the given order.

Perhaps it is conceptually helpful that in the spectral domain, time invariance relates to a particular subset of moment products for which time shifts cancel, making it more natural to think of the measurement of time-invariant spectral moment products as something which does not prejudge the nature of the signal. That is, one does not need to make assumptions about time-dependent moments to measure time-invariant moments, any more than measuring a mean requires an assumption about variance. This view of the matter may serve to unclutter the question of stationarity, making it about the boundedness of the distribution on the one hand, and dependence on particular events on the other.

3 The Bispectrum

The previous section laid out some general context and motivation for appealing to higher-order spectra in signal processing. The remainder of the work is concerned with addressing the two stumbling blocks identified earlier in relation to third-order statistics: computation and interpretation. The appendices contain a more technically oriented overview of the bispectrum in which the question of how to design estimators is given particular attention. To address the problem of interpretation, we broaden the range of estimators beyond what is normally considered. This extension will allow us to link one class of non-standard bispectral estimator to common measures of phase-amplitude coupling. The following section summarizes the main results presented in A and B.

3.1 Definition

Let $X$ be a random realization of a harmonizable time-series process, obeying usual assumptions of boundedness and integrability, with moments defined up to order 3. In most cases, it will be assumed that $X$ is real-valued, unless otherwise specified. Frequency-domain representations will be indicated with a tilde; for example, $\tilde{h}(\omega)$ is the Fourier transform of $h(t)$. Scaling and normalizing constants, such as $1/2\pi$ in the inverse Fourier transform, will also be suppressed for notational economy where doing is not expected to create confusion.

The bispectrum of the process generating $X$ is given by the following expectation in the frequency domain:

[TABLE]

Some insight into the meaning of this quantity comes by introducing the inverse Fourier expansion of $\tilde{X}$ into Eq. (8):

[TABLE]

which integrates a two-dimensional Fourier transform over time. Directly paralleling the ordinary power spectrum [14], the integrand here may be understood as the 3rd-order Wigner-Ville distribution, up to a time-centering change of variables [22, 25, 20]:

[TABLE]

For a third-order stationary process, the expectation of the third-order moment depends only on the relative lags of the times:

[TABLE]

in which case, substituting (11) into Eq. (9)

[TABLE]

In other words, for a third-order stationary process, $B$ is the two-dimensional Fourier transform of the third moment, $B=\tilde{\mu}_{3}(\omega_{1},\omega_{2})$. This relationship parallels the equivalence between the power spectrum and the second moment (i.e. auto-correlation), and extends likewise to higher polyspectra and moments [43, 25]. For a more review of the main properties of the bispectrum, the reader is invited to spend time in A. Matters related to the estimation of the bispectrum and, in particular, the design of bispectral estimators are tackled in B

3.2 The Quick Summary

Of the properties reviewed in 3.1 and A, a particularly useful set for our purposes are those related to convolution. Just as the convolution of two signals in time becomes the multiplication of their spectra in the frequency domain, the bispectrum of the convolution of two signals is the product of the separate bispectra—provided the signals are statistically independent or deterministic. Because this property allows the result to be analyzed in terms of the separate bispectra of the convolved signals, it will be useful for developing a model of recurring transient features whose timing is governed by a point process (described in 5.1 and D). The bispectra of these processes are separately amenable to more detailed analysis, and properties of the bispectrum of a transient feature, in particular, will be used to establish the link with PAC.

B takes a closer look at the problem of estimating bispectra. Just as techniques for estimating power spectra which use windowed Fourier transforms can be understood as smoothing or integrating over the ordinary Wigner-Ville distribution (see Eq. 7), bispectral estimators (as well as higher-order spectral estimators generally) likewise involve integrating and smoothing over higher-order Wigner distributions [22, 20, 25]. Differences between alternative estimators relate to properties of the windows each applies in the smoothing or integration, but all estimators reflect the same underlying quantities related to signal moments of a given order.

Estimators that rely on windowed Fourier transforms can be formally understood as employing a bank of single-sideband filters [9, 3, 36]. Taking this filter-bank perspective opens a wider range of possible estimators to consideration because filter properties might be permitted to vary across bands, whereas they remain fixed in the more conventional windowed-Fourier view. Three different classes of estimators are considered in section B.2, which differ according to the relative bandwidths of the three bands that participate in the bispectral product. The results laid out in that section are used in C to show the close relationship between PAC estimators and one such class of bispectral estimator.

4 Phase-Amplitude Coupling

Phase-amplitude coupling refers to dependence between the analytic envelope of an oscillatory signal component within one band and phase within another. For the envelope of the first component to fluctuate at the scale of the second, the bandwidth of the first component must be at least as great as the center frequency of the second and its center frequency correspondingly higher, for which reason the first component is a “fast oscillation” (FO), while the second is the “slow oscillation” (SO).

Measures of PAC relate the analytic signal in a band encompassing the SO to the analytic envelope from the band of the FO, treating the two as separate signals, typically using a second-order statistic such as coherence, phase-locking or weighted mean vector strength. While it is also most common to use signal amplitude in quantifying PAC, one might use signal power (squared amplitude) in place of amplitude in the same way, giving phase-power coherence (PhPC). C demonstrates the direct relationship between PhPC and the bispectrum. Because conventional measures that use amplitude are related to those that use squared amplitude by a simple scalar transformation of the filtered input signals, they reflect the signal bispectrum to a first approximation and might otherwise be interpreted as bispectral estimates computed on suitably transformed input signals.

4.1 Multi-modal dependence

The second-order measures described above are most suitable when there is a single mode in the distribution of the phase difference between the SO and FO envelope. A number of extensions develop tests that improve sensitivity to multimodal distributions, using alternative distributional measures over phase, such as Kullback-Leibler divergence or other entropy-related measures [56]. Such alternative methods of quantifying dependence can be translated to the bispectrum with relatively little effort; the extension involves a similar treatment of the distribution of phase within the unaveraged terms in the bispectral estimator (Eq. 28). The question is briefly considered again in 6.1. The important take-away message for the present purpose is that the essential argument developed in the following sections remains generally valid for measures of PAC derived from the joint distribution of analytic amplitude and phase in different frequency bands.

5 Bispectral Signatures of PAC

C contains a proof that phase-power coherence is fundamentally a bispectral estimator, which differs from conventional bispectral estimators only in the shape of the associated smoothing kernel. While past authors have noted some similarities between bispectral and PAC measures [37, 28], this formal relationship appears not to have been previously described. The following sections extend this observation with a signal model that will offer some more practical insight into how the bispectrum relates to PAC.

5.1 Transient and Oscillatory Models

A case traditionally used to illustrate PAC is that of a sinusoidally modulated FO added to a sinusoidal SO of the same frequency as the FO modulating window (see Figure 1). Here, the SO is a pure sinusoid, $\cos(\theta t)$, and the FO is the sinusoidally modulated sinusoid, $(1-\cos(\theta t))\cos(\gamma t)$ . It is easy to get a handle on the bispectrum of this signal because we need only worry about 4 quantities (ignoring the negative half of the spectrum, which is symmetric)): peaks at $\omega=\theta$, $\omega=\gamma-\theta$, $\omega=\gamma$ and $\omega=\gamma+\theta$:

[TABLE]

The bispectrum contains two peaks in the principal domain where the product $\tilde{X}(\omega_{1})\tilde{X}(\omega_{2})\tilde{X}^{*}(\omega_{1}+\omega_{2})$ does not vanish, at $(\omega_{1}=\theta,\omega_{2}=\gamma-\theta)$ and $(\omega_{1}=\theta,\omega_{2}=\gamma)$. From this, it is clear that it conveys something about phase-amplitude coupling.

It is natural to try to understand spectral measures in terms of the sinusoidal functions that are elemental to spectral decompositions; but this example offers only modest insight into what the bispectrum means in more realistic settings with signals that may be spectrally broad. Part of the problem is that the bispectrum, by definition, relates to the interaction between the separate components in this example. But this means that one can’t rely on any conceptually helpful decomposition according to the FO and SO considered separately to build up insight into more complicated cases.

A different tack might bring us closer to such a decomposition. We have already noted that the spectrum of the convolution of two signals is the product of their separate spectra, a fact that generalizes to higher order spectra, provided the signals are deterministic or statistically independent at the corresponding order. We might approach the problem by considering the signal as a series of transient features, $f_{i}(t-\tau_{i})$, whose timing is governed by a point process, $N$. This view of the problem decomposes the signal into a series of local features convolved with impulses generated by $N$, which encode global signal properties. Details of this model are worked out in D. For our purpose, we will want to consider transient functions that can be decomposed into two parts, a slowly varying low-frequency part, $f_{\mathrm{SO}}(t)$ and a rapidly modulated high-frequency part $f_{\mathrm{FO}}(t)$, with $f(t)=f_{\mathrm{FO}}+f_{\mathrm{SO}}$ where both SO and FO are transient. 1 As illustrated in Figure 1, the transient model explains the two peaks in the bispectrum of the SO+FO signal as the restriction of the bispectrum of $f$ to peaks in the bispectrum of $N$ through multiplication.

5.1.1 A Caveat

The attentive reader may have noticed a problem with Figure 1, which gives the occasion to point out some limitations of the transient model: it illustrates the case when $f$ is fixed, meaning that the phase of the FO and SO are always in the same relation and by extension that the frequency of the FO must be an integer multiple of the SO, a rather severe restriction which does not arise in the oscillatory view of the same problem. Fortunately, the transient model is not bound to this assumption, and a key result described in section 5.3.1 is that the bispectrum retains energy in the SO $\times$ FO range of frequencies even if the phase of $f_{\mathrm{FO}}$ is random. The limitation of the analysis presented there relates instead to the fact that it takes as a simplifying assumption that N and emitted transients, the $f_{i}$’s, are mutually independent.

Among the results described in 5.3.1, when the phase of the $f_{\mathrm{FO}}$’s are random or vary uniformly, the contribution of the point process to the bispectrum drops out along a given dimension in the FO frequency range, leaving only the contribution of the transient. In such cases, the analysis concludes that the bispectrum will still be confined to peaks in the spectrum of $N$ in the range of the SO but will not be likewise restricted in the range of the FO. To arrive at the fact that there are only two peaks in the bispectrum when $\gamma$ is anything but an integer multiple of $\theta$ (causing the relative phase of the $f_{\mathrm{FO}}$’s to vary), it is necessary to account for the dependence between adjacent $f_{\mathrm{FO}}$’s, which is something the analysis neglects. Nevertheless, as a conceptual tool it gets the general picture right, and as an analytic tool, the assumption of sequential independence is, if not always benign, at least one of the more commonly accepted transgressions in such settings.

5.2 The Transient Bispectrum

The main advantage of the transient model and its relevance for understanding PAC is explained next. If the mean of $f$, $\tilde{f}(0)$, is nonzero, it is easily seen that the bispectrum of $f$ contains its own power spectrum, $\tilde{f}(0)\left|\tilde{f}(\omega)\right|^{2}$, along the axes $\omega_{1}=0$ and $\omega_{2}=0$. This observation can be extended more generally to slow oscillations beyond the DC axis: suppose $f$ is the sum of two components, $f=f_{\mathrm{SO}}+f_{\mathrm{FO}}$, the first of which, $f_{\mathrm{SO}}$, occupies a narrower bandwidth than the second at center frequency, $\xi$, while the amplitude of the second rises and falls transiently at a time scale shorter than the period of the first. This modulation implies that $f_{\mathrm{FO}}$ is effectively windowed at a time scale less than $2\pi\xi^{-1}$. The consequence of such time-domain windowing on the frequency domain mirrors exactly the effect of frequency-domain windowing on the time domain (i.e. filtering), which is to say the result is a smoothed version of the unwindowed spectrum. In particular, windowing at the scale $2\pi\xi^{-1}$ implies that the spectrum of $f_{\mathrm{FO}}$ is smooth at the scale of $\xi$, justifying the following approximation

[TABLE]

where $\tau$ accounts for some arbitrary time delay. These points are illustrated in Figure 2. The term $e^{-i\xi\tau}$ describes the phase of the amplitude modulation relative to the emission time. We may partition the bispectral plane according to the spectral ranges of $f_{\mathrm{SO}}$ and $f_{\mathrm{FO}}$, as shown in Fig. 3. Where energy falls within this partitioning gives the first set of criteria for distinguishing PAC.

5.3 Bispectral Definition of PAC

First, we need a more precise statement of what is meant by PAC. An often-cited example of what is not meant are cases of “spurious” PAC that arise from recurring transient signal features, such as spectrally broad “spikes” or sharp-edged waves, which may or may not be periodic (see Figure 4) [37, 41]. Such false PAC will tend to exhibit a consistent phase relationship across a range of frequencies, in particular, along the harmonics of any fundamental periodicity. “True” PAC is taken here to mean nested oscillations with the following characteristics:

For a fast oscillation (FO) to be “nested,” within a slow oscillation (SO), its amplitude must vary at a time scale around the period of the slow oscillation. 2. 2.

The FO should be concurrent with the SO. 3. 3.

The phase of the FO must either lack any consistent phase relationship with that of the SO or it must fall within a band that is isolated from the SO and any related harmonics.

The first point implies that each burst of the FO occupies a bandwidth wider than the center frequency of the slow oscillation in which it is embedded, giving a characteristically smooth spectrum at the scale of the SO. The second point excludes non-concurrent responses; for example, a fast oscillation followed by a slower one. The third point excludes spectrally broad features associated with sharp transients and ensures that any oscillation nested in the SO with a consistent phase relationship is genuinely oscillatory. These points justify guidelines for identifying PAC in the bispectrum, outlined in the following sections. More generally, they provide a way to characterize recurring signal features through the bispectrum.

5.3.1 Defining PAC: The Outside Criterion

The first criterion for PAC is that the time scale over which the amplitude of the nested fast oscillation, $f_{\mathrm{FO}}$, varies must be less than the period of the slow oscillation, $f_{\mathrm{SO}}$, in which it is embedded. When this holds, the approximation in Eq. (14) implies that the power spectrum of $f_{\mathrm{FO}}$ appears in the bispectrum parallel to the $\omega_{2}$ axis within the support of $\tilde{f}_{\mathrm{SO}}$ because within this range

[TABLE]

and likewise for the $\omega_{1}$ axis under symmetry. The time delay, $\tau$, recovers the relative lag between $f_{\mathrm{SO}}$ and the envelope of $f_{\mathrm{FO}}$. It can be seen that the power spectrum of $f_{\mathrm{FO}}$ is reproduced parallel to the $\omega_{1}$ axis at $\xi$, while in the orthogonal direction, the spectrum of $f_{\mathrm{SO}}$ is reproduced parallel to $\omega_{2}$ within the support of $\tilde{f}_{\mathrm{FO}}(\omega_{1})$. Examples of this effect can be found in Figure 4 (middle panels) for both the nested-oscillation and transient test signal.

Applying a similar analysis to the stochastic case describe in D.2, suppose that the spectrum of $f_{\mathrm{FO}}$ overlaps negligibly with $f_{\mathrm{SO}}$, so that $\tilde{f}_{\mathrm{SO}i}(\omega)\tilde{f}_{\mathrm{FO}j}(\omega)\approx 0$, then within the support of $\tilde{f}_{\mathrm{SO}}(\omega_{1})\tilde{f}_{\mathrm{FO}}(\omega_{2})$ we are left with

[TABLE]

Suppose the phase of $f_{\mathrm{FO}}$ is random such that $\left<\tilde{f}_{\mathrm{FO}}\right>=\left<\tilde{f}_{\mathrm{FO}}\tilde{f}_{\mathrm{SO}}\right>=0$, then only the first two terms of Eq. (16) remain:

[TABLE]

The first term does not vanish when there is a consistent relationship of phase between $f_{\mathrm{SO}}$ and the amplitude modulation of $f_{\mathrm{FO}}$. The second term remains when there is also a consistent lag between both terms and the emission of the point process, which is reflected accordingly in the weighting by the spectrum of the point process, $\tilde{\mu}_{2}(\omega_{1})$. But because we are free to define the emission times of the point process such that $\tau_{\mathrm{SO}}=0$, the second term merely reflects the power spectral contribution of the driving point process. The bispectral estimate within the support of $f_{\mathrm{SO}}(\omega_{1})f_{\mathrm{FO}}(\omega_{2})$ therefore depends on the timing of the amplitude of $f_{\mathrm{FO}}$ relative to the phase of $f_{\mathrm{SO}}$, but not the phase of $f_{\mathrm{FO}}$.

This form of cross-frequency dependence fulfills sensible criteria for phase-amplitude coupling: it implies a correlation between the amplitude of the amplitude-providing signal component, $f_{\mathrm{FO}}$, and an underlying oscillation with a longer period than the time scale of the amplitude modulation of $f_{\mathrm{FO}}$, which is what is meant conventionally by phase-amplitude coupling. Reflecting the scale of its modulation, the spectrum of $f_{\mathrm{FO}}$ must be both relatively broad and smooth, but because the power spectrum discards phase information in $f_{\mathrm{FO}}$, the phase of $f_{\mathrm{FO}}$ does not require any consistent relationship with $f_{\mathrm{SO}}$, which is also a frequent characteristic of nested oscillations in phase-amplitude coupling. This condition does not, however, rule out a contribution from spectrally broad transients (see example in Fig. 4, lower middle panel), a possibility that motivates the second and third criteria, next.

5.3.2 Defining PAC: The Inside Criterion

The central part of the bispectrum covers the regions of support for $\tilde{f}_{\mathrm{SO}}(\omega_{1})\tilde{f}_{\mathrm{SO}}(\omega_{2})$ and $\tilde{f}_{\mathrm{FO}}(\omega_{1})\tilde{f}_{\mathrm{FO}}(\omega_{2})$ (white boxes in Fig. 3). Non-vanishing terms in the bispectrum within these regions relate to the separate bispectra of the components. Energy in the $\mathrm{FO}\times\mathrm{FO}$ “inside” box of Fig. 3 reflects the consistency between the phase of the FO and its own modulating envelope. For example, if the FO is obtained by shifting $\cos(\beta(t+\tau))$ by some random delay, $\tau$ and windowing with an unshifted envelope, $h$:111It is assumed here for the sake of argument that $h$ varies slowly relative to $2\pi\beta^{-1}$ to avoid any complications related to overlapping spread from negative and positive frequencies (i.e. Bedrosian’s theorem applies) but quickly enough to allow for a non-vanishing bispectral product.

[TABLE]

then its bispectrum is given by

[TABLE]

which vanishes if the variability in $\tau$ suffices to make $e^{i\beta\tau}$ uniform in the spectral range of the FO. Energy near the center of the inside box therefore implies a consistent relationship between the modulating envelope of the FO and its own phase.

We observed similarly in 5.3.1 that energy in the outside box has to do with consistency between the envelope of the FO and the phase of the SO. Such consistency may of course, in both cases, be the result of some spectrally broad feature that we do not wish to call a nested oscillation, which will tend to manifest as consistent phase across a broad range of frequencies spanning the transition from inside to outside boxes. One might reason that energy in both boxes implies a consistent phase relationship between the SO and the FO, but this need not always be the case. Variability in the timing of the FO modulating window may happen on the scale of the SO without abolishing energy in the outside box, and this scale might be large relative to the frequency range of the FO. Variability of the FO delay will, however, attenuate the bispectrum around the transition from outside to inside boxes, because any energy in the SO at frequencies above the scale of the variability will be washed out. This will have the effect of accentuating the separation between the inside and outside boxes as the contribution of higher-frequency components of the SO is suppressed by variability in the timing of the FO window.

One may conclude from this that FO’s which are either spectrally isolated, and hence oscillatory, or those that occur with some variability of delay may generate modes in the inside region that clearly stand apart from energy in the outside region.222It can also be pointed out that energy in any of the regions might be the result of entirely unrelated signal components; it is assumed for the present that the signal bispectrum is dominated by one component. The problem of decomposing multi-components signals into their independent parts is an important open question briefly considered in section 8.2. In contrast, energy generated by a sharp-edged transient will tend to exhibit a tight phase relationship among all components, leading to a continuous smearing of energy across both regions. In the presence of periodicity, such energy will cluster at a lattice of harmonically-related peaks along axes given by $m\omega_{1}=n\omega_{2}$ where $\tilde{f}(n\omega)\tilde{f}(m\omega)\tilde{f}^{*}\left((n+m)\omega\right)$ is non-vanishing [51]. This effect is a well-known source of spurious phase-amplitude coupling [41], as well as spurious $n:m$ phase locking [50].

Figure 4 shows a comparison between such a periodic transient and nested oscillations involving a purely random FO. Figure 5 shows a very similar comparison for real data recorded from two neighboring regions of lateral occipital cortex in an epilepsy patient. Frequent inter-ictal spiking in one channel (Fig. 5, channel B, lower panels) results in what looks superficially like PAC between phase around 10 Hz and power in the 50-150 Hz range. In the bispectral measure, this apparent phase-amplitude coupling is accompanied by clearly visible harmonic banding and a wide monotonically diminishing spread of energy, which is much less easy to see in the PAC measure, due to the smoothing and limited field of view of the latter.

In a different nearby channel, both measures reveal an association between phase in the 1-4 Hz range and power in the 80-200 Hz range, in the absence of a similar broad smear of energy—this looks more like the picture expected of phase-amplitude coupling related to nested oscillations. In the bispectral measure, some energy also appears near the center of the inside region (Fig. 5, inset of the upper right panel), implying that the FO exhibits some degree of consistency in its form (that is, phase is not completely random with respect to the modulating window), but the mode clearly separates from the outside region, making the attribution of PAC reasonable. Because PAC measures do not include the inside region within their field of view, this detail is unavailable from them. As this example illustrates, both kinds of phenomena, PAC and sharp waves, occur as part of ongoing activity [58]; discriminating between them is crucial to the interpretation of PAC.

5.3.3 Defining PAC: The Time-Delay Criterion

Eq. (15) implies that the bispectrum contains potentially useful information about the relative timing of SO phase and FO amplitude modulation. Information about timing can be recovered with an inverse Fourier transform along one of the frequency dimensions, as shown in Fig. 6. For the homogenous point process with constant $\lambda$, in the presence of pure PAC:

[TABLE]

where the integral covers the support of $f_{SO}$, excluding the diagonal symmetry region in quadrants II and IV. Because $f_{\mathrm{SO}}$ and $f_{\mathrm{FO}}$ occupy non-overlapping spectral ranges, the result is a two-dimensional function with two bands, the first around the center frequency of $f_{\mathrm{SO}}$, containing the autocorrelation of $f_{\mathrm{FO}}$ scaled by $\tilde{f}_{\mathrm{SO}}(\omega_{2})$ and the second within the support of $\tilde{f}_{\mathrm{FO}}$, containing $f_{\mathrm{SO}}$ scaled by $\left|\tilde{f}_{\mathrm{FO}}(\omega_{2})\right|^{2}$, both shifted by the time delay between them. For a periodic $N$, the terms are reduplicated at the corresponding periodicity.

This time-domain representation is useful for verifying that the FO and SO overlap in time, as should be the case for nested oscillations. More generally it provides a way to relate PAC to specific features of SO. It also illustrates the broader point that because the bispectrum preserves phase information, it encodes details of waveform shape that are lost in standard spectral analyses [6].

6 Extensions

6.1 Multi-modal dependence

As mentioned in section 4.1, a variety of extensions to basic measures of PAC address more complicated forms of dependence between phase and amplitude, as when peaks of amplitude occur at two or more separate phases. It was also noted that similar extensions might be developed for the distribution of phase of within bispectral estimators. A major advantage of standard bispectral estimators is that their capacity to reveal phase-amplitude dependence is not inherently limited by the analysis bandwidths, but the need to construct a distribution over phase in such extensions is bound to increase sensitivity to the parameters of the analysis filter.

An alternative strategy for this problem might appeal to subregions within 4th and higher order spectra. For example, if the modulation of the FO is doubled in frequency compared to the SO, meaning there are two peaks of amplitude at opposite phase over every $2\pi$ cycle of phase in the SO, this might be revealed within the $\omega_{1}=\omega_{3}$ plane of the trispectrum:

[TABLE]

Likewise, a distribution with $k$ peaks in the amplitude of the FO should be revealed within the $\omega_{1}=\omega_{3}=\omega_{4}=\dots=\omega_{k}$ plane of the $(k+1)^{\mathrm{th}}$ order spectrum:

[TABLE]

Because these pertain to a two-dimensional subspace in each case, the full polyspectrum does not need to be estimated to obtain a test of $k$-mode dependence, and estimation should therefore not greatly increase the computational burden over that of the bispectrum.

6.2 The Bispectrum and Phase-Frequency Coupling

So far, our interest in phase-amplitude coupling has led us to focus on smooth bispectra, which gain their smoothness from the windowing of the FO. As shown in C, conventional PAC measures aggressively smooth the bispectrum over the FO frequency range, giving them some specificity for bispectral features expected of PAC. One might view this as an advantage or disadvantage of PAC measures, depending on context, but it raises the question, what should one make of those sharply-varying features they suppress?

To gain some insight on the question, we may turn back to the oscillatory model of section 5.1; in the case of an FO composed of a sinusoidally modulated sinusoid, the bispectrum contained two peaks of the same sign arising from the product of each of the two sidebands in the spectrum of the FO with the central “baseband” peak (See Fig. 1, panel D). It will be useful to think of the FO in this case as the sum of two sinusoids offset in frequency from each other by $\theta$ and both modulated by $\cos(\theta/2\>t)$:

[TABLE]

[TABLE]

The resulting FO is illustrated on the left side of Fig. 7.

We may model the opposite of a smooth bispectrum by inverting one of the two side peaks in the spectrum of the FO within Eq. (13), which has the effect of inverting one of the two peaks in the bispectrum of the FO + SO composite signal:

[TABLE]

Because this change causes the two peaks in the bispectrum to have opposite sign, they will tend to be suppressed by estimators that blur them together. The FO in this case is illustrated on the right side of Fig. 7. As in the former case, this signal may also be written out as the sum of two modulated sinusoids:

[TABLE]

The crucial difference lies in the fact that the modulating envelopes of the two components are out of phase with each other, which has two important consequences; first, the envelope of the composite signal fluctuates over a smaller range than in the previous case:

[TABLE]

resulting in amplitude modulation which varies, in this case, over $\left[1,\sqrt{2}\right]$ rather than $[0,2]$. Second, the instantaneous frequency of the signal gravitates back and forth between between $\gamma+\theta/2$ and $\gamma-\theta/2$ according to which of the two components has the greater amplitude at any moment. One can therefore attach two interpretations to this case; the first is to regard the components as separate AM signals whose modulating windows happen to be out of phase with each other. The second interpretation is that the FO is composed of a single frequency modulated (FM) signal and the bispectrum reflects phase-frequency coupling (PhFC) between the frequency of the FO and the phase of the SO. Estimators tailored to PhFC might replace the smoothing kernel of PAC and conventional bispectral measures with a high-pass “sharpening” kernel. As with PAC, whether the interpretation of PhFC is appropriate in any given case depends on other details.

7 Methods

7.1 Algorithm

General principles of estimation are reviewed in B. The following section briefly describes the implementation with discretely sampled data used in the present examples. The “direct technique” of bispectral estimation starts with a windowed Fourier time-frequency decomposition, summing bispectral products over time. Although the technique is typically described in its application to the short-time Fourier transform, any time-frequency decomposition, $S[k,n]$, that can be expressed as a bank of single sideband filters might be used.333Short-time Fourier transform, complex demodulation and single sideband filtering all belong to a family of techniques that are formally equivalent, up to a linear phase term which is negligible for the present purpose because it cancels from the time-invariant bispectral products [36]. Here $x$ is a discretely sampled signal and $h$ is the time window used by the estimator:

[TABLE]

[TABLE]

where $j,k$ and $l$ index frequency bands $\omega_{1}[j],\omega_{2}[k]$ and $\omega_{3}[l]$ such that

[TABLE]

or more generally

[TABLE]

with $\Delta W_{3}[l]$ denoting the bandwidth of the $S_{3}[l,.]$ band. The use of subindices on $S$ and $h$ here is meant to show that the analysis windows used by the estimator are permitted to vary according to frequency and place in the bispectral product, encompassing the non-standard estimators described in B, as well as those with bandwidths that vary over frequency, such as continuous wavelet decompositions [30].

7.1.1 Normalization

Bispectral statistics are shown in normalized form

[TABLE]

using the normalization described in B.1.1, which takes the sum of bispectral product magnitudes over time as the normalizing term [24, 35], given by

[TABLE]

$A[j,k]$ amounts to the maximum possible value of $|B[j,k]|$, leaving the magnitudes of $S$ unchanged, which is attained when phase remains constant across all terms Eq. (28). With this normalization, $\beta$ may be interpreted as a weighted-average phase locking value, which avoids some of the ambiguity related to the contribution of amplitudes in more standard forms of bicoherence [54, 40]. For a review of a closely-related problem in the context of ordinary coherence, which is directly applicable in this case as well, see Kovach [35].

7.1.2 Bias Correction

This weighted-average interpretation also allows for a straightforward bias correction [35], which compensates for amplitude-dependent biases that may vary by frequency [54]. The biascorrection subtracts an expected bias from the magnitude bispectrum and renormalizes

[TABLE]

where expected bias is given by

[TABLE]

Note that the use of the magnitude symbol on the left side of Eq. (31) is an abuse of notation as random error might cause $\left|\beta_{\mathrm{BC}}\right|$ to become negative. The correction can also be applied in a way that preserves phase information:

[TABLE]

which has the effect of shrinking the positive bias of $\beta$. For further explanation, see Kovach [35].

7.1.3 Software

In the present examples and simulations, frequency-domain complex demodulation [9] (DBT) was employed as described by Kovach and Gander [36], using a software implementation for Matlab (Mathworks, Nattick, MA) developed by the first author and available at https://github.com/ckovach/DBT. Polyspectral estimation with fixed bandwidths across all terms for a single signal is implemented in the function PSPECT, while estimation that accommodates different bandwidths, as well as cross polyspectra, is implemented in PSPECT2. A DBT implementation of phase-power coupling is given in DBTPAC.

For the demonstrations in figures 5 and 4, phase-power coupling plots were generated with a bandwidth of 1 Hz for the SO and 40 Hz for the FO. Bispectra were computed with a bandwidth of 2 Hz.

7.2 Procedure

Example data shown in Fig. 5 were acquired from an epilepsy patient during a period of clinical monitoring with invasive intracranial recordings. Data were obtained over 20 minutes during which the patient passively watched a sitcom episode. All procedures were approved by the Internal Review Board of the University of Iowa and conducted with voluntary informed consent of the patient. For a detailed review of techniques, see [45].

8 Discussion

With respect to the relationship between PAC and the bispectrum, we have shown:

Common measures of PAC based on analytic phase and amplitude are in fact bispectral estimators, and as such provide no unique information beyond what is recovered by standard bispectral estimators. 2. 2.

PAC measures are severely biased with respect to the symmetry of the bispectrum and introduce artificial constraints on the range and resolution of the estimator.

These limitations provide a clear rationale for favoring standard bispectral estimators in the evaluation of phase-amplitude coupling. We have further given a more detailed, albeit qualitative, framework by which to evaluate the presence and nature of PAC from the bispectrum based on the distribution of energy within the bispectrum. The signal model developed in pursuit of this framework lends itself to a particularly straightforward interpretation: it describes some recurring feature embedded at times determined by the driving point process.

8.1 Tuning In to the Bispectrum

Beyond methodological questions, it is worth considering the relevance of this work to theoretical models of the gating and transmission of information in the brain and elsewhere. One of the epoch-defining technologies of the last century, radio transmission, developed from the discovery that information can be efficiently conveyed over narrow slices of the radio spectrum. In essence, a sender and a receiver agree in advance on the power spectrum of the transmitting signal and tune resonators within their respective devices accordingly. Information is therefore transmitted within the space of signals of a given power spectrum, and the degrees of freedom within that space, which encode any to-be-transmitted information, correspond to the Fourier-domain phase of the signal. The crucial problem of overcoming noise is addressed with a “brute force” strategy: broadband noise is overwhelmed within a small frequency interval, along with any competing signal, by concentrating the transmission energy in the narrowest practical region of the spectrum, resulting in a tradeoff between signal-to-noise ratio and maximum transmission rate.

Starting from the model described in 5.1 and D, one may carry a similar transmission scheme over to the bispectrum and in return draw forth some promising new features. The signal space in this scheme is constrained to have a given agreed-upon bispectrum within the $\mathrm{SO\times FO}$ division of the bispectral plane, both with respect to magnitude and phase. Transmitted information is encoded by the phase of the $\mathrm{FO}$, which Eq. (17) allows to vary freely without affecting either the phase- or magnitude-bispectrum. In contrast, the phase of the $\mathrm{SO}$ and the relative delay of the $\mathrm{FO}$-modulating envelope determine the phase bispectrum within this region, meaning that they are not likewise free to vary.

Three noteworthy features emerge from this scheme:

The encoding space is specified within a two-dimensional spectral plane, rather than the one dimension afforded by the power spectrum, increasing the number of ways the space may be divided. Among other things, this means that separate channels may transmit information within spectrally overlapping bands without creating interference. 2. 2.

The phase bispectrum preserves information about timing, which allows for the addition of time division to the encoding scheme. Eq. (15) demonstrates how this might work by spacing the modulating window of the $\mathrm{FO}$ at different delays. Frequency division is determined by the magnitude bispectrum while time division is determined by the phase bispectrum, thus channels which overlap in both dimensions of the magnitude bispectrum may still avoid interference through a division of the phase bispectrum. 3. 3.

Strategies for overcoming Gaussian noise are not confined to brute-force commandeering of the spectrum in a given region, but may take advantage of the fact that Gaussian processes have vanishing bispectra [26]. This holds true in the sense of an expectation, meaning in practice, signal-to-noise ratio might be improved through the averaging of repeated transmissions.

An intuition-friendly summary of this scheme is as follows: it describes the encoding of information within a high-frequency waveform, the $\mathrm{FO}$, which is embedded at a predetermined time relative to a low-frequency carrier waveform the $\mathrm{SO}$; Figure 2 gives an example of how the components of such a signal might appear. The $\mathrm{SO}$ may serve as a timing signal which indicates when a new transmission has arrived. The receiver might identify the transmission by filtering for the $\mathrm{SO}$ and then isolate the $\mathrm{FO}$ with the combination of a second filter and time-window around the interval encoded by the phase bispectrum. In this example, the bispectral equivalent of a “resonator” is a multistep process involving: (1) identification of the $\mathrm{SO}$ through a suitable filter, which can be derived from the bispectrum as shown by Eq. (19). (2) Peaks in the output of the SO filter triggers the isolation of the $\mathrm{FO}$ through a second filter and time window. Finally, if the signal is corrupted by Gaussian noise affecting the output of the $\mathrm{SO}$ filter or the phase of the $\mathrm{FO}$, it may (3) be cleaned up by averaging multiple transmissions, provided the same $\mathrm{FO}$ is also repeated with the $\mathrm{SO}$.

Although this scheme raises the complexity of encoding and decoding considerably beyond the simple resonators of ordinary power-spectral transmission, it greatly increases flexibility in the choice of transmitting signal, accommodating transient as well as oscillatory, spectrally broad as well as narrow signals. This scheme also meshes with recent theoretical accounts of the possible functions of PAC. Such models emphasize the role of PAC in multiplexing and routing of information through temporal gating and frequency division, both properties that emerge from our consideration of the bispectrum [38, 31, 61, 1, 39, 2, 29]. By fixing such recent developments within a more formal spectral framework, it may be possible to better elucidate questions of efficiency and optimality as well as to generate predictions about the behavior that efficient and optimal systems should exhibit.

Common techniques of spectral analysis developed over the first half of the last century in concert with narrowband radio and telegraphic transmission, and it is no accident that they are suited in the first place to the study of spectrally narrow signals. Their availability, familiarity and versatility accounts in part for the special attention given to the neural functions of those periodic and oscillatory phenomena about which classical techniques have the most to say [59, 11], an emphasis carried over into the study of phase-amplitude coupling. Yet spectrally broad, often aperiodic phenomena, exemplified most obviously by the action potential and transient evoked responses, also figure prominently in neural signaling at many scales. The importance of what traditional spectral methods neglect from these signals is becoming increasingly apparent [15]. The tools afforded by higher-order spectra may go a long way towards illuminating the forms and functions of neural responses beyond the confines of the narrow passband.

8.2 Future Directions

We have laid out a qualitative framework for understanding and interpreting bispectral statistics. As one would hope of any useful framework, it brings several interesting questions to the fore, the answering of which must be part of a larger project. Our aim has been less to give a comprehensive set of answers than to provoke interest in how the general picture might be further developed. No attempt has been made here to define quantitative scores that summarize the qualitative features outlined above. The extension of bispectral statistics to the study of phase-frequency coupling was given a cursory and highly schematic treatment and much more might be done to flesh out this application.

One problem that is particularly ripe for further development is that of distinguishing between bispectra which result from a single recurring feature or the co-occurrence of multiple independent features, and developing a suitable decomposition in the latter case. If the underlying processes are independent, the signal bispectrum is a sum of their separate bispectra. An important question is how one might go about decomposing the bispectrum or other higher-order spectra according to distinct sets of independent features. This question is related to one addressed by various blind deconvolution algorithms, a number of which proceed by identifying filters that maximize scalar moments [60, 17] (that is, the value of the $k^{th}$ moment at 0 lag); of particular relevance in the case of the bispectrum are those that rely on skewness [47, 46]. Future work will consider this problem in greater detail.

Acknowledgements

We wish to thank Phil Gander, Alex Billig and Matthew A. Howard, III.

Appendix A Properties of the Bispectrum

Here we review a handful of properties that will be called upon in developing the transient model of section 5.1 and D. Unless otherwise specified, variables $r,s,t$ refer to time and $\tau$ to time lags.When $\delta$ is within the argument of a function it is also a lag variable, while outside of an argument it will represent the Dirac delta function. Frequency variables will be represented with $\omega,\nu,\eta,\xi$ and $\lambda$.

A.1 Symmetry

The third moment is unchanged under a permutation of the time delays, $\tau_{1}$ and $\tau_{2}$, over the three instances of $X$ in Eq. (10), creating a six-fold symmetry, illustrated in Figure 8, with the following relations:

[TABLE]

Because the spectrum of a real-valued signal is conjugate symmetric about the origin, the 6-fold symmetry of the third moment becomes a 12-fold conjugate symmetry in the bispectrum under the following relations:

[TABLE]

For a real signal, the full bispectrum may be recovered from estimates within any one of the symmetry regions. These relations apply as well to the frequency dimensions of the third-order Wigner-Ville distribution, $W_{3}$.

A.2 Convolution

The bispectrum inherits a number of basic spectral properties, among which are those related to convolution. Because the convolution of two signals in time equates to multiplication in the spectral domain, it is directly apparent that the bispectrum of a signal formed from the convolution of two other signals is likewise the product of their separate bispectra, provided the signals are third-order independent or deterministic; that is,

[TABLE]

for third-order independent $X,Y$.

A.3 Multiplication

The reverse relationship, by which the product of two independent signals in time becomes a convolution in the frequency domain holds generally for the bispectrum only if at least one of the signals is also third-order stationary; that is

[TABLE]

if $X$ and $Y$ are third-order independent and if $X$ or $Y$ is third-order stationary. Because of the stationarity requirement, this relation applies only in trivial cases (i.e. $X(t)=c$) when both signals are deterministic, such that $E\left[X(t)\right]=X(t)$.

The relationship does, however, hold generally for the bispectral Wigner-Ville distribution; that is

[TABLE]

for third-order independent $X$ and $Y$, where $\circ_{\omega}$ denotes convolution over the frequency dimension. Then

[TABLE]

For third-order stationary $X$, time drops out of the expectation as $\mathrm{E}\left[W_{3}[X]\right]=\tilde{\mu}_{3}[X]$, so that

[TABLE]

if $X$ is stationary and third-order independent of $Y$.

A.4 Alternative Forms

The lags in Eq. (10) can be defined in different ways while preserving the essential properties of $W_{3}$ [55]; it will at times be useful to apply a change of variables giving:

[TABLE]

where $\nu_{1}\equiv\omega_{1}-\omega_{2}$ and $\nu_{2}\equiv\omega_{1}+\omega_{2}$, $\delta\equiv(\tau_{1}-\tau_{2})/2$, and $\tau\equiv(\tau_{1}+\tau_{2})/2$. Where necessary, the distinction between these forms will be made through a similar abuse of notation in the arguments.

Similarly, in the spectral domain, it will later become useful to consider the change of variables $\omega^{\prime}_{1}=\omega_{1}$ and $\omega^{\prime}_{2}=\omega_{2}+\omega_{1}/2$. This change effectively modifies the definition of the bispectrum to a form that is symmetrical around $\omega_{2}$ :

[TABLE]

Under this change of variables, the bispectrum assumes a radial symmetry with each of the axes shown in Fig. 8 separated by 30 degrees. The corresponding change in the time-domain involves a modification of the lags in the 3rd moment, $\tau^{\prime}_{1}=\tau_{1}-\tau_{2}/2$ and $\tau^{\prime}_{2}=\tau_{2}$, adding $t^{\prime}=t+\tau_{2}/2$ for the sake of symmetry gives

[TABLE]

Appendix B Bispectral Estimators

A central concern here is how to obtain an estimator of $B$ for a third-order stationary signal. The extension to non-stationary estimation parallels that of ordinary spectral estimation and will not be considered separately here. The most straightforward approaches follow the example set by conventional methods for computing Fourier power spectra [44, 10]. The so-called “indirect” method first estimates the third moment of the signal in the time domain and windows it, giving $\hat{\mu}^{H}_{3}=\hat{\mu}_{3}(\tau_{1},\tau_{2})H(\tau_{1},\tau_{2})$, from which the smoothed bispectrum is obtained with a two-dimensional Fourier transform. The “direct” method first computes a series of windowed Fourier spectra, $\chi_{j}(\omega)$, and averages over the product $\chi_{j}(\omega_{1})\chi_{j}(\omega_{2})\chi_{j}^{*}(\omega_{1}+\omega_{2})$. Both approaches generate a largely equivalent family of estimators.

To minimize bias in the estimate and give an estimate with the proper symmetry, it is usually required that $H$ itself exhibit the same symmetry as the third moment, given in Eq. (34). Windows should also be real-valued in time and frequency and non-negative in the frequency domain [48]. Two-dimensional windows fulfilling these requirements may be constructed in either the time or frequency domains from a window defined in one dimension, $h(t)$:

[TABLE]

Although specific windows are shown to be optimal for minimizing bias of bispectral estimators [48], we will consider a more general case, which also includes windows for which symmetry relation may not hold. We will relate PAC measures to bispectral estimators that lack the proper symmetry constraint, which motivates dropping the constraint in the present discussion. It will still be assumed that the windows are non-negative and real-valued in the frequency domain and symmetric under a 180 degree rotation, although the complex conjugate will be indicated where appropriate for notational consistency.

The “direct” estimator most commonly begins with a series of windowed Fourier transforms, evenly spaced in time (here treated as continuous for simplicity), which may be regarded equivalently as the outcome of a bank of filters whose outputs are analytic signals shifted in frequency (i.e. complex demodulates) obtained by bandpass-filtering the original signal [23, 9, 36]:

[TABLE]

An estimate of the bispectrum is obtained from this time-frequency representation as

[TABLE]

The last integrand in (47) defines the smoothing window:

[TABLE]

Comparing this to Eq. (45) makes it clear that the resulting estimator fulfills the symmetry requirement of Eq. (34). Combining these results:

[TABLE]

which makes the smoothing property of the bispectral estimator explicit. With a time-centering change of variables, this relation also gives a time-varying estimate of the bispectrum as a smoothing of the third-order Wigner-Ville distribution, absent the integration over time [20, 25]. As in the second-order case, suppressing cross terms in the time-varying spectrum requires some additional smoothing over time.

B.1 Bicoherence

Each of the estimators considered above can be regarded as an inner product between two signals, the first being the higher analytic band, $\chi_{3}(\omega_{1}+\omega_{2},t)$ and the second obtained by multiplying the two lower bands $\chi_{1}(\omega_{1},t)\chi_{2}(\omega_{2},t)$. The degree of alignment between these two terms is often used as a normalized measure of dependence, which is obtained in the same way as for an ordinary correlation, by dividing the inner product with the magnitudes of the separate terms [34]:

[TABLE]

This is a commonly used definition of bicoherence.

B.1.1 Alternative Definition

It is worth noting that this common definition of bicoherence does not share the full set of symmetry properties with the bispectrum because the normalizing term is not likewise symmetric. Whereas $\beta(\omega_{1},\omega_{2})=\beta(\omega_{2},\omega_{1})$, it is not generally true that $\beta(\omega_{1},\omega_{2})=\beta(-\omega_{1}-\omega_{2},\omega_{2})$. Some alternative definitions of bicoherence do preserve the appropriate symmetry. The following defines an alternative index [24]:

[TABLE]

This index may be interpreted as a weighted-average vector strength over the phase component of the bispectrum, thus it gives a measure of phase locking [35].

B.2 Estimators with Different Analysis Filters

Eq. (47) constrains the analysis window for each of the three bands used in the estimate to be identical, which had the benefit of yielding an estimator with the proper symmetry. While the symmetry conditions in (34) are usually regarded as a criterion of admissibility for windows used in estimating the bispectrum, there is in fact no essential reason why one may not smooth the bispectrum with windows that do not obey these relations. Although the resulting estimator will be biased with respect to the symmetry of the bispectrum, it might conceivably have other useful properties that outweigh this consideration. For instance, one might select estimators that are better tuned to specific features of interest in the signal bispectrum, which we will find illustrated with PAC measures. We therefore consider next the consequence of estimating the bispectrum using different analysis filters for each band, which leads to the following estimation window:

[TABLE]

with separate analysis filters, $h_{1}$, $h_{2}$, and $g$.

In general, because the smoothing occurs over the two-dimensional bispectral plane, the choice of third window is over-constrained. This is easily recognized in the limit as the bandwidth of one filter diverges from that of the other two, in which case the smoothing window factors over the two dimensions. For example, as $\tilde{g}\rightarrow\delta$ in (52), one obtains:

[TABLE]

Whereas in the limit $\tilde{h}_{2}\rightarrow\delta$,

[TABLE]

and likewise for $\tilde{h}_{1}\rightarrow\delta$. In each of these cases the resulting window in one dimension is the product of the two broader filters. The estimator therefore produces the same result as when using two filters: the first given by a narrow filter and the second by the root product of the remaining two. Without loss of generality, the following sections therefore consider estimators for which two of the three analysis filters use the same window. We consider permutations of the bandwidth of the three analysis filters applied to the three bands sorted by frequency: $\omega_{1}\leq\omega_{2}<\nu_{2}$. Shapes of filter kernels resulting from each combination are depicted in Figure 9.

B.2.1 Two Filters: Narrow, Narrow, Broad

When $\tilde{g}$ is wider than $\tilde{h}$ in (52), $\tilde{H}$ becomes symmetric so that $W$ is smoothed equally along the $\nu_{1}$ and $\nu_{2}$ axes. As $\tilde{g}(\nu_{2})\rightarrow 1$ we obtain the approximation:

[TABLE]

An estimate of the bispectrum with radially symmetric smoothing is therefore recovered by:

[TABLE]

Note that the radial symmetry here does not obey the symmetry relations of (34), so the estimate does not exactly recover the symmetry of the bispectrum. On the other hand, (54) suggests a less cumbersome approach to estimation which foregoes the summation over the product of three separate bands. The estimator is obtained instead from the frequency-domain covariance of the time-frequency decomposition after remodulation and weighting by the original signal. In practice, however, the time-frequency decomposition is often downsampled, so it will still be necessary to either apply an anti-aliasing filter to $X$ appropriate for the sampling rate and frequency range under consideration, or oversample $\chi_{h}$. Both options tend to negate any computational advantage that might otherwise be gained from such a simplification.

B.2.2 Two filters: Broad, Broad, Narrow

The estimator in Eq. (54) results in a smoothing window with radial symmetry, but the symmetry came about from greater smoothing in the $\nu_{2}$ direction entailing a sacrifice of resolution. One may consider the converse possibility, allowing the bandwidth of $g$ (in 52) to be narrower than that of $h$. This produces a smoothing window with the limit $\tilde{g}(\nu_{2})\rightarrow\delta(\nu_{2})$. In the approach to this limit one obtains the following approximation:

[TABLE]

Provided that $g$ is narrower than $h$, this estimator allows the smoothing along $\nu_{1}$ and $\nu_{2}$ to be separately controlled.

B.2.3 Two filters: Narrow, Broad, Broad

In the limit $\tilde{h}_{1}\rightarrow\delta$ in (52), we have $g(\nu_{2})\rightarrow g(\omega_{2})$. Letting $h_{2}=g$, without loss of generality, we have

[TABLE]

The same result holds for the BNB estimator swapping the arguments. The smoothing kernel here is related to the previous examples through a 45 degree rotation in the bispectral plane. One important potential drawback of this family of estimators is that the outcome will not reflect the symmetry of the bispectrum about the diagonal, $\omega_{1}=\omega_{2}$. Measures of phase-amplitude coupling relate to bispectral estimators of this type, as shown in C.

Appendix C Phase-Power Coherence as a Smoothed Bispectral Estimator

The following appendix provides a mathematical demonstration of the equivalence between phase-power coherence (PhPC) and the NBB class of bispectral estimators. PhPC is calculated from the cross spectrum between the squared analytic envelope at one frequency and the original signal at another. The goal is normally to isolate the fast oscillation within the first band and the slow oscillation within the second. The squared envelope of the FO band, centered at $\gamma$ filtered with the analysis window, $h$, is given by

[TABLE]

PhPC is calculated from the cross spectrum between $P$ and $X$, obtained with a second analysis filter, $g$, which will be of a narrower band than $h$, targeted to the range of the SO around frequency $\theta$:

[TABLE]

where $g^{(2)}$ denotes the autocorrelation of $g$. Expanding $P$ in (58) according to (57), we have

[TABLE]

Following the foregoing series of substitutions, the integral over $t$ in (59) no longer depends on the signal (which is now indexed by $r$,$s$ and $u$). Pulling out the term that integrates over $t$ yields a smoothing kernel:

[TABLE]

from which, setting $r=s+\tau_{2}$ and $u=s+\tau_{1}$:

[TABLE]

Substituting back into (59) and applying (10), one arrives at

[TABLE]

which with another change of variables, $\omega_{1}\equiv-2\lambda_{2}$ and $\omega_{2}\equiv\lambda_{2}-\lambda_{1}+\gamma$ and applying the symmetry relation from Eq. (35), $W(-\omega_{2}-\omega_{1},\omega_{2},z)=W(\omega_{1},\omega_{2},z)$, this becomes

[TABLE]

Finally, because in measuring phase-amplitude coupling, the SO filter, $g$, is selected with a narrower bandwidth than the FO filter, $h$, approaching the limit $\tilde{g}\rightarrow\delta$, we may use the following good approximation

[TABLE]

C.1 Relationship between PhPC and the Bispectrum

The equivalent smoothing window for PhPC estimators, given by Eq. (64), very nearly describes an NBB smoothing kernel like the one in (56). The only nontrivial difference arises from the fact that the smoothing kernel varies as a function of the SO frequency, $\theta$. The effect is easy to appreciate for the Gaussian window of bandwidth $\sigma_{\omega}$:

[TABLE]

For the Gaussian case, the result is the same as smoothing with an NBB window, only shifted along $\omega_{2}$ by $\theta/2$ and attenuated as a function of the phase-providing frequency by $e^{-\theta^{2}/2}$ without any other change in the size or shape of the smoothing window. More generally, if $h$ is symmetric, then $\tilde{h}^{*}\left(\omega_{2}-\gamma+\theta\right)\tilde{h}\left(\omega_{2}-\gamma\right)$ is symmetric and centered at $\gamma-\theta/2$.

PAC measures require the FO filter to have a bandwidth at least as great as the range of SO frequencies of interest, because the spectral range of fluctuations in the power amplitude is limited to the bandwidth of the filter used to extract them, so that $\theta$ must remain small relative to the bandwidth of $h$ and $\left|h(\omega+\theta)\right|\approx\left|h(\omega)\right|$[4]. Precisely within this range, PAC estimates may be treated directly as NBB estimators of the bispectrum, which use the window:

[TABLE]

This constraint that limits the spectral range of the SO according to the bandwidth of the FO seems at first glance to be a fundamental consideration in measuring PAC, but the preceding analysis shows it to be in fact a dispensable property of the estimator unrelated to anything inherent in the underlying measured quantity.

One approach to circumventing the spectral limitation of PAC estimators modifies the FO filters used according to the SO frequency. For example, the problem is sometimes addressed by scaling the bandwidth of the FO filter, $h$, by the center frequency of the SO filter [8]: $\tilde{h}_{\theta}(\omega)=\tilde{h}\left(\frac{\omega}{\theta}\right)$, giving

[TABLE]

This adjustment removes the attenuation with SO frequency by ensuring that the bandwidth of the FO filter is appropriate for modulation frequencies in the range of $\theta$. It still results in diminishing resolution of the FO with increasing SO, which is not a necessary constraint on bispectral estimators.

An apparent disadvantage of standard bispectral estimators might be noted at this point, which is the shift of the FO spectrum by $\theta/2$. This shift from the FO frequency is however not essential to the bispectral estimator and can be removed through the symmetrizing change of variables given in Eq. (43).

It should be reemphasized that these differences between PhPC and bispectral estimators involve rather arbitrary properties of the respective smoothing kernels and do not reflect any essential difference in the quantity measured. It might be argued that the limit on the combined range and resolution of the estimator imposed by PAC helps tune it to bispectral characteristics of PAC. Such tuning may improve signal-to-noise sensitivity to PAC as well as computational efficiency, but it nevertheless forces some a priori assumptions about the bandwidth and spectral range of PAC. More importantly, smoothing may sacrifice information vital for distinguishing between alternative forms of cross-frequency coupling.

C.2 Interregional PAC and Cross Bicoherence

The foregoing discussion extends naturally to bivariate cross-frequency coupling involving two signals, for example, obtained from two separate recording channels. The bivariate extension of PAC measures describes the dependence between amplitude in one signal with phase in another. Using the same argument as in the single-channel case, measures of interregional PAC may be equated to the cross bispectrum of the form:

[TABLE]

where in the language of PAC, $i$ indexes the “phase-providing” signal and $j$ the “amplitude-providing” signal.

Appendix D Transient Signal Model

The following appendix presents a basic signal model that will serve as a starting point in describing bispectral features associated with phase-amplitude coupling and other third-order signal features. It describes the case when recurring transient features with characteristic but unknown spectra lie embedded at unknown times in the observed record. The generating process can be separated into two parts: a (1) point process, $N$, whose increments determine the times, $\tau_{i}$, at which (2) some random or deterministic transient feature, $f_{i}(t)$ is embedded in the signal. The first process is concerned only with timing and the second only with the emitted waveform. The spectral representation of such a signal is

[TABLE]

The following sections consider the bispectrum that results from this process when $f_{i}$ is independent of time and of N and, in the first case, deterministic, and, in the second, i.i.d. random.

D.1 Deterministic Features

If the $f_{i}$ in (69) are deterministic with $f_{i}=f$, it follows from the convolution property (Eq. 36)

[TABLE]

with

[TABLE]

For a “simple” point process, events coincide with vanishing probability, so that $E[(dN(t))^{k}]=E[dN(t)]$. For a process with the first three moments defined, $\mu_{1}$, $\mu_{2}$, $\mu_{3}$ [7, 16], the expectation in Eq. (71) yields

[TABLE]

The bispectrum for this process is then

[TABLE]

where $\lambda$ is the average rate ($N/T\rightarrow\lambda$).

For a stationary homogenous Poisson driving process (having constant $\mu(t)=\lambda$), Eq. (73) reduces to

[TABLE]

so that, neglecting the probability mass at the origins and along $\omega_{1}+\omega_{2}=0$,

[TABLE]

This quantity is nonzero when $f$ contains harmonics or is otherwise spectrally broad, so that its support covers some $\omega_{1},\omega_{2}$, and $\omega_{1}+\omega_{2}$. In fact, because $f$ is transient by assumption, this condition is already given: $f$ must decay to zero within some finite time window, $\Delta T$, and therefore its spectrum contains at a minimum a main lobe and side lobes spaced at $2\pi/\Delta T$. But if $f$ has zero mean, a high center frequency and is highly oscillatory with an otherwise narrow spectrum, the product may still be negligible. The bispectrum will contain larger peaks when $f$ has a broad spectrum on the order of 2/3 its center frequency, non-zero mean, or contains one or more harmonic complexes.

D.2 Random Features

Extending the preceding analysis to the case when each $f_{i}$ in (69) is itself the realization of a random process independent of $N$,

[TABLE]

If the process generating $f_{i}$ is independent of $N$ and time, so that $\tilde{f}_{i}\sim\tilde{f}_{j}$, for all $i,j$, and $N$ simple, with the first three moments defined, $\mu_{1}$, $\mu_{2}$, $\mu_{3}$ [7, 16]:

[TABLE]

with the brackets denoting expectation, and a normalization by the cumulative event count, $N$, is left implicit. This reduces to $B_{X}=B_{f}B_{N}$ if the error of $f_{i}$ involves independent additive Gaussian noise.

This division of labor allows one to separate the spectral contributions of local features described by $f$ from those of the large-scale driving process; these respective contributions can be distinguished in the bispectrum444The Poisson formulation adopted here for the driving point process is mildly restrictive in that it assumes the event count within a given interval follows a Poisson distribution. Periodicity in the driving process therefore does not entail a train of evenly spaced events, but instead a train of evenly spaced event clusters with Poisson-distributed size. An alternative point-process model that conditions on history, such as a renewal process, might handle periodicity more naturally in this setting, but elucidating the spectra of such processes is more technically involved, making them less useful for the present purpose. One way to address the case of simple uniform periodicity with the Poisson model, which leaves the foregoing spectral analysis unchanged, is to allow $\lambda$ to take the form of transient bursts of infinitesimal duration, each burst generating a cluster with Poisson-distributed size. The resulting signal then contains $f$ at the cluster times modulated by cluster size. The variance of amplitude in $f$ introduced by this weighting with cluster sizes can be made arbitrarily small by scaling $\lambda$ so as to increase cluster size, with the inverse scaling applied $f$ to keep amplitude fixed. . In the case of phase-amplitude coupling, $f$ will be treated as the sum of two components: a single burst of fast oscillations, $f_{\mathrm{FO}}$, whose spectrum is characteristically broad and smooth, and a slow oscillation, $f_{\mathrm{SO}}$, with a narrow spectrum. In constructing the full signal, narrowband large-scale features can be viewed as the consequence of a filter applied to the train of impulses generated by the point process, where $f_{\mathrm{SO}}$ plays the role of filter function.