Permutation Entropy: Too Complex a Measure for EEG Time Series?

Abstract

Permutation entropy (PeEn) is a complexity measure that originated from dynamical systems theory. Specifically engineered to be robustly applicable to real-world data, the quantity has since been utilised for a multitude of time series analysis tasks. In electroencephalogram (EEG) analysis, value changes of PeEn correlate with clinical observations, among them the onset of epileptic seizures or the loss of consciousness induced by anaesthetic agents. Regarding this field of application, the present work suggests a relation between PeEn-based complexity estimation and spectral methods of EEG analysis: for ordinal patterns of three consecutive samples, the PeEn of an epoch of EEG appears to approximate the centroid of its weighted power spectrum. To substantiate this proposition, a systematic approach based on redundancy reduction is introduced and applied to sleep and epileptic seizure EEG. The interrelation demonstrated may aid the interpretation of PeEn in EEG, and may increase its comparability with other techniques of EEG analysis.

1. Introduction

Electroencephalography (EEG) is a neurophysiological technique widely used in research and medical diagnostics. Regardless of its broad spectrum of applications, a common denominator of working with EEG is the assessment of motifs/patterns/features in the signals, and to relate them to observations from other modalities. Analyses are either performed by visual inspection or supported by computer-based parameter extraction. Ultimately, the purpose of EEG analysis is the inference on cortical state changes.

EEG signals are varying electrical potentials that are (in most cases) non-invasively measured across the scalp. These voltage fluctuations represent a very low-dimensional projection of the underlying cerebral processes, and their relation to higher cortical functions is not fully understood. Most techniques of quantitative EEG analysis are therefore intrinsically empirical. Be it the classic frequency band-related measures or more recent additions to the set of analysis tools: the mechanistic background of the vast majority of EEG parameters is not entirely known. In this inherently probabilistic framework, any given signal property that correlates with behavioural observations constitutes a parameter worth considering for EEG analysis.

One quantity that is receiving increasing attention in the field is permutation entropy (PeEn). Introduced by Christoph Bandt and Bernd Pompe as a complexity measure for time series, PeEn is the Shannon entropy of a sequence of ordinal patterns—the latter being discrete symbols that encode how consecutive time series entries relate to one another in terms of position and value [1]. PeEn has been successfully used for EEG analyses in sleep scoring, general anaesthesia monitoring and research on disorders of the central nervous system, most notably epilepsy. Comprehensive overviews regarding the manifold applications of PeEn can be found in [2,3].

New applications of PeEn in EEG processing, and especially ways of extending the method are constantly reported on. For instance, “weighted” [4], “multiscale” [5], “multivariate multi-scale” [6], and “amplitude-aware” [7] variants of PeEn have so far been proposed. As literature suggests, the quantity may well be relevant to any research field that benefits from a widely accepted measure of the complexity of EEG.

Considering the great interest taken in PeEn, and given its virtually universal applicability in EEG analysis, a particular set of questions is surprisingly seldom addressed: what are the dynamic signal characteristics of EEG that PeEn responds to? Can the parameters pattern order and time delay be selected on a phenomenological basis? How should the abstract notion of complexity be interpreted in EEG?

From the outcome-oriented perspective, these may be minor issues. Nevertheless, a careful examination of PeEn in EEG might foster our general understanding of this class of electrophysiological signals. The present work takes a step in this direction, pioneering the following approach: (1) instead of extending PeEn, try to simplify it without compromising its suitability for EEG analysis; and (2) if successful, examine whether the result can more easily be interpreted than the initial chaos-theoretic complexity measure.

Credit for inspiring this idea is due to Bandt [8], who recently questioned the usage of terms like “complexity”, “chaos” and “disorder” in the EEG context and suggested that the PeEn of an epoch of EEG is equivalent to its distance from white noise—a simplification that immediately increases interpretability. In the following article, we shall advance this strategy, systematically deriving an intuitive explanation for the behaviour of PeEn in EEG.

2. An Overview of Permutation Entropy

Both the terminology and mathematical formulation regarding PeEn may vary between authors. To avoid ambiguity, it seems reasonable to shortly recapitulate the basic principles underlying the quantity—using the definitions, notations and nomenclature this work will adhere to.

2.1. Sequences of Ordinal Patterns

The ordinal pattern of an m-tuple of real numbers $(x_1,x_2,\dots ,x_m)$ describes how its elements relate to one another in terms of position and value. For instance, the ordinal pattern of the tuple $(17,7,8)$ is fully specified by: “There are three elements, the first is the greatest, the second is the least.”

Apart from didactic purposes, formal representations of ordinal patterns are usually preferred. Let us hereto agree upon the following simple notation: each element $x_i$ of an m-tuple $(x_1,x_2,\dots ,x_m)$ shall be assigned a distinct index $k_i\in \{1,2,\dots ,m\}$, whereby the greatest element is assigned the index m, the second-greatest is assigned the index $m-1$, and so forth. The ordinal pattern of the aforementioned tuple $(17,7,8)$ thus corresponds to the indices $(k_1,k_2,k_3)=(3,1,2)$. We shall simply write $\mathbf{312}$ when referring to this specific pattern.

Only tuples of pairwise distinct elements can have an unambiguous ordinal pattern. Let us neglect this latent ambiguity, limiting our considerations to tuples $(x_1,x_2,\dots ,x_m)$ with $x_i=x_j$ if and only if $i=j$. For any such tuple, each permutation of its set of values $\left\{x_i\right\}$ results in a different ordinal pattern. Consequently, a set of $m!$ ordinal patterns ${\mathsf{\Omega }}_m=\{{\mathsf{\Pi }}_1,{\mathsf{\Pi }}_2,\dots ,{\mathsf{\Pi }}_{m!}\}$ exists for tuples of length m. The parameter $m\ge 2$ is called the order of the ordinal pattern.

Consider, for example, the tuple $(x_1,x_2,x_3)$ with $x_1\ne x_2\ne x_3$. Its $m!=6$ possible ordinal patterns of order $m=3$ shall be denoted ${\mathsf{\Omega }}_3=\{\mathbf{123},\mathbf{132},\mathbf{213},\mathbf{231},\mathbf{312},\mathbf{321}\}$. Figure 1 displays the shapes of these ordinal patterns.

Any sequence of real numbers $(x_1,x_2,\dots ,x_n)$ can be mapped onto a sequence of ordinal patterns. Given the order m, the particular pattern ${\pi }_t$ at time index t is obtained from the m-tuple of values

$$(x_t,x_{t+\tau },x_{t+2\tau },\dots ,x_{t+(m-1)\tau })\mapsto {\pi }_t\in {\mathsf{\Omega }}_m.$$

(1)

This encoding is frequently called the symbolisation technique and yields an ordinal pattern sequence of the form

$$({\pi }_1,{\pi }_2,\dots ,{\pi }_k)\mathrm{with}k=n-(m-1)\tau ,$$

(2)

where k is the number of patterns obtainable from a sequence of n samples. The parameter $\tau \ge 1$ represents the distance between samples compared and is usually named the time delay.

2.2. Estimating the Permutation Entropy of EEG

PeEn of order m and time delay $\tau$ is the Shannon entropy [9] of a discrete probability distribution of $m!$ ordinal patterns. In EEG analysis, the probability distribution vector

$$\mathbf{p}={\left(\begin{array}{cccc}{p}_1& p_2& \dots & p_{m!}\end{array}\right)}^T\mathrm{and}{\parallel \mathbf{p}\parallel }_1=1$$

(3)

is obtained by first mapping a sequence of EEG samples $(x_1,x_2,\dots ,x_n)$ onto its corresponding sequence of ordinal patterns $({\pi }_1,{\pi }_2,\dots ,{\pi }_k)$ in terms of Equations (1) and (2). The probability distribution $\mathbf{p}$ is then estimated by counting pattern occurrences. For a sequence of k elements, the probability $p_i$ of the pattern ${\mathsf{\Pi }}_i\in {\mathsf{\Omega }}_m$ is hence

$$p_i=\frac{1}{k}\sum _{t=1}^k\left[{\pi }_t={\mathsf{\Pi }}_i\right]\mathrm{for}1\le i\le m!=\left|{\mathsf{\Omega }}_m\right|.$$

(4)

The operator $\left[·\right]$ represents the Iverson bracket [10], a convenient means of notating anonymous indicator functions. Shannon entropy is then computed from the probability distribution $\mathbf{p}$ according to the fundamental formula

$$\mathrm{H}\left(\mathbf{p}\right)=\sum _{i=1}^{m!}{p}_{i}log\frac{1}{p_i}.$$

(5)

For any $p_i=0$, the corresponding summand is replaced by the value of its limit towards zero,

$$\underset{p_i\to 0}{lim}{p}_{i}log\frac{1}{p_i}=0.$$

(6)

PeEn thus converges for any given ordinal pattern distribution. The base of the logarithm can be chosen arbitrarily and determines the unit of measure. Entropy values reported in this work were divided by $logm!$ for normalisation—they are dimensionless quantities on the interval $[0,1]$.

3. Shannon Entropy and Pairwise Probability Balances

No matter the unit of measure: the absolute value of PeEn as obtained from a single epoch of EEG is insignificant. Conversely, the quantity is highly relevant for EEG analysis when comparing entropies across multiple EEG epochs. Considerable contrast in value can then be observed (and be correlated with behavioural observations as is common practice). Let us call these decisive inter-epoch fluctuations the dynamics of PeEn. Their relation with the signal characteristics of EEG are the main topic of the article at hand. To approach the subject, we shall consider some properties of Shannon entropy and derive implications for the dynamics of PeEn therefrom.

Pairwise Probability Balances

From a model-free perspective, Shannon entropy can be interpreted as a statistical parameter on a discrete probability distribution

$$\mathbf{p}={\left(\begin{array}{cccc}{p}_1& p_2& \dots & p_n\end{array}\right)}^T\mathrm{and}{\parallel \mathbf{p}\parallel }_1=1.$$

(7)

One of its core properties is that it increases with any pairwise modification of $\mathbf{p}$ towards the uniform distribution [9]. To utilise this behaviour, let us partition the n elements of $\mathbf{p}$ into a set of $n^2$ probability pairs

$$X=\left\{(p_i,p_j)\mid (i,j)\in \left\{1,\dots ,n\right\}\times \left\{1,\dots ,n\right\}\right\},$$

(8)

and let us further define for each such pair $(p_i,p_j)$ a probability balance coefficient

$${\beta }_{ij}=\left\{\begin{array}{cc}{p}_i/(p_i+p_j)& \mathrm{for}{p}_i+p_j>0,\hfill \\ 1/2& \mathrm{for}{p}_i+p_j=0.\hfill \end{array}\right.$$

(9)

It then holds that ${\beta }_{ij}=1/2$ if and only if $p_i=p_j$. In analogy to the steady continuation of Shannon entropy as per Equation (6), this justifies the stipulation ${\beta }_{ij}=1/2$ for $p_i+p_j=0$. Moreover, it holds that

$${\beta }_{ij}+{\beta }_{ji}=1.$$

(10)

Using its balance coefficient ${\beta }_{ij}$, the contribution of a particular probability $p_i$ to the overall value of Shannon entropy can be expressed as

$$p_{i}log\frac{1}{p_i}=p_{i}log\frac{1}{p_i+p_j}+p_{i}log\frac{p_i+p_j}{p_i}=p_{i}log\frac{1}{p_i+p_j}+p_{i}log\frac{1}{{\beta }_{ij}}.$$

(11)

The summed contribution of a probability pair $(p_i,p_j)$ then is

$$\begin{array}{cc}\hfill p_{i}log\frac{1}{p_i}+p_{j}log\frac{1}{p_j}& =\left(p_{i}log\frac{1}{p_i+p_j}+p_{i}log\frac{1}{{\beta }_{ij}}\right)+\left(p_{j}log\frac{1}{p_i+p_j}+p_{j}log\frac{1}{{\beta }_{ji}}\right)\hfill \\ \hfill & =(p_i+p_j)\left(log\frac{1}{p_i+p_j}+{\beta }_{ij}log\frac{1}{{\beta }_{ij}}+(1-{\beta }_{ij})log\frac{1}{1-{\beta }_{ij}}\right)\hfill \\ \hfill & =(p_i+p_j)\left(log\frac{1}{p_i+p_j}+{\mathrm{H}}_{\mathrm{b}}\left({\beta }_{ij}\right)\right).\hfill \end{array}$$

(12)

The right-hand side of this equation is an affine transformation of the binary entropy function

$${\mathrm{H}}_{\mathrm{b}}\left(x\right)=xlog\frac{1}{x}+(1-x)log\frac{1}{1-x}.$$

(13)

Introducing $p_{ij}=p_i+p_j$, Equation (12) can be reformulated as the bivariate function

$${\mathrm{H}}_{\Delta }(p_{ij},{\beta }_{ij})=p_{ij}log\frac{1}{p_{ij}}+p_{ij}{\mathrm{H}}_{\mathrm{b}}\left({\beta }_{ij}\right),$$

(14)

which behaves as depicted in Figure 2. In particular, if we fix the summed probability $p_{ij}$ at a constant value, ${\mathrm{H}}_{\Delta }$ is a concave function on the domain $0\le {\beta }_{ij}\le 1$. This function is symmetrical with regard to its maximum at ${\beta }_{ij}=1/2$, where $p_i=p_j$. Thus, if we reallocate among a pair $(p_i,p_j)$ their summed probability $p_{ij}$ such that $p_i$ and $p_j$ approach their average (while $p_{ij}$ remains constant), ${\beta }_{ij}$ approaches $1/2$ and Shannon entropy increases. The same principle applies to any pair of probabilities drifting apart, where the entropy decreases. Two well-known corner cases of this generic property are that (1) Shannon entropy is maximal for the uniform distribution, while (2) it is zero if all but one $p_i$ are zero [9].

For a probability distribution $\mathbf{p}$ of n elements, a quadratic matrix of $n^2$ balance coefficients

$$\left(\begin{array}{ccc}{\beta }_{11}& \cdots & {\beta }_{1n}\\ ⋮& \ddots & ⋮\\ {\beta }_{n1}& \cdots & {\beta }_{nn}\end{array}\right)$$

(15)

can be created. We limit our considerations to the entries above its main diagonal, namely the subset

$$B=\left\{{\beta }_{ij}|1\le i<j\le n\right\}\mathrm{such}\mathrm{that}\left|B\right|=\left(\genfrac{}{}{0pt}{}{n}{2}\right)=\frac{n^2+n}{2}.$$

(16)

The remaining coefficients, ${\beta }_{ij}$ with $i\ge j$, are easily obtained via the property of symmetry as given by Equation (10).

In conjunction with the relation ${\parallel \mathbf{p}\parallel }_1=1$, the coefficients ${\beta }_{ij}\in B$ constitute a system of $\left|B\right|+1$ linear equations in the probabilities $\mathbf{p}$. Its augmented system matrix is

$$\left(\begin{array}{ccccccc}1& 1& 1& 1& \dots & 1& 1\hfill \\ {\beta }_{12}-1& {\beta }_{12}& 0& 0& \dots & 0& 0\hfill \\ {\beta }_{13}-1& 0& {\beta }_{13}& 0& \dots & 0& 0\hfill \\ {\beta }_{14}-1& 0& 0& {\beta }_{14}& \dots & 0& 0\hfill \\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\hfill \\ 0& 0& 0& 0& \dots & {\beta }_{n-1,n}& 0\hfill \end{array}\right).$$

(17)

By Gaussian elimination, one confirms this system to be overdetermined but consistent. Any probability distribution vector $\mathbf{p}$ is hence unambiguously determined by its corresponding set of balance coefficients ${\beta }_{ij}\in B$. The same holds true for the value of Shannon entropy $\mathrm{H}\left(\mathbf{p}\right)$. Consequently, the dynamics of PeEn can be studied by analysing variations within the probability balances ${\beta }_{ij}\in B$.

4. Complexity and Pseudo-Complexity

PeEn depends on a set of parameters and surrounding conditions. Besides the order m and time delay $\tau$, also the sampling rate $f_{\mathrm{s}}$, the window length n as well as any filters applied to the EEG may affect the results. Those are a lot of parameters to decide upon, especially when considering the empirical nature of ordinal pattern-based EEG analysis. Existing guidelines on parameter selection are essentially based on computational or statistical feasibility concerns [3], and are rather not motivated by (electro-)physiological considerations. Therefore, any parametrisation of PeEn reported suitable for EEG analysis must be regarded as the result of extensive experimentation.

With that in mind, it is remarkable that independent groups have successfully used the six ordinal patterns of order $m=3$ for EEG encoding [11]. Most likely by experimental optimisation, investigators ended up using the lowest-possible order that still accounts for more than just the patterns “the signal increases” $\left(\mathbf{12}\right)$ versus “the signal decreases” $\left(\mathbf{21}\right)$. Paradoxically, if six ordinal patterns suffice for the analysis, most of what makes PeEn a complexity measure can apparently be omitted.

Accepting this immanent contradiction, one may raise the following questions: can the number of ordinal patterns be even further reduced without impairing the dynamics of PeEn? If so, do the remaining patterns permit a more tangible interpretation of this EEG analysis technique?

4.1. A New Class of Ordinal Patterns?

Regarding these questions, variegating the pattern order m obviously does not permit any further investigation. A finer-grained approach is necessary and shall here be introduced. It is based on the idea that any pattern pair may be the redundant bifurcation of a less specialised, yet more decisive kind of motif.

By way of illustration, let us define an augmented class of ordinal patterns of order $m=3$. For any tuple of distinct values $(x_1,x_2,x_3)$, our extra-ordinal patterns shall not only describe how its elements relate to one another, but additionally encode the value of the Boolean expression

$$(x_{\mathrm{hi}}-x_{\mathrm{mid}})\ge (x_{\mathrm{mid}}-x_{\mathrm{lo}}).$$

(18)

Figuratively, we assess if the centroid of a pattern lies above or below its equatorial axis. This is visualised in Figure 3.

Consequently, we shall call the respective pattern either top-heavy or bottom-heavy. As an example, the tuple $(23,11,14)$ features the bottom-heavy pattern $\mathbf{312}\mathbf{b}$, while the tuple $(28,6,17)$ maps onto a top-heavy $\mathbf{312}\mathbf{t}$. Both patterns are derivates of the same (conventional) ordinal pattern $\mathbf{312}$, though, and it holds for their probabilities that $p_{\mathbf{312}}=p_{\mathbf{312}\mathbf{t}}+p_{\mathbf{312}\mathbf{b}}$. Generally, a pair of extra-ordinal patterns ${\mathsf{\Pi }}_{\mathrm{t}i}$ and ${\mathsf{\Pi }}_{\mathrm{b}i}$ exists for each ordinal pattern ${\mathsf{\Pi }}_i$, and their probabilities relate to each other as in

$$p_i=p_{\mathrm{t}i}+p_{\mathrm{b}i}.$$

(19)

In practice, calculating this extra-ordinary permutation entropy (XPeEn) only marginally differs from the standard approach: probability distribution estimates for $2\times m!=12$ distinct patterns are to be obtained.

4.2. Pseudo-Complexity in Ordinal Pattern Analysis

Apart from nomenclature, XPeEn is theoretically a reasonable extension. Using conventional PeEn of order $m=3$, we implicitly agree to distinguish two upward-peak patterns ($\mathbf{132}$ and $\mathbf{231}$) as well as two downward-peak patterns ($\mathbf{213}$ and $\mathbf{312}$). By implication, the additional distinction between top-heavy and bottom-heavy patterns is then equally justified. Both kinds of partitions stem from the same principle, and the degree of granularity that better matches the signal characteristics of EEG is not evident: ultimately, the relations between ordinal patterns and EEG are not understood.

Understood are, however, the properties of Shannon entropy. According to Equation (14), and in conjunction with Equation (19), it holds for any pair of extra-ordinal patterns ${\mathsf{\Pi }}_{\mathrm{t}i}$ and ${\mathsf{\Pi }}_{\mathrm{b}i}$ that their summed contribution to XPeEn is

$${\mathrm{H}}_{\Delta }(p_i,{\beta }_i)=p_{\mathrm{t}i}log\frac{1}{p_{\mathrm{t}i}}+p_{\mathrm{b}i}log\frac{1}{p_{\mathrm{b}i}}=p_i{\mathrm{H}}_{\mathrm{b}}\left({\beta }_i\right)+p_{i}log\frac{1}{p_i},$$

(20)

where ${\beta }_i$ is the probability balance coefficient

$${\beta }_i=\frac{p_{\mathrm{t}i}}{p_{\mathrm{t}i}+p_{\mathrm{b}i}}=\frac{p_{\mathrm{t}i}}{p_i}.$$

(21)

By estimating the variance of a particular ${\beta }_i$ from a suitably large collection of EEG epochs, we can thus quantify the impact its corresponding pattern pair ${\mathsf{\Pi }}_{\mathrm{t}i}$ and ${\mathsf{\Pi }}_{\mathrm{b}i}$ has on the dynamics of XPeEn. As an example as intuitive as unlikely, imagine that the variance of some balance coefficient ${\beta }_i$ is found to be zero, such that ${\beta }_i$ and ${\mathrm{H}}_{\mathrm{b}}\left({\beta }_i\right)$ remain constant for any epoch of EEG. Equation (20) then implies that ${\mathsf{\Pi }}_{\mathrm{t}i}$ and ${\mathsf{\Pi }}_{\mathrm{b}i}$ can merely contribute to the dynamics (= value changes) of XPeEn via the sum of their probabilities $p_i$, which is the probability of the conventional ordinal pattern ${\mathsf{\Pi }}_i$. In this case, the distinction between ${\mathsf{\Pi }}_{\mathrm{t}i}$ and ${\mathsf{\Pi }}_{\mathrm{b}i}$ is redundant in the first place: without impairing the dynamics of XPeEn, the respective extra-ordinal patterns ${\mathsf{\Pi }}_{\mathrm{t}i}$ and ${\mathsf{\Pi }}_{\mathrm{b}i}$ can be merged into their common ancestor, the traditional ordinal pattern ${\mathsf{\Pi }}_i$.

More realistically, any balance coefficient will most likely feature some variance. Nevertheless, estimating the variances of all balance coefficients ${\beta }_{ij}\in B$ as defined by Equation (16), we can rank the corresponding pattern pairs by their relative contribution to the overall dynamics of XPeEn. Thus, we can separate decisive pattern pairs from redundant ones—and tell apart the complexity of the data from the pseudo-complexity of the method.

Admittedly, if the lengths of the EEG epochs analysed suffice for accurate probability estimation, it is not necessarily detrimental to use more ordinal patterns than the underlying signal characteristics demand for. However, doing so conveys the impression of a more intricate phenomenon than there might actually be, and could hence promote the misinterpretation of analysis results.

5. The Entropy of Peaks

The principle just demonstrated is directly transferable to actual PeEn. To this end, let us choose the six ordinal patterns of order $m=3$ as a starting point. In addition, let us limit our considerations to the time delay $\tau =1$ for now—arbitrary time delays will be discussed near the end of the manuscript. We are thus looking at probability distributions of the form

$$\mathbf{p}={\left(\begin{array}{cccccc}{p}_{\mathbf{123}}& p_{\mathbf{132}}& p_{\mathbf{213}}& p_{\mathbf{231}}& p_{\mathbf{312}}& p_{\mathbf{321}}\end{array}\right)}^T$$

(22)

and at their corresponding balance coefficients in terms of Equations (9) and (15),

$$\left(\begin{array}{cccccc}{\beta }_{\mathbf{123}/\mathbf{123}}& {\beta }_{\mathbf{123}/\mathbf{132}}& {\beta }_{\mathbf{123}/\mathbf{213}}& {\beta }_{\mathbf{123}/\mathbf{231}}& {\beta }_{\mathbf{123}/\mathbf{312}}& {\beta }_{\mathbf{123}/\mathbf{321}}\\ {\beta }_{\mathbf{132}/\mathbf{123}}& {\beta }_{\mathbf{132}/\mathbf{132}}& {\beta }_{\mathbf{132}/\mathbf{213}}& {\beta }_{\mathbf{132}/\mathbf{231}}& {\beta }_{\mathbf{132}/\mathbf{312}}& {\beta }_{\mathbf{132}/\mathbf{321}}\\ {\beta }_{\mathbf{213}/\mathbf{123}}& {\beta }_{\mathbf{213}/\mathbf{132}}& {\beta }_{\mathbf{213}/\mathbf{213}}& {\beta }_{\mathbf{213}/\mathbf{231}}& {\beta }_{\mathbf{213}/\mathbf{312}}& {\beta }_{\mathbf{213}/\mathbf{321}}\\ {\beta }_{\mathbf{231}/\mathbf{123}}& {\beta }_{\mathbf{231}/\mathbf{132}}& {\beta }_{\mathbf{231}/\mathbf{213}}& {\beta }_{\mathbf{231}/\mathbf{231}}& {\beta }_{\mathbf{231}/\mathbf{312}}& {\beta }_{\mathbf{231}/\mathbf{321}}\\ {\beta }_{\mathbf{312}/\mathbf{123}}& {\beta }_{\mathbf{312}/\mathbf{132}}& {\beta }_{\mathbf{312}/\mathbf{213}}& {\beta }_{\mathbf{312}/\mathbf{231}}& {\beta }_{\mathbf{312}/\mathbf{312}}& {\beta }_{\mathbf{312}/\mathbf{321}}\\ {\beta }_{\mathbf{321}/\mathbf{123}}& {\beta }_{\mathbf{321}/\mathbf{132}}& {\beta }_{\mathbf{321}/\mathbf{213}}& {\beta }_{\mathbf{321}/\mathbf{231}}& {\beta }_{\mathbf{321}/\mathbf{312}}& {\beta }_{\mathbf{321}/\mathbf{321}}\end{array}\right).$$

(23)

Due to the symmetry relation of Equation (10), it suffices to consider the matrix entries above the main diagonal, that is, the coefficients ${\beta }_{ij}$ with $1\le i<j\le 6$. We shall aggregate those 15 elements into a vector of balance coefficients

$$\mathit{\beta }={\left(\begin{array}{ccccccc}{\beta }_{\mathbf{123}/\mathbf{132}}& {\beta }_{\mathbf{123}/\mathbf{213}}& {\beta }_{\mathbf{123}/\mathbf{231}}& \dots & {\beta }_{\mathbf{231}/\mathbf{312}}& {\beta }_{\mathbf{231}/\mathbf{321}}& {\beta }_{\mathbf{312}/\mathbf{321}}\end{array}\right)}^T.$$

(24)

Henceforth, we will use linear indexing, that is, the notation ${\beta }_i$ with $1\le i\le 15$, when referring to the individual elements of the vector $\mathit{\beta }$.

5.1. An Open-Source, Open-Data Approach

The statistics of balance coefficients in EEG have to be assessed empirically, that is, from a suitably large collection of data. The analyses presented here were carried out on the CAP Sleep Database [12], a set of 108 polysomnographic recordings conducted at the Ospedale Maggiore di Parma, and kindly dedicated to the public domain. Sleep EEG data are particularly suitable for the task at hand: PeEn reportedly varies for different stages of natural sleep [8,13], and polysomnographic recordings usually encompass hours of contiguous data per subject. In addition, a multitude of pathologies justify polysomnographic clarification, which increases the heterogeneity of the data—another plus for our use case.

To aid reproducibility, using a dataset available to all of the scientific community was considered obligatory. The choice for the CAP Sleep Database was motivated by Bandt’s recent publication on PeEn [8]. The dataset is hosted by PhysioNet [14] and provided free of charge. Once more in the interest of reproducible science [15], the free and open-source numerical computation software GNU Octave 4.2 was utilised for all data analyses. The code is available in the supplementaries—and was tested compatible with MATLAB R2015b (MATLAB is a registered trademark of: The MathWorks, Inc., Natick, MA, USA).

5.2. EEG Segmentation and Processing

The CAP Sleep Database is a collection of 108 files stored in European Data Format (EDF) [16]. Four of them had to be rejected due to apparent inconsistencies in their file headers. From the remaining polysomnograms, all channels containing EEG were selected, irrespective of electrode location and referencing. These signals were converted to a common sampling rate of $200\mathrm{Hz}$. No additional preprocessing nor any artefact correction were applied. The resampled EEG traces were then split into non-overlapping epochs of $20\mathrm{s}$ duration each. To avoid biasing the ordinal pattern distributions, epochs containing ties (= triples of values that violate the constraint $x_1\ne x_2\ne x_3$) were discarded. The procedure yielded $N_{\mathrm{cap}}=1.6\times {10}^6$ individual segments—more than a year of (single-channel) EEG data.

Using the order $m=3$ and time delay $\tau =1$, epochs were then mapped onto sequences of ordinal patterns according to Equations (1) and (2). For each of the resulting $N_{\mathrm{cap}}$ pattern sequences, its distribution vector

$${\mathbf{p}}_i={\left(\begin{array}{cccccc}{p}_{\mathbf{123}i}& p_{\mathbf{132}i}& p_{\mathbf{213}i}& p_{\mathbf{231}i}& p_{\mathbf{312}i}& p_{\mathbf{321}i}\end{array}\right)}^T$$

(25)

was estimated in terms of Equation (4), and PeEn calculated as per Equation (5). Dividing by $logm!$, results were normalised to the interval $[0,1]$. As is to be expected for sleep EEG, the data feature a wide range of PeEn values. Figure 4 depicts a density estimate of the non-trivially distributed entropy vector

$${\mathbf{h}}_{\mathrm{cap}}={\left(\begin{array}{cccc}{h}_1& h_2& \dots & h_{N_{\mathrm{cap}}}\end{array}\right)}^T.$$

(26)

5.3. The Principle Components of Probability Balances

Furthermore, a vector of balance coefficients ${\mathit{\beta }}_i$ as per Equation (24) was computed for each ordinal pattern distribution ${\mathbf{p}}_i$. These vectors were concatenated to form the $N_{\mathrm{cap}}\times 15$ data matrix

$$\mathbf{B}={\left(\begin{array}{cccc}{\mathit{\beta }}_1& {\mathit{\beta }}_2& \cdots & {\mathit{\beta }}_{N_{\mathrm{cap}}}\end{array}\right)}^T.$$

(27)

Interpreting each vector ${\mathit{\beta }}_i$ (each row in $\mathbf{B}$) as one point in a 15-dimensional feature space, principle component analysis (PCA) was then carried out on the data matrix. This method constitutes a perfect match for the problem at hand, given that probability balance coefficients shall be ranked by their contribution to the overall variance in the data matrix $\mathbf{B}$. Let us notate PCA in terms of computing the set of eigenpairs

$$\left\{({\lambda }_i,{\mathbf{w}}_i)|1\le i\le 15\right\}$$

(28)

to the $15\times 15$ covariance matrix ${\mathbf{\Sigma }}_{\mathbf{B}}$ of the data matrix $\mathbf{B}$. The eigenvectors ${\mathbf{w}}_i$ are then the loading vectors of the PCA and map $\mathbf{B}$ onto its principle components ${\mathbf{z}}_i=\mathbf{B}{\mathbf{w}}_i$. Each principle component ${\mathbf{z}}_i$ is a linear combination of the balance coefficients ${\beta }_i$ and explains a fraction of the variance in the dataset $\mathbf{B}$, whereby it holds that ${\lambda }_i=\mathrm{Var}\left({\mathbf{z}}_i\right)$. As is common practice, principle components shall be enumerated such that ${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_{15}$. By convention, the first principle component ${\mathbf{z}}_1$ then contributes most to the variance in $\mathbf{B}$. For the actual CAP dataset, PCA revealed that

$$\frac{1}{{\lambda }_0}\sum _{i=1}^4{\lambda }_i>0.999\mathrm{with}{\lambda }_0=\sum _{i=1}^{15}{\lambda }_i.$$

(29)

Thus, more than $99.9\%$ of variance in $\mathbf{B}$ are due to the principle components ${\mathbf{z}}_1$, ${\mathbf{z}}_2$, ${\mathbf{z}}_3$ and ${\mathbf{z}}_4$ (see 

Table 1. Eigenvalues and explained variations (normalised eigenvalues) of the first four principle components, as well as their Spearman correlation coefficients with the permutation entropy (PeEn) vector ${\mathbf{h}}_{\mathrm{cap}}$.

Component ${\mathbf{z}}_i$	Eigenvalue ${\lambda }_i$	Explained Variation ${\lambda }_i/{\lambda }_0$	Spearman Correlation $\rho ({\mathbf{h}}_{\mathrm{cap}},{\mathbf{z}}_i)$
${\mathbf{z}}_1$	$5.0$ × 10${}^{-2}$	$94.4\%$	$0.99996$
${\mathbf{z}}_2$	$1.5$ × 10${}^{-3}$	$2.8\%$	$-0.006$
${\mathbf{z}}_3$	$1.2$ × 10${}^{-3}$	$2.2\%$	$-0.008$
${\mathbf{z}}_4$	$2.5$ × 10${}^{-4}$	$0.5\%$	$0.02$

). Further considerations were hence limited to the first four principle components. Their density estimates, as shown in Figure 5, are highly instructive.

The density of ${\mathbf{z}}_1$ resembles the PeEn distribution depicted in Figure 4, suggesting some correlation between ${\mathbf{z}}_1$ and the entropy values ${\mathbf{h}}_{\mathrm{cap}}$. In contrast, the density estimates of ${\mathbf{z}}_2$, ${\mathbf{z}}_3$ and ${\mathbf{z}}_4$ seem to be normally distributed around zero—they appear to be random noise.

Testing for monotonic correlation between the PeEn values ${\mathbf{h}}_{\mathrm{cap}}$ and the principle components ${\mathbf{z}}_1$, ${\mathbf{z}}_2$, ${\mathbf{z}}_3$ and ${\mathbf{z}}_4$, respectively, Spearman correlation coefficients as reported in Table 1 were obtained. Consistent with the graphical representations of Figure 4 and Figure 5, the first principle component ${\mathbf{z}}_1$ is almost perfectly correlated with PeEn, while ${\mathbf{z}}_2$, ${\mathbf{z}}_3$ and ${\mathbf{z}}_4$ clearly contain uncorrelated noise. Moreover, ${\mathbf{z}}_1$ alone explains more than $94\%$ of the variance in $\mathbf{B}$. It is therefore safe to assume that the dynamics of PeEn are exclusively driven by ${\mathbf{z}}_1$.

5.4. Eliminating Pseudo-Complexity

The nature of the first principle component’s loading vector ${\mathbf{w}}_1$ is the pivotal result of the manuscript at hand. According to the data presented in

Table 2. The principle component ${\mathbf{z}}_1={\beta }^T{\mathbf{w}}_1$ is a (trivially simple) linear combination of a subset of balance coefficients. In addition, all balance coefficients not included in ${\mathbf{z}}_1$ (listed in the lower part of the table) appear to be symmetrically distributed around $1/2$.

Index $i$	Balance ${\beta }_i$	Weight		Mean $\mathrm{E}\left[{\beta }_i\right]$	Media $\mathrm{Md}\left[{\beta }_i\right]$	Mode $\mathrm{Mo}\left[{\beta }_i\right]$
		$w_{1,i}$	$w_{1,i}/w_{1,1}$
1	${\beta }_{\mathbf{123}/\mathbf{132}}$	$-0.353$	$1.000$	$0.8588$	$0.8842$	$0.9000$
2	${\beta }_{\mathbf{123}/\mathbf{213}}$	$-0.354$	$1.002$	$0.8586$	$0.8842$	$0.9000$
3	${\beta }_{\mathbf{123}/\mathbf{231}}$	$-0.353$	$0.999$	$0.8588$	$0.8842$	$0.9000$
4	${\beta }_{\mathbf{123}/\mathbf{312}}$	$-0.352$	$0.997$	$0.8590$	$0.8842$	$0.9000$
9	${\beta }_{\mathbf{132}/\mathbf{321}}$	$0.354$	$-1.003$	$0.1415$	$0.1157$	$0.0909$
12	${\beta }_{\mathbf{213}/\mathbf{321}}$	$0.355$	$-1.005$	$0.1417$	$0.1158$	$0.0909$
14	${\beta }_{\mathbf{231}/\mathbf{321}}$	$0.354$	$-1.002$	$0.1414$	$0.1158$	$0.1111$
15	${\beta }_{\mathbf{312}/\mathbf{321}}$	$0.353$	$-1.000$	$0.1412$	$0.1158$	$0.0909$
5	${\beta }_{\mathbf{123}/\mathbf{321}}$	$0.002$	$-0.005$	$0.5004$	$0.5004$	$0.5000$
6	${\beta }_{\mathbf{132}/\mathbf{213}}$	$-0.001$	$0.003$	$0.4997$	$0.5000$	$0.5000$
7	${\beta }_{\mathbf{132}/\mathbf{231}}$	$0.001$	$-0.001$	$0.5001$	$0.5000$	$0.5000$
8	${\beta }_{\mathbf{132}/\mathbf{312}}$	$0.001$	$-0.004$	$0.5004$	$0.5000$	$0.5000$
10	${\beta }_{\mathbf{213}/\mathbf{231}}$	$0.001$	$-0.004$	$0.5004$	$0.5000$	$0.5000$
11	${\beta }_{\mathbf{213}/\mathbf{312}}$	$0.002$	$-0.007$	$0.5007$	$0.5008$	$0.5000$
13	${\beta }_{\mathbf{231}/\mathbf{312}}$	$0.001$	$-0.003$	$0.5003$	$0.5000$	$0.5000$

, it holds in very good approximation that

$$\frac{{\mathbf{z}}_1}{w_{1,1}}=\frac{{\mathit{\beta }}^T{\mathbf{w}}_1}{w_{1,1}}\approxeq {\left(\begin{array}{c}{\beta }_{\mathbf{123}/\mathbf{132}}\\ {\beta }_{\mathbf{123}/\mathbf{213}}\\ {\beta }_{\mathbf{123}/\mathbf{231}}\\ {\beta }_{\mathbf{123}/\mathbf{312}}\\ {\beta }_{\mathbf{132}/\mathbf{321}}\\ {\beta }_{\mathbf{213}/\mathbf{321}}\\ {\beta }_{\mathbf{231}/\mathbf{321}}\\ {\beta }_{\mathbf{312}/\mathbf{321}}\end{array}\right)}^T\left(\begin{array}{c}\phantom{-}1\\ \phantom{-}1\\ \phantom{-}1\\ \phantom{-}1\\ -1\\ -1\\ -1\\ -1\end{array}\right)={\left(\begin{array}{c}{\beta }_{\mathbf{123}/\mathbf{132}}\\ {\beta }_{\mathbf{123}/\mathbf{213}}\\ {\beta }_{\mathbf{123}/\mathbf{231}}\\ {\beta }_{\mathbf{123}/\mathbf{312}}\\ {\beta }_{\mathbf{321}/\mathbf{132}}\\ {\beta }_{\mathbf{321}/\mathbf{213}}\\ {\beta }_{\mathbf{321}/\mathbf{231}}\\ {\beta }_{\mathbf{321}/\mathbf{312}}\end{array}\right)}^T\left(\begin{array}{c}1\\ 1\\ 1\\ 1\\ 1\\ 1\\ 1\\ 1\end{array}\right).$$

(30)

Here, the property of symmetry as per Equation (10) and a division by $w_{1,1}$ were applied to aid interpretability.

This is a bubble of pseudo-complexity bursting: while we were trying to isolate some pattern pairs that are irrelevant for the dynamics of PeEn in EEG, Equation (30) implies that the balances ${\beta }_i$ among all pattern pairs

$${\mathsf{\Omega }}_{†}=\left\{(\mathbf{123},\mathbf{321}),(\mathbf{132},\mathbf{213}),(\mathbf{132},\mathbf{231}),(\mathbf{132},\mathbf{312}),(\mathbf{213},\mathbf{231}),(\mathbf{213},\mathbf{312}),(\mathbf{231},\mathbf{312})\right\}$$

(31)

contribute virtually nothing to the variance in ${\mathbf{z}}_1$, and are thus negligible in good approximation. Those are the balances between any two peak patterns, as well as the balance ${\beta }_{\mathbf{123}/\mathbf{321}}$ between the all-rising and all-falling pattern.

Furthermore, estimating from data matrix $\mathbf{B}$ the mean, median and mode of each balance coefficient, results as presented in Table 2 were obtained. In particular, it holds true for all pattern pairs in ${\mathsf{\Omega }}_{†}$ that

$$\mathrm{E}\left[{\beta }_i\right]\approxeq \mathrm{Md}\left[{\beta }_i\right]\approxeq \mathrm{Mo}\left[{\beta }_i\right]\approxeq 1/2,$$

(32)

which strongly suggests that their respective balance coefficients are symmetrically distributed around $1/2$. Under this premise, and in conjunction with the definition as per Equation (9),

$$\mathrm{E}\left[{\mathrm{p}}_{\mathbf{123}}\right]\approxeq \mathrm{E}\left[{\mathrm{p}}_{\mathbf{321}}\right]\mathrm{as}\mathrm{well}\mathrm{as}\mathrm{E}\left[{\mathrm{p}}_{\mathbf{132}}\right]\approxeq \mathrm{E}\left[{\mathrm{p}}_{\mathbf{213}}\right]\approxeq \mathrm{E}\left[{\mathrm{p}}_{\mathbf{231}}\right]\approxeq \mathrm{E}\left[{\mathrm{p}}_{\mathbf{312}}\right]$$

(33)

immediately follow. Moreover, and again due to the linear combination of Equation (30), the balance coefficients for the pattern pairs

$$\begin{array}{cc}\hfill {\mathsf{\Omega }}_{★}=\{& (\mathbf{123},\mathbf{132}),(\mathbf{123},\mathbf{213}),(\mathbf{123},\mathbf{231}),(\mathbf{123},\mathbf{312}),\hfill \\ & (\mathbf{321},\mathbf{132}),(\mathbf{321},\mathbf{213}),(\mathbf{321},\mathbf{231}),(\mathbf{321},\mathbf{312})\}\hfill \end{array}$$

(34)

contribute almost equally to the variance in $\mathbf{B}$. These are the balances between any peak pattern and either the all-rising or all-falling pattern. Defining the peak probability $\widehat{p}$ such that

$$\widehat{p}=p_{\mathbf{132}}+p_{\mathbf{213}}+p_{\mathbf{231}}+p_{\mathbf{312}}\mathrm{and}1-\widehat{p}=p_{\mathbf{123}}+p_{\mathbf{321}},$$

(35)

and using the relations given in Equation (33), each of these balance coefficients can be approximated by

$${\beta }_{★}=\frac{(1-\widehat{p})/2}{(1-\widehat{p})/2+\widehat{p}/4}=\frac{2(\widehat{p}-1)}{\widehat{p}-2}.$$

(36)

The linear combination given by Equation (30) can therefore be reduced to the expression

$$\frac{{\mathbf{z}}_1}{w_{1,1}}\approxeq 8{\beta }_{★}=\frac{16(\widehat{p}-1)}{\widehat{p}-2}.$$

(37)

Let us assume that equality holds in the above relation, such that all variance within the decisive principle component ${\mathbf{z}}_1$ is exclusively induced by the peak probability $\widehat{p}$. Without compromising the dynamics of PeEn, all four peak patterns $\mathbf{132}$, $\mathbf{213}$, $\mathbf{231}$ and $\mathbf{312}$ can then be merged into just one ordinal pattern—the peak pattern. Likewise, the all-rising and all-falling patterns $\mathbf{123}$ and $\mathbf{321}$ can be unified, yielding a single edge pattern. On this basis, PeEn of order $m=3$ and time delay $\tau =1$ can be replaced by a much simpler entropy of peaks function

$$\mathrm{H}\left(\widehat{p}\right)=\widehat{p}log\frac{4}{\widehat{p}}+(1-\widehat{p})log\frac{2}{1-\widehat{p}}={\mathrm{H}}_{\mathrm{b}}\left(\widehat{p}\right)+\widehat{p}log2+log2.$$

(38)

This function is concave with its maximum at $\widehat{p}=2/3$, where $p_i=1/6$ for all i. The latter is consistent with Shannon entropy being maximal for the uniform distribution. The graph of $\mathrm{H}\left(\widehat{p}\right)$ is depicted in Figure 6.

For the purpose of validation, $\mathrm{H}\left(\widehat{p}\right)$ was calculated for all EEG epochs obtained from the CAP Sleep Database. Using the same overall approach as described in Section 5.2, but now computing the entropy of peaks according to Equations (35) and (38), a vector of $N_{\mathrm{cap}}$ results ${\widehat{\mathbf{h}}}_{\mathrm{cap}}$ was obtained. Comparison with the conventionally obtained PeEn values ${\mathbf{h}}_{\mathrm{cap}}$ yielded a Pearson correlation coefficient of $r({\mathbf{h}}_{\mathrm{cap}},{\widehat{\mathbf{h}}}_{\mathrm{cap}})=0.99998$ and a mean relative error of

$$\overline{\eta }({\mathbf{h}}_{\mathrm{cap}},{\widehat{\mathbf{h}}}_{\mathrm{cap}})=\frac{1}{N_{\mathrm{cap}}}\sum _{i=1}^{N_{\mathrm{cap}}}\frac{|h_i-{\widehat{h}}_i|}{h_i}=0.03\%.$$

(39)

Thus, the ordinal patterns obtained from the CAP Sleep Database can essentially be replaced by binary symbols that encode the disjoint properties “three consecutive values form a peak” and “three consecutive values form an edge”.

6. Linearising Permutation Entropy

At this point, we have substantiated the working hypothesis that, for EEG analysis, PeEn of order $m=3$ may be a function of the peak probability $\widehat{p}$ only. This probability can be estimated from an epoch of n EEG samples $(x_1,x_2,\dots ,x_n)$ by computing

$$\widehat{p}=\frac{\widehat{n}}{n}\mathrm{with}\widehat{n}=\sum _{k=2}^{n-1}\left|\left[x_k-x_{k-1}\ge 0\right]-\left[x_{k+1}-x_k\ge 0\right]\right|,$$

(40)

where the number of peaks $\widehat{n}$ is expressed using Iverson bracket notation [10]. Let us look in more detail at the properties of the decisive peak pattern.

6.1. Signal Peaks and Spectral Bandwidth

In mathematical terms, a peak is a local extremum. Any real-valued, twice differentiable function s has a peak at $t_k$ if it holds that $s^{\prime }\left(t_k\right)=0$ and $s^{″}\left(t_k\right)\ne 0$. In other words, a peak in s is a zero crossing in its first derivative $s^{\prime }$. Now, consider a real-valued function s, periodic in $T=1/f_0$ and band limited to frequencies including $N{f}_0$, which is

$$s\left(t\right)=\sum _{k=0}^N{a}_{k}sin(2\pi k{f}_{0}t+{\varphi }_k)\mathrm{with}{a}_k\in \mathbb{R}.$$

(41)

Its first derivative is again a real-valued function, periodic in T and band limited to $N{f}_0$, namely

$$s^{\prime }\left(t\right)=\sum _{k=0}^{N}2\pi a_{k}k{f}_{0}cos(2\pi k{f}_{0}t+{\varphi }_k)=\sum _{k=0}^N{b}_{k}sin(2\pi k{f}_{0}t+{\vartheta }_k).$$

(42)

Being the N-th partial sum of a Fourier series, the function can alternatively be represented using complex exponentials,

$$s^{\prime }\left(t\right)=\sum _{k=-N}^N{c}_k{e}^{\mathrm{j}2\pi k{f}_{0}t}\mathrm{with}\mathrm{j}=\sqrt{-1}\mathrm{and}{c}_k\in \mathbb{C}.$$

(43)

This expression constitutes a trigonometric polynomial of order N and therefore has exactly $2N$ roots per period [17]. As for any real-valued polynomial, its zeroes

$$Z=\left\{t_k|s^{\prime }\left(t_k\right)=0,k\in \left\{1,2,\dots ,2N\right\}\right\}$$

(44)

are not necessarily distinct, and are either real or occur in complex conjugate pairs. Let us focus on the subset of zero crossings exclusively, which is

$$Z_{ø}=\left\{t_k\in \mathbb{R}|s^{″}\left(t_k\right)\ne 0\right\}\subseteq Z.$$

(45)

With $Z_{ø}$ being a subset of Z, it obviously holds that $|Z_{ø}|\le 2N$. Furthermore, the number of peaks $\widehat{n}$ in the function s equals the number of zero crossings $|Z_{ø}|$ in its first derivative $s^{\prime }$. It therefore holds that

$$\widehat{n}=\left|Z_{ø}\right|\le \frac{2f_{\max }}{f_0}\mathrm{with}{f}_{\max }=N{f}_0.$$

(46)

Sampling one full cycle of the function s at a sampling rate of $f_{\mathrm{s}}=n{f}_0$, we obtain a sequence of n values $(x_1,x_2,\dots ,x_n)$. We can then express the maximum number of peaks per period by

$$\widehat{n}\le \frac{2n{f}_{\max }}{f_{\mathrm{s}}}\mathrm{with}{f}_{\max }\le \frac{f_{\mathrm{s}}}{2}.$$

(47)

Consequently, it holds for the peak probability $\widehat{p}$ that

$$\widehat{p}=\frac{\widehat{n}}{n}\le \frac{2f_{\max }}{f_{\mathrm{s}}}.$$

(48)

From a practical point of view, this relation is not limited to periodic signals. To assess $f_{\max }$ in any signal segment, Fourier transform needs to be applied, which by definition constitutes a periodic continuation of the analysis window. Thus, the inequality given by Equation (48) provides an upper bound on the number of peaks in any band limited signal.

This bound may be naïvely conservative, but it is easily derived and suffices our purpose: let us reconsider the entropy of peaks function $\mathrm{H}\left(\widehat{p}\right)$ as defined by Equation (38), which increases monotonically on the interval $0\le \widehat{p}\le 2/3$. In combination with Equation (48), this means that PeEn is monotonic with $\widehat{p}$ for EEG epochs band limited such that

$$f_{\max }\le \frac{f_{\mathrm{s}}}{3}.$$

(49)

As stated earlier, the absolute value of PeEn obtained from a single EEG epoch is insignificant. It is the contrast in value across distinct epochs—its dynamics—that renders PeEn a valuable tool for EEG analysis. However, this contrast persists for any parameter that grows monotonically with PeEn. Therefore, the following constitutes an astounding, yet sensible proposition: given an epoch of EEG, digitised at sampling rate $f_{\mathrm{s}}$ and band limited such that $f_{\max }\le f_{\mathrm{s}}/3$, PeEn of order $m=3$ and time delay $\tau =1$ can be substituted by the peak probability $\widehat{p}$, that is, the number of peaks per unit time.

6.2. Counting the Zigzags of EEG

For an experimental validation of the above, further EEG analyses were performed. To rule out possible peculiarities of the CAP Sleep Database (or sleep EEG in general), the data corpus was expanded by a set of EEG signals from a considerably different clinical setting. The CHB-MIT Scalp EEG Database is a collection of EEG recordings obtained from 22 paediatric patients, all suffering from epileptic seizures [18]. Acquired at Boston Children’s Hospital, the data were generously put into the public domain and are available from PhysioNet [14]. In strict analogy to the procedure described in Section 5.2, a total of $N_{\mathrm{mit}}=4.1\times {10}^6$ EEG epochs of $20\mathrm{s}$ duration were obtained from this dataset—another two and a half years of (single-channel) EEG data. For each of these epochs, PeEn of order $m=3$ and time delay $\tau =1$ was computed according to Equations (1)–(6), yielding a vector of results ${\mathbf{h}}_{\mathrm{mit}}$.

From each EEG epoch of the CAP and CHB-MIT datasets, a peak probability $\widehat{p}$ as per Equation (40) was also obtained. Results were aggregated to form the respective vectors ${\widehat{\mathbf{p}}}_{\mathrm{cap}}$ and ${\widehat{\mathbf{p}}}_{\mathrm{mit}}$. For both datasets, individual Spearman correlation coefficients between peak probabilities $\widehat{\mathbf{p}}$ and PeEn values $\mathbf{h}$ were then computed for (1) the subset of EEG epochs with $\widehat{p}\le 2/3$, (2) the subset of EEG epochs with $\widehat{p}\ge 2/3$, and (3) the entire set of epochs. Results are listed in

Table 3. Spearman correlation coefficients between PeEn values $\mathbf{h}$ and peak probabilities $\widehat{\mathbf{p}}$ for the electroencephalographic data from the CAP Sleep Database and the CHB-MIT Scalp EEG Database. Quantities subscripted “high” and “low” correspond to results obtained from those subsets of signal epochs for which $\widehat{p}\ge 2/3$ or $\widehat{p}\le 2/3$ hold, respectively.

Database	Number of Epochs			Spearman Correlation
	$N$	$N_{\mathrm{low}}$	$N_{\mathrm{high}}$	$\rho (\mathbf{h},\widehat{\mathbf{p}})$	${\rho }_{\mathrm{low}}(\mathbf{h},\widehat{\mathbf{p}})$	${\rho }_{\mathrm{high}}(\mathbf{h},\widehat{\mathbf{p}})$
CAP	$1.6\times {10}^6$	$1.6\times {10}^6$	$6.0\times {10}^2$	$0.99998$	$0.99998$	$-0.923$
CHB-MIT	$4.1\times {10}^6$	$4.0\times {10}^6$	$1.1\times {10}^5$	$0.99937$	$0.99993$	$-0.963$

.

The above results clearly support the hypothesis that PeEn of order $m=3$ and time delay $\tau =1$ increases monotonically with the peak probability for $\widehat{p}\le 2/3$. The correlation is inverted for peak probabilities exceeding this bound. Both results are in line with the function graph of Figure 6. When testing on the full set of peak probabilities, strong positive correlation is maintained: epochs for which $\widehat{p}>2/3$ holds are too seldom to make a difference. This directly relates to the bandwidth constraint of Equation (49). Given our chosen sampling rate $f_{\mathrm{s}}=200\mathrm{Hz}$, the peak probability cannot increase beyond $\widehat{p}=2/3$ for any EEG epoch band limited to frequencies below $f_{\max }\approxeq 66.7\mathrm{Hz}$.

Not only is the relative bandwidth $f_{\max }\le f_{\mathrm{s}}/3$ rarely exceeded in the datasets considered here, but the same presumably holds true for the vast majority of EEG in general. Sampling rates below $100\mathrm{Hz}$ are quite uncommon for EEG, while spectral content above $30\mathrm{Hz}$ is often regarded as insignificant (and actively suppressed using low-pass filtering). Moreover, an investigator actually interested in gamma-band oscillations will choose a considerably higher sampling rate. It is therefore likely that the decrease in PeEn for peak probabilities $\widehat{p}>2/3$ is not a part of the actual effect, but a rarely observed flaw of the analysis technique. In essence, utilising PeEn of order $m=3$ and time delay $\tau =1$ for EEG analysis may constitute an involved way of counting zigzags.

7. Permutation Entropy as a Spectral Estimator

While counting signal peaks may seem a trivial procedure at first glance, a profound mathematical theory supports this approach to signal analysis. It was formulated by Benjamin Kedem in the 1980s and is centred around the notion of higher order crossings. Only a humble outline of his framework—and one limited to aspects that aid the interpretation of PeEn—shall here be provided. For a large-scale, yet low-threshold introduction to the theory of higher order crossings, the reader is referred to Kedem’s comprehensive article “Spectral Analysis and Discrimination by Zero-Crossings” [19], which served as the primary source for the section to follow.

7.1. Zero Crossings and the Dominant Frequency Principle

Signal analysis by means of higher order crossings is a generalisation of the more commonly known zero crossing rate—a classic method of fundamental frequency estimation [20]. Self-descriptively, the zero crossing rate $r_{ø}$ of a signal is its average number of x-intercepts per unit time. It is obtained from a discrete sequence of samples $(x_1,x_2,\dots ,x_n)$ by means of

$$r_{ø}=\frac{n_{ø}}{n-1}\mathrm{with}{n}_{ø}=\sum _{i=2}^n\left|\left[x_i\ge 0\right]-\left[x_{i-1}\ge 0\right]\right|.$$

(50)

Kedem refined the meaning of this quantity, demonstrating that the zero crossing rate of a discrete-time signal constitutes an approximation of the centroid of its power spectrum

$$c_S=\frac{1}{m_S}{\int }_0^{\infty }fS\left(f\right)\mathrm{d}f\mathrm{with}{m}_S={\int }_0^{\infty }S\left(f\right)\mathrm{d}f.$$

(51)

Here, $S\left(f\right)$ denotes the power spectral density of the signal and $m_S$ is its (single-sided) spectral mass. Intuitively, if $S\left(f\right)$ contains a predominant frequency $f_{\mathrm{c}}$, the spectral centroid $c_S$ will gravitate towards this particular frequency. The more $S\left(f_{\mathrm{c}}\right)$ outweighs the rest of the spectrum, the shorter the distance between $c_S$ and $f_{\mathrm{c}}$ will be. Kedem termed this the “dominant frequency principle”. Regarding its relation with the zero crossing rate, he further deduced that

$$r_{ø}=\frac{n_{ø}}{N-1}\approxeq 2c_S.$$

(52)

Maintaining the interdisciplinary scope of the article, we do not discuss the mathematical derivation of this relation—it can be found in [19]. Instead, a very conclusive example from the same publication shall be reproduced. Let our signal be the sinusoidal function $s\left(t\right)=Asin(2\pi f_{\mathrm{c}}t+\phi )$, which has the power spectrum

$$S\left(f\right)=\left\{\begin{array}{cc}{A}^2/4& \mathrm{for}f=\pm f_{\mathrm{c}},\hfill \\ 0& \mathrm{otherwise}.\hfill \end{array}\right.$$

(53)

With $f_{\mathrm{c}}$ being the only frequency in the spectrum, it is obviously the dominant frequency. Calculating the spectral centroid of the signal, we see that the equality $c_S=f_{\mathrm{c}}$ holds true in this case. Hence, this example is a showcase for the dominant frequency principle. Moreover, and just as comprehensible, the relation

$$r_{ø}=2f_{\mathrm{c}}=2c_S$$

(54)

also applies: any sinusoidal signal has exactly two zero crossings per cycle.

7.2. Kedem’s Higher Order Crossings

Within the framework of higher order crossings, the peak probability $\widehat{p}$ of our simplified PeEn is closely related to the zero crossing rate $r_{ø}$. Let $(x_1,x_2,\dots ,x_n)$ be a sequence of values sampled from a continuous-time EEG signal $s\left(t\right)$ at discrete time steps $t=i/f_{\mathrm{s}}=i{T}_{\mathrm{s}}$. Estimating the peak probability $\widehat{p}$ of that sequence as per Equation (40) is equivalent to obtaining the zero crossing rate $r_{ø}$ of the difference sequence

$$(\nabla x_2,\nabla x_3,\dots ,\nabla x_n),$$

(55)

whereby the difference operator ∇ is defined such that

$$\nabla x_i=x_i-x_{i-1}=s\left(i{T}_{\mathrm{s}}\right)-s(i{T}_{\mathrm{s}}-T_{\mathrm{s}}).$$

(56)

In terms of Kedem’s theory, the generalised zero crossings obtained by applying the ∇-operator k times recursively are called: higher order crossings of the ${\mathrm{D}}_{k+1}$ type. Thus, actual zero crossings are higher order crossings of the ${\mathrm{D}}_1$ type, while the peak patterns considered in this manuscript are of type ${\mathrm{D}}_2$.

Kedem’s approach is highly convenient in that the ∇-operator consitutes a linear time-invariant system with the power transfer function

$${\left|H\left(f\right)\right|}^2={\left|(1-e^{-\mathrm{j}2\pi f{T}_{\mathrm{s}}})\right|}^2=2-2cos\left(2\pi f{T}_{\mathrm{s}}\right).$$

(57)

This system acts as a high-pass filter in the digital domain. Its frequency response is depicted in Figure 7. In practical terms, the filter compresses the spectral mass, effectively shifting its centroid further to the right with each recurrent application of the ∇-operator.

7.3. The Spectral Estimation Hypothesis

We have now discussed all aspects necessary to propose an explanation for PeEn of order $m=3$ and time delay $\tau =1$ in EEG. Let us summarise: due to their specific distribution in EEG, it may suffice to distinguish two kinds of ordinal patterns—peaks and edges. Their respective probabilities are $\widehat{p}$ and $1-\widehat{p}$. The peak probability $\widehat{p}$ is the zero crossing rate $r_{ø}$ of a high-pass filtered version of the signal, whereas the zero crossing rate $r_{ø}$ is an estimate of its power spectral centroid. For EEG epochs with $\widehat{p}\le 2/3$, the peak probability $\widehat{p}$ and PeEn are monotone and can therefore be used interchangeably in comparative analyses.

On that basis, it is a rational hypothesis that applying PeEn of order $m=3$ and time delay $\tau =1$ to EEG epochs with peak probability $\widehat{p}\le 2/3$ effectively means: high-pass filtering the signal according to Equation (57) and estimating the centroid of the resulting weighted power spectrum. Both steps are standard procedures of linear signal processing, possibly rendering this EEG parameter a spectral estimator in disguise.

8. Higher Pattern Orders and Time Delays

Throughout the manuscript, considerations and propositions were repeatedly limited to the pattern order $m=3$ and the time delay $\tau =1$. Regarding the latter, a general discussion was postponed for reasons of simplicity only. It will be presented in the following section. Conversely, formulating a universal interpretation for pattern orders $m>3$ is a subject still under investigation—and hopefully one to be reported on in the future. For the time being, the reader is invited to consider the following preliminary results.

8.1. Prospects for Patterns of Higher Order

An ordinal pattern of order $m=3$ is either a peak or an edge. For higher orders, there is no such simple duality. In the general case, an ordinal pattern of order m can contain up to $m-2$ peaks, while the peaks of consecutive ordinal patterns may overlap. Moreover, the marginal probability relation of Equation (33) does not translate to the “sub-patterns” of order $m=3$ that form a pattern of higher order. Conditional probabilities have to be considered instead.

Despite such intricacies, the principle behaviour of PeEn does apparently not change when increasing the order m within statistically reasonable limits. Calculating PeEn of orders $m\in \{3,4,5\}$ for all EEG epochs obtained from the CAP and CHB-MIT datasets, results as depicted in Figure 8 were obtained. Strong monotonic correlation between PeEn and the peak probability $\widehat{p}$ is maintained for higher pattern orders. This hints at a possible generalisation of our spectral estimation hypothesis.

8.2. Sampling, Resampling and Aliasing

In comparison to higher pattern orders m, interpreting arbitrary time delays $\tau$ is straightforward. Under the premise that PeEn acts as a spectral estimator in EEG analysis, the principles of linear systems theory fully apply. Using time delays $\tau >1$ is then a means of downsampling (and thus, of band limiting), but also constitutes a violation of the Nyquist–Shannon sampling theorem. Let us substantiate.

The Nyquist–Shannon sampling theorem [21] states that any time-discrete representation of an analogue signal is implicitly band limited to frequencies below and including half its sampling rate $f_{\mathrm{s}}$. Due to this fundamental constraint of digital signal processing, an analogue low-pass filter is commonly applied during signal acquisition. This filter prevents distortions known as aliasing—spectral content beyond the Nyquist frequency $f_{\mathrm{s}}/2$ folding back into the usable frequency range as given by the mapping

$$f\mapsto \frac{f_{\mathrm{s}}}{2}-|(fmod{f}_{\mathrm{s}})-\frac{f_{\mathrm{s}}}{2}|\mathrm{with}f\ge 0.$$

(58)

This relation describes a triangle waveform with its extrema at multiples of $f_{\mathrm{s}}/2$. Its graph is displayed in Figure 9.

The same constraint applies to sampling rate conversion in the digital domain. Correctly reducing the sampling rate $f_{\mathrm{s}}$ of a digital signal by a factor r therefore involves two steps: limit the bandwidth to $f_{\mathrm{s}}/2r$, then omit all but every r-th sample. If the band limiting is skipped, downsampling is prone to aliasing. Any spectral content exceeding $f_{\mathrm{s}}/2r$ is not compatible with the targeted sampling rate $f_{\mathrm{s}}/r$ and folds back into the representable frequency range, corrupting the resulting signal.

8.3. The Time Delay as a Downsampling Factor

Let us look again at what is called the symbolisation technique. Given a sequence of samples $(x_1,x_2,\dots ,x_n)$, its ordinal pattern ${\pi }_t$ of order m and time delay $\tau$ at time index t is obtained from the m-tuple of values

$$(x_t,x_{t+\tau },x_{t+2\tau },\dots ,x_{t+(m-1)\tau }).$$

(59)

To estimate the ordinal pattern distribution of the sequence, we determine ${\pi }_t$ for $1\le t\le n-\tau (m-1)$ and count the number of occurrences of each pattern ${\mathsf{\Pi }}_i\in {\mathsf{\Omega }}_m$. This is the method described by Equations (1)–(4).

For $\tau >1$, consider an alternative, yet entirely equivalent estimation algorithm: first, partition the sequence of samples $(x_1,x_2,\dots ,x_n)$ into $\tau$ sub-sequences of the form

$$\forall \lambda \in \{1,2,\dots ,\tau \}:(x_{\lambda },x_{\lambda +\tau },x_{\lambda +2\tau },\dots ,x_{\lambda +k\tau })\mathrm{whereby}k=⌊\frac{n-\lambda }{\tau }⌋+1.$$

(60)

Next, estimate from each of the $\tau$ sub-sequences its ordinal pattern distribution ${\mathbf{p}}_{\lambda }$, now using the time delay ${\tau }^{★}=1$. The ordinal pattern distribution $\mathbf{p}$ of the initial sequence can then be obtained by cumulation, that is

$$\mathbf{p}=\sum _{\lambda =1}^{\tau }{\mathbf{p}}_{\lambda }={\mathbf{p}}_1+{\mathbf{p}}_2+\dots +{\mathbf{p}}_{\tau }.$$

(61)

Observe that generating a sub-sequence as per Equation (60) is identical to downsampling the sequence $(x_1,x_2,\dots ,x_n)$ by a factor $\tau$ (whereby the index $\lambda$ merely specifies which sample to keep per block of $\tau$ consecutive samples). This circumstance is both the explanation and the fundamental problem of using $\tau >1$ for ordinal pattern encoding: we implicitly—but no less effectively—downsample the signal without any anti-aliasing measures taken. Inevitably, all spectral content exceeding $f_{\mathrm{s}}/2\tau$ folds back into the representable frequency range. Therefore, calculating PeEn from EEG using $\tau >1$ is equivalent to computing PeEn for the time delay ${\tau }^{★}=1$ from a set of $\tau$ downsampled (and therefore band limited), yet nonlinearly distorted versions of the actual signal.

Notice that this is not a conceptual flaw of PeEn. One has to keep in mind that the quantity was developed as a complexity measure for time series. Working with less generic data, additional constraints and requirements may arise. In particular, while any digital signal can be described as a time series, not all time series are digital signals in terms of linear systems theory. Only if the requirements of the sampling theorem are adhered to can a discrete sequence of voltage measurements be a valid representation of analogue EEG.

8.4. Frequency Aliasing in Permutation Entropy

Computing the entropy of peaks $\mathrm{H}\left(\widehat{p}\right)$ as per Equation (38) from sinusoidal signals of varying frequency $f_{\mathrm{c}}\le f_{\mathrm{s}}/2$, the aliasing introduced by using time delays $\tau >1$ is easily observed. In accordance with the frequency mapping of Equation (58), and depending on the effective sampling rate $f_{\mathrm{s}}/\tau$, the signal that is actually fed into the estimator is a sinusoidal wave of frequency

$$\widetilde{f_{\mathrm{c}}}=\frac{f_{\mathrm{s}}}{2\tau }-|(f_{\mathrm{c}}mod\frac{f_{\mathrm{s}}}{\tau })-\frac{f_{\mathrm{s}}}{2\tau }|.$$

(62)

Given that any sinusoidal signal has exactly two peaks per cycle, it further holds that $\widehat{p}=2\widetilde{f_{\mathrm{c}}}$. In conjunction with Equation (38), the relation

$$\mathrm{H}\left(\widehat{p}\right)=\mathrm{H}(2\widetilde{f_{\mathrm{c}}})={\mathrm{H}}_{\mathrm{b}}(2\widetilde{f_{\mathrm{c}}})+2\widetilde{f_{\mathrm{c}}}log2+log2$$

(63)

is obtained. Variegating the time delay $\tau$, the periodic foldback effect manifests in the graphs of this family of functions. This is depicted and further discussed in Figure 10.

The presented approach is based on the work of Erik Olofsen, Jamie Sleigh and Albert Dahan. To experimentally investigate the frequency dependence of PeEn, the authors calculated PeEn of order $m=3$ and time delays $\tau \le 2$ for sinusoidal signals of increasing frequency [22]. Our here proposed spectral estimation hypothesis conclusively explains the authors’ observations, including the discontinuous frequency dependence for $\tau =2$.

8.5. Anti-Aliased Permutation Entropy

Using time delays $\tau >1$ for EEG encoding may cause unintended signal distortions. Nevertheless, a multitude of applications in EEG analysis depend on higher time delays [11]. In practice, the impact of aliasing may be less detrimental than the graphs of Figure 10 suggest: EEG is anything but sinusoidal, and while PeEn literally locks onto a single sine wave, the discrepancy will be more subtle for actual EEG. In addition, EEG is commonly oversampled ($f_{\mathrm{s}}\gg 2f_{\max }$), which further mitigates the issue (see Section 6).

Nothing is lost by resolving a source of inaccuracy, though. This is easily achieved by means of sample rate conversion. A possible procedure is as follows: instead of using $\tau >1$ for PeEn calculation, correctly convert the EEG to a new sampling rate $f_{\mathrm{s}}^{★}=f_{\mathrm{s}}/\tau$ and compute PeEn therefrom, now choosing ${\tau }^{★}=1$ for the time delay. An apparent disadvantage of this approach is that it reduces the number of ordinal patterns available for probability estimation by a factor $\tau$. However, the sampling theorem implicitly assures that a single rate-reduced sequence suffices to fully describe the band limited EEG—including its signal peaks. Nothing is gained by counting each peak $\tau$ times.

In the anti-aliased case, a very tangible interpretation of the time delay $\tau$ immediately follows: PeEn effectively splits the signal’s single-sided frequency spectrum into $\tau$ spectral bands of equal width, then exclusively responds to the content of the lowest band. Given our spectral estimation hypothesis, PeEn of order $m=3$ and time delay $\tau \ge 1$ thus approximates the weighted power spectral centroid of the EEG’s frequency band ranging from 0 to $f_{\mathrm{s}}/2\tau$. For applications depending on $\tau >1$, it may hence be worthwhile to further investigate the influence of anti-aliasing. This could be a chance of concisely relating PeEn-based techniques with classic band-spectral methods of EEG analysis.

9. Conclusions

EEG has long become an indispensable tool in research and medical diagnostics, even more so as the progress in digital technology has revolutionised our general approach towards data—biomedical data being no exception.

9.1. There Are Only So Many Signal Characteristics

When analysing EEG, one has to accept its predominantly probabilistic nature. Any signal property that statistically correlates with some behavioural observation is a perfectly valid signal property because, ultimately, none of them are fully explainable mechanistically. While a plethora of EEG analysis techniques are available, the application of a new or more sophisticated method does not necessarily yield new information: analysing the EEG data considered for this manuscript, (1) PeEn of order $m=3$ and time delay $\tau =1$; (2) the power spectral centroid of the signal’s first derivative; and (3) Kedem’s higher order crossings of the ${\mathrm{D}}_2$ type can all be used interchangeably.

The signal property apparently underlying all three methods has been known and used for more than 50 years now. In 1964, Neil Burch and colleagues reported on transferring the EEG analysis method of counting “baseline crossing[s]” from analogue circuitry to a digital computer [23]. The authors explicitly pointed out the “well-known accentuation of the higher frequency component in the derivatives”. Back then, the procedure was called period analysis.

9.2. Next Steps in Ordinal Pattern-Based EEG Analysis

It cannot be over-emphasized that EEG analysis is mostly empirical—including the results of this manuscript. What applies for the CAP and CHB-MIT databases does not necessarily hold true for EEG in general. Hence, the work presented does neither falsify the existence of complexity in EEG, nor its observability using PeEn. Still, it demonstrates that PeEn-based EEG processing does not exclusively relate to this phenomenon.

This article thus extends the framework of ordinal pattern-based EEG processing. For a given analysis task, concurrently computing both PeEn and the peak probability $\widehat{p}$ provides a means of rejecting the null hypothesis of PeEn acting as a spectral estimator. This may enable the investigator to narrow down on cases that truly demand for complexity considerations, while otherwise resorting to a simpler and time-tested analysis technique. To this end, the peak probability $\widehat{p}$ as specified by Equation (40) integrates seamlessly into the framework of ordinal pattern analysis. A convenient side effect: if for a particular analysis task PeEn is found to be replaceable by $\widehat{p}$, the issue of pattern order selection does not arise.

In closing, adding the peak probability parameter to your codebase is straightforward. The expression

$$\mathtt{y}\mathtt{=}\mathtt{mean}\mathtt{\left(}\mathtt{abs}\mathtt{\right(}\mathtt{diff}\mathtt{\left(}\mathtt{diff}\mathtt{\right(}\mathtt{x}\mathtt{)}\mathtt{>}\mathtt{=}\mathtt{0}\mathtt{)}\mathtt{\left)}\mathtt{\right)}\mathtt{;}$$

suffices to have the quantity available within a GNU Octave or MATLAB environment. An extended version that supports sliding window analysis is provided in the supplementaries (see the zztop.m function). While not extensively optimised, its execution speed should meet the requirements of major analysis campaigns: the main challenge in data analysis today is not computational feasibility, but the interpretation of the manifold results obtainable.