# Outlier detection and removal in multifractal analysis of electrophysiological brain signals

Merlin Dumeur, J. Matias Palva, Philippe Ciuciu

► To cite this version:

Merlin Dumeur, J. Matias Palva, Philippe Ciuciu. Outlier detection and removal in multifractal analysis of electrophysiological brain signals. EURASIP Journal on Advances in Signal Processing, 2025, 2025, pp.35. 10.1186/s13634-025-01245-4. hal-05379794

HAL Id: hal-05379794

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

Submitted on 24 Nov 2025

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>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH

Open Access

<figure><img src="image_2.png" /></figure>

Check for
updates

# Outlier detection and removal in multifractal analysis of electrophysiological brain signals

Merlin Dumeur<fnref n="1" />,<fnref n="2" />,<fnref n="3" />*, J. Matias Palva<fnref n="3" />,<fnref n="4" />,<fnref n="5" /> and Philippe Ciuciu<fnref n="1" />,<fnref n="2" />

**Abstract**

The analysis of the scale-invariant properties of brain activity has received growing attention in the last 15 years, especially in the context of the brain criticality hypothesis. Multifractal analysis (MFA) extends the ability of the usual tools in brain criticality research (e.g., detrended fluctuation analysis) to characterize scale invariance or scale-free dynamics by investigating the scaling properties of statistical moments beyond second-order analysis. However, MFA is very sensitive to the presence of outliers in the recorded time series, for instance, in the form of impulsive noise. In this work, we propose a novel algorithm for outlier detection and removal, which relies on the wavelet p-leader MFA formalism. The proposed approach consists in temporally segmenting the time-scale representation of the data (time series), and then tagging and removing the outlier segments from the final MFA estimation process. We first compare on synthetic noisy mono- and multifractal time series the performances of our method with those of existing techniques such as the wavelet p-leader MFA formalism and a naive robust cumulant-based approach. As a proof of concept, the proposed method is then evaluated on a single noisy resting-state magnetoencephalography (MEG) recording of a healthy volunteer. Overall, our results show that in the presence of noise impulses, this method removes the existing bias with standard MFA tools (wavelet leader-based analysis) on both synthetic and neurophysiological signals, with consistent findings across the two brain hemispheres in MEG.

**Keywords:** Multifractal analysis, Scale invariance, Outlier detection, MEG

# 1 Introduction

Scale-free time series have their statistical properties preserved across multiple temporal scales. This means that statistics of the time series considered at different time scales are related to one another via a power-law relation. In practice, given a time series $X(t)$, the power-law scaling is related to the concept of *pointwise regularity h*: A given scaling relation at the point $t_0$ bounds the fluctuations of $X$ around $t_0$ by a polynomial of order $h(t_0)$.

On the one hand, time series in which the pointwise regularity $h(t)$ is constant *over time* are monofractal. Those are characterized by a single parameter, the so-called Hurst exponent $H$ with $h(t) = H$ for all $t$. The associated power-law scaling manifests itself in the Fourier power spectrum $\Gamma(f) \propto 1/f^\beta$, which in a log-log representation has the characteristic linear trend with slope $-\beta$. The immediate consequence is that all

<figure><img src="image_3.png" /></figure>

<footnote>*Correspondence:
merlin@dumeur.net</footnote>

<footnote n="1">¹ NeuroSpin, CEA, DRF,
Joliot, Paris-Saclay University,
91191 Gif-sur-Yvette, France</footnote>

<footnote n="2">² MIND Inria, Paris-Saclay
University, 1 Rue Honoré
d'Estienne d'Orves,
91120 Palaiseau, France</footnote>

<footnote n="3">³ NBE, Aalto University, Otakaari
3, 02150 Espoo, Finland</footnote>

<footnote n="4">⁴ Neuroscience Center, Helsinki
University, Haartmaninkatu 3B,
00014 Helsinki, Finland</footnote>

<footnote n="5">⁵ Center for Cognitive
Neuroimaging, Institute
of Neuroscience and Psychology,
University of Glasgow,
Glasgow G12 8QQ, UK</footnote>
scale-invariant properties of the time series can be determined from the second-order statistics of the time series, either in the time domain (e.g., with detrended fluctuation analysis (DFA) [1]), in the spectral domain (e.g. with Welch's periodogram) or using a time-scale representation (e.g. dyadic wavelet decomposition) [2–6]. Typically, for a fractional Gaussian noise, $H = (\beta + 1)/2$ [7]. Fractal methods find a wide spectrum of applications [8–10]. For a fair comparison of the performance of different techniques in the neuroimaging context, the reader can refer to [11].

On the other hand, for the *multifractal* time series, $h(t)$ takes different values for different time points, which means that a single parameter is no longer enough to fully characterize the scaling-invariant properties. Instead, higher-order moments need to be investigated: This is the role of *multifractal analysis*.

Currently in neuroscience, the scale invariance of brain activity has been primarily interpreted within the framework of *brain criticality*: The idea is that the brain can be understood as a dynamical system operating near the critical point of a synchronization phase transition [12–14]. Critical systems typically display scale invariance, which is usually measured using tools such as DFA [15]. However, these methods only measure the self-similarity properties of the time series and are blind to the fluctuations in regularity that are characteristic of multifractality.

Although the neuroscience community working on the brain criticality hypothesis has paid attention to the low-frequency fluctuations of neuronal oscillations [15], the multifractal literature in neuroscience has focused on characterizing the low-frequency fluctuations of the broadband recordings of brain activity [16–20]. In general, MFA as a tool in neuroscience has seen relatively few developments [21]. A possible explanation is that the inherent need to leverage higher-order statistics in MFA makes it very sensitive to outliers in the data [22, 23].

The present work is motivated by multiple and consistent observations in neuroscience, where various noise sources corrupt electrophysiological and functional neuroimaging recordings of brain activity. For instance, in magnetoencephalography (MEG), muscle activation, eye blinks, movement, flux jumps from the SQUID sensors, evoked responses to large stimuli are all known to generate high-powered artifacts in the recorded time series. Those high-power artifacts constitute *impulsive* noise sources that have a small temporal support compared to the length of the recorded signals, but that significantly impede the reliability of MFA on brain signals.

## 1.1 Problem statement

Multifractal analysis as an analytic tool was originally developed and evaluated on "pure" time series, with no impulsive noise contaminating the data. In practice, impulsive noise is not accounted for in MFA frameworks, even when applied to neurological time series analysis. Outliers have been shown to have a large influence on MFA, with examples in finance [24, 25] and on response time series [23]. Denoising techniques (e.g., ICA [26]) specific to each recording modality (e.g., EEG, MEG, fMRI) exist; however, they are not designed with MFA in mind. In electrophysiology (EEG and MEG mainly), the standard analysis pipelines handle noise in the context of *epoched* data: Recordings from repeated trials of the same experimental conditions are averaged together to cancel out measurement noise, and individual noisy trials can easily be discarded by simple thresholding.
However, estimates of multifractality (or scale-free dynamics in the wide sense) require analyzing the statistics of the *whole* time series at multiple large time scales, which means that the input data are the whole uninterrupted recorded signals instead of epoched trials lasting a few seconds.

Visual inspection and delimitation of noisy segments in time series could be consid-
ered, but since noise may be present at many different time scales, it is not feasible for a
human observer to be perfectly accurate, and what may appear upon visual inspection as
barely significant noise may actually result in significant perturbation of the final scale-
free or multifractal estimates. An additional difficulty lies in the fact that multifractal
signals can naturally exhibit fluctuations of the same order of magnitude as impulsive
noise, making the decision whether or not to remove a segment an arbitrary decision.
Furthermore, human-based (i.e., not automated) outlier removal does not scale up to
large data sets, which are becoming increasingly common in neuroscience.

1.2 Related works

Prior contributions for localizing outliers in scale-invariant time series were developed some time ago. For example, [24] uses MFA based on the continuous wavelet transform to identify local singularities in financial time series; [27] finds outliers in chemical data where the box-counting dimension of curves varies heavily with a cross-validation scheme. [28] removes high-frequency noise from a regular function using low-frequency regularity estimates. However, the aforementioned methods were only aiming to detect outliers in time series and do not attempt to correct the scale-free estimates for those outliers. More recently, [29] describes finding outliers but does not describe how they were removed. [30] describes the presence of outliers, but there was no attempt to remove them from the data.

1.3 Outlier detection for multifractal analysis

In a previous contribution [31], we proposed using an iterative Monte Carlo approach to identify outlier wavelet coefficients in the wavelet leader-based MFA [32, 33]. However, that method was limited in its ability to mitigate noise: The final estimates of multifractality were still biased after correction. Additionally, the method was too slow for practical use on large data sets, as Monte Carlo sampling was used at the core of the algorithm to draw a large number of samples in order to discriminate rare events. Apart from [31], two recent contributions proposed a robust estimate of the multifractal parameter: The first [34] assumed a specific Gaussian mixture for the Fourier transform of the multi-scale log-leaders. This assumption was essential to derive a robust EM algorithm for MFA. However, it is not necessarily tenable in the present context, where we are looking for techniques that may identify the impact of white noise impulses (as opposed to a continuous white noise in [34]) on the estimation of multifractal properties of time series. The second [35] uses a projection statistic to compute a threshold, which combined with a blind source separation method allows for the removal of outliers, although this time the impact of the outlier removal procedure on the estimated multifractal exponents is not evaluated.
## 1.4 Approach

In this contribution, we first simulate mono- and multifractal time series with known parameters and evaluate the impact of additive impulsive noise to the multifractal exponent estimates (Fig. 1). Secondly, we compare the performance of a naive approach to robust MFA with the proposed detection algorithm based on the multifractal p-leader framework [36, 37], over multiple simulated scenarios with increasing power in the noise impulses. Thirdly, we apply the detection algorithm to resting-state MEG recordings of brain activity in a single healthy volunteer and show that it regularizes the distribution of multifractal estimates in a consistent manner. In addition, we highlight the reproducible and consistent impact of the outlier detection/removal strategy in well localized brain regions (occipital and parietal areas) over both hemispheres and cortical (medial and lateral) surfaces. Applying this algorithm always induces the same type of correction, namely (lower self-similarity) and stronger multifractality when it exists.

# 2 Methods

## 2.1 Multifractal analysis

Let $X(t)$ be a time series. A multi-scale statistic $T_X(j, k)$ is a description of the statistical properties of $X$ on different temporal scales $2^j$ and temporal shifts $k(t)$. Formally, given a scale-invariant time series $X(t)$, its multiscale statistics considered at different time scales follow a power law [37]: For each generalized moment $q \in \mathbb{R}$:

<figure><img src="image_4.png" /><figcaption>**Fig. 1** Wavelet p-leader multifractal analysis. A Example simulated MRW time series, with theoretical exponents $c_1 = 0.825$, $c_2 = -0.05$. Clean (left) and noisy (right) with added white noise impulses. B Absolute value of the wavelet coefficients (top row) and wavelet 2-leaders (bottom row) for both the clean (left column) and noisy (right column) time series. The light blue rectangle in the top left plot represents the support upon which is computed the 2-leader in the same colored rectangle in the bottom left plot. C Fit of the log-cumulants $c_1$ (top) and $c_2$ (bottom) for both the clean (blue) and noisy (orange) simulations. The top plot has the noisy $c_1$ fit offset by 0.3 for clarity. D Multifractal spectra of the clean (blue) and noisy (orange) simulations</figcaption></figure>
$$ \frac{1}{N_j} \sum_{k=1}^{N_j} |T_X(j, k)|^q \propto 2^{j\zeta(q)}. \qquad (1) $$

For each $q$, the scaling function $\zeta(q)$ thus describes the scaling of $T_X$ as a function of $j$. Regularity is properly defined around the time $t$ as the limit inferior scaling $h$ of $T_X$ as scales get progressively smaller ($j \to -\infty$ represents infinitely small temporal scales—$j = 0$ corresponds to the sampling frequency in observed time series):

$$ h(t) = \liminf_{j \to -\infty} j^{-1} \log T_X(j, k(t)). \qquad (2) $$

Therefore, monofractal time series as introduced above have a function $\zeta(q)$ that depends on a single parameter $H: \zeta(q) = qH$ and $\forall t, h(t) = H$. On the other hand, the function $\zeta(q)$ is concave for multifractal time series: A single parameter is not enough to describe its properties.

In order to characterize the regularities present in a time series, it is not desirable to investigate the regularity near every time point $t$. In practice, we do not have access to infinitely fine time scales, but instead we can aim to characterize the degree of prevalence of each regularity within the time series. A useful characterization is the multifractal spectrum $D(h')$ (also known as the singularity spectrum), which relates to $h'$ the Hausdorff dimension of the set of points for which the pointwise regularity of $X$ is $h'$. We obtain from the Endre transform of $\zeta(q)$ an often sharp upper bound on $D(h')$:

$$ D(h') = \dim_H(\{t: h(t) = h'\}) \le \inf_{q \in \mathbb{R}} (1 + qh - \zeta(q)). $$

## 2.2 Wavelet p-leader multifractal formalism

This formalism forms the basis for our outlier detection algorithm; however, other multifractal formalisms may be chosen [38] depending on the application.

### 2.2.1 Wavelet p-leaders

From (2), it follows that the specific meaning of regularity $h$ depends on the choice of the multiscale quantity $T_X(j, k)$. In this contribution, we use the wavelet $p$-leader formalism [37], which is based on the $p$-exponent operationalization of regularity [36]. This choice is motivated by the presence of a relevant indirect characterization of the properties of the multifractal spectrum via the log-cumulants upon which the method is built. For completeness, this characterization is summarized below.

Let $(d_{j,k}(X))_{j,k}$ be the discrete wavelet coefficient of a signal $X$ at scale $j$ and temporal shift $k$, the wavelet $p$-leader of $X$ at the same $(j, k)$ indices is defined as

$$ \ell_{j,k}^{(p)}(X) := \left( \sum_{j' \in [1, j]} \sum_{k' \in \lambda_{j,j'}(k, r)} |d_{j',k'}(X)|^p 2^{j-j'} \right)^{1/p}, \qquad (3) $$
where $\lambda_{j,j'}(k, r)$ is the ball of radius $r \times 2^{j-j'}$ centered at $k2^{j-j'}$. Note that $r$ is an odd integer, commonly $r = 3$.

### 2.2.2 Multifractal formalism

It is possible to express the scaling function $\zeta(q)$ via its cumulant expansion

$$ \zeta(q) = c_1 q + c_2 \frac{q^2}{2} + c_3 \frac{q^3}{3!} + \sum_{k=4}^{+\infty} c_k \frac{q^k}{k!}, \quad (4) $$

where the log-cumulants $(c_m)_m$ determine the properties of the multifractal spectrum. $c_1$ locates the mode of the spectrum, $c_2$ the spread of the spectrum, $c_3$ the asymmetry of the spectrum, and $c_4$ the tail behavior, and so on for higher $m$ values by analogy with the interpretation of higher-order cumulants.

The $(c_m)_m$ themselves are estimated by linear regression on the cumulants of order $m$, $C_m(j)$, of the log $p$-leaders:

$$ C_m(j) = \text{Cum}_m \left[ \left( \log \ell_{j,k}^{(p)} \right)_k \right] = \hat{c}_m^0 + j \frac{\hat{c}_m}{\log_2 e} + \varepsilon_j. \quad (5) $$

### 2.2.3 Wavelet scaling function

To make multifractal analysis a reliable tool, the analyzed time series needs to be regular enough in order for the above equations to hold, and for the final estimates of the $(c_m)_m$ to be correct. Formally, within the $p$-leader formalism, this requirement is that the wavelet scaling function $\eta(p)$ must be positive:

$$ \eta(p) = \liminf_{j \to -\infty} j^{-1} \log_2 2^j \sum_k |d_{j,k}|^p. \quad (6) $$

If $\eta(p) \le 0$, it is necessary to decrease $p$ or to fractionally integrate the time series with the fractional integration coefficient $\gamma_{int}$ by substituting $d_{j,k} \to 2^{j\gamma_{int}} d_{j,k}$ when computing the $p$-leaders. In brain recordings such as EEG and MEG data, time series are frequently so irregular that a value of $\gamma_{int} > 1$ is often required. Furthermore, it may vary in space between sensors or parcels on the cortical surface, further complicating the finding of a single operational value of $\gamma_{int} > 1$ over throughout the brain.

In what follows, we present two approaches to robustly estimate multifractal parameters: robust cumulants in Section 2.4, and a segmentation approach in Sect. 2.5.

## 2.3 Simulation of noisy time series

We compare the impact of additive impulsive white noise on two scale-invariant signal models: the fractional Brownian motion (fBM), which is monofractal, and the multifractal random walk (MRW) [39]. Other possibilities could include using the Weierstrass function [40].

The fBM is parameterized by its Hurst exponent $H$. In contrast, the MRW is parameterized by $H$ and an additional intermittency parameter $\lambda$, which controls the multifractality of the time series and consequently the $c_2$ log-cumulant.
Monte Carlo simulations with $N = 2^{18}$ samples are produced for both fBM and MRW models. Subsequently, additive white noise is superimposed in 3 impulses that are evenly distributed through the time series. Each noise impulse covers 1% of the duration of the whole time series. The noise level is defined by the signal-to-noise ratio (SNR), which is given here as the signal power at scale 9 over the noise power.

## 2.4 Robust cumulants

A simple approach to robust estimation is to consider replacing the k-statistics used to estimate the cumulants $C_m(j)$, with a robust estimator. The k-statistic is the minimal variance unbiased estimator of the cumulants of a distribution [41], but it is sensitive to the presence of outliers in the data. As the first- and second-order cumulants are identical to the first-order and central second-order moments, respectively, one can replace estimators of the first two cumulants with robust estimators of the mean and variance. Higher-order cumulants naturally place more emphasis on the most outlying values, where it becomes critical to properly distinguish between outliers and inliers, i.e., regular samples.

In the robust cumulant approach, we replace the sample variance estimator with the $Q_n$ estimator of scale [42], which is based on rank-order statistics of the data. For a vector $\mathbf{x}$ of random variables of size $n$, it is defined as the k-th value of the ordered set of pairs of absolute differences:

$$Q_n(\mathbf{x}) = 2.2191 \{ |x_i - x_{i'}| \mid i < i' \le n \}_{(k)}, \quad k = \left\lfloor \frac{n}{2} \right\rfloor + 1. \qquad (7)$$

The sample mean estimator was replaced with the M-estimator of the location parameter $\mu$, using Tukey's biweight [43] as the weighting function $\psi$. Using the median as an initial guess, its objective is to iteratively find the solution to the following problem:

$$\hat{\mu} = \arg \min_{\mu \in \mathbb{R}} \sum_{i} \psi((x_i - \mu)/Q_n(\mathbf{x})). \qquad (8)$$

Tukey's biweight was chosen for the $\psi$ function due to its redescending property, meaning it allocates zero weight to extreme observations, which is useful in scenarios where the noise power is high.

Higher-order cumulants may be replaced with robust estimators of higher-order moments [44] although it is not investigated here.

## 2.5 Outlier detection and removal algorithm

The proposed algorithm is presented in Algorithm 1. It consists in first estimating the multifractal properties of the signal using 1-leaders computed with $r = 1$. Those multifractal estimates provide the basis for a mapping of the 1-leaders at every scale $j \in [j_1, j_2]$ to a common uniform distribution on $[0, 1]$, at a temporal resolution matching the resolution of the finest temporal scale $j_1$. The mapped data are then temporally segmented using an original distance metric which aims to separate high-power 1-leader outliers from the rest of the data. At each temporal scale $j$, each segment is then compared to the rest of the data, and outlier segments are identified based
on the dissimilarity of the distribution of their values to the rest of the time series, as measured by the Wasserstein metric. The output set of vectors indicating outlier 1-leaders at each scale *j* is built from the segments of data designated as outliers and then returned. Consequently, the SEGMENTATION method takes the 1-leaders and returns the segmented mapped data.

Then, the algorithm DETECTIONREMOVAL takes the segmented mapped data as input and returns the set of indices for each scale *j*, for which the associated *p*-leaders are considered outliers. Figure 2 illustrates the successive steps of the overall procedure. The proposed outlier detection approach is presented on a noisy MRW case in the context of an inverse SNR of 0.5, as it comprises the more difficult case to handle due to natural fluctuations in the simulated signal's power.

Two minor helper functions are introduced in Algorithm 2. First, COSTFUNC-TION implements the cost function upon which the segmentation is based. Second,

<figure><img src="image_5.png" /><figcaption>**Fig. 2** Detection algorithm steps. **A** Simulated time series (in blue) with added noise impulses (in orange). **B** Wavelet 1-leader decomposition of the time series, with *r* = 1. **C** Estimation via linear regression of the log-cumulants for the clean (blue) and noisy (orange) version of the time series. **D** Histogram of the 1-leaders at the scale *j* = 8. Leaders which were computed from the noisy segments are colored in orange, while the cleaned ones appear in blue. The black solid line indicates the distribution obtained via cumulant matching based on the $(\hat{c}_m)_{m \text{ estimated in panel C}}$. **E** Probability integral transformed 1-leader matrix **Y**, taking values between 0 and 1. **F** Pairwise distances (Eq. (17)) between the 1-leader vectors for each pair of temporal shifts $k_x, k_y \in [1, N_j]$. **G** Segmentation of the time series using $\beta_{PELT} = 50$. **H** Fit of the log-cumulants both without (orange) and with (green) removal of the detected outliers from the analysis. **I** 2-leader representation of the time series (*r* = 2). The outlier 2-leaders were correctly identified by the procedure, and are outlined in blue.</figcaption></figure>
REJECTLEADERS marks as outlier all the 1-leaders that need to be removed, given a segment that was detected as outlier. MATCHCDF takes the multifractal estimates, the chosen model for the statistical distribution of log-leaders, the scale j, and returns the modeled cumulative distribution of the 1-leaders at that scale j.

In the text below, $Y$ denotes the matrix of data mapped to $[0, 1]$, of which $Y_{j,*}$ is the jth row and $Y_{*,k}$ the kth column. $1_n$ is the vector of ones of size $n$. The notations relevant to each step are introduced and detailed in their corresponding paragraphs hereafter. Additionally, the successive steps of Algorithm 1 are described in a sequential and appearance order in what follows.

### 2.5.1 Initial multifractal estimate

As $p \to +\infty$, the value of the p-leaders depends increasingly on the largest-valued wavelet coefficients irrespective of the scale j. Similarly, increasing values of r smooth out the resulting p-leaders in the time dimension, spreading the influence of a single wavelet coefficients over r different p-leaders. For outlier detection, it is desirable that the majority of the variance of $\ell_{j,k}^{(p)}$ depends on the coefficient $d_{j,k}$. We choose thus $p = 1$ and $r = 1$ for the detection part and keep the usual values of $p = 2$ and $r = 3$ for the final estimation.

We note $\ell_j^{(1)}$ the initial 1-leader representation of the input time series at scale j. The initial fit of the scaling exponents is noted as $(\hat{c}_m)_m$ for the slope and $(\hat{c}_m^0)$ for the intercept.

### 2.5.2 Probability integral transform

Given a real-valued random variable $V$ with cumulative distribution function $F: \mathbb{R} \to [0, 1]$, we have $F(V) \sim U(0, 1)$ [45]. This property means that it is possible to map the variables belonging to different distributions to a common space via this process known as the *Probability Integral Transform*.

Previous literature supports modeling the log p-leaders as normally distributed: $\log \ell_j^{(1)} \sim N(\mu_j, \sigma_j)$ [46]. In this contribution, we compare that approach with the modeling of the log p-leaders according to an exponential power distribution: $\log \ell_j^{(1)} \sim EP(\mu_j, \alpha_j, \beta_j)$.

The proposed approach consists in using the $(\hat{c}_m)_m$ and $(\hat{c}_m^0)_m$ to determine the parameters of the distribution of the log p-leaders at each temporal scale. The fit of the linear regression (5) (Fig. 2C, D) for a given cumulant order $m$ gives a predicted value $\tilde{C}_m(j)$ (as opposed to the observed $\hat{C}_m(j)$) of the cumulants of the 1-leaders at each scale j: $\tilde{C}_m(j) = \frac{\hat{c}_m}{\log_2 e}j + \hat{c}_m^0$. From those cumulants, the parameters of the modeling distribution are determined via cumulant matching.

In the case of the normal model for the log 1-leaders, we thus have:

$$ \mu_j = \tilde{C}_1(j) = \frac{\hat{c}_1}{\log_2 e}j + \hat{c}_1^0, \qquad (9) $$

$$ \sigma_j = \sqrt{\tilde{C}_2(j)} = \sqrt{\frac{\hat{c}_2}{\log_2 e}j + \hat{c}_2^0}, \qquad (10) $$

where the cumulative distribution is defined as $F: x \mapsto 2^{-1}(1 + erf(x))$.
In the case of the exponential power model for the log 1-leaders, $\mu_j$ is the same as in (9), and the other parameters are determined as follows:

$$f : \beta, j \mapsto \frac{\Gamma(5/\beta)\Gamma(1/\beta)}{\Gamma(3/\beta)^2} - 3 - \frac{\hat{c}_4}{\log_2 e}j - \hat{c}_4^0 \quad (11)$$

$$\beta_j = \operatorname*{root}_{\beta \in \mathbb{R}_+} (f(\beta, j)), \quad (12)$$

$$\alpha_j = \sqrt{\left(\frac{\hat{c}_2}{\log_2 e}j + \hat{c}_2^0\right) \frac{\Gamma(1/\beta_j)}{\Gamma(3/\beta_j)}} \quad (13)$$

Since there is no analytic formula for $\beta_j$, it is estimated via finding the unique root of function $f(\beta, j)$ on $\mathbb{R}_+^*$. The cumulative distribution function is then written as:

$$F : x \mapsto \frac{1}{2} \left( 1 + \operatorname{sign}(x - \mu_j) \gamma \left( 1/\beta_j, \left| \frac{x - \mu_j}{\alpha_j} \right|^{\beta_j} \right) \right), \quad (14)$$

where $\gamma(\cdot, \cdot)$ is the regularized lower incomplete gamma function.

Hereafter, we note $\mathbf{Y}_j = F(\log \ell_j^{(1)})$ the probability integral transformed log 1-leaders at scale $j$.

### 2.5.3 Matching temporal scales

In order to merge the information over multiple scales, for $j \in [j_1 + 1, j_2]$, we upsample $(\mathbf{Y}_{j,k})_k$ via linear interpolation. This allows us to match the number of coefficients $N_{j1}$ at scale $j_1$, which has the highest number of coefficients among the scales used. The upsampling operation allows for the fusion of all information at the resolution of the scale $j_1$ and allows both aggregation of information across all scales, as well as fine-grained temporal segmentation later on. Figure 2 shows the $\mathbf{Y}$ matrix.

### 2.5.4 Signal segmentation

There are multiple approaches to performing signal segmentation, as outlined in [47]. Optimal segmentation algorithms solve an optimization problem by minimizing an objective function. In the case of an unknown number of change points as is here, the optimal “pruned exact linear time” (PELT) algorithm [48] suits our needs.

Formally, the PELT algorithm finds the partition $(I_n)_n$ of the interval $[1, N_{j1}]$ into $N$ segments, which minimizes the objective function:

$$\mathrm{PELT}(\mathbf{Y}, \beta) = \arg \min_{(I_n)_n} \beta(N-1) + \sum_{n=1}^{N} c(\mathbf{Y}_{*,I_n}), \quad (15)$$

where $\mathbf{Y}_{*,I_n}$ denotes the temporal slicing associated with the segment $I_n$.

The objective function balances the cost function $c$ that associates a cost to each segment, with the penalization term in $\beta$. The linearity of the penalization term allows the PELT algorithm to find the optimal partition after iterating over all samples in $[1, N_{j1}]$
only once. $\beta$ affects the number of detected change points: Higher values of $\beta$ reduce the number of detected change points, while lower values of $\beta$ increase it.

In general, distance metrics $d(\cdot, \cdot)$ provide a natural cost function $c$:

$$c(\mathbf{Y}_{*,I_n}) = \sum_{k \in I_n} d(\mathbf{Y}_{*,k}, \bar{\mu}_{I_n})^2, \quad \bar{\mu}_{I_n} = \prod_{\alpha \in I_n} Y_{*,\alpha}^{1/2}. \qquad (16)$$

In the following section, we introduce an original distance metric (17) designed to highlight outliers. We used the implementation of PELT in the python package 'ruptures' [47] to perform the segmentation.

### 2.5.5 Distance metric

In order to segment $\mathbf{Y}$ temporally, we need a distance metric between pairs of columns $\mathbf{Y}_{*,k_1}, \mathbf{Y}_{*,k_2}$ that produces high values when one of the samples is an outlier and the other is not, and low values otherwise.

Intuitively, outliers manifest as high-power events that will cause affected coefficients $\mathbf{Y}_{*,k} \in (0, 1)^s$ to tend toward $\mathbf{1}_s$, where $s = j_2 - j_1 + 1$ is the number of temporal scales. Close to the limit, the absolute difference between the samples tends to 0, which renders the usual distance metrics meaningless. Instead, the ratios between the values are meaningful, motivating the use of the Hilbert metric (details in Appendix A), which generalizes hyperbolic geometry to any convex subspace in a finite-dimensional Euclidean space.

The Hilbert metric is defined in a way that does not make the difference between high-power (outliers) and low-power events. However, in practice the outliers in the data tend to push the $\mathbf{Y}_{*,k}$ associated with low-power (high regularity) fluctuations toward zero. This means that in the presence of high-power outliers, the Hilbert metric ends up treating regular data as low-power outliers. Thus, we adjusted the Hilbert metric not to compute the distance from the samples to points on the boundary of the space $(0, 1)^s$, as that affects the low-power events equally as the high-power outliers, but instead to the point $\mathbf{1}_s$, to where only the high-power outliers are close.

The resulting distance metric is defined as follows:

$$d(\mathbf{x}, \mathbf{y}) = \log \left[ \left( 1 + \frac{\|\mathbf{x} - \mathbf{y}\|_2}{\|\mathbf{x} - \mathbf{1}_s\|_2} \right) \left( 1 + \frac{\|\mathbf{x} - \mathbf{y}\|_2}{\|\mathbf{y} - \mathbf{1}_s\|_2} \right) \right], \quad \mathbf{x}, \mathbf{y} \in \mathbb{R}^s \times \mathbb{R}^s. \qquad (17)$$

Further details on the Hilbert metric and the derivation of the above formula are available in Appendix A.

### 2.5.6 Outlier detection and removal

The 1-dimensional Wasserstein distance $W$ provides a way to measure a distance between two distributions and has a fast closed-form implementation for empirical distributions. For each scale $j$, the segments $(\mathbf{Y}_{j,I_n})_n$ are classified as outliers based on the Wasserstein distance of the empirical cumulative distribution $\pi_{j,I_n}$ of $\log[1 - \mathbf{Y}_{j,I_n}]$, to the empirical cumulative distribution $\pi_{j,I_n^C}$ of $\log[1 - \mathbf{Y}_{j,I_n^C}]$ over the complement of the segment, $I_n^C = [1, N_j] \setminus I_n$.

The Wasserstein distance between those two distributions is then expressed as:
$$W(\pi_{j,l_n}, \pi_{j,l_n^C}) = \int_{-\infty}^{0} |\pi_{j,l_n}(x) - \pi_{j,l_n^C}(x)| dx. \quad (18)$$

The user-defined threshold for identifying outliers is termed $\theta$. At each scale $j$, the segments $I_n$ for which $W(\pi_{j,l_n}, \pi_{j,l_n^C}) > \theta$ and median($Y_{j,l_n}$) > 0.75 are considered outlier segments. The outlier segments detected at each $j$ are removed from the subsequent final MFA step. Since the $p$-leaders at a given scale are computed using coefficients at lower scales, when removing $p$-leaders from the analysis at scale $j$, it is necessary to remove coefficients at coarser scales $j' > j$ which depend on the removed $p$-leaders at scale $j$. Further, we decided to also remove the coefficients at finer scales $j' < j$ that build into the outlier $p$-leaders at scale $j$ as it helps to improve the identification of outliers in finer scales, and it carries a comparatively small cost, given the exponentially larger number of coefficients at finer scales.

Coefficients close to the edge of the noise impulses may not have extreme enough values to be detected as outliers and may not be part of the detected segments, but still continue to affect the final estimates. Moreover, since the final multifractal estimate is computed using $r = 3$, for each rejected segment of $p$-leaders we need to additionally reject the two $p$-leaders to the left and to the right of the segment.

**Algorithm 1** Main outlier detection/removal algorithm. The procedure OUTLIERDETECTION takes several input arguments, the time series to be analyzed $x$, the penalty term $\beta_{PELT}$ for signal segmentation, the range $[j_1, j_2]$ of temporal scales which constitutes the support of the linear regression used to estimate the log-cumulants $(c_m)_m$, the threshold $\theta$ used to separate outlier segments from the rest of the data, and finally the model of the probability distribution placed on the log-leaders: either normal or exponential power distribution.
**procedure** SEGMENTATION($\ell^{(1)}$, $\beta_{PELT}$, $j_1$, $j_2$, model)
($\hat{c}_m)_m$, ($\hat{c}_m^0$)$_m$ ← MFA($\ell^{(1)}$, $j_1$, $j_2$)

$$N_{j_1} \leftarrow |\ell_{j_1}^{(1)}|$$

5: $\mathbf{Y} \leftarrow 0 \in \mathcal{M}_{j_2-j_1 \times N_{j_1}}(\mathbb{R})$

for $j \in [j_1, j_2]$ do

$$n_j \leftarrow j - j_1$$
$$N_j \leftarrow |\ell_j^{(1)}|$$

$$F \leftarrow \text{MATCHCDF}(\text{model}, (\hat{c}_m)_m, (\hat{c}_m^0)_m, j)$$

10: $\mathbf{y} \leftarrow F(\log \ell_j^{(1)}) \quad \triangleright \text{ Probability integral transform}$
$\mathbf{Y}_{j,*} \leftarrow \text{Interpolate}(\mathbf{y}, (k + 0.5)_{k \in [0, N_{j-1}]}, (2^{-n_{j-1}} + k2^{-n_j})_{k \in [0, N_{j-1}2^{n_j}]})$

end for

$(I_n)_n \leftarrow \text{PELT}(\mathbf{Y}, \text{COSTFUNCTION}, \beta_{PELT}) \quad \triangleright \text{ Set of intervals}$

return $\mathbf{Y}, (I_n)_n$

15: end procedure

**procedure** DETECTIONREMOVAL($\mathbf{Y}$, $(I_n)_n$, $\theta$)
($\mathbf{o}_j)_j \leftarrow \{\{\text{True}, \text{False}\}^{N_j} | j \in [j_1, j_2]\} \quad \triangleright \text{ Set of boolean vectors}$
for $j \in [j_1, j_2]$ do

20: for $n \in [1, |(I_n)_n|]$ do
$I_n^C \leftarrow [0, N_j] \setminus I_n$
$\pi_{j, I_n} \leftarrow \text{EmpiricalCDF}(\log[1 - \mathbf{Y}_{j, I_n}])$
$\pi_{j, I_n^C} \leftarrow \text{EmpiricalCDF}(\log[1 - \mathbf{Y}_{j, I_n^C}])$
$w \leftarrow W(\pi_{j, I_n}, \pi_{j, I_n^C}) \quad \triangleright \text{ Wasserstein distance}$

25: if $w > \theta$ & median($\mathbf{Y}_{j, I_n}$) > 0.75 then
($\mathbf{o}_j$)$_j \leftarrow \text{REJECTLEADERS}((\mathbf{o}_j)_j, j, I_n, j_1, j_2)$
end if

end for

end for

30: return $(\mathbf{o}_j)_j$

end procedure

**procedure** OUTLIERDETECTION($\mathbf{x}$, $\beta_{PELT}$, $j_1$, $j_2$, $\theta$, model)
$\ell^{(1)} \leftarrow \text{WAVELET p-LEADER}(\mathbf{x}, p = 1, r = 1)$

35: $\mathbf{Y}, (I_n)_n \leftarrow \text{SEGMENTATION}(\ell^{(1)}, \beta_{PELT}, j_1, j_2, \text{model})$
return $(\mathbf{o}_j)_j \leftarrow \text{DETECTIONREMOVAL}(\mathbf{Y}, (I_n)_n, \theta)$

end procedure
<figure><img src="image_6.png" /><figcaption>**Fig. 3** Simulation results. A: fBM simulations. Average $c_1$ (left) and $c_2$ (right) estimates under increasing noise levels. B: MRW simulations. MRW simulation: average $c_1$ (left) and $c_2$ (right) estimates under increasing noise levels. Error bars denote empirical confidence intervals at 90%. Dashed black line: naive estimates at zero noise level; dashed line: robust cumulants approach; dotted line: proposed approach with log-normal model on 1-leaders; dashed-dotted line: proposed approach with log-exponential power model on 1-leaders</figcaption></figure>

simulation-based analyses covering a broader range of scenarios, notably one based on
real and clean MEG data, are detailed in Appendices B and C, respectively.

**3.1 Monofractal fBM**

Benchmarking the multifractal estimates of several competing approaches in a pure self-
similar process such as fBM is informative, as it allows us to realize that both $c_1$ and $c_2$
are affected.

In essence, it is worth noting that the original uncorrected multifractal analysis yields $c_1$ and $c_2$ estimates that diverge for increasing noise levels. Naive $c_1$ estimates remain accurate for low noise levels, while they totally overestimate the amount of self-similarity beyond a certain noise level. Additionally, naive $c_2$ estimates first make the fBM appear multifractal ($c_2 < 0$) for lower noise levels, before turning monofractal ($c_2 \ge 0$) for higher noise levels. In between, there is a point for which the noisy signal has the same $c_2$ exponent as the clean signal. This leads to a situation for which comparing $c_1$ and $c_2$ exponents between a clean and a noisy time series could lead to the same $c_2$ exponents but different $c_1$ exponents, which could be misleading. Strikingly, spurious multifractality can appear from impulsive noise that constitutes only a fraction of the temporal support—which warrants the use of denoising.

The robust cumulant approach improves the estimates, particularly at higher noise lev-
els, but remains biased. In fact, we can observe that it reduces the degree of polynomial
growth of the error: Beyond a certain noise level, naive $c_1$ estimates linearly diverge from the true value with a linear trend, and the robust cumulants yield constant estimates. Similarly for $c_2$, the naive estimates diverge quadratically when the robust cumulants exhibit a linear trend. The proposed approach outperforms the state-of-the-art methods using either probability distribution for the wavelet leaders. More precisely, the confidence intervals of the parameter estimates of $c_1$ and $c_2$ always include the true values regardless of the noise level, illustrating the robustness of the proposed method.

Table 1 presents the average relative error $e_{c_m}$ as defined in (19): the reference values $c_m^{\text{ref}}$ are the naive estimate at zero noise level.

$$e_{c_m} = \frac{\hat{c}_1 - c_m^{\text{ref}}}{c_m^{\text{ref}}} \qquad (19)$$

## 3.2 Multifractal random walk

Compared to the fBM results, the impact of noise appears more detrimental to the estimation of the multifractal properties of MRWs: The noise now strongly positively biases the naive $c_1$ estimates. The relationship between naive and robust cumulant estimates that was observed on fBM simulation also holds here, with a higher bias for the robust cumulants which is likely due to the power fluctuations in MRWs making the estimation process more difficult. In addition and similar to previous observations on fBMs, $c_2$ estimates are more biased in a way that causes the apparent multifractality to be stronger when estimated by the naive or robust cumulant method.

The proposed method aligns with the expected values, for values of inverse SNR above 0.125. The low-noise power setting remains challenging for the detection method, as the natural fluctuations of the underlying signal are of magnitude similar to the noise impulses. The results also show a small negative bias in the zero-noise estimate of $c_1$.

We provide in Table 2 the average relative error, as defined above for Table 1, for the estimation of $c_1$ and $c_2$ on the noisy MRW simulations.

# 4 Experimental analysis

## 4.1 Data set

To showcase the applicability of the proposed approach to electrophysiological recordings, we apply the method to MEG signals collected on the scalp over 102 sensor

<table><caption>**Table 1** Relative error in estimating $c_1$ and $c_2$ on noisy fBM simulations, according to noise power</caption><thead><tr><th rowspan="2">SNR<sup>-1</sup></th><th colspan="3">c<sub>1</sub></th><th colspan="3">c<sub>2</sub></th></tr><tr><th>Naive</th><th>Robust C<sub>m</sub></th><th>Ours</th><th>Naive</th><th>Robust C<sub>m</sub></th><th>Ours</th></tr></thead><tbody><tr><td>0.000</td><td><strong>0.0000</strong></td><td>0.0014</td><td><strong>0.0000</strong></td><td><strong>0.00</strong></td><td>1.51</td><td><strong>0.00</strong></td></tr><tr><td>0.125</td><td>0.0071</td><td><strong>0.0040</strong></td><td>0.0067</td><td>9.00</td><td>2.14</td><td><strong>1.98</strong></td></tr><tr><td>0.250</td><td><strong>0.0052</strong></td><td>0.012</td><td>0.0067</td><td>19.3</td><td>2.59</td><td><strong>1.98</strong></td></tr><tr><td>0.500</td><td><strong>0.0056</strong></td><td>0.011</td><td>0.0067</td><td>25.9</td><td>5.53</td><td><strong>1.96</strong></td></tr><tr><td>1.000</td><td>0.018</td><td>0.0079</td><td><strong>0.0068</strong></td><td>21.3</td><td>9.23</td><td><strong>2.00</strong></td></tr><tr><td>2.000</td><td>0.039</td><td>0.0074</td><td><strong>0.0067</strong></td><td>10.4</td><td>13.4</td><td><strong>1.99</strong></td></tr><tr><td>4.000</td><td>0.065</td><td>0.0091</td><td><strong>0.0064</strong></td><td>40.0</td><td>18.5</td><td><strong>2.25</strong></td></tr><tr><td>8.000</td><td>0.097</td><td>0.010</td><td><strong>0.0068</strong></td><td>104</td><td>23.1</td><td><strong>2.20</strong></td></tr></tbody></table>

Bold values indicate the best estimate for a given noise level
<table><caption>**Table 2** Relative error in estimating c₁ and c₂ on noisy MRW simulations, according to noise power</caption><thead><tr><th rowspan="2">SNR<sup>-1</sup></th><th colspan="3">c<sub>1</sub></th><th colspan="3">c<sub>2</sub></th></tr><tr><th>Naive</th><th>Robust C<sub>m</sub></th><th>Ours</th><th>Naive</th><th>Robust C<sub>m</sub></th><th>Ours</th></tr></thead><tbody><tr><td>0</td><td><strong>0.0000</strong></td><td>0.00395</td><td>0.00619</td><td><strong>0.0000</strong></td><td>0.165</td><td>0.0734</td></tr><tr><td>0.125</td><td><strong>0.00749</strong></td><td>0.00696</td><td>0.0120</td><td>0.308</td><td><strong>0.247</strong></td><td>0.36</td></tr><tr><td>0.25</td><td><strong>0.00557</strong></td><td>0.00669</td><td>0.0163</td><td>0.641</td><td>0.238</td><td><strong>0.225</strong></td></tr><tr><td>0.5</td><td><strong>0.00656</strong></td><td>0.0228</td><td>0.0152</td><td>0.887</td><td>0.229</td><td><strong>0.226</strong></td></tr><tr><td>1</td><td>0.0182</td><td>0.0346</td><td><strong>0.0151</strong></td><td>0.873</td><td>0.511</td><td><strong>0.227</strong></td></tr><tr><td>2</td><td>0.0376</td><td>0.0353</td><td><strong>0.0147</strong></td><td>0.551</td><td>0.806</td><td><strong>0.229</strong></td></tr><tr><td>4</td><td>0.0628</td><td>0.0345</td><td><strong>0.0146</strong></td><td>0.565</td><td>1.15</td><td><strong>0.228</strong></td></tr><tr><td>8</td><td>0.0932</td><td>0.0358</td><td><strong>0.0146</strong></td><td>1.79</td><td>1.56</td><td><strong>0.227</strong></td></tr></tbody></table>

positions (i.e., sensor space), each sensor location aggregating both magnetometer and gradiometer data, giving a total of 306 time series.

The recording is a regular recording that shows noise, from an internal data set [49]. The study protocol for MEG data was approved by the Coordinating Ethical Committee of Helsinki University Central Hospital (ID 290/13/03/2013), written informed consent was obtained from each subject prior to the experiment, and all research was carried out according to the Declaration of Helsinki.

Of note, interested users could apply this algorithm to public data sets (available, e.g., in OpenNeuro [50]). In that case, pre-processing steps must be run to reproduce the source-level analysis.

## 4.2 Data pre-processing

The data are a 10 min resting-state recording from a healthy subject [49], which was preprocessed using the standard pipeline available in MNE-Python [51]: Maxwell filtering, ICA denoising of heartbeat and eye blink artifacts, source space projection. Details of the pre-processing method are described in [49]. Source localization on the cortical surface was then performed with the preprocessed MEG data using the minimum norm solution [51]. Subsequently, the source reconstructed time series were averaged for each of the 200 parcels of the Yeo 2018 parcellation of the cortical surface [52].

## 4.3 Data analysis

The envelope of the alpha band (8–12 Hz) oscillatory activity was extracted in each parcel of the cortical surface via a continuous wavelet transform, using the complex-valued Morlet wavelet with 5 cycles and with its central frequency set at 10 Hz. The outlier detection method and subsequent multifractal analysis were carried out using the absolute value of the resulting time series, following the methodology introduced in [15]. Note that the proposed approach could also be tested on the broadband spectrum of the preprocessed MEG time series as multifractal analysis was introduced and compared to classical scale-free methods in [18, 53, 54].

The parameters used to obtain the results presented here are as follows: $\theta = 1.8$, $\beta_{PELT} = 12$ due to the shorter duration compared to the simulated data. The model used was the exponential power distribution for the log 1-leaders. In order for the final multifractal analysis to be made possible, fractional integration of the time series was
necessary: $\gamma_{int} = 0.5$ was chosen such that $\eta(2) > 0$ for all time series. For the sake of comparison of the results with what is obtained from DFA where $\gamma_{int} = 1$, the $c_1$ values presented have been increased by 0.5. Given their poor performance on simulations, robust cumulant estimates are not reported on the experimental data.

## 4.4 Results

We first present the application of the proposed method on a single noisy time series from the recording. Using $p = 2$ for the final analysis, the outcome is shown in Figure 4: Panel A introduces the noisy time series and its associated alpha band envelope signal that is used for multifractal analysis. After applying the outlier detection procedure, the wavelet 2-leaders identified as outliers are highlighted in panel B. The spectra that are obtained after correction match the multifractal formalism: The function $\zeta(q)$ and the multifractal spectrum are concave, while it was not the case before correction as presented in Panel C: Although it approximately matches the corrected spectrum for low values of $q$ (associated with high regularity), the rest of the spectrum shows a cusp that does not fit the multifractal formalism. In contrast, the corrected time series has a properly formed parabolic concave spectrum. This difference is also reflected in the estimation of the log-cumulants presented in panel D: $c_2$ was positive before correction, whereas it is negative after correction in accordance once again with the multifractal formalism.

The outlier detection and removal method was then applied to the whole 200 parcel ensemble, with results presented in Fig. 5. Corrected and uncorrected estimates are systematically compared: Panel A presents the co-distribution of estimates before and after correction where the outlier data points move toward the core of the distribution. It is also important to note that certain time series remain entirely unaffected by the detection procedure, highlighting the fact that it does not indiscriminately

<figure><img src="image_7.png" /><figcaption>**Fig. 4** Single MEG parcel time series. **A** (top) preprocessed time series which is the average of the source projected time courses in the right hemisphere extrastriate parcel. (bottom) Envelope of the alpha band obtained via the absolute value of the complex Morlet wavelet transform with $N_{cycles} = 5$ and central freq. 10 Hz. **B** Wavelet 2-leader transform of the envelope signal, between scales ($j_1, j_2$) = (10, 13). Blue rectangles denote 2-leaders removed from the corrected analysis. **C** Scaling functions $\zeta(q)$ (left) and multifractal spectra (right) from both the original (blue) and corrected (orange) fits. **D** Log-cumulant fit of the $c_1(d)$ and $c_2(e)$ estimates from the cumulant scaling functions. The fit on the original signal is in blue, while the fit on the corrected signal is in orange</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 5** Brain-wide cortical source level comparative analysis of MEG data contrasting the proposed methodology with the absence of outlier detection. **A** Paired distribution of $c_1$ and $c_2$ estimates, both original (blue crosses) and corrected (colored dots, color matches colormap in (B, C)); each dot represents a brain parcel. In the margins: histogram of the $c_1$ (x-axis), and the $c_2$ (y-axis) estimates, both original(blue) and corrected (orange). The displayed $c_1$ estimates were increased by 0.5 to match the usual range of values obtained from DFA. **B** Paired distribution of the original and corrected estimates: $c_1$ (left) and $c_2$ (right). The color in B and C denotes the average percentage of coefficients rejected across scales $j_1$ to $j_2$. **C** Paired $\eta(2)$ estimates before and after correction (left); paired difference (corrected - original) for $c_1$ and $c_2$ estimates (right). **D** Percentage of 2-leaders removed at each scale $j$ (left). Shift in $c_1$ (middle) and $c_2$ (right) estimates as a function of the average number of 2-leaders rejected over $[j_1, j_2]$</figcaption></figure>

eliminate coefficients in every instance. Panel B shows that correction generally reduces $c_1$ estimates and greatly reduces positive $c_2$ estimates so that they are almost all lower than zero, conforming to the multifractal formalism. Furthermore, panel C shows that without correction, $\eta(2)$ is often overestimated. This is an important issue, as this quantity serves as the basis for evaluating the need for fractional integration before applying multifractal analysis. If it is overestimated, it will lead in some cases to the under-integration of the time series, which means that the analysis would not comply with the multifractal formalism. Finally, the proportion of removed 2-leaders and its impact on the variation in $c_1$ and $c_2$ are presented in panel D: The proportion of coefficients removed increases with the temporal scale, which is expected since there are fewer 2-leaders at higher temporal scales and we propagate rejections at a given scale to all other temporal scales. Averaging the proportion of rejected 2-leaders across all temporal scales, we observe that the large corrections in $c_1$ and $c_2$ occur even when the proportion of rejected coefficients remains below 30%. This confirms that for parcels where the original $c_2$ estimate was positive and far from zero, a minority of the temporal support drove the majority of the deviation in the $c_1$ and $c_2$ outliers.

Figure 6 shows $c_1$ and $c_2$ estimates localized on the cortical surface. This reveals that the corrected time series are co-localized in specific regions. The greatest amount of correction for both $c_1$ and $c_2$ occurs in the medial occipital cortex, consistently in
<figure><img src="image_9.png" /><figcaption>**Fig. 6** Localized multifractal estimates. Maps of the $c_1$ (A) and $c_2$ (B) estimates on the lateral (top half) and medial (bottom half) views of the cortex, of both hemispheres: first the left hemisphere on top, then the right hemisphere below. The results are presented as originally estimated (left), and after correction (right)</figcaption></figure>

both hemispheres. The occipital lateral cortex receives significant corrections for a few parcels in both hemispheres. The parietal lobe (medial and lateral) shows minimal amounts of correction in both $c_1$ and $c_2$. The frontal and temporal lobes show isolated parcels that received small corrections in $c_1$ and $c_2$.

**5 Discussion**

The classic approach to MFA is very sensitive to outliers, because it is designed to pick up variations in local pointwise regularity, and outliers locally disrupt the pointwise regularity of the underlying signal. A first straightforward approach to correct for the presence of outliers in the data consists in replacing the calculation of cumulants of log-leaders with their robust counterpart. However, this basic intuitive approach fails as robust cumulants do not correctly propagate the correction of outliers across multiple scales. In a previous attempt [31], we designed a Monte Carlo-based approach to identify the presence of outliers; however, its sampling step was too slow, and the final estimates of multifractal parameters were still biased for large noise power. The method developed here rectifies these issues and proposes a modality-agnostic corrected estimate of multifractality in the presence of a small number of noise impulses. On simulation results, we showed that in situations where the noise has high amplitude, the estimation error is reduced 6.3- and 7.9-fold for $c_1$ and $c_2$, respectively, making those results usable in practice. When the noise power is of the same amplitude as the signal, we still observe a 3.8-fold improvement in the estimation of $c_2$. Importantly, beyond the
lowest noise amplitude level, there is no relationship between noise power and estimated $c_1$ and $c_2$, after application of our method.

## 5.1 Limits of the algorithm

However, the proposed approach has a few limitations. Firstly, the method remains univariate, whereas electrophysiological data are multivariate recordings where noise is located across multiple time series at specific time points. Gathering information from multiple sensors should, in principle, help with outlier rejection. Secondly, there are user-defined parameters $\theta$ and $\beta_{PELT}$ in the proposed algorithm, as well as a modeling choice with regard to the distribution of the log $p$-leaders. Based on the simulation results, the similarity of the results between the two choices of distribution gives us confidence in the stability of the method with respect to small changes in the modeling distribution. Furthermore, the parameters $\theta$ and $\beta_{PELT}$ are rather easily tuned: Small changes in $\theta$ induce small updates in the number of rejected coefficients; $\beta_{PELT}$ simply needs to be set low enough such that the outliers are segmented accurately, without being too low so that the segments become too small and the evidence needed to reject them becomes too small. A large range of parameter values give qualitatively very similar results, so fine-tuning the parameters is not critical.

Finally, as the removed coefficients at a given scale propagate to the other scales, since the number of coefficients at coarser scales decays with $2^{-j}$—i.e., going from one temporal scale $j$ to $j + 1$ halves the number of coefficients, even a few deletions at finer scales will make it so that the coarser scales are unusable. However, for small enough values of $p$, when the difference in temporal scales is large enough, the central limit theorem dictates that the variance of the contribution of the coefficients at the lower scale is small enough to be neglected and those coefficients could easily be replaced. In addition, methods such as wavelet quantile normalization [55] can be envisaged to replace the holes left by removing coefficients.

## 5.2 Relevance to neuroscience

Historically, the application of multifractal analysis to electrophysiological recordings (EEG, MEG and intracranial EEG) has proved to be a difficult challenge when showing significant effects between experimental conditions, or even on resting-state recordings [18, 54]. Indeed, as a small amount of noise may corrupt the final estimates of multifractal properties of brain signals, controlling for the presence of noise impulses in brain signals is absolutely necessary to obtain a reliable and rich characterization of multifractal brain activity. Furthermore, trivial robust statistical approaches directly applied to the raw time series will fail due to the scale-invariant nature of the data. It is necessary to propagate the rejection of outliers across the temporal scales; otherwise, the rejection algorithm will not be effective. The sensitivity of multifractal analysis to outliers may also provide a sensitive detector of outliers in the data: Values of $c_2$ significantly above zero strongly suggest the presence of outliers, for instance due to noise impulses, and thus, particular attention should be paid to this phenomenon, otherwise leading to wrong interpretations.
## 5.3 Applicability

As a consequence of the high impact of impulsive noise on the final estimates, the ability of the proposed algorithm to reduce its impact should give greater statistical power to studies of neurological functions and opens the door to using MFA to answer neuroscience questions. To that end, the code will be available at the time of publication,<fnref n="1" /> and care will be taken to ensure that it can be used by neuroscientists without the need for specific expertise in MFA. The simplicity of the algorithm makes it very fast: applying the algorithm on a 4 min 20 s-long recording sampled at 1000 Hz takes less than 10 s. This enables it to scale smoothly to large data set sizes that are becoming more and more common in neuroscience due to the reproducibility crisis. Further speed gains in Python are currently limited by the availability of a fast implementation of the PELT algorithm.

## 5.4 Impact of other noise sources

The possible presence of other sources of noise is a potential concern: The proposed method only aims to detect coefficients which fall outside the distribution of the coefficients obtained from the rest of the time series, in the context of scale-invariant signals. If the other noise sources do not affect the stationarity or scale invariance of the time series, then we expect that the method will perform as well in removing the outlier coefficients, but will not mitigate the other noise sources' potential effects on the estimation of the $(c_m)_m$ coefficients. On the other hand, if other noise sources introduce a significant deviation in the scale-invariant behavior of the observed signal, then the method could potentially be affected; although those kinds of noise sources would affect the multifractal estimates even in the absence of the application of the proposed method. This means that, in noisy scenarios, it is important to evaluate the consistency of the outcome when applying the proposed algorithm.

## 5.5 Extension to non-stationary multifractality

Beyond noise impulses, another source of positive $c_2$ exponents is present in time series where the multifractal properties of the signal are not stationary in time, but instead shift from one regime to another at some specific point in the recording. This means that computing the cumulants $C_m(j)$ throughout the recording corresponds to averaging quantities drawn from different statistical distributions, leading to incorrect multifractal estimates. Estimating the change points in multifractal exponents as proposed in the exposed methodology would address this problem and may be thought of as a more general version of the problem we handled in this paper.

# 6 Conclusion

The present contribution identifies the impact of noise impulses on the estimation of scale-free properties of simulated monofractal and multifractal time series. We proposed a new method, which minimizes the impact of the presence of noise impulses on the final estimates, compared to the classical multifractal estimates, as well as a simple technique based on robust cumulants. The proposed approach has been first validated

<footnote n="1">¹ https://github.com/neurospin/pymultifracs</footnote>
on synthetic monofractal and multifractal signals, where its capacity to remove bias from the results has been illustrated. The method has then been assessed on resting-state MEG data, and we have shown that it properly handles the current issues raised by the naive application of the standard MFA, such as the recovery of positive $c_2$ values, which by construction must be negative to reflect the concave shape of the multifractal spectrum.

The high impact of outliers on the final multifractal estimates means that it is important to address their presence in recordings of brain activity. The method proposed here then represents a breakthrough step toward achieving reliable and fast multifractal estimates of brain activity recordings.

## Appendix A: Hilbert metric

Given a bounded open convex set $\Omega$, the Hilbert metric is a natural extension of the notion of a hyperbolic metric. For $A, B \in \Omega^2$, the Hilbert metric is the logarithm of the cross ratio of the chord $X, A, B, Y$ where X and Y are the points where the line that goes through A and B intersects the boundary $\delta\Omega$:

$$d_{\mathcal{H}}(A, B) = \log \frac{\|Y - A\|_2 \cdot \|X - B\|_2}{\|Y - B\|_2 \cdot \|X - A\|_2}, \qquad (20)$$

$$d_{\mathcal{H}}(A, B) = \log \left[ \left( \frac{\|B - A\|_2 + \|Y - B\|_2}{\|Y - B\|_2} \right) \left( \frac{\|A - B\|_2 + \|X - A\|_2}{\|X - A\|_2} \right) \right], \qquad (21)$$

$$d_{\mathcal{H}}(A, B) = \log \left[ \left( 1 + \frac{\|B - A\|_2}{\|Y - B\|_2} \right) \left( 1 + \frac{\|B - A\|_2}{\|X - A\|_2} \right) \right]. \qquad (22)$$

Using the Hilbert metric on the open convex $(0, 1)^s$ is possible, but the distance to the boundary is, in fact, not the most relevant property in our case. The data points which lie close to 0 would contribute greatly to the distance; however, they are not relevant to detecting outliers, which will instead tend to have samples with values close to 1.

Outliers, because they are located at the right tail of their distribution, tend to “push” low-power fluctuations to the point $0_r$, in contrast with the outlier themselves that are close to $1_r$.

## Appendix B: Detecting added noise on clean MEG data

Although the simulations allow us to objectively evaluate the performance of the proposed algorithm, they represent a simple model for multifractal time series. Indeed, MEG data may come from a process that generates multifractality in a way that prevents the algorithm from properly functioning. To alleviate this concern, we simulated a noisy scenario by adding three white noise bursts to a MEG recording, different from the one presented in the body of the article, which has been identified to be free of outliers; the noise bursts were produced following the procedure established for the simulations.
<figure><img src="image_10.png" /><figcaption>**Fig. 7** Proposed method detects added white noise to MEG alpha band envelope (N = 200 brain parcels). The multifractal estimates $c_1$ (left), and $c_2$ (right) are compared before (x-axis) and after (y-axis) addition of white noise bursts. The estimates are represented with (orange) and without (blue) application of the outlier detection method.</figcaption></figure>

Figure 7 shows the multifractal estimates before and after adding white noise bursts. As expected based on the simulations, the impact of the noise bursts is larger on the $c_2$ estimates than it is on the $c_1$ estimates. The outlier removal procedure brings both the $c_1$ and $c_2$ estimates closer to their original values. The corrected estimates are not exactly the same as the uncorrected ones, which is due to the fact that the corrected estimates are obtained from a subset of the data used to estimate the original $c_1$ and $c_2$ values.

## Appendix C: Simulating multiple $H_{noise}$ values

The results presented in Fig. 3 show how varying the power of the noise impulses affects the multifractal estimates on simulated data. Noise may be controlled by more than a single parameter, and those parameters may cause different effects on the multifractal estimates.

To ensure that the analysis method is robust to changes in the nature of the noise, we extended the simulations to include different values of $H$ for the noise process: from an anti-correlated noise ($H = 0.3$) to a very auto-correlated noise ($H = 0.9$). Figure 8 presents the

<figure><img src="image_11.png" /><figcaption>**Fig. 8** Simulation results with noise that has a varying $H_{noise}$ parameter. Results are presented using the same parameters as in Fig. 3, except for $H_{noise}$, which varies from 0.3 to 0.9</figcaption></figure>
results: Auto-correlated noise ($H_{noise} > 0.5$) has a larger impact on the $c_1$ and $c_2$ estimates, but is also easier for the proposed algorithm to handle. Anti-correlated noise ($H_{noise} < 0.5$) has a lower impact at low noise levels, but is consequently more difficult to identify.

## Funding

This work is supported by the "ADI 2020" project funded by the IDEX Paris-Saclay, ANR-11-IDEX-0003-02 to M.D.; by the Instrumentarium Foundation to M.D.; by ANR 19-CE48-0002-04 DARLING to M.D. and P.C.; by the Academy of Finland (SA 1266745, 1296304 to J.M.P.) by the Sigrid Jusélius Foundation to J.M.P. The calculations presented above were performed using computer resources within the Aalto University School of Science "Science-IT" project.

Received: 22 October 2024 Accepted: 28 July 2025

Published online: 25 August 2025

## References

1. C.K. Peng, S.V. Buldyrev, S. Havtin, M. Simons, H.E. Stanley, A.L. Goldbergerz, Mosaic organization of DNA nucleotides. Phys. Rev. E **49**, 1685 (1994)

2. A.A.P. Abry, Semi- and biorthogonal MRA-type wavelet designs and their fast algorithms. SPIE Proc. **2569**, 452–463 (1995). https://doi.org/10.1117/12.217601

3. P. Abry, Wavelet based estimators for self-similar α-stable processes. WCC 2000 - ICSP 2000. 2000 5th International Conference on Signal Processing Proceedings. 16th World Computer Congress 2000 **1**, 369–375 (2000). https://doi.org/10.1109/icsp.2000.894512

4. S.S. Emanuel Guariglia, Fractional-wavelet analysis of positive definite distributions and wavelets on $d'(c)$. Springer Proceedings in Mathematics & Statistics pp. 337–353 (2016). https://doi.org/10.1007/978-3-319-42105-6_16

5. R.C. Guido, Effectively interpreting discrete wavelet transformed signals [lecture notes]. IEEE Signal Process. Mag. **34**(3), 89–100 (2017). https://doi.org/10.1109/bsp.2017.2672759

6. L. Yang, Hyperspectral image classification using wavelet transform-based smooth ordering. Int. J. Wavelets Multiresolut. Inf. Process. **17**(06), 1950050 (2019). https://doi.org/10.1142/s0219691319500504

7. G. Samorodnitsky, M.S. Taqqu, *Stable Non-Gaussian Random Processes* (Routledge, 2017). https://doi.org/10.1201/9780203738818

8. E. Guariglia, Harmonic sierpinski gasket and applications. Entropy **20**(9), 714 (2018). https://doi.org/10.3390/e20090714

9. E. Guariglia, Primality, fractality, and image analysis. Entropy **21**(3), 304 (2019). https://doi.org/10.3390/e21030304

10. R. Achour, A multifractal formalism for new general fractal measures. Chaos, Solitons Fractals **181**, 114655 (2024). https://doi.org/10.1016/j.chaos.2024.114655

11. P. Ciuciu, Interplay between functional connectivity and scale-free dynamics in intrinsic FMRI networks. Neuroimage **95**, 248–263 (2014). https://doi.org/10.1016/j.neuroimage.2014.03.047

12. J.M. Beggs, D. Plenz, Neuronal avalanches in neocortical circuits. J. Neurosci. **23**, 11167–11177 (2003). https://doi.org/10.1523/jneurosci.23-35-11167.2003

13. M. Breakspear, J.A. Roberts, J.R. Terry, S. Rodrigues, N. Mahant, P.A. Robinson, A unifying explanation of primary generalized seizures through nonlinear brain modeling and bifurcation analysis. Cereb. Cortex **16**, 1296–1313 (2006). https://doi.org/10.1093/cercor/bhj072

14. W.L. Shew, D. Plenz, The functional benefits of criticality in the cortex. Neuroscientist **19**, 88–100 (2013). https://doi.org/10.1177/1073858412445487

15. K. Linkenkaer-Hansen, V.V. Nikouline, J.M. Palva, R.J. Ilmoniemi. Long-range temporal correlations and scaling behavior in human brain oscillations. J. Neurosci. **21**(4), 1370–1377 (2001)

16. P. Ciuciu, G. Varoquaux, P. Abry, S. Sadaghiani, A. Kleinschmidt, Scale-free and multifractal time dynamics of FMRI signals during rest and task. Front. Physiol. **3**, 186 (2012). https://doi.org/10.3389/fphys.2012.00186

17. N. Zilber, P. Ciuciu, P. Abry, V.V. Wassenhove, Modulation of scale-free properties of brain activity in meg. Proceedings-International Symposium on Biomedical Imaging pp. 1531–1534 (2012). https://doi.org/10.1109/ISBI.2012.6235864

18. D.L. Rocca, N. Zilber, P. Abry, V. van Wassenhove, P. Ciuciu, Self-similarity and multifractality in human brain activity: a wavelet-based analysis of scale-free brain dynamics. J. Neurosci. Methods **309**, 175–187 (2018). https://doi.org/10.1016/j.neumeth.2018.09.010

19. Z. Fayyaz, M. Bahadorian, J. Doostmohammadi, V. Davoodnia, S. Khodadadian, R. Lashgari, Multifractal detrended fluctuation analysis of continuous neural time series in primate visual cortex. J. Neurosci. Methods **312**, 84–92 (2019). https://doi.org/10.1016/j.neumeth.2018.10.039

20. G. Alamian, T. Lajnef, A. Pascarella, J.M. Lina, L. Knight, J. Walters, K.D. Singh, K. Jerbi, Altered brain criticality in schizophrenia: new insights from magnetoencephalography. Front. Neural Circuits **16**, 630621 (2022). https://doi.org/10.3389/fncir.2022.630621

21. R. Lopes, Multifractal analysis in neuroimaging. Adv. Neurobiol. **36**, 79–93 (2024). https://doi.org/10.1007/978-3-031-47606-8_4

22. J. Ludescher, M.I. Bogachev, J.W. Kantelhardt, A.Y. Schumann, A. Bunde, On spurious and corrupted multifractality: the effects of additive noise, short-term memory and periodic trends. Phys. A **390**, 2480–2490 (2011). https://doi.org/10.1016/j.physa.2011.03.008
23. E.A.F. Ihlen, Multifractal analyses of response time series: a comparative study. *Behav. Res. Methods* **45**, 928–945 (2013). https://doi.org/10.3758/s13428-013-0317-2

24. Z.R. Struzik, A.P. Siebes, Wavelet transform based multifractal formalism in outlier detection and localisation for financial time series. *Phys. A* **309**, 388–402 (2002). https://doi.org/10.1016/S0378-4371(02)00552-6

25. K. Domino, Multivariate cumulants in outlier detection for financial data analysis. *Phys. A* **558**, 124995 (2020). https://doi.org/10.1016/j.physa.2020.124995

26. M. Jas, D.A. Engemann, Y. Bekhti, F. Raimondo, A. Gramfort, Autoreject: automated artifact rejection for MEG and EEG data. *Neuroimage* **159**, 417–429 (2017). https://doi.org/10.1016/j.neuroimage.2017.06.030

27. J.A. Cramer, S.S. Shah, T.M. Battaglia, S.N. Banerji, L.A. Obando, K.S. Booksh, Outlier detection in chemical data by fractal analysis. *J. Chemom.* **18**, 317–326 (2004). https://doi.org/10.1002/cem.875

28. A. Echelard, J.L. Vehel, J. Lévy-Véhel. Wavelet denoising based on local regularity information (2008)

29. S. Wang, D. Mou, X. Qi, Z. Wang, Analysis of multifractal characteristics and detrended cross-correlation of conventional logging data regarding igneous rocks. *Fract. Fract.* **9**(3), 163 (2025). https://doi.org/10.3390/fractalfract9030163

30. C. Bucur, B.G. Tudorică, A. Băra, S.V. Oprea, Multifractal analysis of bitcoin price dynamics. *J. Bus. Econ. Manag.* **26**(1), 21–48 (2025). https://doi.org/10.3846/jbm.2025.23025

31. M. Dumeur, P. Ciuciu, in 29° colloque sur le Traitement du Signal et des Image (GRETSI, 2023)

32. S. Jaffard, B. Lasermes, P. Abry, in *Wavelet Analysis and Applications* (Springer, 2007), pp. 201–246

33. H. Wendt, P. Abry, S. Jaffard, Bootstrap for empirical multifractal analysis. *IEEE Signal Process. Mag.* **24**(4), 38–48 (2007)

34. L. Leon, H. Wendt, J.Y. Tourneret, P. Abry, in *2023 31st European Signal Processing Conference (EUSIPCO)* (IEEE, 2023), pp. 1808–1812

35. B. Hbibi, C. Khiari, K. Wirsing, L. Mili, K. Baccar, A. Mami, Identifying and removing interference and artifacts in multifractal signals with application to EEG signals. *IEEE Access* **11**, 119090–119105 (2023). https://doi.org/10.1109/ACCESS.2023.3325786

36. S. Jaffard, C. Melot, R. Leonarduzzi, H. Wendt, P. Abry, S.G. Roux, M.E. Torres, P-exponent and p-leaders, part i: Negative pointwise regularity. *Phys. A* **448**, 300–318 (2016). https://doi.org/10.1016/j.physa.2015.12.061

37. R. Leonarduzzi, H. Wendt, P. Abry, S. Jaffard, C. Melot, S. Roux, M. Torres, p-exponent and p-leaders, part ii: multifractal analysis relations to detrended fluctuation analysis. *Phys. A* **448**, 319–339 (2016). https://doi.org/10.1016/j.physa.2015.12.035

38. B. Selmi, A review on multifractal analysis of Hewitt-Stromberg measures. *J. Geometric Anal.* **32**(1), 1 (2022). https://doi.org/10.1007/s12220-021-00753-7

39. E. Bacry, J. Delour, J.F. Muzy, Multifractal random walk. *Phys. Rev. E* Stat. Phys. Plasmas Fluids Relat. Interdiscip. Top. **64**, 4 (2001). https://doi.org/10.1103/PhysRevE.64.026103

40. J.F.N. Michael Victor Berry, Z. V. Lewis, On the Weierstrass-Mandelbrot fractal function. *Proc.R. Soc. London A Math. Phys. Sci.* **370**(1743), 459–484 (1980). https://doi.org/10.1098/rspa.1980.0044

41. E.D. Nardo, G. Guarino, D. Senato, A unifying framework for *k*-statistics, polykays and their multivariate generalizations. *Bernoulli* **14**, 440 (2008). https://doi.org/10.3150/07-BEJ6163

42. P.J. Rousseeuw, C. Croux, Alternatives to the median absolute deviation. *J. Am. Stat. Assoc.* **88**, 1273–1283 (1993). https://doi.org/10.1080/01621459.1993.10476408

43. A.E. Beaton, J.W. Tukey, The fitting of power series, meaning polynomials, illustrated on band-spectroscopic data. *Technometrics* **16**, 147 (1974). https://doi.org/10.2307/1267936

44. M. Welling, in *Proceedings of the Tenth International Workshop on Artificial Intelligence and Statistics* (PMLR, 2005)

45. J.E. Angus, The probability integral transform and related results. *SIAM Rev.* **36**, 652–654 (1994). https://doi.org/10.1137/1036146

46. H. Wendt, N. Dobigeon, J.Y. Tourneret, P. Abry, in *IEEE International Conference on Acoustics, Speech and Signal Processing - Proceedings* (IEEE, 2013), pp. 6556–6560. https://doi.org/10.1109/ICASSP.2013.6638929

47. C. Truong, L. Oudre, N. Vayatis, Selective review of offline change point detection methods. *Signal Process.* **167**, 107299 (2020). https://doi.org/10.1016/j sigpro.2019.107299

48. R. Killick, P. Fearnhead, I.A. Eckley, Optimal detection of changepoints with a linear computational cost. *J. Am. Stat. Assoc.* **107**, 1590–1598 (2012). https://doi.org/10.1080/01621459.2012.737745

49. M. Fuscà, F. Siebenhühner, S.H. Wang, V. Myrov, G. Arnulfo, L. Nobili, J.M. Palva, S. Palva, Brain criticality predicts individual levels of inter-areal synchronization in human electrophysiological data. *Nat. Commun.* **14**, 4736 (2023). https://doi.org/10.1038/s41467-023-40056-9

50. C.J. Markiewicz, K.J. Gorgolewski, F. Feingold, R. Blair, Y.O. Halchenko, E. Miller, N. Hardcastle, J. Wexler, O. Esteban, M. Goncavles, A. Jwa, R. Poldrack, The openneuro resource for sharing of neuroscience data. *Elife* **10**, e71774 (2021). https://doi.org/10.7554/eLife.71774

51. A. Gramfort, M. Luessi, E. Larson, D.A. Engemann, D. Strohmeier, C. Brodbeck, R. Goj, M. Jas, T. Brooks, L. Parkkonen et al., MEG and EEG data analysis with MNE-Python. *Front. Neurosci.* **7**, 267 (2013)

52. A. Schaefer, R. Kong, E.M. Gordon, T.O. Laumann, X.N. Zuo, A.J. Holmes, S.B. Eickhoff, B.T. Yeo, Local-global parcellation of the human cerebral cortex from intrinsic functional connectivity MRI. *Cereb. Cortex* **28**(9), 3095–3114 (2018)

53. N. Zilber, P. Ciuciu, P. Abry, V. Van Wassenhove, in *2012 9th IEEE International Symposium on Biomedical Imaging (ISBI)* (IEEE, 2012), pp. 1531–1534

54. N. Zilber, P. Ciuciu, P. Abry, V. Van Wassenhove, in *2013 IEEE 10th International Symposium on Biomedical Imaging* (IEEE, 2013), pp. 998–1001

55. M. Dora, D. Holcman, Adaptive single-channel EEG artifact removal with applications to clinical monitoring. *IEEE Trans. Neural Syst. Rehabil. Eng.* **30**, 286–295 (2022). https://doi.org/10.1109/TNSRE.2022.3147072

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.