# Time-frequency analysis of event-related brain recordings: Connecting power
of evoked potential and inter-trial coherence

Jonas Benhamou, Michel Le Van Quyen, and Guillaume Marrelec All authors are
with the Laboratoire d’imagerie biomédicale, LIB, Sorbonne Université, CNRS,
Inserm, Paris, France. J. Benhamou is also with the École Nationale Supérieure
de Techniques Avancées, ENSTA, Saclay, France. Manuscript received ???;
revised ???. Copyright (c) 2021 IEEE. Personal use of this material is
permitted. However, permission to use this material for any other purposes
must be obtained from the IEEE by sending an email to pubs-
<EMAIL_ADDRESS>
###### Abstract

Objective. In neuroscience, time-frequency analysis has been used to get
insight into brain rhythms from brain recordings. In event-related protocols,
one applies it to investigate how the brain responds to a stimulation repeated
over many trials. In this framework, three measures have been considered: the
amplitude of the transform for each single trial averaged across trials,
avgAMP; inter-trial phase coherence, ITC; and the power of the evoked
potential transform, POWavg. These three measures are sensitive to different
aspects of event-related responses, ITC and POWavg sharing a common
sensitivity to phase resetting phenomena. Methods. In the present manuscript,
we further investigated the connection between ITC and POWavg using
theoretical calculations, a simulation study and analysis of experimental
data. Results. We derived exact expressions for the relationship between
POWavg and ITC in the particular case of the S-transform of an oscillatory
signal. In the more general case, we showed that POWavg and ITC are connected
through a relationship that roughly reads POWavg $\approx$ avgAMP2 $\times$
ITC2. This result was confirmed on simulations. We finally compared the
theoretical prediction with results from real data. Conclusion. We showed that
POWavg and ITC are related through an approximate, simple relationship that
also involves avgAMP. Significance. The presented relationship between POWavg,
ITC, and avgAMP confirms previous empirical evidence and provides a novel
perspective to investigate evoked brain rhythms. It may provide a significant
refinement to the neuroscientific toolbox for studying evoked oscillations.

###### Index Terms:

Time-frequency transform; inter-trial coherence; brain recordings;
electroencephalography (EEG); magnetoencephalography (MEG); High frequency
oscillations (HFOs).

##  I Introduction

In signal processing, time-frequency (TF) analysis is a generic term that
encompasses a variety of methods—such as the short-time Fourier transform [1],
wavelets [2], and the S-transform [3]—for the analysis of signals whose
frequency features vary with time [4, 5]. All methods have in common that they
map a one-dimensional real or complex signal $x(t)$ into a two-dimensional
complex-valued function $T_{x}(t,f)$, called time-frequency transform. This
transform provides, for each time $t$ and frequency $f$, a modulus
$|T_{x}(t,f)|$ (also called amplitude) and an argument $\arg[T_{x}(t,f)]$ or
$\theta_{x}(t,f)$ (also called phase). The terms amplitude and phase are meant
to emphasize the expected relation between the amplitude and phase of an
oscillatory signal and the modulus and argument of its time-frequency
transform[2, Chap. 4]. Time-frequency analysis has been used in a wide variety
of research fields, from geophysics to engineering to medicine [6].

In neuroscience, neuronal assemblies are hypothesized to generate oscillating
electromagnetic fields within specific frequency ranges, also called brain
rhythms. There has been overwhelming evidence that these brain rhythms can be
captured by in-vivo brain recording techniques, from noninvasive methods, such
as electroencephalography (EEG) or magnetoencephalography (MEG), to invasive
methods, such as intracranial EEG (iEEG) or microelectrode recordings
(measuring local field potentials, LFPs) [7, 8, 9, 10, 11, 12]. A common
procedure to investigate brain rhythms in brain recordings is the so-called
event-related protocol, where one records how the brain responds to a given
stimulation over many repetitions, called trials. For $N$ trials, one then has
a collection of $N$ signals $x_{n}(t)$, $n=1,\dots,N$. The brain activity
measured in response to each stimulation is called an event-related response.
The exact interpretation of event-related responses remains an open issue in
neuroscience. It was suggested that they may result from a variable
combination of two phenomena: stimulus-evoked neuronal activity and stimulus-
induced phase resetting of ongoing neuronal dynamics [13, 14, 15, 16]. In the
former, the stimulus leads to the activation of specific neuronal populations,
which contribute to new time-locked oscillations at each trial, translating
into new power contributions to the ongoing brain activity. In the latter, the
stimulus resets the phase of pre-existing oscillations so that a portion of
the ongoing activity becomes transiently phase-locked in a particular
frequency range. To analyze brain recordings from event-related protocols and
be able to discriminate between evoked responses and phase resetting / induced
responses, increasingly sophisticated methods have been used, including time-
frequency analysis [17, 18, 19]. So far, three time-frequency measures have
been considered. A first measure is the average time-frequency amplitude,
i.e., the amplitude $|T_{x_{n}}(t,f)|$ of the transform for each single trial
$n$ averaged across trials, [20, Chap. 9]

$\mathrm{avgAMP}=\frac{1}{N}\sum_{n=1}^{N}\left|T_{x_{n}}(t,f)\right|.$ (1)

A second measure is the so-called inter-trial phase coherence (ITC), which
quantifies the phase consistency of brain responses over experimental trials
[15]

$\mathrm{ITC}=\left|\frac{1}{N}\sum_{n=1}^{N}e^{i\theta_{x_{n}}(t,f)}\right|.$
(2)

A last measure applies the time-frequency transform to the evoked potential,
i.e., the average through trials of all signals, and extracts its power: If
$\overline{x_{n}}(t)$ is the the signal averaged over all $N$ trials, then
[21]

$\mathrm{POWavg}=\left|T_{\overline{x_{n}}}(t,f)\right|^{2}.$ (3)

Importantly, these three measures are sensitive to different aspects of event-
related responses: evoked responses for avgAMP, and phase resetting for ITC
and POWavg [21, 15]. To date, with these three measures, evidence for both
evoked and induced oscillations has been reported in multiple cognitive
protocols generating time-locked responses to specific sensory, cognitive or
motor events [22, 23, 24].

We are here interested in the observed common sensitivity of ITC and POWavg to
induced responses as well as the overall similarity of their outputs [25].
More precisely, we investigate the theoretical underpinning of this empirical
evidence. Intuitively, we would expect a low ITC to correspond to a low
POWavg, as distinct components with phases that are not consistent across
trials (i.e., with low ITC) should cancel out in the averaging process,
leading to an average signal with low POWavg. In the present manuscript, we
further investigate the connection between ITC and POWavg using theoretical
calculations and a simulation study. We show that POWavg is approximately
equal to the product of avgAMP squared by ITC squared,

$\mathrm{POWavg}\approx\mathrm{avgAMP}^{2}\times\mathrm{ITC}^{2},$ (4)

under certain conditions. More precisely, we first derive several
approximations of the relationship between POWavg, ITC and avgAMP in the
particular case of the S-transform of an oscillatory signal with increasing
levels of complexity. We then investigate the more general case and derive
both the exact form of this relationship—see (71)—, as well as the assumptions
leading to the simpler form of (4). We show that this result also holds in the
different setting of synthetic data simulating high frequency oscillations. We
finally investigate the validity of this result for experimental data.

The outline of the manuscript is the following. In Section II, we derive
several variants of (4). In Section III, we perform a time-frequency analysis
of synthetic data originating from a simulation study. The real data are
analyzed in Section IV. Further issues are discussed in Section V.

##  II Theoretical developments

In this section, we provide the theoretical connection between ITC and POWavg,
thus specifying the form of (4). We start by stating some basic results
(Section II-A). We then expose the general framework for time-frequency
analysis (Section II-B) and give some results regarding avgAMP, ITC, and
POWavg (Section II-C). We then provide an intuition for the existence of a
relationship (Section II-E1) and derive the its exact form, both in the
particular case of the S-transform of an oscillatory signal (Section II-D),
and in the general case (Section II-E).

###  II-A Some basic results

We need the following results.

####  II-A1 Complex numbers

First, for any complex number $z$, absolute value $|z|$ and conjugate $z^{*}$
are related by $|z|^{2}=zz^{*}$. We also have $z+z^{*}=2\Re(z)$, where
$\Re(z)$ is the real part of $z$.

####  II-A2 Random variables

If $K$ random variables $X_{1},\dots,X_{K}$ are independent, then the
expectation of their product is equal to the product of their expectations
[26, §2.2.3]

$\mathrm{E}\left(\prod_{k=1}^{K}X_{k}\right)=\prod_{k=1}^{K}\mathrm{E}\left(X_{k}\right).$
(5)

For two complex-valued random variables $X$ and $Y$, we have

$\mathrm{Cov}(X,Y)=\mathrm{E}(XY^{*})-\mathrm{E}(X)\mathrm{E}(Y)^{*}.$ (6)

In particular, this entails that

$\mathrm{Var}(X)=\mathrm{E}\left(|X|^{2}\right)-\left|\mathrm{E}(X)\right|^{2}.$
(7)

A circular random variable is a random variable that, like an angle, is
defined on the unit circle, i.e., whose value is only relevant modulo $2\pi$.
A key quantity for any circular random variable $\theta$ is its circular mean,
defined as

$\mathrm{E}\left(e^{i\theta}\right)=\int\mathrm{p}(\theta)\,e^{i\theta}\,\mathrm{d}\theta.$
(8)

The argument of $\mathrm{E}\left(e^{i\theta}\right)$ is the mean angle or mean
direction, while the modulus of $\mathrm{E}\left(e^{i\theta}\right)$ is the
mean resultant length. It lays between 0 and 1 and is a measure of
concentration of $\mathrm{p}(\theta)$.

A circular variable $\theta$ is said to follow a von Mises distribution with
mean direction $\theta_{0}$ and concentration parameter $\kappa$, denoted
$\theta\sim\mathrm{VonMises}(\theta_{0},\kappa)$, if its distribution is given
by [27, §3.5.4]

$\mathrm{p}(\theta)=\frac{1}{2\pi
I_{0}(\kappa)}e^{\kappa\cos(\theta-\theta_{0})},$

where $I_{0}(\kappa)$ is the modified Bessel function of order 0,

$I_{0}(\kappa)=\frac{1}{2\pi}\int_{-\pi}^{\pi}e^{\kappa\cos(\theta-\theta_{0})}\,\mathrm{d}\theta.$
(9)

The usual Gaussian distribution with mean $\mu$ and variance $\sigma^{2}$ is
denoted by $\mathcal{N}(\mu,\sigma^{2})$.

####  II-A3 Unbiased estimator

For the estimator $\hat{\theta}$ of a quantity $\theta$, the bias is defined
as the difference between the sampling expectation of $\hat{\theta}$ and
$\theta$,

$\mathrm{bias}(\hat{\theta})=\mathrm{E}(\hat{\theta})-\theta.$ (10)

If the bias is equal to 0, then the estimator is said to be unbiased. If the
bias vanishes as the data size tends to infinity, then the estimator is said
to be asymptotically unbiased.

####  II-A4 Bachmann–Landau notation

Finally, a function $f(N)$ is said to be $O(1/N)$ if it is bounded by a
function proportional to $1/N$, that is, for which there exists an $N_{0}$ and
a $k>0$ such that

$\forall N>N_{0}\quad\left|f(N)\right|<\frac{k}{N}.$

###  II-B Time-frequency analysis

We here introduce the general framework for time-frequency analysis. More
specifically, we consider $N$ signals $x_{n}(t)$, $n=1,\dots,N$, each signal
being measured as the response to a repetition of the stimulus. The
$x_{n}(t)$’s are assumed to be $N$ independent and identically distributed
(i.i.d.) realizations of a signal $x(t)$. Time-frequency analysis maps each
$x_{n}(t)$ into a complex time-frequency function $T_{x_{n}}(t,f)$ that we
express as

$T_{x_{n}}(t,f)=|T_{x_{n}}(t,f)|\,e^{\theta_{x_{n}}(t,f)}.$ (11)

$|T_{x_{n}}(t,f)|$ is the amplitude of the time-frequency transform and
$|T_{x_{n}}(t,f)|^{2}$ its power; $\theta_{x_{n}}(t,f)$ is the phase of the
time-frequency transform. Since the $x_{n}(t)$’s are i.i.d., so are their
time-frequency transforms. As a consequence, the $T_{x_{n}}(t,f)$’s can be
considered as i.i.d. realizations of a random time-frequency transform denoted
$T_{x}(t,f)$ with amplitude $|T_{x}(t,f)|$ and phase $\theta_{x}(t,f)$.

###  II-C Three quantities of interest

We can now go back to our three quantities of interest: avgAMP, ITC and
POWavg. To emphasize the fact that they depend on a a series of signals
$x_{1}(t)$, …, $x_{N}(t)$, as well as on the time $t$ and frequency $f$ at
which the time-frequency transform is computed, we write
$\mathrm{avgAMP}_{x_{1:N}}(t,f)$, $\mathrm{ITC}_{x_{1:N}}(t,f)$ and
$\mathrm{POWavg}_{x_{1:N}}(t,f)$.

####  II-C1 avgAMP

avgAMP is defined as in (1). It can be shown (see §1.1 of Supplementary
Material #1) that the expectation of avgAMP2 can be expressed as

$\mathrm{E}\left[\mathrm{avgAMP}_{x_{1:N}}(t,f)^{2}\right]=\mathrm{E}\left[\left|T_{x}(t,f)\right|\right]^{2}+\frac{1}{N}\mathrm{Var}\left[\left|T_{x}(t,f)\right|\right].$
(12)

This first shows that $\mathrm{E}(\mathrm{avgAMP}^{2})$ is only a function of
$|T_{x}(t,f)|$, the amplitude of the time-frequency transform—and therefore
not of the phase $\theta_{x}(t,f)$. Also, from (12), we obtain that
$\mathrm{E}[\mathrm{avgAMP}_{x_{1:N}}(t,f)^{2}]$ tends to
$\mathrm{E}[|T_{x}(t,f)|]^{2}$ as $N\to\infty$, i.e., avgAMP2 is an
asymptotically unbiased estimator of $\mathrm{E}[|T_{x}(t,f)|]^{2}$. For
finite $N$, its bias is given by $\mathrm{Var}[|T_{x}(t,f)|]/N$. A last
expression that can be derived from (12) is the following asymptotic
expression obtained by applying Bachmann–Landau notation:

$\mathrm{E}\left[\mathrm{avgAMP}_{x_{1:N}}(t,f)^{2}\right]=\mathrm{E}\left[\left|T_{x}(t,f)\right|\right]^{2}+O\left(\frac{1}{N}\right).$
(13)

####  II-C2 ITC

ITC is defined as in (2). It can be shown (see §1.2 of Supplementary Material
#1) that the expectation of ITC2 yields

$\mathrm{E}\left[\mathrm{ITC}_{x_{1:N}}(t,f)^{2}\right]=\left|\mathrm{E}\left[e^{i\theta_{x}(t,f)}\right]\right|^{2}+\frac{1}{N}\mathrm{Var}\left[e^{i\theta_{x}(t,f)}\right].$
(14)

This result shows that $\mathrm{E}(\mathrm{ITC}^{2})$ is a only function of
$\theta_{x}(t,f)$, the phase of the time-frequency transform—and therefore not
of the amplitude $|T_{x}(t,f)|$. It also shows that
$\mathrm{E}[\mathrm{ITC}_{x_{1:N}}(t,f)^{2}]$ tends to
$|\mathrm{E}[e^{i\theta_{x}(t,f)}]|^{2}$ when $N\to\infty$, that is, ITC2 it
is an asymptotically unbiased estimator of
$|\mathrm{E}[e^{i\theta_{x}(t,f)}]|^{2}$, the squared mean resultant length of
$\theta_{x}(t,f)$. For finite $N$, the bias is given by
$\mathrm{Var}[e^{i\theta_{x}(t,f)}]/N$. Finally, we obtain from (14) the
asymptotic expansion

$\mathrm{E}\left[\mathrm{ITC}_{x_{1:N}}(t,f)^{2}\right]=\left|\mathrm{E}\left[e^{i\theta_{x}(t,f)}\right]\right|^{2}+O\left(\frac{1}{N}\right).$
(15)

####  II-C3 POWavg

POWavg is defined as in (3), with

$\overline{x_{n}}(t)=\frac{1}{N}\sum_{n=1}^{N}x_{n}(t).$ (16)

Due to the linearity of the time-frequency analysis,
$T_{\overline{x_{n}}}(t,f)$ is equal to the average of the time-frequency
transforms,

$T_{\overline{x_{n}}}(t,f)=\frac{1}{N}\sum_{n=1}^{N}T_{x_{n}}(t,f),$ (17)

so that

$\mathrm{POWavg}_{x_{1:N}}(t,f)=\left|\frac{1}{N}\sum_{n=1}^{N}T_{x_{n}}(t,f)\right|^{2}.$
(18)

It can be shown (see §1.3 of Supplementary Material #1) that the expectation
of POWavg reads

$\mathrm{E}\left[\mathrm{POWavg}_{x_{1:N}}(t,f)\right]=\left|\mathrm{E}\left[T_{x}(t,f)\right]\right|^{2}+\frac{1}{N}\mathrm{Var}\left[T_{x}(t,f)\right].$
(19)

Importantly, this result shows that, unlike avgAMP and ITC, POWavg depends on
both the amplitude and the phase of the time-frequency transform through the
modulus and argument of $T_{x}(t,f)$—recall (11). From (19), we also have that
$\mathrm{E}[\mathrm{POWavg}_{x_{1:N}}(t,f)]$ tends to
$|\mathrm{E}[T_{x}(t,f)]|^{2}$ as $N\to\infty$, leading to the conclusion that
POWavg is an asymptotically unbiased estimator of
$\left|\mathrm{E}\left[T_{x}(t,f)\right]\right|^{2}$. For finite $N$, the bias
is given by $\mathrm{Var}[T_{x}(t,f)]/N$. Similarly to the two other measures,
we asymptotically have

$\mathrm{E}\left[\mathrm{POWavg}_{x_{1:N}}(t,f)\right]=\left|\mathrm{E}\left[T_{x}(t,f)\right]\right|^{2}+O\left(\frac{1}{N}\right).$
(20)

####  II-C4 Summary

Results regarding the general properties of the expectation of avgAMP2, ITC2
and POWavg are summarized in §1.4 of Supplementary Material #1.

###  II-D Investigation of oscillatory model

In this section, we investigate the behaviors of avgAMP, ITC, and POWavg in
the particular case of an oscillatory signal.

####  II-D1 Model

We assume that $x(t)$ is a cosine with amplitude $\Omega$, frequency $\nu$,
and phase $\phi$, i.e., of the form

$x(t)=\Omega\cos(2\pi\nu t+\phi).$ (21)

The parameters $\Omega$, $\nu$, and $\phi$ will either be assumed be constant,
or follow a certain distribution over trials depending on the model (see §2.6
of Supplementary Material #1 for a summary of models and results).

####  II-D2 S-transform

The S-transform is defined as [3]

$T_{x}(t,f)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{+\infty}x(u)\,|f|\,e^{-\frac{f^{2}(u-t)^{2}}{2}}\,e^{-2i\pi
fu}\,\mathrm{d}u.$ (22)

This is of the form

$\frac{1}{\sqrt{2\pi\alpha^{2}}}\int_{-\infty}^{+\infty}x(u)\,e^{-\frac{(u-t)^{2}}{2\alpha^{2}}}\,e^{-2i\pi
fu}\,\mathrm{d}u$

with111Note that $\alpha$ is sometimes defined as $1/(f\sqrt{2\pi})$, so that
$1/\sqrt{2\pi\alpha^{2}}=|f|$ and $1/(2\alpha^{2})=\pi f^{2}$. $\alpha=1/f$.
The S-transform can therefore be seen as a band-pass filter or a windowed
Fourier transform with a window whose width decreases with increasing
frequency.

####  II-D3 time-frequency transform of signal

Direct calculation shows that the S-transform of $x(t)$ is given by (see §2.2
of Supplementary Material #1)

$\displaystyle T_{x}(t,f)$ $\displaystyle=$
$\displaystyle\frac{\Omega}{2}e^{-\frac{1}{2}(2\pi)^{2}\left(1-\frac{\nu}{f}\right)^{2}}e^{i[\phi-2\pi(f-\nu)t]}$
$\displaystyle\qquad+\frac{\Omega}{2}e^{-\frac{1}{2}(2\pi)^{2}\left(1+\frac{\nu}{f}\right)^{2}}e^{i[-\phi-2\pi(f+\nu)t]}.$

The amplitude and phase of this complex number do not have simple exact forms.
Still, it can be shown that a very good approximation of $T_{x}(t,f)$ can be
obtained by only keeping the first term of the right-hand side of
(LABEL:eq:osc:gen:tf). More precisely, we have (see §2.3 of Supplementary
Material #1)

$T_{x}(t,f)=T_{x}^{\dagger}(t,f)\left[1+\epsilon(t,f)\right]$ (24)

with

$T_{x}^{\dagger}(t,f)=\frac{\Omega}{2}e^{-\frac{1}{2}(2\pi)^{2}\left(1-\frac{\nu}{f}\right)^{2}}e^{i[\phi-2\pi(f-\nu)t]}$
(25)

and

$\left|\epsilon(t,f)\right|=e^{-8\pi^{2}\frac{\nu}{f}}.$ (26)

The modulus and argument of $T_{x}^{\dagger}(t,f)$ are given by

$\left|T_{x}^{\dagger}(t,f)\right|=\frac{\Omega}{2}e^{-\frac{1}{2}(2\pi)^{2}\left(1-\frac{\nu}{f}\right)^{2}}$
(27)

and

$\arg\left[T_{x}^{\dagger}(t,f)\right]=\phi-2\pi(f-\nu)t,$ (28)

respectively. Since we are mostly interested in frequencies in the range
$0<f<10\nu$, the relative error $\epsilon(t,f)$ is bounded by
$|\epsilon(t,f)|<3.8\times 10^{-4}$ (see §2.3 of Supplementary Material #1).
The amplitude of $T_{x}(t,f)$ is given by

$\left|T_{x}(t,f)\right|=\frac{\Omega}{2}e^{-\frac{1}{2}(2\pi)^{2}\left(1-\frac{\nu}{f}\right)^{2}}\left[1+\epsilon_{m}(t,f)\right]$
(29)

and the phase of $T_{x}(t,f)$ by

$\theta_{x}(t,f)\approx\phi-2\pi(f-\nu)t+\epsilon_{a}(t,f).$ (30)

For more detail on the errors $\epsilon_{m}(t,f)$ and $\epsilon_{a}(t,f)$, see
§2.3 of Supplementary Material #1.

As expected, there is a strong connection between the amplitude and phase of
the signal and those of the time-frequency transform: the amplitude of the
time-frequency transform, (29), only depends on the amplitude (and frequency)
of the signal, while the phase of the time-frequency transform, (30), only
depends on the phase (and frequency) of the signal.

Note that the amplitude of the time-frequency transform is not a function of
$t$. As a function of $f$, it is maximal for $f=\nu$ (and equal to $\Omega/2$)
and decreases when $f$ moves away from $\nu$. The phase of the time-frequency
transform is a (bilinear) function of both $t$ and $f$ that is equal to $\phi$
for pairs $(t,f)$ of the form either $(0,f)$ for any $f$ or $(t,\nu)$ for any
$t$.

####  II-D4 Simplest model of repetition

We now consider $N$ repetitions of the oscillatory signal of (21) with
constant amplitude ($\Omega=\Omega_{0}$) and frequency ($\nu=\nu_{0}$) but
varying phases $\phi_{n}$, $n=1,\dots,N$,

$x_{n}(t)=\Omega_{0}\cos(2\pi\nu_{0}t+\phi_{n}),\qquad n=1,\dots,N.$ (31)

We assume that the $\phi_{n}$’s are i.i.d. repetitions of
$\phi\sim\mathrm{VonMises}(\phi_{0},\kappa)$. We define the circular mean as
in (8)

$\mathrm{E}(e^{i\phi})=\rho e^{i\phi_{0}},$ (32)

with $\phi_{0}$ the mean direction and $\rho$ the mean resultant length.
Applying the approximation of (25), the time-frequency transform of each
signal $x_{n}(t)$ is given by

$T_{x_{n}}^{\dagger}(t,f)=\frac{\Omega_{0}}{2}e^{-\frac{1}{2}(2\pi)^{2}\left(1-\frac{\nu_{0}}{f}\right)^{2}}e^{i[\phi_{n}-2\pi(f-\nu_{0})t]}.$
(33)

The amplitude is equal to

$\left|T_{x_{n}}^{\dagger}(t,f)\right|=\frac{\Omega_{0}}{2}e^{-\frac{1}{2}(2\pi)^{2}\left(1-\frac{\nu_{0}}{f}\right)^{2}},$
(34)

and the phase to

$\arg\left[T_{x_{n}}^{\dagger}(t,f)\right]=\phi_{n}-2\pi(f-\nu_{0})t.$ (35)

In this model, the phase depends on $n$ and is random (as it depends on
$\phi_{n}$), whereas the amplitude does not depend on $n$ and is
deterministic. From (1) and (34), we obtain

$\mathrm{avgAMP}_{x_{1:N}}^{\dagger}(t,f)=\frac{\Omega_{0}}{2}e^{-\frac{1}{2}(2\pi)^{2}\left(1-\frac{\nu_{0}}{f}\right)^{2}},$
(36)

and, from (2) and (35),

$\mathrm{ITC}_{x_{1:N}}^{\dagger}(t,f)=\left|\frac{1}{N}\sum_{n=1}^{N}e^{i\phi_{n}}\right|.$
(37)

Besides, the time-frequency transform of the average signal can be calculed
using (17) and (33), yielding

$T_{\overline{x_{n}}}^{\dagger}(t,f)=\frac{\Omega_{0}}{2}e^{-\frac{1}{2}(2\pi)^{2}\left(1-\frac{\nu_{0}}{f}\right)^{2}}\,e^{-2i\pi(f-\nu_{0})t}\,\frac{1}{N}\sum_{n=1}^{N}e^{i\phi_{n}},$
(38)

which has amplitude

$\left|T_{\overline{x_{n}}}^{\dagger}(t,f)\right|=\frac{\Omega_{0}}{2}e^{-\frac{1}{2}(2\pi)^{2}\left(1-\frac{\nu_{0}}{f}\right)^{2}}\left|\frac{1}{N}\sum_{n=1}^{N}e^{i\phi_{n}}\right|$
(39)

(the phase is irrelevant for our calculation). This, together with (3), leads
to

$\mathrm{POWavg}_{x_{1:N}}^{\dagger}(t,f)=\left[\frac{\Omega_{0}}{2}e^{-\frac{1}{2}(2\pi)^{2}\left(1-\frac{\nu_{0}}{f}\right)^{2}}\right]^{2}\left|\frac{1}{N}\sum_{n=1}^{N}e^{i\phi_{n}}\right|^{2}.$
(40)

From (36), (37), and (40), it is straightforward to see that we have the
following relationship

$\mathrm{POWavg}^{\dagger}=\left(\mathrm{avgAMP}^{\dagger}\right)^{2}\times\left(\mathrm{ITC}^{\dagger}\right)^{2}.$
(41)

####  II-D5 Model with varying amplitude and phase

We now consider a more complex model, in which we observe repetitions of an
oscillatory signal with constant frequency ($\nu=\nu_{0}$) but varying
amplitude and phase,

$x_{n}(t)=\Omega_{n}\cos(2\pi\nu_{0}t+\phi_{n}),\qquad n=1,\dots,N.$ (42)

We assume that the $\Omega_{n}$’s are i.i.d. repetitions of
$\Omega\sim\mathcal{N}(\Omega_{0},\tau_{\Omega}^{2})$, while the $\phi_{n}$’s
are i.i.d repetitions of $\phi\sim\mathrm{VonMises}(\phi_{0},\kappa)$.
Furthermore, $\Omega_{n}$ and $\phi_{n}$ are assumed to be independent from
each other for every $n$.

According to (25), the (approximate) S-transform of each signal is given by

$T_{x_{n}}^{\dagger}(t,f)=\frac{\Omega_{n}}{2}e^{-\frac{1}{2}(2\pi)^{2}\left(1-\frac{\nu_{0}}{f}\right)^{2}}e^{i[\phi_{n}-2\pi(f-\nu_{0})t]},$
(43)

which has amplitude

$\left|T_{x_{n}}^{\dagger}(t,f)\right|=\frac{\Omega_{n}}{2}e^{-\frac{1}{2}(2\pi)^{2}\left(1-\frac{\nu_{0}}{f}\right)^{2}}$
(44)

and phase

$\arg\left[T_{x_{n}}^{\dagger}(t,f)\right]=\phi_{n}-2\pi(f-\nu_{0})t.$ (45)

In this case, both amplitude and phase vary with $n$ and are random (through
their dependence in $\Omega_{n}$ for amplitude, and $\phi_{n}$ for phase).
Also, since $\Omega_{n}$ and $\phi_{n}$ are independent for every $n$, so are
$|T_{x_{n}}^{\dagger}(t,f)|$ and $\arg[T_{x_{n}}^{\dagger}(t,f)]$.

From (1) and (44), we have

$\mathrm{avgAMP}_{x_{1:N}}^{\dagger}(t,f)=\frac{1}{2}e^{-\frac{1}{2}(2\pi)^{2}\left(1-\frac{\nu_{0}}{f}\right)^{2}}\frac{1}{N}\sum_{n=1}^{N}\Omega_{n},$
(46)

and, from (2) and (45),

$\mathrm{ITC}_{x_{1:N}}^{\dagger}(t,f)=\left|\frac{1}{N}\sum_{n=1}^{N}e^{i\phi_{n}}\right|,$
(47)

so that

$\displaystyle\mathrm{avgAMP}_{x_{1:N}}^{\dagger}(t,f)^{2}\times\mathrm{ITC}_{x_{1:N}}^{\dagger}(t,f)^{2}$
$\displaystyle=$
$\displaystyle\left[\frac{1}{2}e^{-\frac{1}{2}(2\pi)^{2}\left(1-\frac{\nu_{0}}{f}\right)^{2}}\right]^{2}\left(\frac{1}{N}\sum_{n=1}^{N}\Omega_{n}\right)^{2}\left|\frac{1}{N}\sum_{n=1}^{N}e^{i\phi_{n}}\right|^{2}.$

Besides, the time-frequency transform of the average signal is equal to the
average of the $T_{x_{n}}(t,f)$’s from Equation (43), that is,

$T_{\overline{x_{n}}}^{\dagger}(t,f)=\frac{1}{2}e^{-\frac{1}{2}(2\pi)^{2}\left(1-\frac{\nu_{0}}{f}\right)^{2}}e^{-2i\pi(f-\nu_{0})t}\,\frac{1}{N}\sum_{n=1}^{N}\Omega_{n}\,e^{i\phi_{n}}.$
(49)

The amplitude of this quantity is given by

$\left|T_{\overline{x_{n}}}^{\dagger}(t,f)\right|=\frac{1}{2}e^{-\frac{1}{2}(2\pi)^{2}\left(1-\frac{\nu_{0}}{f}\right)^{2}}\left|\frac{1}{N}\sum_{n=1}^{N}\Omega_{n}e^{i\phi_{n}}\right|.$
(50)

Incorporating this result into (3), we are led to

$\mathrm{POWavg}_{x_{1:N}}^{\dagger}(t,f)=\left[\frac{1}{2}e^{-\frac{1}{2}(2\pi)^{2}\left(1-\frac{\nu_{0}}{f}\right)^{2}}\right]^{2}\left|\frac{1}{N}\sum_{n=1}^{N}\Omega_{n}e^{i\phi_{n}}\right|^{2}.$
(51)

Comparing (LABEL:eq:osc:m2:ITC-avgAMP) and (51), we see that, in this
particular case, the relationship of (41) does not hold. Indeed, the
influences of the signals’ amplitudes (the $\Omega_{n}$’s) and phases (the
$\phi_{n}$’s) are separated in (LABEL:eq:osc:m2:ITC-avgAMP) but appear jointly
in (51).

Still, one can take advantage of the statistical independence of the
$\Omega_{n}$’s and the $\phi_{n}$’s to derive the expectations of the various
quantities. We start by calculating the expectation of amplitude from (44),
leading to

$\mathrm{E}\left[\left|T_{x_{n}}^{\dagger}(t,f)\right|\right]=\frac{\Omega_{0}}{2}e^{-\frac{1}{2}(2\pi)^{2}\left(1-\frac{\nu_{0}}{f}\right)^{2}}.$
(52)

An asymptotic approximation of $\mathrm{E}[(\mathrm{avgAMP}^{\dagger})^{2}]$
is then available from (13) and (52)

$\mathrm{E}\left[\mathrm{avgAMP}_{x_{1:N}}^{\dagger}(t,f)^{2}\right]=\left[\frac{\Omega_{0}}{2}e^{-\frac{1}{2}(2\pi)^{2}\left(1-\frac{\nu_{0}}{f}\right)^{2}}\right]^{2}+O\left(\frac{1}{N}\right).$
(53)

One can also compute the circular mean by using (8) and (45),

$\mathrm{E}\left[e^{i\arg\left[T_{x_{n}}^{\dagger}(t,f)\right]}\right]=\rho
e^{i[\phi_{0}-2\pi(f-\nu_{0})t]},$ (54)

and, from (15),

$\mathrm{E}\left[\mathrm{ITC}_{x_{1:N}}^{\dagger}(t,f)^{2}\right]=\rho^{2}+O\left(\frac{1}{N}\right).$
(55)

The product $\mathrm{E}[(\mathrm{avgAMPL}^{\dagger})^{2}]$ by
$\mathrm{E}[(\mathrm{ITC}^{\dagger})^{2}]$ can then be calculated from (53)
and (55), yielding

$\displaystyle\mathrm{E}\left[\mathrm{avgAMP}_{x_{1:N}}^{\dagger}(t,f)^{2}\right]\times\mathrm{E}\left[\mathrm{ITC}_{x_{1:N}}^{\dagger}(t,f)^{2}\right]$
$\displaystyle=$
$\displaystyle\left[\frac{\Omega_{0}\rho}{2}e^{-\frac{1}{2}(2\pi)^{2}\left(1-\frac{\nu_{0}}{f}\right)^{2}}\right]^{2}+O\left(\frac{1}{N}\right).$

Since $|T_{x}(t,f)|$ and $\arg[T_{x}^{\dagger}(t,f)]$ are independent, so are
avgAMP† and ITC†, so that, according to (5),

$\displaystyle\mathrm{E}\left[\mathrm{ITC}_{x_{1:N}}^{\dagger}(t,f)^{2}\times\mathrm{avgAMP}_{x_{1:N}}^{\dagger}(t,f)^{2}\right]$
(56) $\displaystyle=$
$\displaystyle\left[\frac{\Omega_{0}\rho}{2}e^{-\frac{1}{2}(2\pi)^{2}\left(1-\frac{\nu_{0}}{f}\right)^{2}}\right]^{2}+O\left(\frac{1}{N}\right).$

Deriving an asymptotic form for $\mathrm{E}(\mathrm{POWavg}^{\dagger})$ from
(20) and (49) is lengthier, but it can be done and leads to (see §2.4 of
Supplementary Material #1)

$\mathrm{E}[\mathrm{POWavg}_{x_{1:N}}^{\dagger}(t,f)]=\left[\frac{\Omega_{0}\rho}{2}e^{-\frac{1}{2}(2\pi)^{2}\left(1-\frac{\nu_{0}}{f}\right)^{2}}\right]^{2}+O\left(\frac{1}{N}\right).$
(57)

Equality of (56) and (57) shows that we have

$\mathrm{E}\left(\mathrm{POWavg}^{\dagger}\right)=\mathrm{E}\left[\left(\mathrm{avgAMP}^{\dagger}\right)^{2}\times\left(\mathrm{ITC}^{\dagger}\right)^{2}\right]+O\left(\frac{1}{N}\right).$
(58)

This result is weaker than in the simplest model leading to (41), in that (i)
it involves an equality in expectation, and (ii) it is only true
asymptotically.

####  II-D6 Model with varying amplitude, frequency, and phase

We finally consider repetitions of the oscillatory signal with varying
amplitude, phase, and frequency

$x_{n}(t)=\Omega_{n}\cos(2\pi\nu_{n}t+\phi_{n}),\qquad n=1,\dots,N.$ (59)

We assume that we have
$\Omega_{n}\sim\mathcal{N}(\Omega_{0},\tau_{\Omega}^{2})$,
$\nu_{n}\sim\mathcal{N}(\nu_{0},\tau_{\nu}^{2})$, and
$\phi_{n}\sim\mathrm{VonMises}(\phi_{0},\kappa)$.

The approximate time-frequency transform of each signal is given by

$T_{x_{n}}^{\dagger}(t,f)=\frac{\Omega_{n}}{2}e^{-\frac{1}{2}(2\pi)^{2}\left(1-\frac{\nu_{n}}{f}\right)^{2}}e^{i[\phi_{n}-2\pi(f-\nu_{n})t]},$
(60)

with amplitude

$\left|T_{x_{n}}^{\dagger}(t,f)\right|=\frac{\Omega_{n}}{2}e^{-\frac{1}{2}(2\pi)^{2}\left(1-\frac{\nu_{n}}{f}\right)^{2}}$
(61)

and phase

$\arg\left[T_{x_{n}}^{\dagger}(t,f)\right]=\phi_{n}-2\pi(f-\nu_{n})t.$ (62)

Note that, in this case, amplitude and phase are _not_ independent, as they
both depend on $\nu_{n}$. As in the previous model, we will need to calculate
expectations. From Equation (61), we obtain (see §2.5.2 of Supplementary
Material #1)

$\mathrm{E}\left[\left|T_{x_{n}}^{\dagger}(t,f)\right|\right]=\frac{\Omega_{0}}{2}\frac{1}{\sqrt{\left(\frac{2\pi\tau_{\nu}}{f}\right)^{2}+1}}e^{-\frac{1}{2}\frac{(2\pi)^{2}}{\left(\frac{2\pi\tau_{\nu}}{f}\right)^{2}+1}\left(1-\frac{\nu_{0}}{f}\right)^{2}}$
(63)

and, from (13),

$\displaystyle\mathrm{E}\left[\mathrm{avgAMP}_{x_{1:N}}^{\dagger}(t,f)^{2}\right]$
(65) $\displaystyle=$
$\displaystyle\left[\frac{\Omega_{0}}{2}\frac{1}{\sqrt{\left(\frac{2\pi\tau_{\nu}}{f}\right)^{2}+1}}e^{-\frac{1}{2}\frac{(2\pi)^{2}}{\left(\frac{2\pi\tau_{\nu}}{f}\right)^{2}+1}\left(1-\frac{\nu_{0}}{f}\right)^{2}}\right]^{2}$
$\displaystyle\quad+O\left(\frac{1}{N}\right).$

From (62), we obtain (see §2.5.3 of Supplementary Material #1)

$\mathrm{E}\left\\{e^{i\arg\left[T_{x_{n}}^{\dagger}(t,f)\right]}\right\\}=\rho
e^{-\frac{1}{2}(2\pi
t)^{2}\tau_{\nu}^{2}}e^{i\left[\phi_{0}-2\pi(f-\nu_{0})t\right]},$ (66)

so that

$\left|\mathrm{E}\left\\{e^{i\arg\left[T_{x_{n}}^{\dagger}(t,f)\right]}\right\\}\right|=\rho
e^{-\frac{1}{2}(2\pi t)^{2}\tau_{\nu}^{2}},$ (67)

and, from (15),

$\mathrm{E}\left[\mathrm{ITC}_{x_{1:N}}^{\dagger}(t,f)^{2}\right]=\rho^{2}\left[e^{-\frac{1}{2}(2\pi
t)^{2}\tau_{\nu}^{2}}\right]^{2}+O\left(\frac{1}{N}\right).$ (68)

To calculate POWavg†, we need the following quantity (see §2.5.4 of
Supplementary Material #1):

$\displaystyle\mathrm{E}\left[T_{x_{n}}^{\dagger}(t,f)\right]$
$\displaystyle=$
$\displaystyle\frac{1}{2}\mathrm{E}(\Omega_{n})\mathrm{E}\left(e^{i\phi_{n}}\right)$
$\displaystyle\quad\times\mathrm{E}\left[e^{-\frac{1}{2}(2\pi)^{2}\left(1-\frac{\nu_{n}}{f}\right)^{2}-2i\pi(f-\nu_{n})t}\right]$
$\displaystyle=$
$\displaystyle\frac{\Omega_{0}\rho}{2}\frac{1}{\sqrt{\left(\frac{2\pi\tau_{\nu}}{f}\right)^{2}+1}}e^{-2i\pi(f-\nu_{0})t}$
$\displaystyle\qquad\times
e^{-\frac{1}{2}\frac{(2\pi)^{2}}{\left(\frac{2\pi\tau_{\nu}}{f}\right)^{2}+1}\left[\left(1-\frac{\nu_{0}}{f}\right)^{2}+\tau_{\nu}^{2}t^{2}\right]}.$

Finally, from (20),

$\displaystyle\mathrm{E}\left[\mathrm{POWavg}_{x_{1:N}}^{\dagger}(t,f)\right]$
(69) $\displaystyle=$
$\displaystyle\left[\frac{\Omega_{0}}{2}\frac{1}{\sqrt{\left(\frac{2\pi\tau_{\nu}}{f}\right)^{2}+1}}e^{-\frac{1}{2}\frac{(2\pi)^{2}}{\left(\frac{2\pi\tau_{\nu}}{f}\right)^{2}+1}\left(1-\frac{\nu_{0}}{f}\right)^{2}}\right]^{2}$
$\displaystyle\quad\times\rho^{2}\left[e^{-\frac{1}{2}\frac{(2\pi)^{2}}{\left(\frac{2\pi\tau_{\nu}}{f}\right)^{2}+1}\tau_{\nu}^{2}t^{2}}\right]^{2}+O\left(\frac{1}{N}\right).$

From (65), (68), and (69), we are now in position to calculate

$\displaystyle\mathrm{E}\left[\mathrm{POWavg}_{x_{1:N}}^{\dagger}(t,f)\right]$
(70)
$\displaystyle\qquad-\mathrm{E}\left[\mathrm{ITC}_{x_{1:N}}^{\dagger}(t,f)^{2}\times\mathrm{avgAMP}_{x_{1:N}}^{\dagger}(t,f)^{2}\right]$
$\displaystyle=$
$\displaystyle\left[\frac{\Omega_{0}\rho}{2}\frac{1}{\sqrt{\left(\frac{2\pi\tau_{\nu}}{f}\right)^{2}+1}}e^{-\frac{1}{2}\frac{(2\pi)^{2}}{\left(\frac{2\pi\tau_{\nu}}{f}\right)^{2}+1}\left(1-\frac{\nu_{0}}{f}\right)^{2}}\right]^{2}$
$\displaystyle\
\times\left\\{\left[e^{-\frac{1}{2}\frac{(2\pi)^{2}}{\left(\frac{2\pi\tau_{\nu}}{f}\right)^{2}+1}\tau_{\nu}^{2}t^{2}}\right]^{2}-1\right\\}+O\left(\frac{1}{N}\right).$

This difference, which is negative, does not tend to 0 as $N\to\infty$.
However, it becomes increasingly smaller when $\tau_{\nu}$ becomes smaller and
vanishes for $\tau_{\nu}=0$, which corresponds to a situation where the signal
frequency is fixed, leading to independent amplitude and phase for the time-
frequency transform.

###  II-E General expression

In this section, we confirm the results of the previous section and show in a
more general setting that

$\mathrm{E}\left[\mathrm{POWavg}-\mathrm{ITC}^{2}\times\mathrm{avgAMP}^{2}\right]=O\left(\frac{1}{N}\right),$
(71)

provided that a certain condition is met, which includes the case where the
covariance between the amplitude $|T_{x}(t,f)|$ and phase $\theta_{x}(t,f)$ of
$T_{x}(t,f)$ is equal to 0. Since independence entails a covariance of 0,
independence is a sufficient condition.

####  II-E1 Case of independent amplitude and phase

Given the asymptotic expansions of avgAMP, ITC, and POWavg, we can expect the
following behaviors for large $N$. On the one hand, from (13) and (15), we
obtain

$\displaystyle\mathrm{avgAMP}_{x_{1:N}}(t,f)^{2}\times\mathrm{ITC}_{x_{1:N}}(t,f)^{2}$
(72) $\displaystyle=$
$\displaystyle\mathrm{E}\left[\left|T_{x}(t,f)\right|\right]^{2}\,\left|\mathrm{E}\left[e^{i\theta_{x}(t,f)}\right]\right|^{2}+O\left(\frac{1}{N}\right)$

On the other hand, starting from (20) and, remembering the expression of
$T_{x}(t,f)$ in terms of $|T_{x}(t,f)|$ and $e^{i\theta_{x}(t,f)}$, (11), we
have

$\displaystyle\mathrm{E}\left[\mathrm{POWavg}_{x_{1:N}}(t,f)\right]$ (73)
$\displaystyle=$
$\displaystyle\left|\mathrm{E}\left[|T_{x}(t,f)|e^{i\theta_{x}(t,f)}\right]\right|^{2}+O\left(\frac{1}{N}\right).$

Now, if we assume independence of $|T_{x}(t,f)|$ and $e^{i\theta_{x}(t,f)}$,
we can apply (5), leading to

$\displaystyle\mathrm{E}\left[\mathrm{POWavg}_{x_{1:N}}(t,f)\right]$ (74)
$\displaystyle=$
$\displaystyle\mathrm{E}\left[\left|T_{x}(t,f)\right|\right]^{2}\,\left|\mathrm{E}\left[e^{i\theta_{x}(t,f)}\right]\right|^{2}+O\left(\frac{1}{N}\right).$

Since (72) and (74) are identical, we showed that

$\mathrm{E}(\mathrm{POWavg})=\mathrm{E}(\mathrm{avgAMP}^{2}\times\mathrm{ITC}^{2})+O\left(\frac{1}{N}\right)$
(75)

in the case of independent amplitude and phase of the time-frequency
transform.

####  II-E2 General case

We now proceed to a more general proof, showing in the process that having a
covariance of 0 is a sufficient condition. To this end, we first calculate
avgAMP2 $\times$ ITC2 from (1) and (2), expanding the product by brute force,
and then calculate its expectation term by term. We then calculate POWavg and
its expectation using (20), (11), and (6). Calculating the difference between
the two terms leads to (see §3 of Supplementary Material #1):

$\displaystyle\mathrm{E}(\mathrm{POWavg}-\mathrm{ITC}^{2}\times\mathrm{avgAMP}^{2})$
(76) $\displaystyle=$
$\displaystyle\left|\mathrm{E}\left[e^{i\theta_{x}(t,f)}\right]\mathrm{E}\left[\left|T_{x}(t,f)\right|\right]+\mathrm{Cov}\left[e^{i\theta_{x}(t,f)},\left|T_{x}(t,f)\right|\right]\right|^{2}$
$\displaystyle\quad-\left|\mathrm{E}\left[e^{i\theta_{x}(t,f)}\right]\right|^{2}\mathrm{E}\left[\left|T_{x}(t,f)\right|\right]^{2}+O\left(\frac{1}{N}\right).$

If $\mathrm{Cov}\left[e^{i\theta_{x}(t,f)},\left|T_{x}(t,f)\right|\right]=0$,
then the right-hand side of the equation is equal to $O(1/N)$. A sufficient
condition for this covariance to be equal to 0 is for $e^{i\theta_{x}(t,f)}$
and $\left|T_{x}(t,f)\right|$ to be statistically independent. But this is not
a necessary condition, as we can have non-independent (non-Gaussian) variables
whose covariance is equal to 0. If the covariance is different from 0, it is
still possible for the right-hand side of (76) to be $O(1/N)$, although the
corresponding solutions cannot be associated with simple cases.

##  III Simulation study

In this section, we investigated the connection between avgAMP, ITC, and
POWavg by generating synthetic data in the context of high-frequency
oscillations (HFOs). Application of time-frequency analysis to brain
recordings from event-related protocols involving sensory stimulation have
indicated that phase locking of ongoing EEG activity strongly contributes to
HFOs (range 400–800 Hz) which are superimposed on the somatosensory evoked
potential (SEP) [28, 29, 30, 25]. The corresponding detailed temporal
components of short durations and small amplitudes were generally identified
by measuring latencies of the averaged evoked potential in time domain [30].
In addition, ITC analysis was found to mimic the time-frequency power
distribution of the HFOs, while the average of the single-sweep time-frequency
power transforms did not show any significant increment of power [25].

###  III-A Data generation

We generated signals on a time window of $[-100,100]$ ms at a sampling rate of
$f_{s}=2$ kHz (corresponding to a recording every $\delta t=0.5$ ms). Signals
corresponding to the different trials were generated independently. For each
trial $n$, we simulated an induced response in the $[20,30]$ ms time window,
and ongoing activity the rest of the time. Both the ongoing and the induced
activities were generated using (21) with the same amplitude $\Omega_{n}$ and
frequency $\nu_{n}$, but with different phase: $\phi_{n}^{(o)}$ for the
ongoing activity, and $\phi_{n}^{(i)}$ for the induced activity. The exact
values of parameters were sampled according to specific distributions (see
Table I). To limit the effect of spectral leakage on ITC, we added Gaussian
white noise with a standard deviation of 0.01 (corresponding to 1% of the
expected signal amplitude). The signals were analyzed using the $S$-transform,
as well as two wavelet transforms (Morse wavelet and analytic Morlet wavelet).
Only results of the S-transform applied to the signals with $\kappa^{(i)}=10$
are reported here; see §1 of Supplementary Material #2 for all results.

TABLE I: Simulation study. Values of parameters for data generation. Parameter | Distribution | Parameters
---|---|---
$\Omega_{n}$ | $\mathcal{N}(\Omega_{0},\tau_{\Omega}^{2})$ | $\Omega_{0}=1$ | $\tau_{\Omega}=0.1$
$\nu_{n}$ | $\mathcal{N}(\nu_{0},\tau_{\nu}^{2})$ | $\nu_{0}=500$ | $\tau_{\nu}\in\\{0,5\\}$
$\phi_{n}^{(o)}$ | $\mathrm{vonMises}[\phi_{0},\kappa^{(o)}]$ | $\phi_{0}=0$ | $\kappa^{(o)}=0$
$\phi_{n}^{(i)}$ | $\mathrm{vonMises}[\phi_{0},\kappa^{(i)}]$ | $\phi_{0}=0$ | $\kappa^{(i)}\in\\{1,10\\}$

###  III-B Assessment analysis

Based on the simulated data, we computed avgAMP2, ITC2, POWavg, as well as the
difference between POWavg and $\mathrm{avgAMP}^{2}\times\mathrm{ITC}^{2}$. To
assess the validity of the assumption of a covariance of 0 between the
amplitude and phase of the time-frequency transform, we also computed the
module of the sample covariance between $e^{i\theta_{x}(t,f)}$ and
$\left|T_{x}(t,f)\right|$.

###  III-C Results

Results are summarized in Fig. 1 and in §1 of Supplementary Material #2.
avgAMP2 evidenced a (mostly) constant signal power around 500 Hz. For an
oscillatory signal of the form given by (21), the maximum of the power for the
S-transform is reached around $f=\nu$, with value close to
$\Omega_{0}^{2}/4=0.25$. The values of avgAMP2 were in agreement with this
result.

ITC behaved as expected. It was able to detect the increase in phase synchrony
between 20 and 30 ms, with very low values (close to 0) outside the 20–30 ms
window, and large values (close to 1) within this range. As predicted by our
calculations, high values in the 20–30 ms window were not limited to
frequencies close to 500 Hz.

POWavg was found to well fit the values predicted by
$\mathrm{avgAMP}^{2}\times\mathrm{ITC}^{2}$. As a consequence, it shared the
localization properties of both avgAMP and ITC. Like avgAMP, it was larger for
frequencies close to 500 Hz; and, like ITC, it was larger in the 20–30 ms
window. Both features put together made it so that it was able to precisely
localize both the frequency of the signal and the time of increased phase
synchronization.

Looking more specifically at the error made by the relationship of (4), we
observed that it had a maximum value of about 10% of the maximum value of
POWavg. Its distribution in terms of time and frequency was very specific, as
it was lower than 0.005 (corresponding to 2% of the maximum value of POWavg)
for most values of $t$ and $f$, with a relative exception around 20 and 30 ms.
These two times corresponded to discontinuities in all measures; see below.
Also, largest error values were observed for small values of POWavg: for
values of POWavg larger than 0.1 (corresponding to at least 40% of the
maximum), the error was smaller than 0.0025 (corresponding to 1% of the
maximum value of POWavg). The absolute value of the covariance between
$e^{i\theta_{x}(t,f)}$ and $|T_{x}(t,f)|$ was found to follow a distribution
in terms of $t$ and $f$ that was similar to that of the error, with a maximal
value around 0.08.

These results bring evidence in favor of the validity of the approximations
made in the theoretical calculations: (i) approximating measures by their
expectations; (ii) using an approximation of the S-transform, (25); (iii)
using asymptotic approximations in $O(1/N)$.

We also noticed a discontinuity in avgAMP2 around 20 ms and 30 ms, as well as
artifacts of ITC2 at lower frequencies. These are probably artifacts of the
time-frequency transform are due to the sudden change in phase at 20 ms (from
$\phi_{n}^{(o)}$ to $\phi_{n}^{(i)}$) and at 30 ms (back from $\phi_{n}^{(i)}$
to $\phi_{n}^{(o)}$) and limited amount of noise in the data, i.e, artifacts
due to the simple model we used. In particular, the artifacts observed for
ITC2 were likely a consequence of spectral leakage and were reduced by
increasing the standard deviation of the noise.

POWavg | avgAMP2
---|---
|
ITC2 |  $\mathrm{avgAMP}^{2}\times\mathrm{ITC}^{2}$ as a function of POWavg
|
Figure 1: Simulation study. Results from time-frequency transform with the
S-transform for signals with $\kappa^{(i)}=10$ and $\tau_{\nu}=0$. Around the
expected signal frequency, $f=\nu_{0}$, the expected value of the power should
be close to $\Omega_{0}^{2}/4=0.25$.

##  IV Application to experimental data

We rely on [25], who previously investigated the connection between the
behaviors of POWavg, avgAMP, and ITC in the context of high-frequency
oscillations (HFOs) in early cortical somatosensory evoked potentials.

###  IV-A Data

Somatosensory evoked potentials following median nerve stimulations were
recorded in a healthy subject. Brain responses were acquired using
multichannel EEG with a sampling frequency of 3 kHz. Electrical median nerve
stimulation of 1 ms duration was applied to median nerve at the wrist level to
elicit a burst of high-frequency oscillations (HFOs) in the 400–800 frequency
range superimposed onto the cortical N20 potential [28, 29, 30]. The stimulus
was applied 300 times, with a 500-ms inter-trial interval. Following previous
recommendations [30], we studied the fronto-central channels (CP3–Fz). Data
acquisition was performed at the Center for Neuroimaging Research (CENIR) of
the Brain and Spine Institute (ICM, Paris, France). The experimental protocol
was approved by the CNRS Ethics Committee and by the national ethical
authorities (CPP Île-de-France, Paris 6 – Pitié-Salpêtrière and ANSM).

###  IV-B Analysis

The data was analyzed in the window 10–100 ms after stimulation. We started at
10 ms after stimulus to avoid consequences of the artifacts generated by the
stimulation. For the time-frequency transform, we used the S-transform. Akin
to what was done for the simulation study, we computed the three measures
POWavg, avgAMP2, and ITC2; the error
$\left|\mathrm{POWavg}-\mathrm{avgAMP}^{2}\times\mathrm{ITC}^{2}\right|$; and
the absolute value of the covariance between $e^{i\theta_{x}(t,f)}$ and
$|T_{x}(t,f)|$.

###  IV-C Results

Results are summarized in Fig. 2 and in §2 of Supplementary Material #2.
avgAMP2 did not exhibit a clear frequency band with larger values. POWavg and
ITC2 both exhibited an increase around 20 ms in the range of 800–850 Hz. We
found evidence for a linear relationship between POWavg and
$\mathrm{avgAMP}^{2}\times\mathrm{ITC}^{2}$, even though the proportionality
factor was found to be rather in the range 0.6–0.7 than close the predicted
value of 1. The error made by the approximation of (4) was observed to be
roughly proportional to POWavg, corresponding to relative error of the order
of 30%, consequently exhibiting larger values where POWavg was larger. The
covariance between $e^{i\theta_{x}(t,f)}$ and $\left|T_{x}(t,f)\right|$ was
also observed to be larger in the same time-frequency region.

POWavg | avgAMP2
---|---
|
ITC2 |  $\mathrm{avgAMP}^{2}\times\mathrm{ITC}^{2}$ as a function of POWavg
|
Figure 2: Real data. Results from time-frequency transform with the
S-transform: POWavg, avgAMP2, ITC2,
$\left|\mathrm{POWavg}-\mathrm{avgAMP}^{2}\times\mathrm{ITC}^{2}\right|$, and
relation between POWavg and $\mathrm{avgAMP}^{2}\times\mathrm{ITC}^{2}$.
Horizontal solid lines were added at 800 and 850 Hz, and vertical dashed lines
were added at 19 and 23 ms.

##  V Discussion

In neuroscience, time frequency analysis is a key method for the investigation
of brain rhythms from brain recordings, in particular when considering event-
related protocols. In this framework, three measures have been proposed, which
we defined in (1)–(3) and coined avgAMP, ITC, and POWavg. While these three
measures are sensitive to different features of brain responses, similarities
between ITC and POWavg have been observed. In the present manuscript, we
further investigated the connection between ITC and POWavg using theoretical
calculations and a simulation study. With the theoretical calculations, we
derived several variants of a relationship that roughly read POWavg $=$
avgAMP2 $\times$ ITC2 and whose exact forms depended on the specifics of the
underlying models and assumptions. The simulation study nicely confirmed the
global form of the relationship, thereby bringing evidence in favor of the
validity of the approximations made in the theoretical calculations.

The relationship between POWavg, ITC and avgAMP might seem obvious and highly
anticipated, as (i) these three measures are computed from the same brain
signals, and (ii) empirical results have showed strong similarities between
POWavg and ITC2, e.g., [25]. Yet, we argue that the results developed in the
present manuscript are novel in at least two ways. First, we started with the
observed common similarity of results produced using POWavg and ITC and
investigated its theoretical underpinning. We then proved that the
relationship between these _two_ measures actually also involves the _third_
measure, avgAMP. We believe that the involvement of avgAMP was not
anticipated. The second novelty is that the relationship takes a very simple
form, one of proportionality. We do not think that this result was largely
expected either.

We advocate that both results are an unlikely coincidence that is the very
consequence of the choice of POWavg, ITC and avgAMP as measures of interest.
For instance, changing ITC to another measure of angular dispersion, e.g., the
circular mean deviation [27, §2.3.4]

$d_{0}(\tilde{\phi})=\frac{1}{N}\sum_{n=1}^{N}\left[\pi-\left|\pi-\left|\theta_{x_{n}}(t,f)-\tilde{\phi}\right|\right|\right],$
(77)

where $\tilde{\phi}$ is the sample median of the $\theta_{x_{n}}(t,f)$’s,
would render the relationship invalid. Similarly, changing avgAMP (a measure
of amplitude) to a measure of power,

$\mathrm{avgPOW}=\frac{1}{N}\sum_{n=1}^{N}\left|T_{x_{n}}(t,f)\right|^{2},$
(78)

would make for a more complex relationship, as we have

$\mathrm{E}(\mathrm{avgPOW})=\mathrm{E}\left[\left|T_{x}(t,f)\right|\right]^{2}+\mathrm{Var}\left[\left|T_{x}(t,f)\right|\right],$
(79)

to be compared with (12). In the expression of $\mathrm{E}(\mathrm{avgPOW})$,
the variance term does not vanish with increasing $N$, while it does for
$\mathrm{E}(\mathrm{avgAMP}^{2})$—see (13). As a consequence, we have
$\mathrm{POWavg}\neq\mathrm{avgPOW}\times\mathrm{ITC}^{2}$. So it is a subtle
interplay between the three measures selected that make the existence of the
relationship possible.

As a consequence, the observed similarity between the _two_ measures POWavg
and ITC actually emphasizes the fact that the _three_ measures POWavg, ITC and
avgAMP are “naturally” connected in ways that remain to be fully understood,
and provide three views of the same core time-frequency feature of the signal.

With the calculations, we first investigated the particular case of
oscillatory signals (Section II-D). Recordings of brain activity are
hypothesized to be interpretable in terms of brain rhythms, i.e., of signals
of specific frequency ranges. This is a major reason why time-frequency
transform is applied to this kind of signals. It is also evidence that
oscillatory signals are natural candidates for realist models of brain
activity.

To perform the time-frequency transform, we applied the S-transform. There are
two reason for this choice. First, it is a method that is commonly applied in
the analysis of brain recordings. Furthermore, it had the advantage of
rendering our calculations tractable. Its main drawbacks are its (relative)
computational burden and the fact that it is not analytic, i.e., its Fourier
transform has nonzero power in negative frequencies. As a consequence, it
splits the power of the signal on positive and negative frequency components.
Our calculations with the S-transform were made possible by replacing the true
time-frequency transform $T_{x}(t,f)$ of (LABEL:eq:osc:gen:tf) with the
approximation $T_{x}^{\dagger}(t,f)$ of (25). We showed that such an
approximation was valid over a wide range of frequencies.

Using theoretical calculations, we also tackled the general case (Section
II-E). There, we showed that the relationship was of the form of (71),
provided that the covariance between the amplitude and phase of the time-
frequency transform was equal to 0. To reach this conclusion, we relied on the
assumption that avgAMP, ITC, and POWavg, which are averages computed over many
repetitions $N$ of the same experiment (typically, $N$ is equal to one to
several hundreds), could be approximated by expectations. This is equivalent
to neglecting the variability of the measures, i.e., equivalent to assuming
low variance.

We finally used a simulation study (Section III) to generate synthetic HFOs as
observed in event-related protocols involving sensory stimulation. The
generative model allowed for a certain level of variability that could not be
taken into account in the calculations. Note that we did not simulate the main
SEP (N20), as this was not necessary: Since the time-frequency transform is
linear, the time-frequency transform of the SEP would only be superimposed to
the time-frequency transform of th HFOs. Since frequencies do not overlap, the
SEP is expected to have very limited influence on our analysis.

The synthetic data were then analyzed using the S-transform. Despite the
above-mentioned variability, the global relationship of (4) still seemed to
hold relatively well. In particular, the results showed that both ITC and
POWavg were sensitive to the simulated phase resetting phenomenon.

Unexpectedly, the results also hinted that ITC might not be able to determine
the frequency at which the phase resetting occurs, while POWavg seemed to be
able to. This is in agreement with our theoretical calculations—compare, on
the one hand, (37), (55), and (68), which show that these expressions of ITC
(or ITC2, or its expectation) are not a function of $f$, and, on the other
hand, (40), (51), and (69), which show that POWavg (or its expectation) are.
However, to the best of our knowledge, this is contrary to empirical evidence
on real data, where ITC is usually localized in frequency, making possible to
narrow the range where the induced oscillations take place [21, 25]. We aim to
find an explanation to this discrepancy.

Besides the S-transform, another common ways to perform time-frequency
transform is to use the wavelet transform. Wavelets transforms are fast and,
when using analytic wavelets, only decompose a signal in the range of positive
frequencies. For instance, the expected power at $f=\nu_{0}$ of our
oscillatory signal is 1 instead of 0.25. In our simulation study, we found
that the practical difference between S-transform and wavelet transform was
very limited; see, e.g., §3 of Supplementary Material #2 for analysis of the
synthetic data using wavelet transforms (Morse and analytic Morlet wavelets).

Note that the simulation study used a crude model of brain response. As a
consequence, we observed on the simulated data some results that are not
typical of time-frequency diagrams of real brain recordings. Some of these
atypical results corresponded to theoretical predictions, while others did
not. For instance, as mentioned above, the fact that ITC on simulated HFOs was
not localized in a limited frequency band around the true oscillation
frequency was to be expected from the theoretical calculations but not from
the literature. Other features of the simulated time-frequency diagrams did
not correspond to what is typically observed on time-frequency diagrams of
real brain recordings. This includes a discontinuity in avgAMP2 at 20 ms and
30 ms, as well as artifacts of ITC2 at lower frequencies. These are probably
artifacts of the time-frequency transform might be due to the sudden change in
phase at 20 ms (from $\phi_{n}^{(o)}$ to $\phi_{n}^{(i)}$) and at 30 ms (back
from $\phi_{n}^{(i)}$ to $\phi_{n}^{(o)}$) and limited amount of noise in the
data. As a consequence, the possibility that we had in our simulations to
detect phase resetting in avgAMP is merely an artifact resulting from an
oversimplified simulation model rather than an empirically confirmed result.

While we focused on HFOs for the simulation study, nothing in the theoretical
developments depends on the frequency of the underlying signal. As a
consequence, we expect the general relationship to hold for different brain
frequency ranges as well. As a confirmation, we ran a new simulation with
$\nu_{0}=40$ Hz. As expected, the results were very similar to the ones
obtained for the simulated HFOs (for details, see §4 of Supplementary Material
#2). This is particularly relevant in the context of auditory evoked steady-
state responses (ASSR) which are increasingly being used to probe the ability
of auditory circuits to produce synchronous activity at 40 Hz, in response to
repetitive external stimulation [31]. Perturbation in both power or phase-
locking has been also observed in various neuropsychiatric disorders [32, 33].
Because the ASSR is composed of complex time-dependent frequency responses
[34], the present study suggests that the combination of different methods
(avgAMP, ITC, and POWavg) can be used to avoid a misleading or wrong
interpretation of the neuronal modulations, especially as markers of brain
dysfunctions in psychiatric disorders.

Also, the present work focused on the case of a signal composed of a single
frequency. We expect the presence of multiple frequencies to have a limited
impact on the calculation of avgAMP beyond mere additivity, but to strongly
interfere with the calculation of ITC, potentially reducing the frequency
range over which a phase resetting phenomenon is associated with a larger
value of ITC. Further simulation studies could help assess the validity of
such an assumption. It would also be interesting to investigate the particular
case where there exists dependencies between frequencies, e.g., in the form of
cross-frequency coupling.

From the theoretical calculations, we found that the relationship (4) between
POWavg, avgAMP, and ITC requires that amplitude and phase have 0 covariance.
This requirement is slightly more general than independence, as there exist
variables that have 0 covariance but are not independent. On the simulated
data, we quantified the level of validity of this assumption and found a
rather good fit, at the exception of around 20 and 30 ms, where we observed
spectral leakage. On the experimental data, we found the presence of a linear
relationship between quantities, albeit with a proportionality factor that was
lower than the expected value of 1. Also, the values of both
$|\mathrm{POWavg}-\mathrm{avgAMP}^{2}\times\mathrm{ITC}^{2}|$ and
$|\mathrm{Cov}[e^{i\theta_{x}(t,f)},|T_{x}(t,f)|]|$ were roughly linearly
related to the values of POWavg. We hypothesize that the departure from
theoretical expectations might be related to the presence of noise.

Regarding noise, it has to be mentioned that the present work is based on the
assumption that there is no noise and that the totality of the signal measured
is relevant. This is usually not the case, as brain recordings are subject to
noise. The influence of this noise on time-frequency transform needs to be
further investigated. In particular, we need to investigate whether it could
be the source of the “dampening” effect on the linear relationship observed in
the real data. Note that, by contrast, [35, §19.3] expects noise to increase
ITC. We hope to be able to clarify this point in the near future. In this
endeavor, the relationship of (4) might prove to be a valuable tool.

We can now conclude regarding the implication of the present development to
neuroscientific practice. Time-frequency analyses have been widely used in
electrophysiological research for the characterization of oscillatory
responses to stimuli in a large range of brain processes, from simple sensory
processing to higher order cognition [22, 17, 15]. Multiple studies have
reported evidence that evoked oscillations are not exclusively composed of
either a stimulus phasic-related component or some stimulus-induced phase
resetting, but usually include a combination of both in variable proportions
[36, 24]. As a consequence, the respective contributions of both types of
components have to be clarified before attempting to draw physiological
conclusions. Unfortunately, distinguishing between the two phenomena remains
quite a challenge, in particular when multiple oscillations are present in
neighboring frequency bands [37, 38]. One potentially fruitful approach for
future research could be to develop measures quantifying to what extent
changes in ITC and in avgAMP are correlated in the time and frequency domains.
The exact bearing of the relationship will become increasingly clearer as real
data are analyzed with this result in mind.

##  VI Conclusion

In this manuscript, we showed that avgAMP, ITC, and POWavg were related
through the approximate, simple relationship given by (4). This result, based
on theoretical calculations and a simulation study, confirms previous
empirical evidence and provides a novel perspective to investigate evoked
brain rhythms. In particular, looking at ITC and avgAMP as two distinct
measures quantifying different aspects of rhythmic brain activity, and taking
into consideration their level of correlation may help disentangle
contributions from additive evoked activity and phase-reset oscillatory
activity. In this context, the presented relationship between POWavg, ITC, and
avgAMP may provide a significant refinement to the neuroscientific toolbox for
studying evoked oscillations.

As a final note, time-frequency analysis is a method that is commonly used in
many fields of scientific research, including fluids mechanics, engineering,
other fields of biomedical imaging and medicine, and geophysics [6]. We hope
that the present development will find applications in these fields as well.

## Acknowledgment

The authors would like to thank the Center for Neuroimaging Research (CENIR)
of the Brain and Spine Institute (ICM, Paris, France) for the acquisition of
the experimental data, and Véronique Marchand-Pauvert for providing them with
the data.

## References

  * [1] K. Gröchenig, _Foundation of Time-Frequency Analysis_. Birkhäuser, Boston, 2001.
  * [2] S. Mallat, _A Wavelet Tour of Signal Processing_ , 2nd ed., ser. Wavelet Analysis & Its Applications. Academic Press, Amsterdam, 1999.
  * [3] R. G. Stockwell et al, “Localization of the complex spectrum: the $S$ transform,” _IEEE Transactions on Signal Processing_ , vol. 44, no. 4, pp. 998–1001, 1996.
  * [4] L. Cohen, _Time-Frequency Analysis_. Prentice-Hall, New York, 1995.
  * [5] P. Flandrin, _Time-Frequency/Time-Scale Analysis_ , ser. Wavelet Analysis and Its Applications. Academic Press, San Diego, 1999, vol. 10.
  * [6] P. S. Addison, _The Illustrated Wavelet Transform Handbook. Introductory Theory and Applications in Science, Engineering, Medicine and Finance_. Institute of Physics Publishing, Bristol, 2002\.
  * [7] E. Başar et al, “Are cognitive processes manifested in event-related gamma, alpha, theta and delta oscillations in the EEG?” _Neuroscience Letters_ , vol. 259, pp. 165–168, 1999.
  * [8] W. Singer, “Neuronal synchrony: a versatile code for the definition of relations?” _Neuron_ , vol. 24, pp. 49–65, 1999.
  * [9] L. Glass, “Synchronization and rhythmic processes in physiology,” _Nature_ , vol. 410, pp. 277–284, 2001.
  * [10] G. Buzsáki and A. Draguhn, “Neuronal oscillations in cortical networks,” _Science_ , vol. 304, pp. 1926–1929, 2004.
  * [11] H. Laufs, “Brain rhythms,” in _EEG-fMRI. Physiological Basis, Technique, and Applications_ , C. Mulert and L. Lemieux, Eds. Springer, Berlin, 2010, pp. 262–277.
  * [12] P. Fries, “Rhythms for cognition: Communication through coherence,” _Neuron_ , vol. 88, no. 1, pp. 220–235, 2015.
  * [13] W. Penny et al, “Event-related brain dynamics,” _Trends in Neurosciences_ , vol. 25, no. 8, pp. 387–389, 2002.
  * [14] S. Makeig et al, “Dynamic brain sources of visual evoked responses,” _Science_ , vol. 295, pp. 690–694, 2002.
  * [15] S. Makeig et al, “Mining event-related brain dynamics,” _Trends in Cognitive Sciences_ , vol. 8, no. 5, pp. 204–210, 2004\.
  * [16] O. David et al, “Mechanisms of evoked and induced responses in MEG/EEG,” _NeuroImage_ , vol. 31, no. 4, pp. 1580–1591, 2006.
  * [17] F. J. Varela et al, “The brainweb: phase synchronization and large-scale intergration,” _Nature Reviews. Neuroscience_ , vol. 2, no. 4, pp. 229–239, 2001.
  * [18] A. Delorme and S. Makeig, “EEGLAB: an open source toolbox for analysis of single-trial EEG dynamics including independent component analysis,” _Journal of Neuroscience Methods_ , vol. 134, no. 1, pp. 9–21, 2004.
  * [19] M. Le Van Quyen and A. Bragin, “Analysis of dynamic brain oscillations: methodological advances,” _Trends in Neurosciences_ , vol. 30, no. 7, pp. 365–373, 2007.
  * [20] R. Hari and A. Puce, _MEG-EEG Primer_. Oxford University Press, Oxford, 2017.
  * [21] C. Tallon-Baudry et al, “Stimulus specificity of phase-locked and non-phase-locked 40 Hz visual responses in human,” _The Journal of Neuroscience_ , vol. 16, no. 13, pp. 4240–4249, 1996\.
  * [22] C. Tallon-Baudry and O. Bertrand, “Oscillatory gamma activity in human and its role in object representation,” _Trends in Cognitive Sciences_ , vol. 3, no. 4, pp. 151–162, 1999.
  * [23] L. Fuentemilla et al, “Modulation of spectral power and of phase resetting of EEG contributes differentially to the generation of auditory event-related potentials,” _NeuroImage_ , vol. 30, no. 3, pp. 909–916, 2006.
  * [24] S. Hanslmayr et al, “Alpha phase reset contributes to the generation of ERPs,” _Cerebral Cortex_ , vol. 17, no. 1, pp. 1–8, 2007.
  * [25] M. Valencia et al, “High frequency oscillations in the somatosensory evoked potentials (SSEP’s) are mainly due to phase-resetting phenomena,” _Journal of Neuroscience Methods_ , vol. 154, pp. 142–148, 2006.
  * [26] T. W. Anderson, _An Introduction to Multivariate Statistical Analysis_ , ser. Wiley Publications in Statistics. John Wiley and Sons, New York, 1958.
  * [27] K. V. Mardia and P. E. Jupp, _Directional Statistics_ , ser. Wiley Series in Probability and Statistics. Wiley, Chichester, 2000.
  * [28] G. Curio et al, “High-frequency activity (600 Hz) evoked in the human primary somatosensory cortex: a survey of electric and magnetic recordings,” in _Oscillatory Event-Related Brain Dynamics_ , ser. NATO ASI Series (Series A: Life Sciences), C. Pantev, T. Elbert, and B. Lütkenhöner, Eds. Springer, Boston, MA, 1994, vol. 271, pp. 205–218.
  * [29] G. Curio et al, “Somatotopic source arrangement of 600 Hz oscillatory magnetic fields at the human primary somatosensory hand cortex,” _Neuroscience Letters_ , vol. 234, no. 2–3, pp. 131–134, 1997.
  * [30] I. Ozaki et al, “High frequency oscillations in early cortical somatosensory evoked potentials,” _Electroencephalography and Clinical Neurophysiology_ , vol. 108, no. 6, pp. 536–542, 1998.
  * [31] R. Galambos et al, “A 40 Hz auditory potential recorded from the human scalp,” _Proceedings of the National Academy of Sciences of the U.S.A._ , vol. 78, no. 4, pp. 2643–2647, 1981.
  * [32] J. S. Kwon et al, “Gamma frequency–range abnormalities to auditory stimulation in schizophrenia,” _Archives of General Psychiatry_ , vol. 56, no. 11, pp. 1001–1005, 1999.
  * [33] T. W. Wilson et al, “Children and adolescents with autism exhibit reduced MEG steady-state gamma responses,” _Biological Psychiatry_ , vol. 62, no. 3, pp. 192–197, 2007.
  * [34] M. Tada et al, “Global and parallel cortical processing based on auditory gamma oscillatory responses in humans,” _Cerebral Cortex_ , vol. 31, no. 10, pp. 4518–4532, 2021.
  * [35] M. X. Cohen, _Analyzing Neural Time Series Data: Theory and Practice_ , ser. Issues in Clinical and Cognitive Neuropsychology. The MIT Press, Cambridge, MA, 2014.
  * [36] A. S. Shah et al, “Neural dynamics and the fundamental mechanisms of event-related brain potentials,” _Cerebral Cortex_ , vol. 14, no. 5, pp. 476–483, 2004.
  * [37] N. Yeung et al, “Detection of synchronized oscillations in the electroencephalogram: an evaluation of methods,” _Psychophysiology_ , vol. 41, no. 6, pp. 822–832, 2004.
  * [38] R. M. van Diepen and A. Mazaheri, “The caveats of observing inter-trial phase-coherence in cognitive neuroscience,” _Scientific Reports_ , vol. 8, p. 2990, 2018.