# Uncertainty and information in physiological signals: explicit physical trade-off with log-normal wavelets

Alexandre Guillet, Françoise Argoul

► To cite this version:

Alexandre Guillet, Françoise Argoul. Uncertainty and information in physiological signals: explicit physical trade-off with log-normal wavelets. 2023. hal-04260114

HAL Id: hal-04260114

https://hal.science/hal-04260114v1

Preprint submitted on 26 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Uncertainty and information in physiological signals: explicit physical trade-off with log-normal wavelets

Alexandre Guillet<sup>1,2,‡</sup>, Françoise Argoul<sup>1</sup>

1. LOMA, Laboratoire Ondes et Matière d'Aquitaine, CNRS, Université de Bordeaux, Talence, France.

2. Division Biological Physics, Max Planck Institute for the Physics of Complex Systems, Dresden, Germany
‡ to whom correspondence should be sent

## Abstract

Physiological recordings contain a great deal of information about the underlying dynamics of life. The practical statistical treatment of these single-trial measurements is often hampered by the inadequacy of overly strong assumptions. Heisenberg's uncertainty principle allows for more parsimony, trading off statistical significance for localization. By decomposing signals into time-frequency atoms and recomposing them into local and flexible estimators, we propose a concise and expressive implementation of these fundamental concepts based on the choice of a geometric paradigm and two physical parameters. Starting from the spectrogram based on two fixed timescales and Gabor normal window, we then build its scale-invariant analogue, the scalogram based on two quality factors and Grossmann log-normal wavelet. These canonical estimators provide a minimal and flexible framework for single trial time-frequency statistics, which we apply to polysomnographic signals: EEG representations, coherence and mutual information between ECG-derived heart rate and respiration, and their spurious statistics.

## Keywords

Time-frequency analysis, Log-normal wavelet, Uncertainty principle, Polysomnography, Physiological signals, Coherence, Statistical significance, Mutual information, Single trial time-frequency statistics.

## 1. Introduction

Physiological signals have been extensively studied over the last two decades [1, 2], for their complex behaviours, ranging from periodic to stochastic and spanning a wide range of timescales. The formal identification and quantification of distant organ interactions (e.g. between respiration, heart and brain) has since been integrated into the field of “Network Physiology” [3, 4, 5, 6]. Pairwise interactions have been approached from different angles [7], such as synchronization, coherence, time delay, mutual information or transfer entropy, and then extended to higher-dimensional structures [6, 8]. Such developments require dedicated analytic tools to deal with simultaneous signals from different sources, with non-stationary dynamics spanning a broad range of timescales. The probabilistic approach is a key element to model the variability of these signals, considered as a realization of a stochastic process. The considerable difficulty of conditioning living systems into well-controlled and reproducible initial states limits their repeatability, so that in practice we are often forced to consider a physiological recording as a single trial experiment. The assumption of stationarity, at least
locally within well-chosen timescales, is a common simplification for the statistical treatment of such unique signals. However, their study is complicated by the multiple uncertainties associated with the variety of recorded events, in particular their timescales (of relaxation, correlation, evolution) and their possible separation. Time-frequency representation is another key component to investigate complex dynamics [9]. Many types of kernels, windows, filters, wavelets and other atomic shapes have been developed over the years to decipher rich signals. Dynamical information is unfolded in the time and frequency domains as more or less local features. Again, several uncertainty relations are available to characterize localization properties [10]. The diversity of these techniques and their lack of integration with the statistical approach requires expert choices that we intend to simplify in this paper. With this article, we pay tribute to Alain Arneodo for his pioneering work on the intersection between wavelet theory and physics [11, 12]. His late interest in biomedical images [13] has been a great stimulus in initiating the research developed here.

How can we make statistical estimates from individual non-stationary signals without knowing the underlying timescales? Is it possible to manage our uncertainty in terms of statistical significance and localization in time and scales? In this article, we present an elegant framework that succinctly answers these questions. It embeds statistical estimators in a time-frequency or time-scale analysis, to jointly handle the mix of dynamical and statistical information contained in biophysical signals. In particular, we design a flexible and versatile scalogram estimator based on log-normal wavelets, which we apply to a variety of signals. Our uncertainty and assumptions are given a precise shape by the choice of two quality factors, that control both the time-frequency resolution and the statistical significance of the analysis. The classical idea of the uncertainty principle or scale separation, between correlation and evolution timescales, is here revisited in a scale-invariant formalism constructed from time-frequency uncertainty atoms, available as statistical samples in a time-frequency analysis. Exact in the proposed formalism, this straightforward uncertainty relation unlocks the practical management of probabilistic and dynamical uncertainty (or information) in the vicinity of Heisenberg's uncertainty bound.

This article is divided into seven sections. Section 2 introduces concepts arising from a joint sta-
tistical and time-frequency analysis. We discuss, illustrate and formalize the relationships between
the assumption of timescale (or quality factor) separation, Heisenberg's uncertainty principle and the
statistical approach to unique signals. The continuous numbering of effective statistical degrees of
freedom or time-frequency uncertainty atoms is given an information-theoretic interpretation. Section
3 summarizes these concepts in an explicit parametric framework. Two paradigmatic time-frequency
estimators are introduced, one for the spectrogram and the other for the scalogram, called canonical
for their elementary expressions in terms of a pair of physical parameters, respectively timescales or
quality factors. In particular, their ratio gives the number of uncertainty atoms. The first paradigm
is based on the Gabor transform and the normal window, and the second on the wavelet transform
and the (Altes-Grossmann) log-normal wavelet. Sections 4, 5 and 6 contain the application of this an-
alytical framework to human physiological recordings of the neural, cardiac and respiratory activity.
These non-stationary and multiscale signals are part of the same polysomnography (shhs2-200901),
downloaded from the National Sleep Research Resource (Sleep Heart Health Study) [14, 15] and
used after ethics review board approval. In section 4, the intermittent and noisy neural dynamics
captured by an EEG signal illustrate the comparison of the two time-frequency paradigms. Their
limiting cases and their physical parameters are discussed. In contrast, the ECG signal introduced in
section 5 contains the highly rhythmic activity of the heart. The standard procedure for extracting
its fluctuating rate is translated into the log-normal wavelet framework; timescales and quality factors
are specified at each step and optimized against the standard estimation. In section 6, we analyse
the time-frequency coherence between the resulting heart rate signal and the respiratory oscillations.
Its spurious statistics is characterized in terms of the number of atoms and its significance is inter-
preted in terms of a time-frequency mutual information density. This completes the picture of the
trade-off between statistical significance and time-frequency localization in the analysis of signal trial
recordings. We conclude in section 7.

## 2. Single trial statistics in an uncertain dynamical context

In order to obtain valuable statistical information from single trial (unique) measurements of biophysical dynamics, it is often necessary to make assumption about their stationarity and ergodicity. The first moments of their distribution can be estimated at different times if the living system is assumed to be stationary in the short term. However, in the absence of prior dynamical information, the distinction between short and long term is uncertain. In discussing this practical issue, we consider two types of temporal scale that appear in the Heisenberg relation between duration and frequency uncertainties, and the need to separate them in order to create local statistical degrees of freedom. We then present the possibility of designing statistical estimators based on such temporal scale separation, which do not fix any arbitrary timescale parameter. Such a scale-free statistical approach is proposed as a parsimonious way to handle dynamical assumptions in biophysical contexts. Finally, we introduce the formalism associated with non-negative quadratic operators and their time-frequency and statistical interpretations, upon which the framework proposed in the next section is built.

2.1. *Statistics and timescales of evolution and correlation*

In a statistical approach, the data is divided into samples assumed to have the same probability
distribution, and statistical moments are then estimated from an empirical mean over these samples.
For a signal of duration *T* sampled at the frequency *f*<sub>s</sub>, we can access at most *N* = *Tf*<sub>s</sub> samples at
the cost of the strongest dynamical assumption of stationarity. The ratio of the macro and micro
timescales *T* and *dt* = *f*<sub>s</sub><sup>-1</sup> represents the number of data points *N* and the maximum number of
degrees of freedom available for a statistical estimation, if the assumption of a stationary ergodic
process is correct.

In order to reveal a possible time evolution of the statistics, the system can be assumed to be
stationary in the short term, i.e. for durations shorter than a meso-timescale $\tau_0$, and non-stationary
for longer durations. This weaker dynamical assumption comes at the cost of a smaller number of
statistical samples. Correlated dynamics can also be demonstrated by taking signal blocks of duration
$\tau_0$ as samples if the process is stationary. A time delay $\tau < \tau_0$ is introduced into the estimator of
the correlation (second order moment), or equivalently, a power spectrum with frequency $f > \tau_0^{-1}$ is
estimated (e.g. Welch's method). The spectral domain also facilitates the identification of periodic
oscillations in the data.

The analysis of both the temporal evolution on periods longer than $\tau_0$ and the correlation on
intervals shorter than $\tau_0$ have incompatible dynamical requirements. Statistical estimation is indeed
performed either within or across these intervals, based on $n = \tau_0/dt$ and $n = T/\tau_0$ statistical
samples, respectively. Fortunately, the probabilistic treatment of a signal interpreted as a single
realization of an evolving and correlated stochastic process is still possible. This is achieved by
further restricting the range of dynamical timescales and reducing the number of statistical degrees
of freedom: $n = \tau_+/\tau_- \ge 1$ samples are obtained as signal blocks of duration $\tau_-$, assumed to evolve
on timescales longer than $\tau_+$. Statistical estimates on unique signals with complex dynamics are only
possible by separating evolution and correlation timescales.
<figure><img src="image_2.png" /><figcaption>Figure 1: Schematic representation of information quanta as atoms of time-frequency uncertainty (Gabor windows or wavelets), even dispositions in the time-frequency plane, and different ways to compose them. (A) Different atoms, shaping the resolution trade-off between time and frequency. Shapes of order zero (blue) are more compact than shapes of higher orthogonal order (other colours). (B) Even sampling of the time-frequency plane with Gabor windows (left) and wavelets (right), endowing it with distinct geometries (Euclidean and hyperbolic respectively). The logarithmic frequency axis improves the flat visualization of wavelets at multiple scales, although their apparent area is no more representative of their atomic (unit) resolution. (C) A single-trial statistical estimate in time and frequency is based on a local average over several atoms, here $n = 5$, corresponding to both to the number of statistical samples and the area of the resulting time-frequency uncertainty. Comparable estimates can be obtained from different compositions and shapes of atoms, from left to right: contiguous along time or frequency, or superimposed orthogonal shapes.</figcaption></figure>

## 2.2. Time and frequency uncertainty versus statistics

The smallest distinguishable evolution timescale $\tau_+$ represents a time resolution or uncertainty $\Delta t = \tau_+$ in the analysis, and the longest correlation timescale $\tau_-$ corresponds to a spectral bandwidth or frequency resolution $\Delta f = \tau_-^{-1}$. Here, the lower bound of Heisenberg's uncertainty principle, $\Delta t \Delta f \ge 1$, appears as the condition for a statistical time-frequency analysis of the signal. The case of a unique timescale $\tau_+ = \tau_-$, thus $\Delta t \Delta f = 1$, has been described as a decomposition of the signal into information quanta called “logons” [16], or uncertainty atoms. These atoms are sketched in Fig. 1: in (A) for the balance between the time and frequency extent of the atomic uncertainty, in (B, left) for a fixed scale time-frequency decomposition, and in (C) for the uncertainty associated with several atoms forming a statistical estimate.

The notion of statistical degrees of freedom has been introduced into a spectral or a time-frequency analysis with orthogonal “tapers” (orthogonal atoms), yielding multi-taper estimators [17, 18, 19]. Single-taper spectrograms play the role of statistical samples for the multi-taper estimation of the spectrogram. Schematized in Fig. 1 (C, right), the multi-taper approach is implicitly related to the uncertainty relation, since the number of tapers is optimally chosen from the uncertainty product
$n = \Delta t \Delta f$. The zeroth order taper (a simple window) gives the best time-frequency resolution, but a single sample is insufficient for a statistical estimation.

The probabilistic approach to evolution and correlation in single complex biophysical recordings relies on the decomposition of the signals into time-frequency atoms. The separation of timescales to collect enough samples, $n = \Delta t \Delta f = \frac{\tau_+}{\tau_-}$, is a sine qua non. These uncertainty atoms can be further recomposed into a statistical estimator. The existence of such a range of scale separation cannot be checked a priori in living systems, but it appears to be a necessary trade-off between statistical information and dynamical uncertainty. The choice of relevant timescales in the analysis is therefore an important practical issue if we want to measure them. We then discuss how to avoid an arbitrary choice of timescales or resolutions as control parameters.

### 2.3. Scale-free time-frequency statistics with constant quality factors

By choosing two arbitrary timescales for a spectrogram estimation, the information about long-term correlation and short-term evolution is lost. This would allow the analysis of complex biological systems, with the crucial assumption of local stationarity on selected timescales. This goal is achieved with an affine time-frequency (or time-scale) representation, also called constant $Q$ representation, because the relative bandwidth $\Delta f/f$ is inversely proportional to the quality factor $Q$ and the dimensionless duration $f \Delta t \propto Q$ are both fixed. The wavelet transform and the scalogram treat the timescale as a variable (inverse frequency) whose relative precision increases with the parameter $Q$. Therefore, the resulting analysis has no arbitrary reference timescale. A hyperbolic geometry is known to be associated with the resulting time-frequency decomposition [20, 21], see Fig. 1 (B, right).

Analogous to the multi-taper spectrogram, a multi-wavelet scalogram has been formulated from orthogonal analytic wavelets [18] to control the degrees of freedom in the statistical estimation. In practice, the simpler smoothing approach is performed as an empirical statistical mean [22, 23, 24]. This local average in time or frequency provides a clearer control of the time-frequency uncertainty of the analysis.

In this article, we demonstrate that the statistical interpretation of a scalogram estimate is based on the distinction of two quality factors, $Q_+$, which controls durations $\Delta t$, and $Q_-$, which controls bandwidths $\Delta f$. Most importantly, we also relate their ratio $Q_+/Q_-$ to the number of uncertainty atoms or statistical degrees of freedom, playing the role of a scale separation in a scale-free analytic framework.

### 2.4. Continuous numbering of time-frequency atoms

The consistency of the probabilistic perspective on individual trajectories of complex dynamics derives from the assumption of an incompressible separation between correlation and evolutionary timescales. This leaves room for degrees of freedom for a statistical estimation, also interpreted as atoms of time-frequency uncertainty.

We reformulate these concepts by relying on a result in [19] that relates quadratic Hermitian (bilinear) operators on signals to the multi-taper approach of spectral statistics. Let $x$ be a (possibly multidimensional) signal or spectrum $[x_1 x_2 \dots]^T$, such as a discrete data matrix or a set of time series. The global energy (matrix) of $x$, possibly divergent for infinite signals, can be written as a scalar product along time or frequency: $xx^* = \int x(t)x^*(t)dt = \int \hat{x}(f)\hat{x}^*(f)df$, where $x^*$ is the conjugate transpose of $x$, see the definition of the Fourier transform $\hat{x}(f)$ in Eq. (A.1) in Appendix A.1.

A time-frequency statistical estimator operates on the data $x$. In particular, quadratic estimators
can be defined by an operator K which can take different forms of temporal or spectral kernels:

$$xKx^* = \iint x(t_1)x^*(t_2)\tilde{K}(t_1, t_2)dt_1dt_2 = \iint \hat{x}(f_1)\hat{x}^*(f_2)\hat{K}(f_1, f_2)df_1df_2 \quad . \quad (1)$$

The first notation refers to the matrix product, which applies to a discrete time or frequency variable. The following expressions are explicit integrals over the time or spectral continuum. When $x$ is a multidimensional signal (such as a multichannel simultaneous recording), the estimate is a matrix: each entry is a cross-spectrum (the diagonal terms being auto-spectra). If the operator $K$ is non-negative, the estimate can be interpreted as a $K$-localized energy. In the discrete and finite context of multi-taper estimators, the kernel is a Hermitian symmetric matrix (chosen to be positive-definite), whose eigenvectors are the multiple orthogonal tapers and the number of degrees of freedom is readily computed from its (positive) eigenvalues. This computation is equivalently expressed, directly from the kernel, as an effective number of eigenstates (orthogonal atoms), given by a beautiful trace formula [19]:

$$n = \frac{(\mathrm{Tr} K)^2}{\mathrm{Tr} KK^*} = \frac{\left(\int \tilde{K}(t, t) dt\right)^2}{\iint |\tilde{K}(t_1, t_2)|^2 dt_1 dt_2} = \frac{\left(\int \hat{K}(f, f) df\right)^2}{\iint |\hat{K}(f_1, f_2)|^2 df_1 df_2} \quad . \quad (2)$$

Again, the second and third equalities clarify the correspondence between discretized matrix notation and continuous integrals. In time, the Hermitian symmetry writes $\tilde{K}^*(t_2, t_1) = \tilde{K}(t_1, t_2)$, and as for the matrix product, the linear algebraic trace $\mathrm{Tr} K$ corresponds to a continuous integration. In the context of a statistical estimation, $n$ takes the interpretation of the equivalent number of complex degrees of freedom in the operator $K$.

The time-frequency interpretation of Eq. (1) is introduced with the Wigner-Ville representation:

$$W_{xx}(t, f) = \int x(t + \frac{\tau}{2})x^*(t - \frac{\tau}{2})e^{-i2\pi f\tau}d\tau = \int x(f + \frac{\eta}{2})x^*(f - \frac{\eta}{2})e^{i2\pi\eta t}d\eta \quad (3)$$

$$xKx^* = \iint W_{xx}(t, f)K(t, f)dtdf \quad . \quad (4)$$

A Hermitian operator $K$ translates into a real kernel $K(t, f) = K^*(t, f)$, which we choose to be non-negative. The spectrogram and the scalogram of $x$ can be estimated from the Wigner-Ville representation in terms of certain convolution operators $K$. In the case of the square of the Gabor or wavelet transform (introduced in the next section), the kernel covers a single uncertainty atom, $n = 1$, which is insufficient for a statistical interpretation. The Wigner-Ville representation itself corresponds to the degenerate case $n = 0$: no smoothing covers no atom of uncertainty, which is related to its infinite variance and possible negative values [25, 26].

In Wigner-Ville's quadratic time-frequency representation, the logarithm of the number of atoms:

$$\log n = -\log \frac{\iint K(t, f)^2 dt df}{\left(\iint K(t, f) dt df\right)^2} = H_2[K] \quad , \quad (5)$$

turns out to be the order 2 Rényi entropy $H_2[K]$ of the non-negative time-frequency kernel. Relations to temporal and spectral representations are given in Eqs. (B.1, B.2). Since an exponential entropy
counts an effective number of “states” in the distribution $K(t, f)$ [27], we arrive at the following interpretation: to gain statistical significance, $n$ quanta of information (Gabor’s logon) are averaged in the smoothing operation, but the distinct time-frequency localization of these $n$ atoms of uncertainty is lost.

The Rényi entropy have previously been applied to time-frequency representations of a signal [28, 29, 30, 31, 32], as a way of quantifying their concentration. Interestingly, the case of the Rényi entropy of order 2 was rejected because of its unsuitable handling of cross terms. This is irrelevant here since the kernel $K(t, f)$ is non-negative; we may interpret $n$ as the time-frequency area effectively covered by the kernel.

Although the continuous time-frequency integral expression for statistical degrees of freedom has already appeared in the literature [26], its practical usefulness remained unrecognized so that explicit eigendecompositions seemed inevitable. In the following sections, we present two time-frequency statistical estimators, whose number of atoms is simply derived from resolution parameters, shedding a new light on the formulation and the application of the uncertainty principle.

### 3. Canonical quadratic time-frequency estimators

In this section, we present a concise perspective on the prolific and mature field of continuous time-frequency representations by focusing on their single trial statistics. In what follows, we revisit two quadratic time-frequency estimators, the spectrogram and the scalogram, respectively in the Cohen and the affine classes, respectively [33]. They are built from a decomposition of the signal into paradigmatic time-frequency atoms, Gabor’s normal oscillatory window and Grossmann’s log-normal wavelet respectively. We call these representations canonical, each in its own paradigm, for their elementary and exact parameterization of the statistical and time-frequency properties discussed in section 2 in terms of two essential physical quantities. The time-frequency resolution is controlled by either a pair of timescales or a pair of quality factors, and their scale separation (ratio) corresponds exactly to the effective number of statistical degrees of freedom. Heisenberg’s uncertainty principle is strikingly reformulated as a special case of Eqs. (2, 5) and a requirement for a statistical approach of single-trial signals.

#### 3.1. Normal window paradigm

We start with the Gaussian or normal window of width (duration) $\tau$:

$$w_{\tau}(t' - t) = e^{-\pi(t'/\tau)^2}, \qquad (6)$$

whose spectrum $\hat{w}_{\tau} = \tau w_{\frac{1}{\tau}}$ has bandwidth $\frac{1}{\tau}$. The short-time Fourier transform of a signal $x$ under this normal window decomposes the spectrum $\hat{x}$ in the time domain, by shifting the window to any $t$. This is also known as the Gabor transform:

$$G_x(t, f; \tau) = \int x(t')w_{\tau}(t' - t)e^{-i2\pi ft'}dt' = \int \hat{x}(f')w_{\frac{1}{\tau}}(f' - f)\tau e^{i2\pi(f'-f)t}df' . \qquad (7)$$

The pure wave multiplied by the normal windows constitutes Gabor’s oscillating atom of uncertainty, or logon [16], shifted at any time-frequency location $(t, f)$. From this decomposition, schematized in Fig. 1 (B, left), the spectrogram of a stochastic process $x$ can be defined as a second order moment:

$$S_{xx}(t, f; \tau) = \frac{\sqrt{2}}{\tau} \mathbb{E}[G_x G_x^*], \qquad (8)$$

normalized by $\frac{\tau}{\sqrt{2}}$, the energy of the atomic window $w_{\tau}$. It is interpreted as a time-dependent power spectral density, in other words the time-frequency density of energy in $x$.
A statistical estimate of the spectrogram is obtained by smoothing the product $G_x G_x^*$, either along time or along frequency, assuming its probability distribution locally invariant. This recomposition of atoms from a local averaging operation again involves a normal window of different duration. The choice of a Gaussian smoothing leads to two equivalent definitions of the canonical spectrogram estimator:

$$ S_{xx}(t, f; \tau_+, \tau_-) = \frac{2}{\tau_+ \tau_-} \int G_x(t', f; \tau_-) G_x^*(t', f; \tau_-) w_{\tau_+}^2(t' - t) dt' , \quad \tau_+^2 = \tau_+^2 - \tau_-^2 \quad (9) $$

$$ = \frac{2\tau_-}{\tau_+} \int G_x(t, f'; \tau_+) G_x^*(t, f'; \tau_+) w_{\tau_-}^2 (f' - f) df' , \quad \tau_-^2 = \tau_-^2 - \tau_+^2 . \quad (10) $$

The corresponding compositions of atoms are schematized in Fig. 1 (C), the left and middle ones respectively. A time-smoothing extent $\tau_+$ is required to obtain a composition of duration $\tau_+$ from atoms of duration $\tau_-$ (similarly for $\frac{1}{\tau_+}$ in the frequency domain). Two different timescales $\tau_+ \ge \tau_-$ are thus introduced in these decomposition and recomposition operations, as summarized in the time-frequency kernel (Wigner-Ville representation) of this estimator $xK_w x^*$:

$$ K_w(t', f'; t, f) = 2(w_{\tau_+}(t' - t)w_{\tau_-}(f' - f))^2 . \quad (11) $$

This density function provides the exact shape of the composition of atoms in Gabor's normal window paradigm. See Appendix Appendix B.2 for purely temporal or spectral explicit expressions. This averaging operator covers a time-frequency box located at $(t, f)$ with widths $\Delta t = \tau_+$ and $\Delta f = \tau_-^{-1}$, which contains:

$$ n = \frac{\tau_+}{\tau_-} \quad (12) $$

uncertainty atoms, or equivalently $n$ statistical degrees of freedom. The scale order $\tau_+ \ge \tau_-$ is thus equivalent to Heisenberg's uncertainty relation, here $n = \Delta t \Delta f \ge 1$.

## 3.2. Log-normal wavelet paradigm

To change of paradigm, it could be sufficient to let the duration of the atom adapt to the frequency, $\tau = Q|f|^{-1}$, and instead fix a number of oscillations or a quality factor $Q$. The resulting time-frequency decomposition, proportional to $G_x(t, f; Q|f|^{-1})$, is known as the wavelet transform with Morlet's Gabor-like atom [34] or the closely related S-transform [35]. The relationships between these conventions are discussed further in the Appendix Appendix A. However, the previous development would fail at Eq. (10). Moreover, this normal (Gabor-Morlet) wavelet would lack important properties, useful in a constant-Q signal analysis and synthesis: it is neither analytic (only approximately so for $Q \gg 1$), nor admissible (without adjustment).

Instead, like others before us, we propose to use a “particularly good” wavelet in this respect [36], already anticipated in [37]:

$$ \hat{\psi}_Q(f') = e^{-\pi(Q \log(f'/f))^2} \Theta(f') , \quad (13) $$

where the Heaviside step function $\Theta$ handles the sign of the frequency domain (relevant for complex-valued signals). This is the log-normal (Altes-Grossmann) wavelet with quality factor $Q$, a symmetric limit case in the 2-parameter Morse family of analytic wavelets [38], also known as the log-Gabor filter (its radial part) in image processing [39, 40, 41]. Defined from the frequency domain, this uncertainty atom is scaled at frequency $f$ and shifted at time $t$ once multiplied by $e^{-i2\pi f't}$. There is no closed-form
temporal expression for $\psi_Q$, but an effective number of oscillations $Q = |f|\Delta t$ is expected [42]. The log-normal and the normal shapes are similar for $Q \gg 1$: in this regime, the Altes-Grossmann wavelet is indistinguishable from the Gabor-Morlet wavelet, but it retains its analyticity and admissibility with a small quality factor. It has other interesting properties, such as invariance with respect to derivation, orthogonality to any polynomial trend in the signal and a faster than polynomial decay in time and in frequency on both sides, which makes it particularly regular and localized [42].

The log-normal wavelet has a relative bandwidth $\Delta \log f = Q^{-1}$ around the peak frequency $f' = f$, so that we define the corresponding wavelet transform as a time-frequency decomposition:

$$X(t, f; Q) = \int \hat{x}(f')\hat{\psi}_Q\left(\frac{f'}{f}\right)e^{i2\pi f't}d f' = \int x(t')\psi_Q^*(f(t' - t))|f|dt' \quad . \quad (14)$$

Alternative spectral interpretations of the scale variable differ by a factor $\lambda^\alpha = \exp \frac{\alpha}{2\pi Q^2}$ [42], for instance $\alpha = \frac{3}{2}$ for the instantaneous frequency (phase derivative) of $\psi_Q$. From the atomic decomposition Eq. (14) schematized in Fig. 1 (B, right), the scalogram can be defined as a second order moment:

$$S_{xx}(t, f; Q)|f| = \sqrt{2Q} \mathbb{E}[XX^*] \quad (15)$$

normalized by $(\sqrt{2Q})^{-1}$, the admissibility coefficient of the wavelet $\psi_Q$. The scalogram is interpreted as a time-dependent power log-frequency density: the explicit factor $|f|$ is compensated when integrated with respect to $\frac{df}{|f|}$ (equal to $d \log f$ for $f > 0$). This notation as a product $S_{xx}|f|$ conserves the physical unit of the evolutionary spectrum $S_{xx}$ in both paradigms.

A statistical estimate of the scalogram is obtained by smoothing the product $XX^*$. The use of the log-normal wavelet with a different quality factor to compute the local average along log-frequencies (or scales) leads to the canonical scalogram estimator:

$$S_{xx}(t, f; Q_+, Q_-)|f| = 2Q_\downarrow Q_+ \int X(t, f'; Q_+)X^*(t, f'; Q_+) \hat{\psi}_{Q_\downarrow}^2\left(\frac{f'}{f}\right) \frac{df'}{|f'|}, \quad Q_\downarrow^{-2} = Q_-^{-2} - Q_+^{-2}. \quad (16)$$

This smoothing is formally analogous to Eq. (9), where time is replaced by log-frequency, but conceptually similar to Eq. (10) and Fig. 1 (C, middle), where the recomposition step is performed along the vertical (spectral) direction. The change in bandwidth $Q_\downarrow^{-1}$ is the log-frequency smoothing amount to achieve a relative bandwidth $Q_-^{-1}$ for the composition, from atoms with relative bandwidth $Q_+^{-1}$. Again, the two different quality factors $Q_+ \ge Q_-$ introduced here are visible in the kernel of the canonical scalogram estimator, when written in the spectral domain:

$$\hat{K}_{\psi}(f_1, f_2; t, f) = \sqrt{2Q_{-}} \left( \hat{\psi}_{Q_{+}}\left(\frac{f_1}{f_2}\right) \hat{\psi}_{Q_{-}}\left(\frac{f_1 f_2}{f^2}\right) \right)^{\frac{1}{2}} e^{-i2\pi(f_1 - f_2)t} \Theta\left(\frac{f_1}{f}\right) \Theta\left(\frac{f_2}{f}\right) \quad . \quad (17)$$

The time and time-frequency kernels for this operator are intractable integral forms, see Appendix Appendix B.1, as for the log-normal wavelet $\psi_Q$ in the time domain. Nevertheless, we numerically compute the time-frequency kernel $K_\psi(t', f'; t, f)$ from Eqs. (17, B.1) to confirm the extension of the time-frequency uncertainty, and note its non-negativity (up to numerical errors). The kernels are shown for different quality factors and numbers of atoms in Fig. 2. In contrast to the Gabor window paradigm, Eq. (11), they have different shapes in the $(t', f')$ (A, C) and $(t, f)$ (B, D) planes. They represent the area of the time-frequency uncertainty, well summarized by an uncertainty box of widths $\Delta t = Q_+|f'|^{-1}$ and $\Delta \log f = Q_-^{-1}$, thus confirming the interpretation of the quality factor, even for small values. The case $Q_+ = Q_-$ corresponds to the Wigner-Ville distribution of the log-normal wavelet, which is a quadratic time-frequency representation of the uncertainty atom.
<figure><img src="image_3.png" /><figcaption>Figure 2: Log-normal kernels: (A, C) composed of a fixed number of atoms $n = \frac{Q_+}{Q_-} = 4$ with different mean quality factors $Q_0 = \sqrt{Q_+Q_-} = \frac{1}{2}, 2, 8$ (from left to right), (B, D) with a fixed mean quality factor $Q_0 = 2$ and composed of an increasing number of atoms $n = 1, 4, 16$ (from left to right). The kernel is represented as a dimensionless function $K_\psi(t', f'; t, f) \Delta t \Delta \log f$, (A, B) in Wigner-Ville's time-frequency plane $(t', f')$ with fixed wavelet location $(t, f)$, and (C, D) in the wavelet's (target) time-frequency plane $(t, f)$ with fixed Wigner-Ville location $(t', f')$. Each kernel is drawn in 3 copies at different locations (marked with a plus sign). The white boxes materialize the uncertainty $\Delta t = Q_+ |f'|^{-1}$, $\Delta \log f = Q_-^{-1}$, with a time-frequency area equal to $n$.</figcaption></figure>

The kernel of the canonical scalogram estimator contains a number of uncertainty atoms:

$$n = \frac{Q_+}{Q_-}, \qquad (18)$$

which corresponds precisely to the time-frequency area of the illustrated uncertainty boxes, fixed in Fig. 2 (A, C) or increasing from $n = 1$ (single atoms, left) to $n = 16$ (right) in Fig. 2 (B, D). For a uniform visualization of these various time-frequency densities on linear-logarithmic axes, we have plotted their dimensionless form obtained by multiplying $K_\psi$ by $\Delta t \Delta \log f = n |f'|^{-1}$.

As before, there is also a temporal smoothing derived from the log-normal wavelet that corresponds to the horizontal composition of atoms schematized in Fig. 1 (C, left). Instead of acting on a single time variable in the product of wavelet transforms, it acts on two time variables:

$$S_{xx}(t, f; Q_+, Q_-)|f| = 2Q_+Q_- \iint X(t_1, f; Q_-)X^*(t_2, f; Q_-)\chi_{Q_\uparrow}(t_1 - t, t_2 - t)dt_1dt_2 \\ \text{where } \chi_{Q_\uparrow}(t, t') = \int \psi_{Q_\uparrow}^*(f't)\psi_{Q_\uparrow}(f't')|f'|df', \quad Q_\uparrow^2 = Q_+^2 - Q_-^2. \qquad (19)$$

This remarkable expression is specific to the log-normal wavelet paradigm.
Finally, the convenient intuition of the time-frequency statistics in the normal window paradigm is largely preserved in this log-normal wavelet paradigm. The separation of the correlation scale (frequency) and the evolution scale (duration) is achieved for each frequency of interest $f$ by the distinction $Q_+ \ge Q_-$, which corresponds to an uncertainty relation $n = |f|\Delta t\Delta \log f \ge 1$, allowing a statistical interpretation of the scalogram estimation.

### 3.3. Algorithmic considerations

A time-frequency decomposition $X(t, f)$ based on the Gabor or wavelet transform is expressed in terms of convolution integrals, which can be computed numerically from a pair of fast Fourier transforms (FFT). In practice, the continuity of the time-frequency representations is rendered by oversampling with respect to the time-frequency resolution $\Delta t$ and $\Delta f$ or $\Delta \log f$. In both paradigms, we implement reflective boundary conditions by extending the time-frequency domain in order to reduce the often unrealistic effects of periodic boundary conditions. Zero padding can sometimes be more realistic, especially in the spectral direction for band-limited experimental recordings. Boundary effects can remain at a distance of $k = 1$ to 3 times the resolution from the initial and final times and spectral boundaries.

The Gabor transform is naturally represented as a matrix that linearly samples the time and frequency domains on a 2-dimensional grid. The spectrogram $S_{xx}(t, f; \tau_+, \tau_-)$ can thus be efficiently estimated, from Eq. (9) or Eq. (10), by a pair of FFT along one or the other dimension.

Since $\Delta \log f$ is constant for wavelets, frequencies are sampled geometrically (regularly on a log-scale). However, the resolution $\Delta t = Q|f|^{-1}$ suggests a time sampling rate that increases with the frequency. For this reason, a single matrix (sometimes called stationary wavelet transform) is not a parsimonious discrete representation. This problem can be mitigated, at the cost of splitting the grid into multiple matrices, for example one for each octave, downsampling time by a factor of 2 when going down an octave. Details of this sampling approach in the similar context of the S-transform can be found in [43].

From the wavelet decomposition, the canonical scalogram estimator is more easily obtained by a smoothing with respect to log-frequencies, Eq. (16) (one pair of FFT), as compared to the double temporal convolution with a non-separable kernel, Eq. (19). Therefore, the discretization on a single grid is convenient to compute this second step, although it relies on the non-parsimonious matrix representation. Unfortunately, this leads to an unreasonable memory consumption when working on a large data set, such as overnight polysomnographic signals.

In a previous paper [44], we have proposed an alternative scalogram estimator, whose numerical implementation is better suited to long time series. The recomposition step performs a smoothing in the time domain instead of in the log-frequency domain, with a scaled normal window of duration $\tau = mQ|f|^{-1}$, see definition in Eq. (C.1). This scalogram estimator is compatible with a more parsimonious discretization since it implements the horizontal composition of atoms schematized in Fig. 1 (C, left). However, the resulting estimator can only be considered as an approximation of the canonical scalogram estimator ($Q_+ \approx mQ, Q_- = Q$): the exact expression for the number of atoms $n$ in terms of the quality factor of the wavelet and the scaling of the smoothing window is intractable; an approximation $n \approx (1 + m^2 \exp \frac{1}{2\pi Q^2})^{\frac{1}{2}}$ is computed in Appendix C using Lagrange's method.

For the canonical scalogram estimator, the smoothing step along log-frequency complicates the search for a more parsimonious numerical implementation since it has to be repeated for each matrix with a sufficient overlap (several times the resolution $\Delta \log f = Q_-^{-1}$). There is still much room left for improving parsimony and efficiency with an arbitrary degree of continuity in the numerical implementation of the canonical scalogram (or time-scale representations in general). A discrete wavelet
transform (such as the one in [45]) on the dyadic lattice followed by a switch to a continuous (arbitrarily oversampled) wavelet basis, or an analogous algorithm that produces local matrices (rather than coefficients) at each stage of the binary tree, may be potential ways of achievement beyond the FFT algorithmic strategy.

# 4. Application to an electroencephalogram

Time-frequency statistical estimators such as the spectrogram and the scalogram are ideal tools for approaching and deciphering the complex mixture of ever-changing fluctuations and oscillations that make up the brain activity. An EEG recording provides a rich and natural illustration for our discussion.

During sleep, the physiological interaction between organs is orchestrated by the autonomic nervous system; the EEG recording of potential fluctuations provides valuable information about these interactions, their regulation or their imbalance [4, 46, 6]. Polysomnographic signals (including EEG signals) are now used in routine clinical practice to identify not only sleep disorders but also their interaction with other metabolic, organic or cognitive functions [47].

The EEG shown in Fig. 3 is part of the polysomnography *shhs2-200901* from the Sleep Heart Health Study [14, 15] that were recorded on a 54-year-old sleeping woman. The cardiac (ECG) and respiratory signals from this polysomnography are also presented in the next section, in the same selected time interval: from the fourth hour after sleep onset to awakening.

## 4.1. Characteristics of the EEG recording

Electroencephalography (EEG) is a recording technique that uses electrodes in contact with the scalp to record the potential variations produced by neurons in the brain. Each electrode locally integrates a collective and spatial neural activity. The resulting signal, shown in Fig. 3 (B), exhibits combinations of erratic (scale-free) and rhythmic dynamics over a wide range of frequencies (more than two decades). Unlike cardiac or respiratory signals, no single fundamental mode can be identified in this neural signal. This EEG signal measures the electrical potential between the points C4-A1 with a resolution of 1 µV, a sampling frequency $f_s = 125$ s⁻¹, and instrumental filters (high-pass at 0.15 s⁻¹ and notch at 60 s⁻¹).

The complexity of an EEG can be represented in the time-frequency plane to distinguish both its evolutionary and spectral aspects, see Fig. 3 (D, F). The traditional way of dealing with this complexity is to divide time into discrete stages and frequency into discrete bands. Wake-sleep stages have been subdivided into wake (W), rapid eye movement (REM) and non-REM (or NREM) stages, and the NREM stage has been further decomposed into N1 (lightest), N2 and N3 (deep sleep) stages. This classification is based on a set of criteria about the relative EEG power in each frequency band (possibly supplemented by the EOG and EMG), within 30-second epochs. The sequence of these stages, available in most polysomnographic databases as a clinician annotation every 30 s, is a simplified representation of sleep called a hypnogram, see Fig. 3 (A). With the aim of empirically distinguishing neural waves, the frequency bands are roughly defined in slices of 4 s⁻¹ and denoted by Greek letters: $\delta$ up to 4 s⁻¹, $\theta$ from 4 to 8 s⁻¹, $\alpha$ from 8 to 12 s⁻¹, $\sigma$ from 12 to 16 s⁻¹, $\beta$ from 16 to 20 s⁻¹, and $\gamma$ above 20 s⁻¹. Thus, $\gamma$ waves are intense during wakefulness whereas they disappear during NREM sleep, which is characterized by intense $\delta$ waves (especially during the N3 stage). In contrast, during REM sleep, the EEG power is low in all bands.

Time-frequency statistics are particularly well suited to the analysis of such “unique” (in the sense of strictly non-reproducible) signals, which intermittently switch from highly stochastic (scale invariant) to more regular dynamics during the night. To illustrate brain activity, we use the canonical
spectrogram and scalogram estimators, and we discuss the parameters of these continuous atomic decomposition and recomposition operations, in both the window and wavelet atom paradigms.

<figure><img src="image_4.png" /><figcaption>Figure 3: Electroencephalogram (EEG) of a sleeping person, comparison of different representations, in particular, the spectrogram and the scalogram. (A) Hypnogram: EEG scoring of the sleep stages “wake”, “rapid eye movement” (REM) and “non-REM” stages 1, 2 and 3 (light to deep sleep). (B) EEG signal (µV), from electrodes C4-A1. (C) Fourier amplitude spectrum $|\hat{x}(f)|/T$. (D) Spectrogram estimation $S_{xx}(t, f; \tau_+, \tau_-)$ with $(\tau_+, \tau_-) = (20, 1)$ s, and (E) its persistence spectrum (grey) and power spectral density (blue). (F) Scalogram estimation $S_{xx}(t, f; Q_+, Q_-)|f|$ with $(Q_+, Q_-) = (100, 5)$, and (G) its persistence spectrum (grey) and power log-frequency density (red). The colour bars for (D, F) are aligned with the corresponding axes in (E, G), with a perceptually uniform contrast between quantiles [0.1, 0.99]. (H) Time-dependent power estimated from (D, blue) and (F, red), and mean power (EEG variance, grey).</figcaption></figure>

## 4.2. Single trial time-frequency statistics without orthogonal atoms

The use of a time-frequency statistical analysis to study EEG recordings has already been advocated, based on the multi-taper spectrogram [48]. We intend to refine this perspective, starting with a criticism: the explicit use of orthogonal atoms such as multi-tapers is unnecessary, and there is no explicit control of the time-frequency resolution. The continuous numbering of atoms Eq. (2) is straightforward for any (non-negative) quadratic operator, and those we introduce in this article are explicitly parameterized in terms of the resolutions. Some general considerations are worth mentioning when choosing these analytical tools. As we have already pointed out in previous sections, this choice builds on compromises.

On the one hand, it is generally accepted that the best quadratic time-frequency representations are those with the highest resolution. On the other hand, an averaging based on orthogonal windows (tapers) is often pushed forward as optimal for estimating statistics from these representations. This implies a systematic trade-off between time-frequency resolution and the number of statistical degrees
of freedom. The discrete orthogonal approach explicitly details the shape of each independent atom, but the resulting resolution is implicit. The converse is true for the continuous (smoothing) approach: the resolution is explicitly parameterized, but the exact shape of the uncertainty atoms is implicit. We further argue that the latter is contingent: the effective number of atoms is known exactly.

Indeed, for the same resulting composition, alternative continuous decompositions (into atomic windows or wavelets) and recompositions (into a statistical estimator) are available; see Eqs. (9, 10) for the spectrogram or Eqs. (16, 19) for the scalogram. Thus, the duration of a window can be thought of as belonging to $[\tau_-, \tau_+]$, or the number of wavelet oscillations to $[Q_-, Q_+]$, but only the final composition of atoms matters (characterized by the bounds of these intervals). Moreover, we have shown that both the resulting time-frequency resolution and the effective number of statistical degrees of freedom can be parameterized exactly as two timescales or quality factors and their ratio respectively.

Therefore, the spectrogram of the EEG $S_{xx}(t, f; \tau_+, \tau_-)$ in Fig. 3 (D) provides the same quality of time-frequency statistical estimates as the multi-taper spectrogram [48], without relying on orthogonal windows and with an explicit control of the resolution: here $\tau_+ = 20$ s = $\Delta t$ and $\tau_- = 1$ s = $1/\Delta f$. The use of a perceptual colour coding [49] (on a logarithmic axis) is also an important aspect to prevent visual artefacts from hindering a quantitative assessment of the time-frequency representation. The choice between the spectrogram and the scalogram is discussed next.

### 4.3. Resolutions in the spectrogram and the scalogram

The choice of a time-frequency paradigm is another critical aspect of the EEG analysis. In Fig. 3 (D, F), we try to compare the spectrogram (D) and the scalogram (F) as fairly as possible: both estimates are composed of $n = 20$ uncertainty atoms, the spectral ranges are chosen to have the same apparent resolution (along their respective linear and logarithmic axes), and the colour scale is bounded by magnitudes chosen as the same quantiles, at 10% and 99% of the magnitude distribution in the images.

The scalogram $S_{xx}(t, f; Q_+, Q_-)|f|$ in Fig. 3 (F) is similar to the spectrogram in Fig. 3 (D) once the factor $|f|$ is omitted: $S_{xx}(t, f; Q_+, Q_-) \sim S_{xx}(t, f; \tau_+, \tau_-)$, up to their different time-frequency resolution. The factor $|f|$ recalls the main specificity of the wavelet approach: the spectral resolution is equal on a log-frequency axis ($\frac{\mathrm{d}f}{|f|} = \mathrm{d}\log f$ for $f > 0$). The spectrogram and the scalogram are time-frequency densities of the signal energy, and the factor $|f|$ appears to keep this interpretation on a log-frequency axis.

The frequency dependence of the time resolution in the scalogram is well visible for impulsive events. The large upper quality factor chosen in Fig. 3 (F), $Q_+ = 100$, is not the most suitable to visualize single oscillations (as it averages over 100 of them), especially at low frequencies where single events spread easily over several minutes. Nevertheless, the time resolution at high frequencies is much better in the scalogram (F) than in the spectrogram (D): at $f = 20$ s⁻¹, singular events spread over $\Delta t = 5$ s in the former compared to $\Delta t = 20$ s in the latter. The inversion of the temporal resolution between the scalogram and the spectrogram occurs at the frequency $Q_+/\tau_+ = 5$ s⁻¹, and the same inversion occurs for the spectral resolution, $Q_-/\tau_- = 5$ s⁻¹.

Focusing now on the spectral resolution, the spectrogram representation appears to be a natural extension of the traditional discretization of the EEG into (linearly spaced) frequency bands, approximately every 4 s⁻¹ [48]. But the existence of a neural activity on multiple timescales [50], supported by the similarity between the EEG and a pink noise (self-similar process), clearly points to the scale-invariant wavelet paradigm and its scalogram. Indeed, the observation of the EEG fluctuations shows that the spectral range is only limited by the instrumental characteristics, especially at low frequencies, where the effect of the high-pass filter is particularly marked in Fig. 3 (F, G).
This limitation is hidden in the spectrogram, since the spectral range below its frequency resolution $\Delta f = 1 \text{ s}^{-1}$ is not resolved.

The brain waves $\beta-\gamma$ (fast) and $\delta$ (slow) are particularly irregular compared to the intermediate brain waves ($\theta$ and especially $\alpha$ and $\sigma$), which are more rhythmic. In the scalogram, the irregularity of brain waves can be observed independently of their frequency as their relative bandwidth, visible as their spread in log-frequency. We cannot see this aspect from the spectrogram, which relies heavily on its fixed timescales, in particular $\tau_- = 1/\Delta f$. The difference is striking when comparing the extension of weaker high-frequency activity with the narrowing of the stronger low-frequency activity in the spectrogram. This picture is more balanced in the scalogram.

The frequency-localized structures in the intermediate band (between 8 and 16 s⁻¹) are composed of $\sigma$ waves (sleep spindles, in the non-REM phase) and $\alpha$ waves (strongest in the wake phase). In Fig. 3, the spectrogram benefits from a better spectral resolution in this band than the scalogram. The latter can, at best, distinguish between frequencies 9 and 11 s⁻¹ with the quality factor $Q_- = 5$.

Our conclusion is nuanced: the scalogram appears to be the preferred representation for a broad perspective on the diversity of timescales and on the regularity of the recorded events, whereas the spectrogram is useful for focusing on neural phenomena around the reference timescales.

### 4.4. Marginals and limits: retrieval of temporal or spectral statistical estimators

In the following, we relate the canonical spectrogram and scalogram estimators to simpler temporal, spectral or global estimators and non-statistical local representations, and comment on their application to the EEG.

#### 4.4.1. Maximal statistics: no resolution

The time-frequency integration of the spectrogram and the scalogram yields the same global energy: $xx^* = \iint S_{xx}(t, f)dtdf$. For the spectrogram, the integrator $\frac{df}{|f|}$ cancels the factor $|f|$. Various average densities can be obtained dividing this energy by the bandwidth, the duration or the time-frequency domain. The global variance of a centred signal $x(t)$ corresponds to the mean power $xx^*/T$. For the EEG in Fig. 3, the mean power is about 500 $\mu\text{V}^2$ (grey line in H), proportional to the electrical power and an assumed constant or slowly varying impedance.

The average over all time and frequency corresponds to the limits $\tau_+, Q_+ \to \infty$ and $\tau_-, Q_- \to 0$ (in fact, the averages over positive and negative frequencies are different for $Q_- \to 0$). Although the resulting number of degrees of freedom available for the statistical analysis is maximal, the joint assumption of stationarity and self-similarity is the strongest. In this limit, time-frequency representations are not relevant: there are no degrees of freedom left to obtain localization information.

#### 4.4.2. Spectrally resolved statistics

A spectral marginal can be obtained by averaging the spectrogram or the scalogram over time. In the window paradigm, it is known as the power spectral density or the periodogram and can be written as $S_{xx}(f; \tau)$. In the wavelet paradigm, the corresponding spectrum $S_{xx}(f; Q)|f|$ can be interpreted as the power log-frequency density. This purely spectral statistical estimate corresponds to $\tau_+ \to \infty$ in the window paradigm or $Q_+ \to \infty$ in the wavelet paradigm. This limit corresponds to the assumption of a stationary signal. The spectral resolution $\Delta f$ is still given by either $\tau_-^{-1}$ or $Q_-^{-1}|f|$, and the loss of temporal resolution increases the number of statistical degrees of freedom, bounded by the duration of the signal $T$: $n \approx T\Delta f$.

In practice, we compute the spectral marginal directly from the squared time-frequency decomposition, by replacing the ensemble average in Eqs. (8, 15) by the time average. This is equivalent to smoothing $|\hat{x}(f)|^2/T$ along the (linear or logarithmic) frequency axis with a normal kernel.
The spectral marginals of the EEG are shown in Fig. 3 (E, G). The EEG spectrogram (D) is aligned with its power spectral density (E) and the EEG scalogram (F, blue) with its power log-frequency density (G, red). Persistence spectra are also shown (time-sampled spectrogram and scalogram), giving a glimpse of the temporal variability associated with sleep stages and more room to explore higher order moments, or other statistical quantities such as quantiles (used for colour-coding).

We can see that the $1/f$ trend of the EEG in the spectrogram and the power spectral density is compensated in the scalogram and in the power log-frequency density. Switching to the wavelet scale-invariant paradigm therefore provides a better contrast to EEG representations. For instance, the narrowband activity of the $\alpha - \sigma$ waves is clearer Fig. 3 (G) than in (E), in spite of a lower resolution: $\Delta f \approx 2$ to $3$ s$^{-1}$ versus $1$ s$^{-1}$.

### 4.4.3. Temporally resolved statistics

The converse temporal marginal consists of averaging or integrating evolutionary spectra over frequencies. Called the power of the signal $P_{xx}(t; \tau)$ or $P_{xx}(t; Q)$, the integral over all frequencies of $S_{xx}$ estimates $\mathbb{E}[x(t)x^*(t)]$, a time-dependent estimator of the variance for a centred stochastic process. The temporal marginal is associated to a self-similarity assumption, analogous to the stationarity assumption [51] for the spectral marginal. More precisely, the average over $f$ of the spectrogram estimates an intensity per s$^{-1}$, assuming no correlation (white noise); the average over log $f$ of the scalogram estimates the intensity per octave or decade (depending on the base of the logarithm), assuming it constant (pink noise). Regardless of their interpretation, both averages are proportional to the same variance estimator $P_{xx}$.

This broadband case corresponds to the limit $\tau_- \to 0$ in the normal window paradigm, and to $Q_- \to 0$ in the log-normal wavelet paradigm, except for the following subtlety. For the scalogram, this limit does not average positive and negative frequency components, but keeps them separate. The resulting powers are those of the so-called analytic and anti-analytic signals $x^{\pm}(t)$, which carry distinct information when $x(t)$ is complex-valued. Their sum is the total power.

The time-dependent EEG power estimates are compared in Fig. 3 (H), where $P_{xx}(t; \tau)$ (blue) boils down to a simple smoothing of $|x(t)|^2$ with a normal window, while the wavelet estimation $P_{xx}(t; Q)$ (red) has an unspecified resolution that varies with the timescale of the fluctuation under consideration. Both estimates may differ locally, but their trend is very similar. The switch between NREM and REM (see the hypnogram, A) has a clear effect on the time-dependent variance, which can drop from 1 mV during deep sleep (N3 stage) to almost 100 µV during REM.

### 4.4.4. Maximal resolution: no statistics

One could argue that the maximum time-frequency resolution is achieved with the “skeleton” of quadratic representations, namely the Wigner-Ville distribution, Eqs. (4, 3). This extreme case with unspecified atomic uncertainty, $n = 0$, leads to degenerate properties, such as a divergent variance or negative cross terms.

In the present framework, we consider that the highest resolution is achieved by the atomic time-frequency decomposition, $n = 1$. Without recomposition, the corresponding quadratic representations are the non-averaged spectrogram $S_{xx}(t, f; \tau, \tau) = \frac{\sqrt{2}}{\tau}G_x(t, f; \tau)G_x^*(t, f; \tau)$ and the non-averaged scalogram $S_{xx}(t, f; Q, Q) = \sqrt{2}QX(t, f; Q)X^*(t, f; Q)$, for which $\tau_+ = \tau_-$ and $Q_+ = Q_-$. This is the common conception of the spectrogram and the scalogram, but these transforms of the stochastic trajectory $x$ cannot be considered as statistical estimators (at best trivial ones).

The Gabor and wavelet transforms $G_x(t, f; \tau)$ and $X(t, f; Q)$ contain information about both the amplitude and the phase of oscillations. Working directly with linear time-frequency decomposition can be of interest for studying the micro-architecture of sleep and its rhythmic brain waves, such
as sleep spindles ($\sigma$) and slow oscillations ($\delta$) in deep sleep (N3), or $\alpha$ waves. In addition from the amplitude, the phase information enhances the analysis of rhythmic behaviour, as applied to the cardio-respiratory activity in the next section. Although this approach does not typically rely on statistics, a few statistical degrees of freedom are shown to be beneficial.

# 5. Decoding the heart rhythm and its variable rate

As the central conductor of the dual circulatory system, the contractile activity of the heart adapts in real time to the needs of the entire organ system and is finely regulated by the sympathetic and parasympathetic branches of the autonomic nervous system. The heart rate variability (HRV) has been proposed to maintain health more than three decades ago [52]; HRV is still the subject of a long-lasting interest [53, 54]. The further improvement and comparison of HRV identification methods remains an important issue. Non-invasive techniques have been developed to monitor the rhythmic activity of the heart, such as electrocardiography, photoplethysmography, and mechanocardiography. We focus here on an electrocardiogram (ECG), selected from the same polysomnography (shhs2-200901) as before. The ECG waveform has three main components: a P wave (atrial depolarization), a QRS complex (ventricular depolarization), and a T wave (ventricular repolarization). The detection of RR time intervals from successive R-peaks is a standard and discrete method of measuring the cardiac period from the ECG [55]. This discrete method is here translated into the continuous log-normal wavelet framework. We propose a statistical estimator of a heart rate signal, derived from the phase of the rhythmic heart activity.

## 5.1. Identification of QRS complexes from ECG signals

As a reference conventional method and for comparison with the method proposed here, we use the *jqrs* algorithm from the PhysioNet Cardiovascular Signal Toolbox [56], which first computes the QRS intensity using a matched filter (band-pass) and a sequence of signal processing operations [57, 58]. The intensity maxima in the resulting positive pulse train are then detected by a thresholding procedure. To validate the detection method or to reject unreliable RR intervals, a signal quality index can be calculated, for example from a comparison with another detector [59]. Automatic error correction, which interpolates missing beats and removes false beats based on physiological assumptions, can improve the estimation of the heart period (NN intervals for “normal”). Adjustment of sufficient detection parameters (such as thresholds) can give satisfactory results without the need for expert (but less reproducible) modifications.

The electrical activity of the heart, as recorded by the ECG, is shown in Fig. 4. The sharp and prominent QRS oscillations in this rhythm are a clear marker of a heartbeat, visible in the time domain, either in the signal (A, zoom in E top) or in time-frequency representations (wavelet transforms B and C), but not in the spectral domain (D). The QRS complex refers to a sequence of 3 local extrema in the ECG waveform, in the time domain, including the sharpest peak. Its detection (red circles in E, F) is routinely used to extract the heart rate (black dots in B, C, E, F) from the time intervals between beats. Due to the frequent polarity reversal of the ECG signal (depending on the placement of the electrodes on the chest), a proxy for the QRS amplitude is often preferred to the actual voltage peak, after band-pass filtering baseline drift and high frequency noise.

This “preprocessing” operation is done simply as the modulus of the wavelet transform $y(t) = |X(t, f_{qrs}; Q_{qrs})|$ at the characteristic scale of the QRS complex. Note that $y(t)$ is specified by both a frequency location $f_{qrs}$ and a resolution: the quality factor $Q_{qrs}$. The modulus of two wavelet transforms with different quality factors $Q = 1$ and $Q_{qrs} = 3$ are shown in Fig. 4 (B and C) in the time-frequency domain, and a time slice of the latter is shown at the selected frequency $f_{qrs} = 14 \text{ s}^{-1}$ (F top). The non-negative signal $y(t)$ represents the QRS amplitude, whose local maxima are a key
<figure><img src="image_5.png" /><figcaption>Figure 4: (A) ECG signal $x(t)$ (μV). (B) $|X(t, f; Q)|$ with $Q = 1$. (C) $|X(t, f; Q_{qrs})|$ with $Q_{qrs} = 3$. The colour scale (μV) is aligned with (D). (D) Amplitude spectra: $|\hat{x}(f)|/T$ over one RR interval (black triangles) or the selected time interval (grey line), median (thick lines), and root-mean-square (thin lines) over time of (B, blue) and (C, red). (E) $x(t)$ and $X(t, f; Q)$ (zoom in the boxes in A, B), with modulus-angle colour coding near the heart rate. (F) Amplitude of the QRS complexes $y(t) = |X(t, f_{qrs}; Q_{qrs})|$ at $f_{qrs} = 14$ s⁻¹ (grey dashed line in C, D) and $Y(t, f; Q)$ with $Q = 1$. Heart rate (dots in B, C, E, F) estimated from detected QRS complexes (red circles in E, F).</figcaption></figure>

element in the detection of the QRS oscillations and the estimation of the heart rate (inverse of the RR time intervals). This standard discrete estimation is shown on all time-frequency representations in Fig. 4 (B, C): it clearly follows the frequency modulations of the fundamental mode of the ECG.

The angular components of the wavelet transforms $X$ (ECG) and $Y$ (QRS amplitude) are compared in Fig. 4 (E and F bottom). It can be seen that the cardiac cycle, as detected by the discrete reference method, follows the angle at the cardiac frequency. This is particularly robust in the phase of the QRS amplitude (F) as compared to the ECG (E), which can also encode occasional oscillations of the baseline (motion artefacts). This suggests an alternative and continuous description of the heart rate, as the derivative of this cycling heart phase, which will be presented in the following sections.

## 5.2. Important considerations for the choice of quality factors

Defining a rate estimator from a single fluctuating rhythmic signal $x(t)$ is generally difficult without prior contextual information or preprocessing. We propose to translate these decisions in terms of the time-frequency location $(t, f)$ plus two quality factors $(Q_+, Q_-)$ as resolution parameters. Let us recall their dynamical and statistical meaning in the context of a rhythmic or oscillatory behaviour, before working out a definition.

The dynamical uncertainty balance is first established by fixing a quality factor $Q$ ($Q_+ = Q_-$ for the moment) which allows distinguishing between a temporal and a spectral location $(t, f)$. In
the time domain, $Q$ is a dimensionless time resolution, $f\Delta t$, interpreted as the effective number of oscillations in the analysing wavelet. Oscillations are also localized in the frequency domain, where $Q^{-1}$ is a log-frequency resolution $\Delta \log f$ (alternative expression for relative bandwidth $\frac{\Delta f}{f}$). The quality factor $Q$ is smaller than that defined from the bandwidth in [60] by a factor of about $\sqrt{\pi}$. The quantities $Q$ and $Q^{-1}$ are very close to the full width at half maximum of the wavelet [42], respectively in time (for $f = 1$) and in log-frequency.

The specificity of the waveform of a regular rhythm can be distinguished from a purely circular (sine) oscillation in the frequency domain by the presence of multiple harmonic modes with frequencies $f_k = k f_1$ for positive integer orders $k$. In the wavelet paradigm, the order $k + 1$ can be distinguished from the order $k$ whenever $\Delta \log f = Q^{-1} < \log \frac{k+1}{k}$. Furthermore, wavelets with $Q$ oscillations can only analyse a full rhythmic period $\Delta t = f_1^{-1}$ up to the frequency $Qf_1$. As a consequence of the following logarithmic inequality [61]:

$$ \frac{1}{Q + \frac{1}{2}} < \log \frac{Q + 1}{Q} < \frac{1}{Q}, \qquad (20) $$

the spectral discrimination of $Q > 0$ partials requires a wavelet with $Q$ to $Q + \frac{1}{2}$ oscillations. In other words, a quality factor $Q$ allows the distinction of $Q - \frac{1}{2}$ to $Q$ harmonic orders. This inequality illustrates two intertwined and interchangeable interpretations of the quality factor, in terms of time-frequency resolution and the distinction of partials in a rhythm. We can visualize their relevance in Fig. 4, in the wavelet transforms of the ECG signal (A) with quality factor $Q = 1$ (B) or 3 (C).

The introduction of an effective number of atoms $n = \frac{Q_{\pm}}{Q_{-}} \ge 1$ as degrees of freedom of the quadratic statistical estimation actually splits the temporal and the spectral interpretations. The upper quality factor $Q_{+} = f\Delta t$ controls the new time uncertainty, while the lower one $Q_{-} = \Delta \log f$ controls the new frequency uncertainty. This composition of uncertainty atoms complies with Heisenberg's principle, whose uncertainty area $\Delta f\Delta t$ of at least one atom is rewritten:

$$ \frac{Q_{+}}{Q_{-}} \ge 1. \qquad (21) $$

This inequality relates the possibility of a statistical estimation (away from the lower bound) to the extent of a finite scale separation between slow temporal modulations and fast carrier frequencies. Remarkably, in the wavelet paradigm, these timescales do not need to be fixed; only their minimum separation is fixed as the ratio of the quality factors.

These very practical and natural interpretations of the quality factors as the number of wavelet oscillations, the number of distinguishable harmonic orders, and the number of uncertainty atoms (for their ratio), are silently orchestrated by the presence of the circle constant in the definition of the log-normal wavelet Eq. (13). These realizations have motivated a rescaling of the quality factor, and hence a slight redefinition of the log-normal wavelet, as compared to our previous definition in [44, 42, 62].

### 5.3. Fluctuating rates: from the analytical signal to a time-frequency statistical estimation

Given a single-channel oscillatory signal $x$, its amplitude and frequency modulations (AM and FM) can be measured from the modulus and angle of the analytic version of the signal. Denoted $x^{+}$, it is here defined (together with the anti-analytic version $x^{-}$) as the complex helical signal corresponding to the part of the Fourier spectrum $\hat{x}$ that is supported on positive (respectively negative) frequencies:

$$ x^{\pm}(t) = \int \Theta(\pm f) \hat{x}(f) e^{i2\pi ft} df. \qquad (22) $$
For any real signal $x(t)$, a complex rate can be defined from the analytic version of the signal as the following time derivative ($D_t$):

$$F_x(t) = \frac{1}{i2\pi} D_t \log x^\dagger(t) \quad , \qquad (23)$$

whose real part estimates the fluctuating frequency of the oscillation (the imaginary part is a magnitude rate).

In practice, the resulting modulated rate signal is corrupted by random noise and non-linearities (non-circular waveforms), which can be mitigated by band-pass filtering the original signal. This is precisely the purpose of a time-frequency decomposition such as the wavelet transform: the band is determined by a fixed frequency $f$ and its bandwidth, absolute ($\Delta f$) or relative ($\Delta \log f$). The log-normal wavelet representation $X(t, f; Q)$ has the advantage of being an analytic signal for any fixed $f$, for which the analytic (and anti-analytic) signal $x^\pm$ is a limit (as well as a marginal) case:

$$x^{\pm}(t) = X(t, \pm|f|; Q \to 0) = Q \int \Theta(\pm f) X(t, f; Q) \frac{df}{|f|} \quad , \qquad (24)$$

where $\Theta$ is the Heaviside step function. For this reason, $F_x(t)$ can be considered as the limit when $Q \to 0$ of the more general partial time derivative ($\partial_t$):

$$F_x(t, f; Q) = \frac{1}{i2\pi} \partial_t \log X(t, f; Q) = \lambda^{\frac{1}{2}} f \frac{X(t, \lambda f; Q)}{X(t, f; Q)} \quad , \quad \lambda = \exp \frac{1}{2\pi Q^2} \quad . \qquad (25)$$

The second equality highlights the practical way to derive with respect to time, as a slight frequency scaling specific to the log-normal wavelet framework. Note that the limit $Q \to \infty$ corresponds to $F_x \to f$. The complex rate inherits the time-frequency resolution of the wavelet transform. Here, the quality factor counts the number of oscillations of period $|f|^{-1}$ used to estimate the complex rate.

If the oscillatory signal $x$ has a modulated frequency $F(t)$, then the real part $\Re\{F_x(t, f; Q)\}$ can be used to estimate this varying frequency, provided that the frequency band selected with $f$ and $Q$ contains only this mode. The imaginary part captures the modulation rate of the magnitude of the oscillation (instead of the derivative of its phase).

Considering the scale continuum $f$ at a fixed time $t$, there exist basins of attraction that yield the same frequency estimate. This “squeezing” phenomenon can be used to reassign time-frequency atoms (or their intensity) in the same basin of attraction to their estimated frequency: the real part of Eq. (25). This process is known as frequency reassignment for quadratic evolutionary spectra or as synchrosqueezed transform for time-frequency decompositions such as the wavelet and Gabor transforms [63, 64].

In practice, we intend to use $F_x(t, f; Q)$ as a signal, after selecting a convenient band $(f, Q)$. For natural recordings $x$, however, this complex rate both vanishes and diverges at certain locations related to zeros (phase defects) of the time-frequency decomposition. This singular behaviour has been described in terms of a determinantal point process in the case of the white Gaussian noise [21]. Their perturbation on reassigned or synchrosqueezed representations is mitigated by an additional parameter: these regions with modulus below a certain threshold are excluded [64, 65].

We propose here to average the complex rate Eq. (25) over more than one time-frequency atom, in order to preserve the continuity of the estimator in situations with a noisy oscillator, quite common in natural signals. Instead of applying a local mean directly to $F_x(t, f; Q)$, we introduce a quadratic weight $|X(t, f; Q)|^2$ to favour higher amplitudes, and we normalize it with the average local intensity, which is the scalogram estimator. This leads to the generalized definition of the complex rate as
a statistical estimator, involving both quality factors ($Q_+$, $Q_-$), the signal $x(t)$ and its derivative $\dot{x}(t) = D_t x(t)$:

$$F_x(t, f; Q_+, Q_-) = \frac{S_{\dot{x}x}(t, f; Q_+, Q_-)}{i2\pi S_{xx}(t, f; Q_+, Q_-)} \quad (26)$$

The newly introduced control parameter corresponds to the number of atoms, $\frac{Q_+}{Q_-} \ge 1$, taken as statistical samples. The statistical estimator $F_x(t, f; Q_+, Q_-)$ is free from the poles and zeros of $F_x(t, f; Q)$, thus stabilizing the complex rate estimation. A multi-taper reassignment procedure has been proposed in [66]; i.e. the frequency reassignment of the canonical scalogram estimator based on Eq. (26) is expected to achieve a similar goal of improved concentration and stability.

The broadband limit $Q_- \to 0$, which can be denoted $F_x(t; Q_+)$, was introduced in [44] as a more stable alternative to the $F_x(t)$ limit (reached for $Q_+ = Q_- \to 0$). When several distinct modes are contained within a bandwidth, for example with frequencies $f_1, f_2$ satisfying $Q_+^{-1} < \log \frac{f_2}{f_1} < Q_-^{-1}$, then Eq. (25) behaves as the intensity-weighted version of the arithmetic mean $\frac{f_1+f_2}{2}$. For indistinguishable modes such as $\log \frac{f_2}{f_1} < Q_+^{-1}$, the rate estimation is affected by the beating phenomenon. Alternative rate estimates can be developed based on geometric or harmonic averages; they all agree for a sufficiently large $Q_-$ (narrow band). Finally, the degree of squeezing can be expressed in terms of the coherence between the rhythmic signal and its derivative [42]. Such a measure of the accuracy of rate estimation could take over the role of the signal quality index used in practice.

We believe that the explicit balance between the dynamical uncertainty and the number of statistical degrees of freedom in this type of rate estimator makes it quite generic and flexible. It can introduce a statistical perspective into more sophisticated rate estimation methods [67, 68, 69, 70]. In the following, we apply it to the case of the heart rate in the most simple and direct way, which consists in fixing the band frequency $f$.

### 5.4. Statistical estimation of the heart rate signal from the QRS amplitude

We now compute the complex rate estimator Eq. (26) of the selected heart signal, the QRS amplitude $y(t)$ (extracted from the ECG) plotted in Fig. 5 (A top). In this statistical version of the rate estimator, the squeezing phenomenon is also observed at the cardiac frequency in Fig. 5 (B). In a neighbouring range of the order of the bandwidth $Q^{-1}$, the wavelet transform has a nearly constant angle, as can be observed in Fig. 4 (E, F), so that the real part of the complex rate estimator is also nearly constant in this band. We use this flexibility to fix the frequency variable in this range at $f_1 = 1 \text{ s}^{-1}$ for any time, as a first estimate of the heart rate.

The pair of quality factors remains to be determined: we want a stable estimator ($Q_- < Q_+$) as close as possible to achieve a time resolution of one heartbeat ($Q_+ \sim 1$). The presence of a second harmonic mode in addition to the fundamental mode also implies the constraint $Q_+^{-1} < \log 2$, or $Q_+ > 1.5$ using Eq. (20), to avoid its interference (beating). The bandwidth should be chosen wide enough to contain the heart rate variations, but narrow enough ($Q_- \sim 1$) to separate the fundamental mode from other influences (such as the next partial). Such conflicting requirements cannot be satisfied simultaneously when the heart rate displays large frequency modulations. This limitation suggests a refinement of this direct estimation, from fixed $f$ to variable $f(t)$ as for a ridge detection approach [71, 67]. Remarkably, such curve extraction was not necessary in this application at a low quality factor.

Although the choice of the quality factors is already quite limited for this application, we propose to optimize their values by minimizing the mean squared relative error of the estimated heart rate as compared to the reference discrete estimator. In fact, there are two ways of doing this: either we
<figure><img src="image_6.png" /><figcaption>Figure 5: Estimation of the heart rate from the real part of the complex rate Eq. (26) of the heart signal $y(t)$ (QRS amplitude, top). (A) Rate signals $F_j(t) = \Re\{F_y(t, f_j; Q_+, Q_-)\}$ for fixed frequencies $f_j=1,2,...$, see colour bar, with $(Q_+, Q_-) = (2.3, 1.6)$. (B) Continuous heart rate estimation for fixed $f = 1.0 \text{ s}^{-1}$ (green line) and discrete reference estimation from RR intervals (dots). (C) Heart rate estimates and hypnogram (bottom); the time interval selected in (A, B) and in Fig. 4 corresponds to the grey box.</figcaption></figure>

interpolate the reference discrete heart rate as a continuous variable, or we discretize the continuous heart rate signal at discrete times (midpoints of the RR intervals). In both cases, the error is about the same: 3% on average over the selected time interval in Fig. 5 this signal), or even 1% if we exclude the time intervals for which the reference estimation is considered unreliable (using NN intervals). The error is almost minimal for all these different methods for $Q_+ = 2.3$ and $Q_- = 1.6$. A noticeable difference is the observation of a bias (about -0.2%) with the discrete method as compared to no bias with the continuous method. This is caused by a slight oscillation in the heart rate signal at the cardiac frequency, interpreted as a residual beat from unresolved higher harmonic orders. This phenomenon is amplified for smaller quality factors, especially for the broadband marginal considered in [44]. The frequency variable $f_1$ and the band ($f_{qrs}, Q_{qrs}$) defining the QRS amplitude signal from the ECG recording can also be optimally selected. The error is very sensitive to $f_1$, nearly optimal at $1.0 \text{ s}^{-1}$, but less sensitive to the QRS band, fixed to $f_{qrs} = 14 \text{ s}^{-1}$ and $Q_{qrs} = 3$ as in Fig. 4 (C, F).

In Fig. 5 (C) we can see how the heart rate variability evolves with the sleep stages. Non-REM sleep is characterized by a noticeably lower rate, with particularly low variability during deep sleep (N3). The heart rate and its variability increase during REM sleep, with bursts during brief awakenings. Cardiac activity, as measured by heart rate variability, is compared with the respiratory activity in the next section.

In conclusion, rate estimation over several atoms degrades the time-frequency resolution, but a single atom is not sufficient for stability: this problem has also been addressed in [72, 70], where only
n = 2 orthogonal tapers are used, and the averaging is performed over numerous representations each based on a random linear combination of these two atoms. Here, our continuous approach reduces the effective number of atoms to an optimum $n = \frac{Q_{\pm}}{Q_{-}} \approx \sqrt{2}$. As an extension of this approach, the number n could be made adaptive, increasing as the regularity of the rhythm is lost. Finally, an alternative statistical approach based on maximum likelihood estimation of smoothly deformed Gaussian stationary processes has been developed in [73, 74, 75]. From an iterative procedure, the fluctuation rate of a rhythmic signal can also be obtained in this way.

# 6. Coherence and mutual information of cardio-respiratory signals

Cardio-respiratory coupling is one of the most important factors in the complex dynamics of the cardiovascular system and has aroused great interest of both clinicians and theoreticians. The interaction of the heart with other organs ensures the normal functioning of the body [76], while a deterioration of this coordination may indicate health impairment [77]. Some authors have even argued that the irregularity of breathing and the nonlinear nature of heart-respiration coupling could be potential sources of temporal variability in the heart rate [78, 79]. This interaction results in a frequency modulation of the heart rate at the breathing frequency, one of the component of the HRV known as respiratory sinus arrhythmia (RSA) [80]. The distinction between cardio-respiratory synchronization and coordination has been the subject of a debate [81]. The complexity of heart rate variability and cardio-respiratory coupling has previously been estimated from RR intervals using correlation dimension, Lyapunov exponents [82, 83], spectral coherence [84, 85, 86] or information-theoretic measures [87]. We investigate here the time-frequency coherence of cardio-respiratory signals and we to relate its statistical significance to the mutual information and its time-frequency distribution. We elaborate on the adequate choice of the quality factors for a localized and significant analysis of the cardio-respiratory interaction.

## 6.1. Time-frequency coherence: generic definition and single trial estimation

Our aim here is to characterize the common information shared between two distinct and simultaneous signals. Sitting somewhere between pairwise correlations and mutual information, we revisit the concept of time-frequency coherence [88, 22, 89, 90], an approach that offers both dynamical and statistical insights.

Given two processes y, z, an evolutionary cross-spectrum $S_{yz}(t, f)$ can be defined as in Eqs. (8) or (15), once a time-frequency decomposition is specified. This cross-spectrum can be normalized using the auto-spectra, resulting in a dimensionless coefficient called the time-frequency coherence between the processes [91, 92]:

$$ \gamma_{yz}(t, f) = \frac{S_{yz}(t, f)}{\sqrt{S_{yy}(t, f)S_{zz}^*(t, f)}} \quad (27) $$

Like the cross-spectrum, the coherence has a Hermitian symmetry: $\gamma_{zy} = \gamma_{yz}^*$. It characterizes the similarities between two signals over time and frequency, in terms of the strength (modulus) and phase difference (angle) of their linear relationship. Coherence reaches a maximum modulus of 1 in the case of a perfectly linear relationship. In other words, it is a (complex) correlation coefficient between the statistics of the oscillating atoms at the location $(t, f)$ in each signal.

The above definition of the time-frequency coherence is generic, and applies to various definitions of the evolutionary spectra, such as the spectrogram or the scalogram as statistical expectations, or as single trial estimators. When applied to the spectral, temporal or global marginals of the evolutionary spectrum, namely the power spectral density, the time-dependent or the global variance (power),
Eq. (27) reduces to the spectral coherence, the time-dependent or the global Pearson correlation coefficient, respectively.

The estimation of the time-frequency spectra and coherence should require several statistical samples, such as independent and identically distributed realizations of the signal pairs. However, single trial experiments provide unique realizations of the processes (which cannot be identically repeated). In this case, nearby atoms from the time-frequency decomposition of the signals are used as distinct statistical samples. Such an estimation comes at the cost of a loss of time-frequency resolution: around the location $(t, f)$, the original atomic uncertainty is extended to a continuous composition of $n$ atoms, as defined in Eqs. (2) or (5). In the case of the wavelet-based estimator, $n = \frac{Q_+}{Q_-}$, as represented by the Wigner-Ville distribution of the estimation kernel, Fig. 2.

The frequency marginal is equivalent to assuming stationarity, while the time marginal is equivalent to self-similarity (the correspondence between these concepts is discussed in [51]).

The number of atoms plays a crucial role in the interpretation of the coherence value: the larger the number $n$, the better the statistical estimate of the coherence, but the lower the time-frequency resolution and the stronger the assumption of local stationarity or self-similarity. The limit case of the coherence with the maximum number of statistical degrees of freedom and no resolution is the (global) Pearson correlation coefficient. Conversely, no statistical averaging corresponds to the limit case $n = 1$ with maximum resolution (first described in [88]): whatever the pair of signals, they always look perfectly linearly related (unit modulus), which is always true atom-wise. The interference of several oscillating atoms is therefore crucial to obtain a statistical estimate of the coherence [22]. The case $n = 0$ would correspond to an ideal time-frequency localization based directly on the Wigner-Ville distribution of the signals. Its lack of non-negativity leads to pathological properties [91, 93], consequences of the unspecified shape of the uncertainty atom or kernel, see Eq. (4).

The trade-off between dynamical uncertainty and statistical significance is a crucial aspect of the single trial estimation of the time-frequency coherence. The control of the significance of the estimate as $n$ is varied is examined next.

## 6.2. Spurious coherence and statistical significance

To assess the significance of a time-frequency coherence estimate when statistical samples are limited, we propose to compare the value of the coherence modulus with what would be obtained from completely independent (i.e. incoherent) signals. The null hypothesis can be constructed from the residual value of the coherence estimated between two independent stochastic trajectories. We characterize the so-called spurious coherence from extensive simulations of the wavelet-based canonical estimator of the time-frequency coherence in Fig. 6. A simple but very useful model of spurious coherence is provided by two independent stationary self-similar Gaussian processes. In particular, the spurious coherence shown in Fig. 6 (A) is computed from two pink noises.

In the case of Gaussian stationary processes, the spurious (spectral) coherence squared $|\gamma_{sp}|^2$ is known to follow the Goodman distribution [94, 95, 96]. For independent Gaussian stationary processes, this probability distribution simplifies to a certain beta law with a single parameter $B(1, \beta)$:

$$ \mathbb{P}[|\gamma_{sp}|^2 \le \gamma^2] = 1 - (1 - \gamma^2)^\beta, \qquad (28) $$

which can be easily transformed into a significant coherence threshold $\gamma^2(p) = 1 - p^{\frac{1}{\beta}}$ given the complementary probability $p = \mathbb{P}[|\gamma_{sp}|^2 > \gamma^2]$ or $p$-value. Fig. 6 (B) shows that this distribution accurately matches the empirical distribution of the spurious coherence squared $|\gamma_{sp}|^2 = |\gamma_{yz}(t, f; Q_+, Q_-)|^2$ obtained from the simulation (A).

The number $\beta + 1$ is interpreted as the number of (effective and complex) degrees of freedom of the statistical distribution, and it is identified in [96] with Eq. (2), which we call the (effective)
<figure><img src="image_7.png" /><figcaption>Figure 6: Statistics of the spurious coherence for the canonical wavelet estimator, obtained from a realization of a pair of independent self-similar stationary real Gaussian processes. (A) Spurious coherence squared $|\gamma_{sp}|^2 = |\gamma_{yz}(t, f; Q_+, Q_-)|^2$ with quality factors $Q_+ = 20$ and $Q_- = 2$, computed from pink noises $y, z$ of length $N = 2^{19}$ (top, arbitrary unit). To avoid edge effects, the outer region (faint colours) at a distance from the boundaries less than twice the resolution $\Delta t = f^{-1}Q_+$, $\Delta \log f = Q_-^{-1}$ is excluded from further calculations. (B) Cumulative distribution of $|\gamma_{sp}|^2$ (black solid line) compared to the beta distribution $B(1, \beta)$ with the same mean (red dashed line), aligned with the colour bar. Inset: quantile-quantile plot showing close agreement. (C) Parameter difference $\beta(n) - n$ (symbols for two different estimators, lines for linear fits with respect to $n^{-1}$), which seems to converge to $-\frac{1}{3}$ as $n$ increases, despite increasing uncertainty (95% confidence interval, filled area).</figcaption></figure>

number of atoms $n$. For the canonical scalogram and the resulting coherence estimator, the number of atoms corresponds exactly to the ratio of the quality factors, $n = \frac{Q_+}{Q_-}$. Assuming a beta law for $|\gamma_{sp}|^2$, we point out in Fig. 6 (C) a small difference (less than one) between the actual number of atoms $n$ and the estimated statistical degrees of freedom $\beta + 1$. By varying $n$ in the spurious coherence simulation, we estimate the parameter $\beta$ from the mean of the beta-distributed squared spurious coherence $\langle |\gamma_{sp}|^2 \rangle = \frac{1}{\beta+1}$, averaged over the time-frequency domain:

$$ \langle |\gamma_{yz}|^2 \rangle = \|\Omega\|^{-1} \iint_{\Omega} |\gamma_{yz}(t, f; Q_+, Q_-)|^2 dt df . \qquad (29) $$

Here, $\Omega$ is the inner region in Fig. 6 (A), free of boundary effects, with area $\|\Omega\| = \iint_{\Omega} dt df = a - b(1 + \log \frac{a}{b})$, $a = \frac{N}{2}e^{-\frac{k}{Q_-}}$, $b = kQ_+e^{\frac{k}{Q_-}}$, where $N = \frac{T}{dt}$ is the length of the time series, and $k$ is the distance from the boundaries (in resolution units) that is excluded from the inner region $\Omega$. Even a very long single simulation of the spurious coherence is not sufficient to obtain an accurate estimate of $\beta$, especially for a high value of $n = \frac{Q_+}{Q_-}$. Therefore, we average Eq. (29) over numerous simulations $i = 1, 2, \dots$, so that the cumulative averaged area $\sum_i \|\Omega_i\|$, exceeds $2^{21} \times n$ for each value of $n$.
The difference $\beta - n$ is plotted in Fig. 6 (C) (black plus markers) and linearly fitted with respect to $n^{-1}$ (black line). The parameter $\beta(n)$ grows from zero (for $n = 1$) to an asymptotic value compatible with $n - \frac{1}{3} - \frac{1}{2n}$ for large $n$. We have tested the robustness of this asymptotic behaviour by varying different parameters, such as the Hurst exponent of the independent signals, the quality factors $(Q_+, Q_-)$ or numerical parameters, without observing discrepancies beyond the uncertainty. We have tested white and pink noises (independent self-similar stationary Gaussian processes), with Hurst exponents $-\frac{1}{2}$ and 0, respectively. Only the ratio of the quality factors matters, their product (form factor of the time-frequency kernel) does not change the estimate. Any location $(t, f)$ in the inner time-frequency region $\Omega$ (away from the boundaries) yields the same estimate: the spurious coherence is homogeneously distributed. A distance of twice the unit resolution, $k = 2$, is sufficient to separate the excluded outer region from $\Omega$. We could decrease the sampling along log-frequency down to two samples per resolution unit without observing alteration of the estimation.

To assess the relevance and the accuracy of assuming a beta law $B(1, \beta)$ for the spurious coherence squared, we also estimate the parameter $\beta$ using a different estimator. As described in the next section, it exploits the relationship between beta and exponential random variables. The result is marked with blue circles in Fig. 6 (C) and shows an asymptotic behaviour compatible with $\beta \sim n - \frac{1}{3} - \frac{1}{4n}$. Small values ($n < 1.5$) have been discarded due to numerical imprecision. The term $O(n^{-1})$ differs from the first estimate: this discrepancy between the two estimates of $\beta$ highlights the inaccuracy of the beta law to model the distribution of the spurious coherence squared for small values of $n$.

In practice, the case $n = 10$ in Fig. 6 (A, B) yields $\beta = 9.618 \pm 0.004$ (black) and $9.647 \pm 0.005$ (blue), or a very small difference in the significance test. The beta law is not the exact distribution of the spurious coherence squared, but it may be considered as a very good approximation for a sufficiently large $n$. Eq. (28) provides a simple way to assess the significance of the time-frequency coherence estimator $\gamma_{yz}$, as the p-value $p = (1 - |\gamma_{yz}|^2)^\beta$.

### 6.3. Mutual information and its time-frequency distribution

In the case of two Gaussian processes $y$ and $z$, the modulus of their coherence gives the time-frequency distribution of their mutual information in the following form:

$$i_{yz}(t, f) = -\frac{1}{2} \log (1 - |\gamma_{yz}(t, f)|^2) \quad . \qquad (30)$$

The relationship between coherence and the rate of mutual information, $\dot{I}_{yz} = \int i_{yz}(f)df$, is indeed proved in the case of stationary processes [97]. In the simplest case (for white noises), the mutual information of two Gaussian random variables is obtained directly by applying Eq. (30) to the Pearson correlation coefficient (as a non-localized version of the coherence).

This provides a natural non-stationary estimator of the mutual information between two Gaussian processes $y$ and $z$, possibly restricted to a time-frequency domain $\Omega$:

$$I_{yz}[\Omega] = \iint_{\Omega} i_{yz}(t, f; Q_{+}, Q_{-}) dtdf \quad , \qquad (31)$$

in the wavelet paradigm (and similarly in the window paradigm). Mutual information, defined by the natural logarithm, is measured in natural units (nat). As for the coherence, in situations with limited time-frequency atoms ($n$ small), it is necessary to assess the significance and the spurious expected value of this statistical estimator.

A beta law $B(1, \beta)$ for the spurious coherence squared $|\gamma_{sp}|^2$ implies an exponential law $\text{Exp}(2\beta)$ for the spurious density of mutual information, expected for independent Gaussian processes:

$$\langle i_{\text{sp}} \rangle = \frac{1}{2\beta} \quad , \quad \langle I_{\text{sp}}[\Omega] \rangle = \frac{\|\Omega\|}{2\beta} \quad . \qquad (32)$$
This provides the alternative way of estimating the parameter $\beta$ in Fig. 6 (C, blue circles), a small correction away from the number of atoms $n$. Furthermore, it allows us to interpret the magnitude of significance of the estimated coherence as the ratio of the estimated and spurious density of mutual information: $-\log p = i_{yz}/\langle i_{sp} \rangle$. Evenly sampled, this quantity is useful to colour significant values of a coherence estimate (see Fig. 7).

The difference $I_{yz} - \langle I_{sp} \rangle$ can be used as a “corrected” estimator of the mutual information [98], which corresponds to the “debiased” coherence estimator $((|\gamma_{yz}|^2 - \langle |\gamma_{sp}|^2 \rangle)/(1 - \langle |\gamma_{sp}|^2 \rangle))$ [99]. Here, $\langle |\gamma_{sp}|^2 \rangle = \frac{1}{\beta+1}$, but the direct use of $\frac{1}{n}$ (without correction) was shown to be accurate for an effective number of atoms as low as $n = 8$ [99].

In the next section, these time-frequency coherence and mutual information estimators based on Gaussian processes prove useful in practice. We apply them to physiological signals that are assumed to be locally stationary over $Q_+$ oscillations and locally self-similar over a relative bandwidth $Q_-^{-1}$. Due to the non-stationarity of these signals, we believe that this time-frequency estimator of the mutual information based on the Gaussian hypothesis may be more informative than alternative non-parametric estimators based on a stronger assumption of stationarity.

<figure><img src="image_8.png" /><figcaption>Figure 7: Coherence and mutual information between heart rate and breathing effort for the sleeping subject shown in Fig. 3. (A) Signals of the heart rate $y$ (red) and the respiratory effort $z$ (blue), with arbitrary units. (B) Amplitude spectra $y$ and $z$, with frequency in beats per minute (min$^{-1}$). (C) Time-frequency coherence $\gamma_{yz}(t, f; Q_+, Q_-)$ for $Q_+ = 7$ and $Q_- = 1$. Colour coding: the cyclic hue corresponds to the phase, and the 4 levels of saturation are delimited by 3 significance thresholds for the modulus: $p = 10^{-a}$, $a = 1, 2, 3$. Regions affected by boundary effects, recording saturation or a higher than average error (3%) in the heart rate estimation are shaded. Hypnogram (black stairs). (D) Spectral density of mutual information (nat·min) and expected spurious value (black dashed line). (E) Zoom of (A, C) at the end of sleep (grey box in A). (F) Rate (temporal density) of mutual information (nat·min$^{-1}$), with expected spurious value (black dashed line). Significance levels (grey) are the same as in (C).</figcaption></figure>
### 6.4. Different regimes of cardio-respiratory coherence

To study the modulation of the heart rate at the breathing frequency, called respiratory sinus arrhythmia (RSA), and other forms of cardio-respiratory coherence, we apply the above coherence and mutual information estimators to the heart rate and respiratory signals. The former is the rate estimate previously extracted from the QRS amplitude in the ECG, and the latter is a breathing effort signal defined as the sum of recordings from thoracic and abdominal belt transducers. The sign of the respiratory oscillation is set to positive for inhalation. The recordings are taken from the polysomnography of the same subject as in Fig. 3. The pair of signals, together with their time-frequency coherence and information estimates, are shown in Fig. 7.

The oscillatory content of the heart rate signal $y(t)$ (red) and of the respiratory signal $z(t)$ (blue), Fig. 7 (A), is summarized in their respective amplitude spectra (B). Quite regular, the respiratory cycle is concentrated in a relatively narrow frequency band around $10 \text{ min}^{-1}$. In contrast, the heart rate variability covers a wider band centred around $2-3 \text{ min}^{-1}$. The spectrum of the respiratory recording is attenuated at frequencies below 1 oscillation per minute, as a result of instrumental filtering. There is no such high-pass filtering for the heart rate signal, which is extracted from modulations in the ECG recording. Instead, most of the amplitude at the cardiac (carrier) frequency $(60 \text{ min}^{-1})$ has disappeared from the estimated heart rate signal. Note that the slow modulations of the respiratory frequency and amplitude can also be estimated as a complex rate signal. Cardio-respiratory coherence with both types of respiratory signals, oscillations or fluctuating rate, have been analysed extensively for all subjects in the shhs2 dataset in [42]. Rate signals allow the study of very low modulation frequencies. However, that would come at the cost of losing amplitude and changing the phase at the carrier (breathing) frequency, thus affecting the observation of the RSA coherence.

The time-frequency coherence $\gamma_{yz}$ between these signals (A) is then estimated in Fig. 7 (C, zoom in E), with quality factors $Q_+ = 7$ and $Q_- = 1$. Significance thresholds are colour-coded as changes in saturation, in order to visualize only the relevant (density of mutual) information. We can distinguish three bands (timescales) of coherent cardio-respiratory activity, as summarized in the spectral marginal density (time integral) of the mutual information (D). Overall, the most significant is the intermediate band close to and above the respiratory frequency: this is the RSA. Patches of coherence are observed in a lower and wider band, associated to slow intermittent and coordinated cardio-respiratory events. Finally, a cardiac mode (at $60 \text{ min}^{-1}$ and above) indicates that the respiratory recording can capture heart pulses. This coherent mode is absent when the reference (interpolated discrete) estimator of the heart rate is used, which does not contain residual cardiac oscillations (present in the continuous estimator).

The coherence of the respiratory (RSA) mode often reaches high levels of significance, $p < 10^{-3}$, as shown in Fig. 7 (C, E). This is the case not only around the fundamental respiratory frequency $(10 \text{ min}^{-1})$ but also beyond $20 \text{ min}^{-1}$, implying that the specific shape of the respiratory oscillation (carried by higher order partials, here indistinguishable with $Q_- = 1$) is coherent in both signals. Part of the intermittent loss of RSA coherence coincides with regions potentially affected by a poor heart rate estimation (shaded), either caused by motion artefacts or by strong heart rate variability. This coherence is otherwise quite persistent.

The phase difference between the RSA mode in the heart rate and the oscillation of the breathing effort can vary significantly. Mostly in-phase (angle of coherence 0 in green: the heart rate increases when breathing in), the RSA mode can be lagged of $-\frac{\pi}{3}$ (yellow); this tends to happen during deep sleep (N3 stage, see the hypnogram). In contrast, the RSA phase tends to increase during rapid-eye-movement (REM) stage, occasionally reaching positive phase quadrature $(+\frac{\pi}{2}$ in blue) or even phase opposition $(\pi$ in magenta). The role of sleep stages in RSA changes would require further statistical confirmation across many subjects. However, these observations are consistent with a weaker control
of the cardio-respiratory inter-regulation during REM, than during non-REM sleep.

The spectral range is sufficient to observe longer timescales of cardio-respiratory coherence, inter-
mittently reaching high levels of significance ($p < 10^{-3}$) in a broad band below 5 min$^{-1}$. The angle
of coherence facilitate the distinction of two types of coherent cardio-respiratory events in this band,
that were previously hard to separate in a low and a very low-frequency band [86].

Micro-awakenings are easily detected in the respiratory signal, as a deep inhalation, followed by
a long exhalation (or more). The empty lungs then need several breathing cycles to refill, while the
heart rate remains high to pump the oxygenated blood, hence the phase opposition (π in magenta) in
this slow mode. Since the rather singular events that generate it are averaged over $Q_+ = 7$ oscillations,
they can interfere with (or even dominate) less intense coherent events. Lower quality factors would
be better suited to focus on these broadband time-localized structures.

Coherent events of comparatively lower amplitude, such as hyperpnea and hypnopnea, are associ-
ated with patches of coherence close to negative phase quadrature ($-\frac{\pi}{2}$ to $-\frac{2\pi}{3}$ in orange-red). They
disappear during the long episode of deep (N3) sleep, and may be associated with a decrease in RSA
coherence. A strong and persistent occurrence of this coherent mode, as a rhythm near 1 to 2 min$^{-1}$,
is associated to severe sleep apnea [44], hence called the apneic mode.

This time-frequency map of the cardio-respiratory coherence allows a concise visualization of local
and significant statistical similarities between physiological fluctuations and oscillations, together with
their phase. It should be remembered that the time-frequency plane has a hyperbolic metric in the
wavelet paradigm, so that its flat projection as an image should be interpreted with caution. In fact,
low frequency atoms look more important than their actual spread, and vice versa at high frequencies,
as schematized in Fig. 1 (B, right). We can think of it as a perspective view: close-by (low frequency)
wavelets look larger than the (high frequency) ones near the horizon, as can be seen in the spurious
coherence Fig. 6 (A).

6.5. *Quality factors and trade-off between localization and significance of mutual information*

The pictures of the significant estimated coherence in Fig. 7 (C, E) combine the angle of coherence to the empirical time-frequency distribution of mutual information $I_{yz}(t, f; Q_+, Q_-)$. In Fig. 7 (D and F), we plot the spectral and temporal marginals of the mutual information distribution $I_{yz}$, by integrating it over time or frequency respectively. The corresponding significance thresholds are shown with grey levels. From the visible region in (D), we can estimate a corrected value of the total mutual information $I_{yz}[\Omega] - \langle I_{sp}[\Omega] \rangle = 4.2$ knat, within a range $||\Omega|| = 2f_{max}T = 3.9 \times 10^3$ (the factor 2 accounts for negative frequencies), or on average just above 1 nat per unit area.

This mutual information is divided in three visible frequency bands as follows: 2.0 knat in the
cardiac band, 2.0 knat in the respiratory/RSA band, 0.2 knat in the low frequency band. These
bands are bounded by the frequencies $44 \text{ min}^{-1}$ and $5.5 \text{ min}^{-1}$ of two minima close to the spurious
value (black dashed line) of the spectral marginal of mutual information Fig. 7 (D). Note that the
high-frequency (cardiac) band, accumulates much more mutual information than the low-frequency
band, despite a lower spectral significance (D is not a log-frequency density).

For this reason, the estimated mutual information rate (temporal marginal) in Fig. 7 (F) is
dominated by the high-frequency coherence, even though the corresponding modulations have a very
small amplitude in the signals. The expected spurious value (black dashed line) correctly describes
the baseline mutual information rate. However, its value and significance is very sensitive to the
choice of $f_{max}$, which here includes the cardiac band.

Again, we try to refine the choice of quality factors. Essentially, we want to maximize the informa-
tion extracted from the coherence of the signals in the selected time interval by tuning the parameters
of our local stationarity and self-similarity assumptions. In practice, it is quite difficult to achieve
for several reasons. First, we have to deal with the limitations of the information estimators and the
resulting artefacts, namely: boundary effects and the imprecision of the statistical model for small
$n$. The significance decreases near the boundary in a way that is not modelled here, and is therefore
excluded from the inner region of integration $\Omega$. The upper frequency limit $f_{\max} < \frac{1}{2}f_s e^{-k/Q_{-}}$ is
particularly critical. The corrected estimate of the mutual information $I_{yz}[\Omega] - \langle I_{\text{sp}}[\Omega] \rangle$ increases
rapidly when $n$ approaches 1, as a consequence of the imprecision of the beta and exponential models
of significance (Fig. 6 C). We thus consider only situations with at least $n = \frac{Q_{+}}{Q_{-}} \ge 2$ atoms.

The corrected estimate of the mutual information has no clear global maximum in the valid
parameter space: for $f_{\max} = 100 \text{ min}^{-1}$, the estimate is nearly optimal for $Q_{-} \approx 1$, but keeps
increasing for $Q_{+} \to Q_{-}$. Although corrected, the mutual information has very little significance in
this limit, since the dynamical (time-frequency) uncertainty is minimal. The opposite situation occurs
when maximizing the mean magnitude of significance $\langle -\log p \rangle = I_{yz}[\Omega]/\langle I_{\text{sp}}[\Omega] \rangle$: the dynamical
uncertainty is also maximized, $n \to +\infty$. An optimal bandwidth emerges from these observations,
but the number of atoms varies widely depending on our preference for a high corrected value or a
high significance of the mutual information.

We measure a trend $\langle i_{yz} \rangle = I_{yz}[\Omega]/\|\Omega\| \approx \alpha/\log n$, which holds for an optimal value $\alpha = 0.36$ over
a wide range $n \in [3.5, 70]$, for $f_{\max} = 100 \text{ min}^{-1}$ and $Q_{-} \approx 1$. We learn that the mutual information
is here inversely proportional to the loss of localization information. As $n$ decreases, the spurious
contribution increases. In general, $\alpha$ varies with the quality factors and $f_{\max}$; we believe that these
variations are informative about the statistical and dynamical properties of the different coherent
components in the analysed signals.

The low quality factor $Q_{-} = 1$ provides an almost optimal bandwidth in the spectral range of
Fig. 7, sufficient to distinguish the cardiac, respiratory and apneic activity bands. The upper quality
factor $Q_{+} = 7$ represents a compromise between a precise localization and a high level of significance,
highlighting the richness of the interplay between the cardio-respiratory signals.

Thanks to a careful choice of the quality factors, the time-frequency coherence and the mutual
information provide local dynamical insights with a high statistical significance on different forms of
cardio-respiratory coherence. The mutual information and its directed extension, the transfer entropy,
have also been estimated on short timescales for physiological point processes [98]: the control of their
spurious statistics is not straightforward and has been managed from surrogate data. We therefore
believe that the detailed discussion of localization and significance in this type of analysis will be
a valuable aid for future investigations. Mutual information between heart rate and respiration is
also used in other contexts, such as emotional characterization [100], which could constitute another
immediate application of this work.

# 7. Conclusion

From the interplay between respiratory oscillations and heart rhythm fluctuations, to the multi-scale stochastic and intermittent waves of a sleeping brain, a single recording of biophysical observables can contain a wealth of information about the underlying living system dynamics. However, the practical exploitation of this information is hampered by the inadequacy of standard dynamical assumptions: a stationary, self-similar or scale-separated approach wastes some of the richness of the recorded information. The parsimony of our assumptions is limited only by the uncertainty principle, expressed here in the time-frequency domain.

In the context of a single trial experiment, we have shown how the uncertainty principle binds
the localization of the analysis to the significance of the resulting statistical estimate. On Heisenberg's lower bound, signals are continuously decomposed into individual time-frequency atoms, but
the significance of the statistical approach is brought to the brink of collapse. Time-frequency de-
composition is a fertile ground for recomposing uncertainty atoms into flexible quadratic statistical
estimators, whose localization lies away from the bound.

More specifically, the numbering of uncertainty atoms (eigenstates) composing a quadratic statisti-
cal estimator is given a continuous time-frequency interpretation from the Wigner-Ville representation
of the operator. This general property has been put into practice. We have proposed a canonical
implementation of these generic concepts in a concise analytical framework, based on the choice of a
geometric paradigm and two essential physical parameters.

First, a canonical spectrogram estimator was built from oscillating normal windows of Gabor. The translation invariance of the atom in the time-frequency plane was associated with a flat (Euclidean) representation with fixed timescales, one for the decomposition, and a different one for the recompo- sition into a statistical estimate. Their ratio gives the exact number of uncertainty atoms that make up the statistical estimate, while matching the classical expression of the uncertainty principle with its product of resolutions.

We then constructed the scale-invariant analogue, the canonical scalogram estimator, from the
log-normal wavelet of Altes and Grossmann. The translation and scaling of the atom gives the time-
frequency plane with a hyperbolic geometry (hence the delicacy of its representation as a flat image);
quality factors are fixed parameters in this paradigm, rather than timescales. This leads to a scale-free
formulation of a scale separation and of the uncertainty principle.

The EEG of a sleeping subject provides a practical discussion of the two geometric paradigms
and the choice of the localization parameters. The spectrogram approach can be convenient for
studying a neural phenomenon around a fixed timescale that separates well its correlation time from
its evolution, whereas the scalogram allows a broader and clearer study of narrowband (rhythmic)
and broadband (noisy) components across scales.

In the application to heart rhythm, we have developed a continuous statistical estimator of the
heart rate from the ECG recording by translating the standard procedure into the log-normal wavelet
framework. We have also related the quality factor to the highest distinguishable harmonic order in
the time-frequency decomposition of a rhythmic signal.

We have measured the coherence between heart rate modulations and respiratory oscillations
over time and scales, and we have interpreted its significance as the ratio of the estimated mutual
information to its spurious value. The spurious statistics of the coherence and mutual information
estimators is precisely characterized in terms of the number of atoms, and the relation between mutual
information and localization information has been investigated.

This tour of physiological signals completes the picture of the uncertainty trade-off linking the
significance of a single trial estimate and its time-frequency localization. Few physical parameters are
thus left to tune our dynamical and statistical assumptions about the living system, as applied here
to the complexity of polysomnographic information. The proposed canonical framework provides a
practical summary of the fundamental concepts at stake in single trial time-frequency statistics, and
offers a generic and explicit way of dealing with rich but uncertain dynamics.

Acknowledgements

We are grateful to Alain Arneodo for stimulating discussions. AG acknowledges LOMA and MPI-PKS for their hospitality and EURLight S&T for funding. This work was supported by Agence Nationale de la Recherche (ANR) through contract: ANR-18-CE45-0012-01.
References

[1] M. Costa, A. L. Goldberger, C.-K. Peng, Multiscale Entropy Analysis of Complex Physiologic Time Series, Physical Review Letters 89 (6) (2002) 068102. doi:10/dwvzzn.

[2] J. Zschocke, R. P. Bartsch, M. Glos, T. Penzel, R. Mikolajczyk, J. W. Kantelhardt, Long- and short-term fluctuations compared for several organ systems across sleep stages, Frontiers in Network Physiology 2 (2022). doi:10/gsntpq.

[3] A. Bashan, R. P. Bartsch, Jan. W. Kantelhardt, S. Havlin, P. C. Ivanov, Network physiology reveals relations between network topology and physiological function, Nature Communications 3 (1) (2012) 702. doi:10/f4f8v7.

[4] L. Faes, D. Marinazzo, F. Jurysta, G. Nollo, Linear and non-linear brain-heart and brain-brain interactions during sleep, Physiological Measurement 36 (4) (2015) 683. doi:10/grs9hk.

[5] G. D'Agostino, A. Scala (Eds.), Networks of Networks: The Last Frontier of Complexity, 1st Edition, Understanding Complex Systems, Springer International Publishing : Imprint: Springer, Cham, 2014. doi:10/krsz.

[6] R. P. Bartsch, K. K. Liu, A. Bashan, P. C. Ivanov, Network physiology: How organ systems dynamically interact, PloS one 10 (11) (2015) e0142143. doi:10/f8b4vm.

[7] P. Dayan, L. F. Abbott, Theoretical Neuroscience: Computational and Mathematical Modeling of Neural Systems, Computational Neuroscience, MIT Press, Cambridge (Mass.) London, 2001.

[8] F. Battiston, G. Cencetti, I. Iacopini, V. Latora, M. Lucas, A. Patania, J.-G. Young, G. Petri, Networks beyond pairwise interactions: Structure and dynamics, Physics Reports 874 (2020) 1-92. doi:10/ghqt4w.

[9] R. Carmona, W.-L. Hwang, B. Torrésani, Practical Time-Frequency Analysis: Gabor and Wavelet Transforms, with an Implementation in S, Academic Press, 1998.

[10] B. Ricaud, B. Torrésani, A survey of uncertainty principles and some signal processing applications, Advances in Computational Mathematics 40 (3) (2014) 629-650. doi:10/gsvpdj.

[11] A. Arneodo, E. Barcy, F. Argoul, J. Elezgaray, J. F. Muzy, Ondelettes, Multifractales et Turbulences (de l'ADN Aux Crossances Cristallines), Diderot éditeur, arts et sciences, Paris, 1995.

[12] A. Arneodo, B. Audit, E. Bacry, S. Manneville, JF. Muzy, SG. Roux, Scale Invariance and Beyond: What Can We Learn from Wavelet Analysis?, in: Scale Invariance and Beyond, Springer, 1997, pp. 37-51.

[13] E. Gerasimova, B. Audit, S. G. Roux, A. Khalil, O. Gileva, F. Argoul, O. Naimark, A. Arneodo, Wavelet-based multifractal analysis of dynamic infrared thermograms to assist in early breast cancer diagnosis, Frontiers in Physiology 5 (May 2014). doi:10/gjvnrh.

[14] G.-Q. Zhang, L. Cui, R. Mueller, S. Tao, M. Kim, M. Rueschman, S. Mariani, D. Mobley, S. Redline, The National Sleep Research Resource: Towards a sleep data commons, Journal of the American Medical Informatics Association 25 (10) (2018) 1351-1358. doi:10/gdntnm.
[15] S. F. Quan, B. V. Howard, C. Iber, J. P. Kiley, F. J. Nieto, G. T. O'Connor, D. M. Rapoport, S. Redline, J. Robbins, J. M. Samet, The sleep heart health study: Design, rationale, and methods, Sleep 20 (12) (1997) 1077-1085.

[16] D. Gabor, Theory of communication. Part 1: The analysis of information, Journal of the Institution of Electrical Engineers - Part III: Radio and Communication Engineering 93 (26) (1946) 429-441. doi:10/2rz.

[17] D. B. Percival, A. T. Walden, Spectral Analysis for Physical Applications: Multitaper and Conventional Univariate Techniques, Cambridge University Press, Cambridge ; New York, NY, USA, 1993.

[18] M. Bayram, R. Baraniuk, Multiple window time-frequency analysis, in: Proceedings of Third International Symposium on Time-Frequency and Time-Scale Analysis (TFTS-96), IEEE, Paris, France, 1996, pp. 173–176. doi:10/dzsqxh.

[19] A. T. Walden, A unified view of multitaper multivariate spectral estimation, Biometrika 87 (4) (2000) 767–788. doi:10/d8dcbd.

[20] J.-M. Combes, A. Grossmann, P. Tchamitchian, Wavelets: Time-Frequency Methods and Phase Space Proceedings of the International Conference, Marseille, France, December 14–18, 1987, Inverse Problems and Theoretical Imaging, Springer, Berlin Heidelberg, 1989.

[21] R. Bardenet, A. Hardy, Time-frequency transforms of white noises and Gaussian analytic functions, Applied and Computational Harmonic Analysis 50 (2021) 73–104. doi:10/gj8mts.

[22] C. Torrence, G. P. Compo, A Practical Guide to Wavelet Analysis, Bulletin of the American Meteorological Society 79 (1) (1998) 61–78. doi:10/bhbwhf.

[23] B. Cazelles, M. Chavez, A. J. McMichael, S. Hales, Nonstationary influence of El Nino on the synchronous dengue epidemics in Thailand, PLoS medicine 2 (4) (2005) e106. doi:10/ftw9z8.

[24] E. A. K. Cohen, A. T. Walden, A Statistical Study of Temporally Smoothed Wavelet Coherence, IEEE Transactions on Signal Processing 58 (6) (2010) 2964–2973. doi:10/fcngnq.

[25] L. Stanković, S. Stanković, Wigner distribution of noisy signals, IEEE Transactions on Signal Processing 41 (2) (1993) 956–960. doi:10/brfrr8.

[26] J. W. Pitton, The statistics of time-frequency analysis, Journal of the Franklin Institute 337 (4) (2000) 379–388. doi:10/cpzg52.

[27] V. Sucic, N. Saulig, B. Boashash, Estimating the number of components of a multicomponent nonstationary signal using the short-term time-frequency Rényi entropy, EURASIP Journal on Advances in Signal Processing 2011 (1) (2011) 125. doi:10/dj9v3b.

[28] D. L. Jones, T. W. Parks, A high resolution data-adaptive time-frequency representation, IEEE Transactions on Acoustics, Speech, and Signal Processing 38 (12) (1990) 2127–2135. doi:10/bq9xjv.

[29] W. J. Williams, M. L. Brown, A. O. Hero III, Uncertainty, information, and time-frequency distributions, in: F. T. Luk (Ed.), San Diego, '91, San Diego, CA, San Diego, CA, 1991, pp. 144–156. doi:10/dt8pp9.
[30] P. Flandrin, R. Baraniuk, O. Michel, Time-frequency complexity and information, in: Proceedings of ICASSP '94. IEEE International Conference on Acoustics, Speech and Signal Processing, Vol. iii, IEEE, Adelaide, SA, Australia, 1994, pp. III/329–III/332. doi:10/bqmd3g.

[31] L. Stanković, A measure of some time-frequency distributions concentration, Signal Processing 81 (3) (2001) 621–631. doi:10/dmvq77.

[32] R. G. Baraniuk, P. Flandrin, A. J. Janssen, O. J. Michel, Measuring time-frequency information content using the Rényi entropies, IEEE Transactions on Information theory 47 (4) (2001) 1391–1409. doi:10/dcsvwh.

[33] A. Papandreou-Suppappola, F. Hlawatsch, G. F. Boudreaux-Bartels, Quadratic time-frequency representations with scale covariance and generalized time-shift covariance: A unified framework for the affine, hyperbolic, and power classes, Digital Signal Processing 8 (1) (1998) 3–48. doi:10/d2twxd.

[34] R. Kronland-Martinet, J. Morlet, A. Grossmann, Analysis of sound patterns through wavelet transforms, International Journal of Pattern Recognition and Artificial Intelligence 01 (02) (1987) 273–302. doi:10/bfs8dp.

[35] R. Stockwell, L. Mansinha, R. Lowe, Localization of the complex spectrum: The S transform, IEEE Transactions on Signal Processing 44 (4) (1996) 998–1001. doi:10/dj3sx5.

[36] A. Grossmann, J. Morlet, Decomposition of Hardy Functions into Square Integrable Wavelets of Constant Shape, SIAM Journal on Mathematical Analysis 15 (4) (1984) 723–736. doi:10/b9ddn3.

[37] R. A. Altes, Sonar for generalized target description and its similarity to animal echolocation systems, The Journal of the Acoustical Society of America 59 (1) (1976) 97–105. doi:10/d2c3bb.

[38] J. M. Lilly, S. C. Olhede, Generalized Morse Wavelets as a Superfamily of Analytic Wavelets, IEEE Transactions on Signal Processing 60 (11) (2012) 6036–6041. doi:10/gjvnrq.

[39] H. Knutsson, C.-F. Westin, G. Granlund, Local multiscale frequency and bandwidth estimation, in: Proceedings of 1st International Conference on Image Processing, Vol. 1, IEEE Comput. Soc. Press, Austin, TX, USA, 1994, pp. 36–40. doi:10/fj4x5k.

[40] D. J. Field, Relations between the statistics of natural images and the response properties of cortical cells, Journal of the Optical Society of America A 4 (12) (1987) 2379. doi:10/fr7v9c.

[41] P. Kovesi, Invariant measures of image features from phase information, Ph.D. thesis, University of Western Australia Crawley, Australia (1996).

[42] A. Guillet, Revisiting statistical physics tools for the complex time-frequency characterization of physiological signals, These de doctorat, Bordeaux (May 2022).

[43] R. A. Brown, R. Frayne, A fast discrete S-transform for biomedical signal processing, in: 2008 30th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, IEEE, Vancouver, BC, 2008, pp. 2586–2589. doi:10/cf7pkc.

[44] A. Guillet, A. Arneodo, F. Argoul, Tracking Rhythms Coherence From Polysomnographic Records: A Time-Frequency Approach, Frontiers in Applied Mathematics and Statistics 7 (2021) 624456. doi:10/gkdfvx.
[45] I. Selesnick, R. Baraniuk, N. Kingsbury, The dual-tree complex wavelet transform, IEEE Signal Processing Magazine 22 (6) (2005) 123-151. doi:10/fwhmtz.

[46] P. C. Ivanov, R. P. Bartsch, Network Physiology: Mapping Interactions Between Networks of Physiologic Networks, in: G. D'Agostino, A. Scala (Eds.), Networks of Networks: The Last Frontier of Complexity, Understanding Complex Systems, Springer International Publishing, Cham, 2014, pp. 203-222. doi:10/krst.

[47] S. Chokroverty (Ed.), Sleep Disorders Medicine: Basic Science, Technical Considerations and Clinical Aspects, Springer New York, New York, NY, 2017. doi:10/kzrq.

[48] M. J. Prerau, R. E. Brown, M. T. Bianchi, J. M. Ellenbogen, P. L. Purdon, Sleep Neurophysiological Dynamics Through the Lens of Multitaper Spectral Analysis, Physiology 32 (1) (2017) 60-92. doi:10/f9gvdh.

[49] P. Kovesi, Good Colour Maps: How to Design Them (Sep. 2015). arXiv:1509.03700.

[50] G. Buzsáki, A. Draguhn, Neuronal oscillations in cortical networks, science 304 (5679) (2004) 1926-1929. doi:10/bwtf9q.

[51] P. Flandrin, P. Borgnat, P.-O. Amblard, From Stationarity to Self-similarity, and Back: Variations on the Lamperti Transformation, in: R. Beig, B. G. Englert, U. Frisch, P. Hänggi, K. Hepp, W. Hillebrandt, D. Imboden, R. L. Jaffe, R. Lipowsky, H. v. Löhneysen, I. Ojima, D. Sornette, S. Theisen, W. Weise, J. Wess, J. Zittartz, G. Rangarajan, M. Ding (Eds.), Processes with Long-Range Correlations, Vol. 621, Springer Berlin Heidelberg, Berlin, Heidelberg, 2003, pp. 88-117. doi:10/ctsgz8.

[52] R. Pool, Is It Healthy to Be Chaotic?, Science 243 (4891) (1989) 604-607. doi:10/c3bdpm.

[53] E. Tobaldini, L. Nobili, S. Strada, K. R. Casali, A. Braghiroli, N. Montano, Heart rate variability in normal and pathological sleep, Frontiers in Physiology 4 (2013). doi:10/gk42cc.

[54] E. J. C. de Geus, P. J. Gianaros, R. C. Brindle, J. R. Jennings, G. G. Berntson, Should heart rate variability be "corrected" for heart rate? Biological, quantitative, and interpretive considerations, Psychophysiology 56 (2) (2019) e13287. doi:10/gfhfcc.

[55] J. Pan, W. J. Tompkins, A Real-Time QRS Detection Algorithm, IEEE Transactions on Biomedical Engineering BME-32 (3) (1985) 230-236. doi:10/ddj8r2.

[56] A. N. Vest, G. D. Poian, Q. Li, C. Liu, S. Nemati, A. Shah, G. D. Clifford, Cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox: PhysioNet-Cardiovascular-Signal-Toolbox 1.0.2, Zenodo (Sep. 2019). doi:10/krsx.

[57] J. Behar, A. Johnson, G. D. Clifford, J. Oster, A comparison of single channel fetal ECG extraction methods, Annals of biomedical engineering 42 (6) (2014) 1340-1353. doi:10/f524wq.

[58] A. E. Johnson, J. Behar, F. Andreotti, G. D. Clifford, J. Oster, R-peak estimation using multimodal lead switching, in: Computing in Cardiology 2014, IEEE, 2014, pp. 281-284.

[59] A. N. Vest, G. Da Poian, Q. Li, C. Liu, S. Nemati, A. J. Shah, G. D. Clifford, An open source benchmarked toolbox for cardiovascular waveform and interval analysis, Physiological Measurement 39 (10) (2018) 105004. doi:10/gmdh98.
[60] T.-P. Le, P. Argoul, Continuous wavelet transform for modal identification using free decay response, Journal of Sound and Vibration 277 (1-2) (2004) 73–100. doi:10/dkdr47.

[61] F. Topsøe, Some bounds for the logarithmic function, Inequality theory and applications 4 (2007) 137.

[62] A. Guillet, A. Arneodo‡, P. Argoul, F. Argoul, Quantifying the Rationality of Rhythmic Signals, in: P. Flandrin, S. Jaffard, T. Paul, B. Torresani (Eds.), Theoretical Physics, Wavelets, Analysis, Genomics: An Indisciplinary Tribute to Alex Grossmann, Applied and Numerical Harmonic Analysis, Springer International Publishing, Cham, 2023, pp. 565–593. doi:10/krsv.

[63] F. Auger, P. Flandrin, Improving the readability of time-frequency and time-scale representations by the reassignment method, IEEE Transactions on Signal Processing 43 (5) (1995) 1068–1089. doi:10/dpmrhp.

[64] I. Daubechies, J. Lu, H.-T. Wu, Synchrosqueezed wavelet transforms: An empirical mode decomposition-like tool, Applied and Computational Harmonic Analysis 30 (2) (2011) 243–261. doi:10/dnnqvt.

[65] G. Thakur, E. Brevdo, N. S. Fučkar, H.-T. Wu, The Synchrosqueezing algorithm for time-varying spectral analysis: Robustness properties and new paleoclimate applications, Signal Processing 93 (5) (2013) 1079–1094. doi:10/f4v6wh.

[66] J. Xiao, P. Flandrin, Multitaper Time-Frequency Reassignment for Nonstationary Spectrum Estimation and Chirp Enhancement, IEEE Transactions on Signal Processing 55 (6) (2007) 2851–2860. doi:10/chkh7j.

[67] D. Iatsenko, P. McClintock, A. Stefanovska, Extraction of instantaneous frequencies from ridges in time-frequency representations of signals, Signal Processing 125 (2016) 290–303. doi:10/gfzgx2.

[68] D. Fourer, F. Auger, K. Czarnecki, S. Meignen, P. Flandrin, Chirp Rate and Instantaneous Frequency Estimation: Application to Recursive Vertical Synchrosqueezing, IEEE Signal Processing Letters 24 (11) (2017) 1724–1728. doi:10/ggnrtb.

[69] C.-Y. Lin, L. Su, H.-T. Wu, Wave-shape function analysis, Journal of Fourier Analysis and Applications 24 (2) (2018) 451–505. doi:10/gc8dvx.

[70] H.-T. Wu, E. Z. Soliman, A new approach for analysis of heart rate variability and QT variability in long-term ECG recording, BioMedical Engineering OnLine 17 (1) (2018) 54. doi:10/gjvpbj.

[71] R. Carmona, W. Hwang, B. Torrésani, Characterization of signals by the ridges of their wavelet transforms, IEEE Transactions on Signal Processing 45 (10) (1997) 2586–2590. doi:10/cdnsc8.

[72] I. Daubechies, Y. Wang, H.-t. Wu, ConceFT: Concentration of frequency and time via a multi-tapered synchrosqueezed transform, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 374 (2065) (2016) 20150193. doi:10/gmv3x6.

[73] H. Omer, B. Torrésani, Time-frequency and time-scale analysis of deformed stationary processes, with application to non-stationary sound modeling, Applied and Computational Harmonic Analysis 43 (1) (2017) 1–22. doi:10/f97d6p.
[74] A. Meynard, B. Torrésani, Spectral Analysis for Nonstationary Audio, IEEE/ACM Transactions on Audio, Speech, and Language Processing 26 (12) (2018) 2371–2380. doi:10/gh2r36.

[75] A. Meynard, B. Torrésani, Time-Scale Synthesis for Locally Stationary Signals, in: ICASSP 2020 - 2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2020, pp. 5820–5824. doi:10/gsv5gh.

[76] P. C. Ivanov, L. A. N. Amaral, A. L. Goldberger, S. Havlin, M. G. Rosenblum, Z. R. Struzik, H. E. Stanley, Multifractality in human heartbeat dynamics, Nature 399 (6735) (1999) 461–465. doi:10/bk6r24.

[77] V. Ponomarenko, M. Prokhorov, A. Karavaev, A. Kiselev, V. Gridnev, B. Bezruchko, Syn-chronization of low-frequency oscillations in the cardiovascular system: Application to medical diagnostics and treatment, The European Physical Journal Special Topics 222 (10) (2013) 2687–2696. doi:10/f5ffpf.

[78] N. Wessel, M. Riedl, J. Kurths, Is the normal heart rate “chaotic” due to respiration?, Chaos: An Interdisciplinary Journal of Nonlinear Science 19 (2) (2009) 028508. doi:10/cgshg3.

[79] Y. Shiogai, A. Stefanovska, P. V. E. McClintock, Nonlinear dynamics of cardiovascular ageing, Physics Reports 488 (2) (2010) 51–110. doi:10/dmwhvb.

[80] J. A. Hirsch, B. Bishop, Respiratory sinus arrhythmia in humans: How breathing pattern modulates heart rate, American Journal of Physiology-Heart and Circulatory Physiology 241 (4) (1981) H620–H629. doi:10/gjvpcd.

[81] H. Krause, J. F. Kraemer, T. Penzel, J. Kurths, N. Wessel, On the difference of cardiorespiratory synchronisation and coordination, Chaos: An Interdisciplinary Journal of Nonlinear Science 27 (9) (2017) 093933. doi:10/gbz8v7.

[82] T. Henriques, M. Ribeiro, A. Teixeira, L. Castro, L. Antunes, C. Costa-Santos, Nonlinear Methods Most Applied to Heart-Rate Time Series: A Review, Entropy 22 (3) (2020) 309. doi:10/ggtqpg.

[83] Y. M. Ishbulatov, T. S. Bibicheva, V. I. Gridnev, M. D. Prokhorov, M. V. Ogneva, A. R. Kiselev, A. S. Karavaev, Contribution of Cardiorespiratory Coupling to the Irregular Dynamics of the Human Cardiovascular System, Mathematics 10 (7) (2022) 1088. doi:10/gswzfs.

[84] J. P. Saul, R. D. Berger, M. H. Chen, R. J. Cohen, Transfer function analysis of autonomic regulation. II. Respiratory sinus arrhythmia, American Journal of Physiology-Heart and Circulatory Physiology 256 (1) (1989) H153–H161. doi:10/gnm4tp.

[85] L. Faes, G. Pinna, A. Porta, R. Maestri, G. Nollo, Surrogate Data Analysis for Assessing the Significance of the Coherence Function, IEEE Transactions on Biomedical Engineering 51 (7) (2004) 1156–1166. doi:10/dtrzcj.

[86] R. J. Thomas, J. E. Mietus, C.-K. Peng, A. L. Goldberger, An electrocardiogram-based technique to assess cardiopulmonary coupling during sleep, Sleep 28 (9) (2005) 1151–1161. doi:10/gf8nss.
[87] M. Javorka, J. Krohova, B. Czippelova, Z. Turianikova, Z. Lazarova, R. Wiszt, L. Faes, Towards understanding the complexity of cardiovascular oscillations: Insights from information theory, Computers in Biology and Medicine 98 (2018) 48–57. doi:10/gdv4dd.

[88] P. C. Liu, Wavelet Spectrum Analysis and Ocean Wind Waves, in: Wavelet Analysis and Its Applications, Vol. 4, Elsevier, 1994, pp. 151–166. doi:10/krsw.

[89] K. Gurley, T. Kijewski, A. Kareem, First- and Higher-Order Correlation Detection Using Wavelet Transforms, Journal of Engineering Mechanics 129 (2) (2003) 188–201. doi:10/crqrgt.

[90] A. Grinsted, J. C. Moore, S. Jevrejeva, Application of the cross wavelet transform and wavelet coherence to geophysical time series, Nonlinear Processes in Geophysics 11 (5/6) (2004) 561–566. doi:10/d6dhzv.

[91] L. White, B. Boashash, Cross spectral analysis of nonstationary processes, IEEE Transactions on Information Theory 36 (4) (1990) 830–835. doi:10/bhxr6c.

[92] G. Matz, F. Hlawatsch, Time-frequency coherence analysis of nonstationary random processes, in: Proceedings of the Tenth IEEE Workshop on Statistical Signal and Array Processing (Cat. No. 00TH8496), IEEE, 2000, pp. 554–558.

[93] M. Orini, R. Bailon, L. T. Mainardi, P. Laguna, P. Flandrin, Characterization of Dynamic Interactions Between Cardiovascular Signals by Time-Frequency Coherence, IEEE Transactions on Biomedical Engineering 59 (3) (2012) 663–673. doi:10/b4rwmn.

[94] N. R. Goodman, Statistical Analysis Based on a Certain Multivariate Complex Gaussian Distribution (An Introduction), The Annals of Mathematical Statistics 34 (1) (1963) 152–177. doi:10/cgm7kq.

[95] G. Carter, A. Nuttall, Statistics of the estimate of coherence, Proceedings of the IEEE 60 (4) (1972) 465–466. doi:10/ddpd4f.

[96] E. Cohen, A. Walden, A Statistical Analysis of Morse Wavelet Coherence, IEEE Transactions on Signal Processing 58 (3) (2010) 980–989. doi:10/bt572x.

[97] I. M. Gel'fand, A. M. Yaglom, Computation of the amount of information about a stochastic function contained in another such function, Uspekhi Matematicheskikh Nauk 12 (1) (1957) 3–52.

[98] G. Mijatovic, R. Pernice, A. Perinelli, Y. Antonacci, A. Busacca, M. Javorka, L. Ricci, L. Faes, Measuring the Rate of Information Exchange in Point-Process Data With Application to Cardiovascular Variability, Frontiers in Network Physiology 1 (2022) 765332. doi:10/gsmrmc.

[99] A. T. Walden, E. A. K. Cohen, Statistical Properties for Coherence Estimators From Evolutionary Spectra, IEEE Transactions on Signal Processing 60 (9) (2012) 4586–4597. doi:10/gjvpmj.

[100] M. T. Valderas, J. Bolea, P. Laguna, Fellow IEEE, R. Bailón, M. Vallverdú, Mutual information between heart rate variability and respiration for emotion characterization, Physiological Measurement 40 (8) (2019) 084001. doi:10/gktkg3.
# Appendix A. Linear representations: list of relations

The signal and spectrum representations, as related by the Fourier transform, are considered here as limit cases of time-frequency representations, that are purely localized in one domain and delocalized in the other. We recall their definition here to clarify the convention used, and we provide explicit relations between limit and intermediate cases (as defined in Eqs. 14 and 7) with finite quality factor or timescale.

## Appendix A.1. Signal and spectrum conventions

The Fourier transform and its inverse relate the signal to the spectrum, as follows:

$$x(t) = \int \hat{x}(f)e^{i2\pi ft}df \quad ; \quad \hat{x}(f) = \int x(t)e^{-i2\pi ft}dt \quad . \qquad (A.1)$$

The wavelet transform $X(t, f; Q)$ as defined in Eq. (14) is a band-pass filtering, i.e. a spectral decomposition of the signal $x(t)$ (with the same physical unit) into different bands $(f, Q)$. In contrast, the Gabor transform $G_x(t, f; \tau)$, Eq. (7), is a Fourier transform under a sliding window, i.e. a temporal decomposition of the spectrum $\hat{x}(f)$ in different time periods $(t, \tau)$. These signal-like or spectrum-like behaviours are best noticed by recovering the actual signal and spectrum in the limit cases:

$$x(t) = \lim_{Q \to 0} X(t, f; Q) + X(t, -f; Q) = \lim_{\tau \to 0} \tau^{-1} G_x(t, f; \tau) e^{i2\pi ft} \qquad (A.2)$$

$$\hat{x}(f) = \lim_{Q \to +\infty} Q|f|^{-1} X(t, f; Q) e^{-i2\pi ft} = \lim_{\tau \to +\infty} G_x(t, f; \tau) \quad . \qquad (A.3)$$

The Gabor transform can thus be converted into a spectral decomposition of the signal:

$$X(t, f; \tau) = \tau^{-1} G_x(t, f; \tau) e^{i2\pi ft} \quad . \qquad (A.4)$$

Conversely, the wavelet transform becomes a temporal decomposition of the spectrum by reversing this equality with $\tau = Q|f|^{-1}$. The S-transform [35], which would be denoted here $Q|f|G_x(t, f; Q|f|^{-1})$, is a mixture of these two conventions with the phase of the spectrum and the unit of the signal. Based on a Gabor-Morlet wavelet (a scaled normal window), the S-transform behaves as $X(t, f; Q)e^{-i2\pi ft}$ in the log-normal wavelet framework.

## Appendix A.2. Change of localization

The localization parameter of a time-frequency decomposition can be decreased or increased as follows:

$$X(t, f; Q_{-}) = \frac{Q_{+}Q_{-}}{Q_{-}} \int X(t, f'; Q_{+}) \hat{\psi}_{Q_{-}}^{*}(f'/f) \frac{df'}{|f'|} , \qquad Q_{-}^{-2} = Q_{-}^{-2} - Q_{+}^{-2} \qquad (A.5)$$

$$X(t, f; Q_{+}) = \int X(t', f; Q_{-}) \psi_{Q_{+}}^{*}(f(t' - t)) |f| dt' , \qquad Q_{+}^{2} = Q_{+}^{2} - Q_{-}^{2} \qquad (A.6)$$

$$G_x(t, f; \tau_{-}) = \int G_x(t, f'; \tau_{+}) w_{\frac{1}{\tau_{-}}} (f' - f) \tau_{-} e^{i2\pi(f'-f)t} dt' , \qquad \tau_{-}^{-2} = \tau_{-}^{-2} - \tau_{+}^{-2} \qquad (A.7)$$

$$G_x(t, f; \tau_{+}) = \frac{\tau_{+}}{\tau_{-}\tau_{+}} \int G_x(t', f; \tau_{-}) w_{\tau_{+}}(t' - t) dt' , \qquad \tau_{+}^{2} = \tau_{+}^{2} - \tau_{-}^{2} \quad . \qquad (A.8)$$

The signal and the spectrum are also recovered by integration over one domain:

$$x(t) = Q \int X(t, f; Q) \frac{df}{|f|} = \int G_x(t, f; \tau) e^{i2\pi ft} df \qquad (A.9)$$

$$\hat{x}(f) = \int X(t, f; Q) e^{-i2\pi ft} dt = \frac{1}{\tau} \int G_x(t, f; \tau) dt \quad . \qquad (A.10)$$
# Appendix B. Quadratic representations in time and/or in frequency

The relationships between time, frequency and time-frequency kernel representations of a quadratic operator are given here, using the canonical spectrogram estimator as an example.

## Appendix B.1. Wigner-Ville representation as a time-frequency distribution

The Wigner-Ville distribution defined in Eq. (3) can be applied to the temporal or spectral kernel of a quadratic operator, to obtain its time-frequency representation:

$$K(t, f) = \int \check{K}(t + \frac{\tau}{2}, t - \frac{\tau}{2})e^{-i2\pi f\tau}d\tau = \int \hat{K}(f + \frac{\eta}{2}, f - \frac{\eta}{2})e^{i2\pi\eta t}d\eta \quad . \quad (B.1)$$

The temporal and the spectral representations (related by Fourier transforms) of the quadratic operator are obtained from the time-frequency kernel as follows:

$$\check{K}(t_1, t_2) = \int K(\frac{t_1+t_2}{2}, f)e^{i2\pi f(t_1-t_2)}df \quad ; \quad \hat{K}(f_1, f_2) = \int K(t, \frac{f_1+f_2}{2})e^{-i2\pi(f_1-f_2)t}dt \quad . \quad (B.2)$$

## Appendix B.2. Kernel of the canonical spectrogram estimator

By denoting the Gaussian function $\phi(u) = e^{-\pi u^2}$ which satisfies $\hat{\phi} = \phi$, we can obtain symmetric expressions for the operator $K_w$ represented in the time and frequency domains:

$$\check{K}_w(t_1, t_2; t, f) = \frac{\sqrt{2}}{\tau_+} \left( \phi\left(\frac{t_1-t_2}{\tau_-}\right) \phi\left(\frac{t_1+t_2-2t}{\tau_+}\right) \right)^{\frac{1}{2}} e^{i2\pi f(t_1-t_2)} \quad (B.3)$$

$$\hat{K}_w(f_1, f_2; t, f) = \frac{\sqrt{2}}{\eta_+} \left( \phi\left(\frac{f_1-f_2}{\eta_-}\right) \phi\left(\frac{f_1+f_2-2f}{\eta_+}\right) \right)^{\frac{1}{2}} e^{-i2\pi(f_1-f_2)t}, \quad \eta_{\pm} = \tau_{\mp}^{-1} \quad . \quad (B.4)$$

The time-frequency representation of this operator is tractable:

$$K_w(t', f'; t, f) = 2 \left( \phi\left(\frac{t'-t}{\tau_+}\right) \phi\left(\frac{f'-f}{\eta_+}\right) \right)^2, \quad (B.5)$$

it covers $n = \frac{\tau_+}{\tau_-} = \tau_+\eta_+ = \Delta t \Delta f$ atoms and its trace is Tr $K_w = 1$.

# Appendix C. Laplace's method for an unsolvable number of atoms

Given the operator $K$ of a quadratic estimator $xKx^*$, the associated effective number of atoms (or eigenstates) $n$, as defined in Eq. (2), is not always as easy to compute as for the canonical spectrogram ($n = \frac{\tau_+}{\tau_-}$) and canonical scalogram ($n = \frac{Q_+}{Q_-}$) estimators. The discussion in section 3.3 motivates the definition of an alternative scalogram estimator whose numerical implementation is more convenient. The trace formula for its number of atom is intractable, but it can be approximated with Laplace's method.

## Appendix C.1. Definition of the practical time-smoothing scalogram estimator

The practical scalogram estimator $S_{xx}(t, f; Q, m)|f|$, introduced in [44], is based on a time-smoothing over $m$ times the wavelet duration. It is defined from the product of the wavelet decomposition $XX^*$ similarly to the canonical spectrogram estimator Eq. (9):

$$S_{xx}(t, f; Q, m) = \frac{2}{m} \int X(t', f; Q) X^*(t', f; Q) w_{\tau}^2(t' - t) dt', \quad \tau = mQ|f|^{-1} \quad . \quad (C.1)$$
This hybrid estimator is characterized by a kernel that includes both a log-normal wavelet with quality factor $Q$ and the normal window of duration $\tau = mQ|f|^{-1}$. The kernel of the practical scalogram estimator $S_{xx}(t, f; Q, m)|f|$ takes the following form in the frequency domain:

$$ \hat{K}(f_1, f_2; f, t) = \sqrt{2Q} \hat{\psi}_Q(f_1/f) \hat{\psi}_Q^*(f_2/f) w_{\frac{\sqrt{2}|f|}{mQ}}(f_1 - f_2) e^{i2\pi(f_1 - f_2)t} \quad (C.2) $$

The corresponding number of atoms $n$, designed to be approximated by $m$, is defined in Eq. (2) from the ratio of the integral expressions:

$$ \left( \int \hat{K}(v, v) dv \right)^2 = 2 \left( Q \int |\hat{\psi}_Q(v)|^2 dv \right)^2 = \lambda^{\frac{1}{2}}, \quad \lambda = e^{\frac{1}{2\pi Q^2}} \quad (C.3) $$

$$ \iint |\hat{K}(v, v')|^2 dv dv' = 2Q^2 \iint e^{-2\pi Q^2 (\log v)^2 - 2\pi Q^2 (\log v')^2 - \pi(mQ)^2 (v-v')^2} \Theta(v) \Theta(v') dv dv' , \quad (C.4) $$

where we have set $v = f_1/f$, $v' = f_2/f$ and $t = 0$ for clarity, and without loss of generality. Since the second double integration is intractable, we compute its Laplace approximation.

## Appendix C.2. Laplace’s approximation

Laplace’s method aims to approximate such an integral of a positive complicated function by the integral of a Gaussian which has the same global maximum up to the second order (same peak value and curvature). Provided that the integrated function is doubly differentiable at a unique global maximum, the precision of the approximation depends on the existence of a sufficiently large exponent $\alpha$, which ensures that significant contributions to the integral come from the neighbourhood of the maximum value of the integrated function. We apply the bivariate form of this approximation:

$$ G(\alpha) = \iint e^{\pi \alpha g(x, y)} dxdy \approx \frac{2}{\alpha} \det[-\partial_x \partial_y g(x_0, y_0)]^{-\frac{1}{2}} e^{\pi \alpha g(x_0, y_0)}, \quad (C.5) $$

where $(x_0, y_0)$ is the position of the global maximum of the function $g(x, y)$, $\partial_x \partial_y g(x_0, y_0)$ is its Hessian matrix evaluated at the maximum, and det refers to its determinant.

Rewriting Eq. (C.4), we can identify the following quantities:

$$ \iint e^{\frac{1}{4\pi Q^2} - 2\pi Q^2 (\log v - \frac{1}{4\pi Q^2})^2 - 2\pi Q^2 (\log v' - \frac{1}{4\pi Q^2})^2 - \pi(mQ)^2 (e^{\log v} - e^{\log v'})^2} d\log v d\log v' = \lambda^{\frac{1}{2}} \beta^2 \iint e^{\pi \alpha g(x, y)} dxdy \\ \alpha = (mQ)^2 \lambda, \quad \beta = \frac{m}{\sqrt{2}} \lambda^{\frac{1}{2}}, \quad (x, y) = (\log(v, v') - \frac{1}{4\pi Q^2})/\beta, \quad g(x, y; \beta) = -x^2 - y^2 - (e^{\beta x} - e^{\beta y})^2. $$

The maximum is 0 at $(x_0, y_0) = (0, 0)$ and the Hessian matrix at the maximum is:

$$ -\partial_x \partial_y g(x_0, y_0; \beta) = \begin{bmatrix} 2 + 2\beta^2 & -2\beta^2 \\ -2\beta^2 & 2 + 2\beta^2 \end{bmatrix}, \quad (C.6) $$

with determinant $4 + 8\beta^2$. Therefore, we obtain:

$$ \iint |\hat{K}(v, v')|^2 dv dv' \approx \lambda^{\frac{1}{2}} \beta^2 G(\alpha; \beta) = \lambda^{\frac{1}{2}} (1 + m^2 \lambda)^{-\frac{1}{2}} , \quad (C.7) $$

so that Laplace’s approximation of the number of atoms in the time-smoothing scalogram estimator $S_{xx}(t, f; Q, m)$ is:

$$ n \approx (1 + m^2 \lambda)^{\frac{1}{2}} , \quad (C.8) $$

for a sufficiently large value of $\alpha = (mQ)^2\lambda$. Note the dependence of $n$ on the quality factor $Q$ (through $\lambda$).