Brain Topogr (2016) 29:13–26

DOI 10.1007/s10548-015-0448-0

ORIGINAL PAPER

Separating Fractal and Oscillatory Components in the Power

Spectrum of Neurophysiological Signal
Haiguang Wen2 • Zhongming Liu1,2

Received: 8 October 2014 / Accepted: 19 August 2015 / Published online: 29 August 2015
Ó Springer Science+Business Media New York 2015

Abstract Neurophysiological field-potential signals con- demonstrate the efficacy and potential applications of this
sist of both arrhythmic and rhythmic patterns indicative of method in analyzing electrophysiological signatures of
the fractal and oscillatory dynamics arising from likely large-scale neural circuit activity. We expect that the pro-
distinct mechanisms. Here, we present a new method, posed method or its future variations would potentially
namely the irregular-resampling auto-spectral analysis allow for more specific characterization of the differential
(IRASA), to separate fractal and oscillatory components in contributions of oscillatory and fractal dynamics to dis-
the power spectrum of neurophysiological signal according tributed neural processes underlying various brain
to their distinct temporal and spectral characteristics. In functions.
this method, we irregularly resampled the neural signal by
a set of non-integer factors, and statistically summarized
the auto-power spectra of the resampled signals to separate Keywords Fractal  Oscillation  Power-law  Self-
the fractal component from the oscillatory component in affinity  Scale-free
the frequency domain. We tested this method on simulated
data and demonstrated that IRASA could robustly separate
the fractal component from the oscillatory component. In Introduction
addition, applications of IRASA to macaque electrocor-
ticography and human magnetoencephalography data RICH temporal dynamics of neural activity is commonly
revealed a greater power-law exponent of fractal dynamics captured by electrophysiological recordings, e.g. local field
during sleep compared to wakefulness. The temporal potential (LFP), electrocorticogram (ECoG), electroen-
fluctuation in the broadband power of the fractal compo- cephalography (EEG) and magnetoencephalography (MEG)
nent revealed characteristic dynamics within and across the (Nunez and Srinivasan 2006). Such neural signals exhibit a
eyes-closed, eyes-open and sleep states. These results variable mixture of rhythmic and arrhythmic patterns
(Buzsaki et al. 2012). The former arises from oscillatory
network activity with a characteristic time scale (Buzsáki
Electronic supplementary material The online version of this
and Draguhn 2004), whereas the latter is not confined to any
article (doi:10.1007/s10548-015-0448-0) contains supplementary specific scale and reflects the so-called fractal (or scale-free)
material, which is available to authorized users. dynamics (He et al. 2010). When these signals are trans-
formed to corresponding power spectra, the oscillatory
& Zhongming Liu
zmliu@[Link] component results in multiple discrete peaks at specific fre-
quencies, and the fractal component manifests itself as a
1
Weldon School of Biomedical Engineering, College of descending straight line on the log–log plot suggesting a
Engineering, Purdue University, West Lafayette, IN 47907, power-law relationship (Miller et al. 2009).
USA
Oscillatory and fractal activities are likely generated
2
School of Electrical and Computer Engineering, College of through different mechanisms (Buzsáki and Draguhn 2004;
Engineering, Purdue University, West Lafayette, IN 47907,
USA He et al. 2010) and report on distinct features of large-

123
14 Brain Topogr (2016) 29:13–26

scale neural networks (Bullmore and Sporns 2009; Engel series. Such interactions are difficult to remove based on
et al. 2013; Siegel et al. 2012). In particular, neuronal the cross-spectral analysis, bringing into question the
oscillations have been of long-standing interest to neuro- practical efficacy of CGSA for separating the spectral
physiologists, and in practice are extracted by applying components attributed to the fractal and oscillatory
band-pass filtering to neural signals. However, the result- dynamics.
ing band-limited signals, although appearing oscillatory, Here we propose a new method, referred to as the
can result from filtering even in the absence of any bio- Irregular-Resampling Auto-Spectral Analysis (IRASA), to
logical oscillator, and thus may sometimes reflect artifi- separate fractal and oscillatory components in the power
cially extracted fractions of the brain’s arrhythmic activity spectrum of neurophysiological signals. In IRASA, we
(He 2014). resample a neural signal by multiple non-integer pairwise
Analysis of fractal dynamics has also attracted growing factors (positive numbers and their reciprocals), and then
interest with the current emphasis directed upon its fre- compute the geometric mean of the auto-power spectra of
quency scaling property (Ciuciu et al. 2012; El Boustani every pair of the resampled signals. In the resulting spec-
et al. 2009; Fransson et al. 2013; Freeman 2007; He et al. trum, the power associated with the oscillatory component
2010; Manning et al. 2009). To characterize this property, is redistributed away from its original (fundamental and
power spectral density (PSD) of neural signals is often harmonic) frequencies by a frequency offset that varies
fitted by a power-law function, or equivalently a linear with the resampling factor, whereas the power solely
function in double-logarithmic coordinates, where the attributed to the fractal component remains the same
slope provides an estimate of the underlying power-law (or power-law statistical distribution independent of the
frequency-scaling) exponent (Manning et al. 2009; Miller resampling factor. It follows that taking the median of the
et al. 2009). The precision of such estimation, however, is mean auto-power spectra of the variously resampled sig-
compromised in the presence of prominent oscillations that nals can extract the power spectrum of the fractal compo-
deviate the PSD from a power-law distribution (He et al. nent, and the difference between the original power
2010; Miller et al. 2009). spectrum and the extracted fractal spectrum offers an
Therefore, it is desirable to separate fractal and approximate estimate of the power spectrum of the oscil-
oscillatory components in order to analyze their indi- latory component.
vidual characteristics without being concerned about In what follows, we will show that IRASA is robust
their mutual interference. Toward this end, a method against the presence of the complex interaction between
known as the coarse graining spectral analysis (CGSA) fractal and oscillatory components, being useful to extract
has been developed to extract and analyze the fractal the broadband activity that conforms to the scale-free
dynamics (Yamamoto and Hughson 1991, 1993) of dynamics. To demonstrate this point, we have evaluated
human ECoG signals (He et al. 2010). Central to CGSA IRASA with both simulation and experimental (MEG and
is the so-called self-affinity property of any fractal time ECoG) data. We have noted that the fractal dynamics in
series: the statistical distribution remains unchanged different behavioral states gives rise to distinct broadband
when resampled at different time scales (Mandelbrot and features, suggesting that IRASA can be used to potentially
Van Ness 1968). It follows that after resampling a time- reveal and differentiate the varying fractal signatures
series signal (e.g. by a factor of 2 or ‘), the cross power across brain states. Other implications and future applica-
spectrum of the original and resampled signals is tions are also discussed. In the online supplementary
expected to follow a power-law distribution if the signal materials, we will further compare IRASA with CGSA in
is self-affine or fractal, but close to zero if the signal is theory and with simulation data.
simply periodic or oscillatory (Yamamoto and Hughson
1991, 1993). For this reason, CGSA has been presumed
to be able to extract the fractal component from the Methods and Materials
power spectrum of a time series that mixes both fractal
and oscillatory processes (González et al. 1999; He et al. The problem addressed here is to separate the power
2010; Pereda et al. 1998; Yamamoto and Hughson 1991, spectra of the unknown fractal activity, f(t), and the
1993). Although it sounds intuitively reasonable, this unknown oscillatory activity, x(t), given the measured
presumption is questionable in theory. When a signal is neural signal, y(t), which is assumed to consist of both f(t)
composed of both fractal and oscillatory time series, the and x(t) through a simple additive model without noise.
cross power spectrum of this signal and its resampled yðtÞ ¼ f ðtÞ þ xðtÞ ð1Þ
version includes non-negligible contributions from the
interactions of the original (or resampled) fractal time To solve this problem, we have developed the IRASA
series and the resampled (or original) oscillatory time method that utilizes the self-affine property of fractal time

123
Brain Topogr (2016) 29:13–26 15

oscillation x(t)
PSD
frequency

PSD
log PSD
total y(t)
frequency
log PSD
time t log frequency IRASA
fractal f(t)

log PSD
log frequency
time t
log frequency

time t

h=1.2
h=1.4
h=1.6
h=1.8
h=1.2, 1.4, 1.6, 1.8
resample median

Fig. 1 Overall schematic illustration of IRASA. The input for IRASA is a mixed time series composed of both fractal and oscillatory signals and
the output is the separated power spectra of these two components

series and the frequency-specific nature of oscillatory time Fh ðxÞ ¼ hH FðxÞ ð3Þ
series (Fig. 1). For IRASA, the spectral component due to
the fractal activity was first extracted from the spectrum of It means that the amplitude spectrum of the resampled
the mixed neural signal, and then taking the difference fractal time series is the same as that of the original time
between the mixed spectrum and the extracted fractal series scaled by a factor of hH. Note that Eq. 3 holds true
spectrum served to approximate the spectral component only if F(x) follows a power-law distribution shown as a
due to the oscillatory activity. In the following, we will first line in log–log coordinates.
introduce the temporal and spectral characteristics of
fractal and oscillatory signals; then we will elaborate the Narrow-Band Nature of Oscillations
theory and algorithm for IRASA; lastly we will describe
the simulation and experimental data used for performance Let x(t) stand for an oscillatory time series and xh(t) be the
evaluation. The theoretical basis and limitations of CGSA time series obtained by resampling x(t) by a factor of h. Let
are elaborated in the online supplementary materials in XðxÞejbðxÞ and Xh ðxÞejbh ðxÞ be the Fourier representations
addition to (Yamamoto and Hughson 1991, 1993). of x(t) and xh(t), respectively. By definition x(t) is a peri-
odic and narrow-banded signal. The power spectra of both
Self-Affinity of Fractal Time Series x(t) and its resampled version xh(t) are non-zero only at
specific frequencies and close to zero elsewhere. The
A fractal time series, f(t), satisfies a self-affine relationship ‘‘sparse’’ spectral distribution is characteristic of oscilla-
expressed as (2). tions in contrast to the fractal dynamics, which has a
‘‘continuous’’ broadband distribution.
fh ðtÞ , hH f ðtÞ ð2Þ
It indicates that when the fractal time series is resampled Irregular-Resampling Auto-Spectral Analysis
by a factor of h (h [ 0), the resampled time series, denoted (IRASA)
as fh(t), has the same statistical distribution as the original
one scaled by a factor of hH, where H is called the Hurst In IRASA, we extract the fractal power spectrum through
exponent (Mandelbrot and Van Ness 1968). Applying the computing the auto-power spectra of the resampled signals
Fourier transform to f(t) and fh(t) yields FðxÞejaðxÞ and that result from downsampling and upsampling the mea-
Fh ðxÞejah ðxÞ respectively, where F(x) and a(x) are the sured signal by a set of non-integer factors. As such, the
amplitude and phase at specific angular frequencies x for resulting time series signals are ‘‘irregularly’’ resampled
f(t), and Fh(x) and ah(x) are the counterparts for fh(t). This from the original signal.
self-affine relationship is equivalent to the frequency The auto-power spectra of yh(t) and y1/h(t) are as (4) and
scaling property expressed as (3). (5), respectively.

123
16 Brain Topogr (2016) 29:13–26

h i
estimates, viewed at every individual frequency, are
Syh yh ðxÞ ¼ Fh ðxÞejah ðxÞ þ Xh ðxÞejbh ðxÞ
h i expected to mostly center about the true power of the
Fh ðxÞejah ðxÞ þ Xh ðxÞejbh ðxÞ ð4Þ underlying fractal component with typically one deviation
due to the interference from the underlying oscillatory
¼ h2H F 2 ðxÞjj1 þ Wh ðxÞejhh ðxÞ jj2 component. This deviation is an ‘‘outlier’’ within a set of
h i otherwise unbiased fractal power estimates at each fre-
Sy1=h y1=h ðxÞ ¼ F1=h ðxÞeja1=h ðxÞ þ X1=h ðxÞejb1=h ðxÞ
quency (as illustrated in Fig. S2). We can eliminate the bias
h i
F1=h ðxÞeja1=h ðxÞ þ X1=h ðxÞejb1=h ðxÞ ð5Þ from the outlier by taking the median of these fractal power
estimates at every frequency.
¼ h2H F 2 ðxÞjj1 þ W1=h ðxÞejh1=h ðxÞ jj2
F 2 ðxÞ ¼ medianh fSh ðxÞg; for each x ð8Þ
where Wh ðxÞ ¼ Xh ðxÞ=Fh ðxÞ; hh ðxÞ ¼ ah ðxÞ  bh ðxÞ;
Since the median is a very robust sampling statistic
W1=h ðxÞ ¼ X1=h ðxÞ=F1=h ðxÞ; h1=h ðxÞ
against the presence of outliers with a breakdown point of
¼ a1=h ðxÞ  b1=h ðxÞ. Note that W and h indicate the 50 % (Bassett 1991), Eq. 8 is expected to yield an unbiased
relationship between the oscillatory and fractal components estimate of the power of the fractal component at any
in terms of their ratio in magnitude and their difference in specific frequency as long as the number of outliers is less
phase, respectively. than half at this frequency. See Figs. 2 and S3 for more
By computing the geometric mean of these two auto- detailed illustrations and explanations.
power spectra, we can obtain an initial estimate of the To further estimate the PSD of the oscillatory compo-
fractal power spectrum, which is denoted as Sh ðxÞ and nent, we derive Eq. (9) from Eq. (1).
expressed as (6).
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Y 2 ðxÞ ¼ F 2 ðxÞ þ X 2 ðxÞ þ 2F ðxÞX ðxÞ cosðaðxÞ  bðxÞÞ
Sh ðxÞ ¼ Syh yh ðxÞSy1=h y1=h ðxÞ ð9Þ
¼ F 2 ðxÞjj1 þ Wh ðxÞejhh ðxÞ jj jj1 þ W1=h ðxÞejh1=h ðxÞ jj It shows that the PSD of the oscillatory component is
ð6Þ related not only to the difference between the PSD of the
mixed signal and the PSD of the fractal component, but
If yðtÞ only contains the fractal component, Wh ðxÞ and also to another confounding term that depends on the
W1=h ðxÞ are equal to zero for any x, and Sh ðxÞ is then relative phase difference between the fractal and oscilla-
independent of h and always equals to the fractal power tory components. Assuming no phase–phase coupling
spectrum. If yðtÞ only contains the oscillatory component between the fractal and oscillatory components, the
with a single harmonic frequency x0 , Sh ðxÞ is equal to expectation of this confounding term is close to zero. To
zero for any frequency. approximate the expectation, we can select multiple time
If both fractal and oscillatory components are present in segments within the given period during which the signal
yðtÞ, Sh ðxÞ is as (7). is assumed to be stationary. After averaging the PSD of
Sh ðxÞ¼ the segmented signal across segments, the confounding
8 term reduces to zero. Therefore, the difference of the
>
> F 2 ðxÞ; x6¼ x0 h and x6¼ x0 =h averaged PSD of the raw signal and the estimated PSD of
<
2 jhh ðxÞ
F ðxÞjj1þWh ðxÞe jj; x¼x0 h the fractal component can serve as an unbiased estimate
>
> jh ð xÞ
: F 2 ðxÞjj1þW1 ðxÞe h jj; x¼x =h
1 for the expectation of the PSD of the oscillatory
h
0
component.
ð7Þ
It suggests that Sh ðxÞ is an unbiased estimate of the Algorithm
fractal power spectrum at the oscillatory frequency as well
as other frequencies except x0 h and x0 =h. Importantly, (a) From a given period of the time series signal, we
since the estimation errors are expected to occur only at chose ten time segments. Every time segment was
frequencies that are dependent on the resampling factor h, 90 % of the total length. These segments were
setting different non-integer values (e.g. between 1.1 and evenly distributed within the total period and were
1.9) to h relocates the oscillation-associated power and partially overlapping.
yields a set of estimates for the fractal power spectrum (b) For each segment defined in (a), the auto-power
fSh ðxÞg. In every one of these estimates, the oscillatory spectrum was estimated from the original signal by
power is relocated to different frequencies (as illustrated in using fast Fourier transform (FFT) tapered with a
Fig. 2). In other words, these h-dependent spectral Hanning window. The number of frequency points

123
Brain Topogr (2016) 29:13–26 17

Downsampled Upsampled Geometric mean

4 4 4
10 10 10
3 3 3
10 10 10
2 2 2
10 10 10
h =1.2 101 101 101
100 100 100
−1 −1 −1
10 10 10
0 1 2 0 1 2 0 1 2
10 10 10 10 10 10 10 10 10 4
4 10
4 4 10 4 3
10 10

Log PSD
10 3 10
3 3 10 3
Log PSD

10 10 10 2
2 2 10
2
2 median 10
10 h =1.4 10 10 101
101
101 101 101 100
100
100 100 −1 100 −1
10
−1 −1 10 −1 0 1 2
10 10 0 1 2 0 1 2 10 0 1 2 10 10 10
0 1 2 10 10 10 10 10 10 10 10 10
10 10 10
4 4 Log Frequency
Log Frequency 10 4 10
3 10 3
10 3 10
2 10 2
10 2 10
10
101 101
h =1.6 101
100 100
−1 100 −1
10 −1 10
10 0 1 2 0 1 2
0 1 2 10 10 10 10 10 10
10 10 10
4 4 4
10 10 10
3 3 3
10 10 10
2 2 2
10 10 10
101 101 101
h =1.8
100 100 100
−1 −1
10−1 10 10
0 1 2 0 1 2
0 1 2 10 10 10 10 10 10
10 10 10

Fig. 2 Procedures to estimate the power spectrum of the fractal component by statistically summarizing the geometric means of the auto-power
spectra of the signals up-sampled and down-sampled by a set of non-integer factors

(NFFT) was chosen to be two times the smallest The above algorithm was implemented in Matlab and will soon
power of two that was greater than the number of be made publicly available with example datasets on our
time points in each segment. Given h \ 2 as used in website ([Link]
(c), NFFT is guaranteed to be greater than the The estimated power spectrum of the fractal component
number of time points in the original signal or its was transformed to the log–log scale only for the display
down-sampled or up-sampled version, ensuring and quantification of the power-law distribution. In the log
equal sampling in the angular frequency without scale, there are increasingly more spectral sampling points
losing any overall spectral energy. for higher frequencies. Therefore, we resampled the fre-
(c) For each segment defined in (a), we resampled the quencies to be evenly spaced in the log scale, in order to
signal by h and 1/h, where h ranged from 1.1 to 1.9 prevent the higher frequency range from dominating the
with the increment of 0.05. In order to irregularly estimation of the power-law properties. Then we used the
upsample the signal by h, we interpolated the signal least squares estimation to obtain a linear function that
by using a cubic spline method. To irregularly down best fitted the estimated fractal power spectrum. The slope
sample the signal, we first applied an anti-aliasing of the linear function was the estimated power-law expo-
low-pass filter and then used the cubic spline nent, b, and the intercept of the linear function was
interpolation. Then we estimated the auto-power referred to as the log-power intercept. The mean power in
spectra of the resampled signals using the same FFT- the log scale was the estimated broadband power (Miller
based method and the same number of frequency et al. 2009).
points as in (b). Note that a sliding window can be used to progressively
(d) For each segment defined in (a), the geometric mean select shorter periods out of a relatively longer period, as
of the auto-power spectra of the up-sampled and commonly used in short-time Fourier transform. Within
down-sampled signals was calculated for each each sliding window, the above algorithm and quantifica-
h value. The median of the results with all h-values tion can be used to separate the power spectra of the fractal
was obtained to estimate the power spectrum of the and oscillatory components of the signal, resulting in time
fractal component. frequency representations of the fractal and oscillatory
(e) The estimated PSD of the fractal component and the components.
calculated PSD of the signal were both averaged
across all time segments defined in (a). Computer Simulation
(f) The average power spectrum of the fractal compo-
nent was subtracted from the average power spec- For the initial tests of the proposed IRASA method, we
trum of the original signal to yield the estimated simulated fractal and oscillatory time series without con-
power spectrum of the oscillatory component. sidering noise. The fractal time series was generated by

123
18 Brain Topogr (2016) 29:13–26

inverse FFT applied to a power-law Fourier spectrum as For both ECoG and MEG data, we applied IRASA to a
described in (Yamamoto and Hughson 1993). Briefly, the sliding time-window to obtain a time-varying fractal and
power distribution followed a 1/f trend with a preset oscillatory spectra (i.e. time–frequency representation or
power-law exponent b. The phase distribution was random spectrogram). The sliding window was set to be 3 s with
and followed a uniform distribution between 0 and hmax , 1 s increments. We fitted a power-law function to the
where hmax  2p. estimated fractal spectrum (1–30 Hz) at every sliding time
The oscillatory component was simulated as a single window, and obtained the time courses of the power-law
sinusoidal signal or the sum of a varying number of sinu- exponent, the intercept of the power-law function in log–
soidal signals, and then added to the simulated fractal log coordinates, and the broadband power. We further
component. The amplitudes of the additive sinusoidal averaged the time-resolved spectral components within
signals were systematically varied from 10 to 400 % of distinct behavioral states and compared the results across
those of the fractal component at the same frequencies. All states. The purpose was to characterize and differentiate
of the simulated time series signals included 8500 time different behavioral states based on their broadband scale-
points with a 1000-Hz sampling rate. With these noise-free free signatures as opposed to narrowband oscillations.
simulation data, we evaluated the performance of separat-
ing and characterizing the power spectrum of the fractal
component as a function of the number of oscillations, the Result
phase distribution of the fractal time series, and the relative
amplitude of the oscillatory component to the fractal We initially evaluated the efficacy of IRASA in compar-
component. In addition, we also applied the IRASA to ison with CGSA by using simulated time-series signals
estimate the PSD of the simulated multifractal time-series with various mixtures of fractal and oscillatory compo-
component that exhibited power-law distributions in mul- nents. Figure 3a, b show the PSD of the simulated signals
tiple separate frequency ranges with different power-law and the PSD of the fractal components estimated by using
exponents. CGSA and IRASA. For CGSA, the extracted fractal PSD
Although the aforementioned theoretical derivation is generally followed a power-law (1/f) trend with ‘‘noisy’’
based on a noise-free model, we also simulated data with deviations especially at the oscillation frequencies where
additive white noises. For an initial testing purpose, we first the spectral peaks remained obvious despite the reduced
mixed one 10 Hz oscillatory signal with the simulated power. Note that CGSA could not lead to the complete
scale-free signal as aforementioned; then we added ran- separation of the fractal and oscillatory components even if
domly generated noise to the mixed signal such that the the oscillation only had a single frequency (as predicted by
signal to noise ratio, defined as the ratio of the variance of its theoretical limitations elaborated in the online supple-
the mixed signal to the variance of the additive noise, mentary material), and its performance became worse
varied at 1, 10, and 100. With these data, we evaluated the when an increasing number of oscillations were included.
robustness of IRASA against noise contamination. In contrast, IRASA could be used to successfully remove
oscillations and obtain an accurate fractal PSD that strictly
Experiment Data and Analyses obeyed the power-law distribution. The residual oscilla-
tion-related spectral peaks were negligible when the signal
We also used the proposed IRASA method to analyze included a small number of oscillations. Even if as many as
in vivo ECoG and MEG data. The ECoG data were 50 oscillations were included, IRASA was still able to
downloaded from the website of the Project Tycho (http:// exclude most oscillations and resulted in a clean power-law
[Link]). The ECoG data were originally collected PSD that well characterized the underlying fractal com-
with a sampling rate of 1 kHz from macaque brain in the ponent with minor interference from the oscillatory com-
eye-closed wakeful resting state and the sleep state by ponent. In addition, IRASA could also be used to estimate
using 128 sensors covering the entire lateral surface of the the PSD of a multifractal time-series component, albeit
left macaque cortex. In this study, we analyzed data from blurred breakpoints between these frequency ranges
two macaques, 4 sessions for each macaque, and each (Fig. 3c).
session included two segments of 5-min time series with With various levels of additive white noises, we evalu-
one in wakefulness and the other in sleep. The MEG data ated the robustness of IRASA against noise contamination.
were obtained from a previously published study (Liu Figure 4 shows that the performance of IRASA in
et al. 2010). It was recorded in three different behavioral extracting the fractal PSD deteriorated with decreasing
states, the eyes-closed, eyes-open and sleep states for a SNR. When the overall SNR was as low as 1, the noise
total of 40 min. We analyzed the MEG data from 5 affected the high frequency range in which the measured
subjects. signal was dominated by the noise; however, the extracted

123
Brain Topogr (2016) 29:13–26 19

Fig. 3 Comparison between the CGSA IRASA IRASA

estimated power spectra of the 10 5 10 5 10 5
fractal component based on Total Total Total

Log PSD (W/Hz)

CGSA (a) and IRASA (b), 10 4 Fractal 10 4 Fractal 10 4 Fractal
given different numbers (1, 2,
50) of oscillation frequencies, 10 3 10 3 10 3
and the result of IRASA in
multiple fractal components (c). 10 2 10Hz 10 2 10Hz 10 2 β= 1
Fractal time series has random
phase uniformly distributed 10 1 1 oscillation 10 1 1 oscillation 10 1
between 0 and 2p and power- 0
10 0 10 0 10
law exponent is set to 1.5 (a, b)
10 0 10
1
10
2
10
0
10
1
10
2
10 0 10 1 10 2
Log Frequency (Hz)
10 5 10 5 10
5
Total Total Total
Log PSD (W/Hz)

10 4 Fractal 10 4 Fractal
Fractal
10 3 10 3
10
4
10Hz 10Hz
10 2 10 2 β = 0.3
20Hz 20Hz
10 1 10 1 2 oscillations β= 1
2 oscillations
10 0 10 0 10
3
10 0 10 1 10 2 10 0 10 1 10 2 10 0 10 1 10 2
Log Frequency (Hz) 5
10 5 10 5 10
Total Total Total
10 4
Log PSD (W/Hz)

10 4 Fractal Fractal Fractal

10 3 10 3
4
10 β = 0.3
10 2 10 2
β= 1
10 1 50 oscillations 10 1 50 oscillations
3 β=2
10 0 10 0 10
10 0 10 1 10 2 10 0 10 1 10 2 10 0 10 1 10 2
Log Frequency (Hz)
(A) (B) (C)

fractal component still followed a clear power-law distri- conditions, respectively. The IRASA-extracted oscillatory
bution in the low frequency range. In general, IRASA component included peaks at the alpha and beta bands, as
performed well with modest to high SNR in a simulation well as the 60 Hz power-line noise. The eyes-open con-
setting. dition had much reduced alpha power and diminished beta
To further test the practical utility of IRASA, we used power relative to the eyes-closed condition. Figures 6 and
IRASA to differentiate fractal and oscillatory neuroelec- 7 show the results of the separated fractal dynamics in
trical activities observed with ECoG and MEG data. Fig- different brain regions as observed with macaque ECoG
ure 5a shows the IRASA-extracted fractal and oscillatory and human MEG data, respectively. These results
PSD of the example ECoG signals recorded from one demonstrate the efficacy of the proposed method to extract
sensor located in the frontal lobe of a macaque brain. The the fractal component from measured electrophysiological
fractal PSD depicted the broadband activity of the ECoG signals.
spectrum. This broadband activity was approximately an Moreover, one can also use IRASA to obtain the spec-
exponential function in the regular coordinate or equiva- trograms showing the temporal dynamics of the fractal and
lently a linear function in the double-logarithmic coordi- oscillatory spectra by using short sliding time windows. To
nate. The PSD of the oscillatory component, obtained by demonstrate this point, we separately extracted the fractal
subtracting the extracted fractal PSD from the raw ECoG and oscillatory spectrograms and obtained the time courses
PSD, provided a clear depiction of alpha (8-13 Hz) of the alpha power (i.e. the average power in 8–13 Hz), the
rhythm. Similarly, Fig. 5b, c show the separated fractal power-law exponent, the log-power intercept, as well as the
and oscillatory components of one MEG sensor over the broadband power from the macaque ECoG and the human
occipital lobe under the eyes-closed and eyes-open MEG signals.

123
20 Brain Topogr (2016) 29:13–26

Fig. 4 Robustness evaluation

of IRASA on extracting the
SNR = inf SNR = 100
4
PSD of fractal component under 10
different levels of additive white Total
noise contamination

log PSD (W/Hz)

Fractal
3
10

2
10

1
10

0
10

SNR = 10 SNR = 1
4
10
log PSD (W/Hz)

3
10

2
10

1
10

0
10 0
10
2 0
10
2
10 10
log frequency (Hz) log frequency (Hz)

Similar to previous studies (Olbrich et al. 2009), the behavioral states may be characterized not only by fre-
distinct behavioral states exhibited different patterns of quency-specific oscillations but also by broadband fractal
oscillations, especially in the alpha frequency range activities.
(Figs. 8a, 9a). We observed much stronger alpha oscilla-
tions during eyes-closed wakefulness relative to the eyes-
open and sleep states (Fig. 9a). Interestingly, we also Discussions
found that the fractal component differed considerably
across these behavioral and arousal states. In the macaque The signal processing technique (IRASA) presented here
ECoG data (Fig. 8), we found that the power-law expo- allows the separation of fractal and oscillatory components
nent was larger in the sleep state, with a mean and the of neurophysiological signals in the spectral domain, with
standard deviation of 1.26 ± 0.49, than the eyes-closed the assumption that neural recordings are formed by the
wakeful state (0.95 ±0.40) (Fig. 8b, c). For the human sum of scale-free and frequency-specific components.
MEG data (Fig. 9), the power-law exponents increased Arguably, such separated oscillatory and scale-free
from 0.81 ± 0.24 to 1.14 ±0.27 as the subjects fell into dynamics may report useful information about rhythmic
sleep from eyes-closed wakefulness. We also observed and arrhythmic dynamics emerging from interacting neural
greater variance in the power-law exponent during the circuits. When applied to various electrophysiological
sleep state than the wakeful states. Moreover, the broad- signals, the values and dynamics of the power-law expo-
band power also differed across these behavioral and nent and broadband power may provide complementary
arousal states, being lower when the eyes were open than perspectives to those that report on neuronal oscillations.
that in the eyes-closed condition and being higher as the However, it awaits future studies to investigate whether the
subjects fell asleep. These results suggest that distinct scale-free broadband component indeed originates from

123
Brain Topogr (2016) 29:13–26 21

−8 −8
ECoG(awake) MEG (eye-closed) x 10 MEG (eye-open)
250 8 x 10 3.5
Total Total

Log PSD
Total

Log PSD
−7
7 −7 3 10
2 Fractal Fractal Fractal

Log PSD
10

(W/Hz)
200 10

(W/Hz)
(W/Hz)

Oscillation 6 Oscillation −8 Oscillation

−8 2.5 10
0
5 10
150 2
10 4 −9
−9
10 10
10 2

PSD
0 101 1.5

PSD
2
PSD

100 10 3 100 101 10

Log Frequency 100 101 102 1 Log Frequency
2 Log Frequency
50 0.5
1
0
0 0
0 10 20 30 40 50 60 0 10 20 30 40 50 60 70 80 0 10 20 30 40 50 60 70 80
Frequency (Hz) Frequency (Hz) Frequency (Hz)
(A) (B) (C)

Fig. 5 Separated power spectra of the fractal and oscillatory components underlying the eyes-closed awake state macaque ECoG (a), as well as
human MEG signals in the eyes-closed (b) and eyes-open (c) conditions

neuronal population activities with fractal dynamics, as described in (Yamamoto and Hughson 1991, 1993; He
well as the specific functional role of such activities. et al. 2010). Combined with the two features mentioned
In the following, we discuss multiple methodological above, this feature utilizes the fact that changing the
considerations as well as speculative implications and resampling factor in IRASA leads to relocation of the
applications to neuroscience research. oscillation-related energy to different frequencies. It fol-
lows that taking the median of the results obtained with
Comparison with CGSA different resampling factors enables IRASA to almost
completely exclude the contribution and interference from
This work is inspired by the earlier CGSA methods (Ya- the oscillatory component in estimating the spectrum of the
mamoto and Hughson 1991, 1993). However, the cross- fractal component. However, this similar ‘trick’ does not
spectral analysis in CGSA suffers from theoretical limita- benefit CGSA in the same way (Figs. 3, S2).
tions (see the online supplementary material). The differ-
ence between IRASA and CGSA has the following aspects. Choosing the Resampling Factors
Firstly, IRASA is based on the auto spectra of the up-
sampled and down-sampled signals, whereas CGSA is We choose the resampling factor (h and 1/h) as non-integer
based on the magnitudes of the cross spectra between the numbers for irregular resampling of the original time ser-
original signal and its up-sampled and down-sampled ies. For any h [ 0, either h or 1/h is greater or equal to 1.
versions. This difference is critical in that it enables Without loss of generality, we choose h [ 1. A large
IRASA to relocate the spectral energy associated with h value would lead to greater dispersion of the spectra of
oscillations away from their occurring frequencies, while the resampled oscillatory component (Fig. 2). However, it
part of such energy always remains at the original oscil- would limit the bandwidth of the final estimate of the
lation frequencies for CGSA. fractal spectrum, because the bandwidth is reduced by a
Secondly, IRASA uses non-integer (irregular) resam- factor of h following the resampling. For example, for a
pling, whereas CGSA, as in (Yamamoto and Hughson typical EEG sampling frequency of 1000 Hz, the band-
1991, 1993; He et al. 2010), uses integer (regular) resam- width is 500 Hz. A resampling factor of 2 would reduce the
pling. Although it requires time-series interpolation with bandwidth to 250 Hz but still covers the frequency range of
slightly higher computation demand, the irregular resam- interest. If the original sampling rate is overly high, it is
pling is more advantageous than regular resampling in that more preferable to use a wider range for h. Within any
the spectra of the resampled oscillations do not overlap range of choice, one may further choose more h values
with the original spectrum when the oscillatory component such that taking the median would likely yield a more
occurs at multiple harmonic frequencies, which is often the accurate and robust estimate of the fractal PSD yet at the
case in realistic neural signals. Although it has not been cost of increasing computational time. Balancing the above
used with CGSA in published studies (to the best of our considerations, we choose the value of h between 1.1 and
knowledge), the irregular resampling is expected to also 1.9 with a 0.05 increment. Although such a choice is
benefit CGSA in addressing issues related to harmonic pragmatic and not entirely principle-driven, it has given
frequencies. rise to robust results in extracting the scale-free broadband
Thirdly, IRASA uses a number of non-integer resam- activity in both simulation and experimental data included
pling factors, whereas CGSA uses only one, at least as in this study.

123
22 Brain Topogr (2016) 29:13–26

Total PSD
Fractal PSD

4 frontal parietal 4 frontal parietal

10 10

log PSD (W/Hz)

2 2
10 10
0 0

log PSD
10 0 1 2 0 1 2 10 0 1 2 0 1 2
10 10 10 10 10 10 10 10 10 10 10 10
4 temporal occipital 4 temporal occipital
10 10
2 2
10 10
frontal temporal 0 0
10 0 10
parietal occipital 1 2 0 1 2
10
0
10
1
10
2
10
0
10
1
10
2
10 10 10 10 10 10
log frequency (Hz)
eye-closed wakefulness sleep
(A) (B) (C)
Fig. 6 Extracted PSD of fractal dynamics from macaque ECoG data regions in the eyes-closed awake condition. c Averaged fractal power
of 128 sensors. a The colored dots are the 128 sensor in different spectra in the sleep condition (Color figure online)
brain regions. b Averaged fractal power spectra in different brain

For any choice of the h values, one can zero-pad all activity. Therefore, another important feature of the scale-
resampled time series to the same length equal to the power free activity has to lie in its phase distribution. A random
of 2. This will permit FFT with the same number of fre- fractal process should possess a random phase distribution
quency points for all resampled signals and allow for within [0, 2p] (Yamamoto and Hughson 1993), whereas
subsequent computation in the frequency domain. fractal processes are not all characterized by random phase
distribution. The IRASA method serves to extract the PSD
Multi-Scale Fractal Dynamics of the signal components that comply with the scale-free
property, as expressed by Eqs. (2) and (3). Such extracted
One assumption for both IRASA and CGSA is that the spectral components contain no phase information, and
fractal dynamics obeys a power law. It should be noted that may or may not reflect a fractal process of biological
this assumption does not always hold true in practice when significance.
a very wide frequency range is of interest. Some neuronal
signals have different frequency scaling properties for Magnitude Versus Phase
different frequency ranges (Hwa and Ferree 2002; Miller
et al. 2009; Robinson 2003). We found that IRASA per- Of note, IRASA, as well as CGSA, can only separate the
formed well even for such multi-scale fractal signals, power spectra of the fractal and oscillatory components of
although it blurred the ‘breakpoints’ where the frequency neurophysiological signals. It is currently difficult, if pos-
ranges of distinct fractal components intersected (Fig. 3c). sible at all, to separate these components in the time
The practical utility of IRASA for multifractal time series domain or uncover their phase differences or relationship.
awaits further investigation and should be taken with This imposes a technical barrier for investigating phase–
caution as of now. phase and phase-amplitude couplings of fractal dynamics
between different frequency points or between different
Fractal Versus Scale-Free spatial locations. Future technical development is highly
desirable to resolve this barrier.
The terms of fractal and scale-free have often been used or
discussed interchangeably. Although these two concepts Other Alternative Solutions
are closely related, there are some conceptual distinctions
perhaps worth noting. For an extreme example, a deter- In addition, it is also of future interest to compare the
ministic signal (e.g. a linear drift) has the scale-free prop- spectral analysis method with other alternative time-series
erty and obeys a 1/f power law. But it is obviously of little analysis methods, e.g. DFA (Peng et al. 1995) or its recent
interest for brain research and may likely be mistaken as variations such as adaptive fractal analysis (AFA) (Gao
the spectral characteristic of the brain’s arrhythmic et al. 2011). Note that although IRASA is based on the

123
Brain Topogr (2016) 29:13–26 23

LF ZF RF

LF ZF RF
LC ZC RC

LT RT

log PSD
LC ZC RC LP ZP RP

LT RT log frequency

LO ZO RO
LP ZP RP

LO ZO
sleep
RO Total PSD
Fractal PSD LF ZF RF
LF ZF RF

LC ZC RC
LC ZC RC

RT LT RT
LT
log PSD
log PSD

LP ZP RP
LP ZP RP
log frequency
log frequency

LO ZO RO LO ZO RO

eye-closed wakefulness eye-open wakefulness

Fig. 7 Extracted PSD of fractal dynamics in human MEG data of 274 regions are separately shown for the eyes-closed, eyes-open and sleep
sensors. The black dots in the upper left figure are the MEG sensor conditions
locations. The fractal power spectra averaged within different brain

Fourier transform, similar concepts might be useful for specific broadband modulation. Arguably, separating the
other time–frequency analysis methods, such as those fractal component may provide more accurate estimation of
based on wavelet transform. the modulation of neuronal oscillations. This is beneficial
when one is interested in the amplitude (not phase) rela-
Implication to Neuronal Oscillations tionship between distinct frequency bands, or in the con-
tributions of various frequency bands to other brain signals,
One motivation of this study is to better extract neural such as neuronal spikes or hemodynamic signals.
oscillations from electrophysiological signals. Here we Unlike the ideal oscillations assumed in our model
assume that oscillations co-occur with the broadband Eq. (1), realistic oscillatory activity is not always localized
activity that approximately follows a power-law distribu- in frequency, but manifests itself in a frequency band. In
tion. As such, the temporal modulation in power or ampli- such cases, the nominal fractal and oscillatory spectra,
tude of oscillatory activity is mixed with the modulation of extracted by using IRASA (and likewise CGSA), may still
the underlying broadband activity. Both modulations con- contain residuals that comprise their corresponding spectral
tribute to the amplitude or power envelops of band-limited behaviors from being strictly scale-free or frequency spe-
signals, or the power fluctuation observed with the time– cific. To partly overcome this, one could choose a larger
frequency analysis. Without separating oscillatory and range of resampling values as permitted by the original
fractal activities, the fluctuation of band-limited signals may sampling frequency and the desired overall frequency
not be entirely specific to the frequency band of interest, but range of interest, as aforementioned for the practical choice
contains a variable fraction owing to the frequency non- of the resampling factors.

123
24 Brain Topogr (2016) 29:13–26

wakefulness sleep

broadband power intercept power-law exponent

spectrogram of fractal component
2
5
frequency (Hz)

80
1
60
sleep wakefulness
wakefulness 0
40 4.5 sleep
0 1 2 3 4 5 6 7 8 9 10
20
log W/Hz
6 4 N(0.95,0.40)
1 2 3 4 5 6 7 8 9 10
time (min)

intersect
4 N (1.26, 0.49)
spectrogram of oscillatory component
90
2
3.5
alpha modulation
frequency (Hz)

0 1 2 3 4 5 6 7 8 9 10
70
log W/Hz
50 4 3

30 3 0 0.5 1 1.5 2 2.5 3

power-law exponent
10 2.5
2 0.5 1 1.5 2 2.5 3
0 1 2 3 4 5 6 7 8 9 10
1 2 3 4 5 6 7 8 9 10
time (min) time (min) power-law exponent

(A) (B) (C)

Fig. 8 a Separated spectrograms of the fractal and oscillatory exponent, intercept and broadband powers derived from the power-
components averaged across all sensors underlying the macaque law fitting of the fractal spectrogram. c Scatter plots and histograms
ECoG signals in the eyes-closed awake (from 0 to 5 min) and sleep of the power-law exponent and the log-power intercept in wakeful-
(from 5 to 10 min) conditions. b Time courses of the power-law ness and sleep

spectrogram of fractal component eyes-closed eyes-open sleep

intercept power-law exponent

2
70
frequency (Hz)

60 1.5 −6
AwakeEyeClosed
50
1 AwakeEyeOpen
40 EC EC EC sleep −6.2
Sleep
30 EO EO 0.5
0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 −6.4
20
log W/Hz
10
−5
−6.6

intersect
0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 −6.8
−6
time (min) N(0.81,0.24)
N(0.85,0.23)
−7 −7
spectrogram of oscillatory component
−8
70 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 −7.2
broadband power

alpha modulation
frequency (Hz)

60 N(1.14,0.27)

50 log W/Hz −7.4

−7.5
40
−7.6
30
−8
20 −7.8 0 0.5 1 1.5 2 2.5
power-law exponent
10
−8.5 −8
0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 0.5 1 1.5 2 2.5 3
0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40
time (min) time (min) power-law exponent

(A) (B) (C)

Fig. 9 a Separated spectrograms of the fractal and oscillatory exponent, the log-power intercept and the broadband power derived
components averaged across all sensors underlying the human MEG from the power-law fitting of the fractal spectrogram. c Scatter plots
signals in the eyes-closed (from 0 to 2 and 4 to 6 min periods), eyes- and histograms the power-law exponent and log-power intercept in
open (from 2 to 4 and 6 to 8 min periods) and sleep (from around the three different states
12 min to the end) conditions. b Time courses of the power-law

Implication to Neuronal Fractal Activity in the temporal fluctuation of fractal dynamics, and to
compare the resulting correlational patterns against those
The fractal dynamics, as the broadband activity underneath based on oscillatory dynamics or the structural networks
neural oscillations, fluctuates over time in a way that seems observed with diffusion MRI. This would help address the
functionally relevant. For example, it varies by voluntary network origin of fractal versus oscillatory dynamics. It is
eyes opening and closing in terms of the broadband power also of potential interest to investigate the spatiotemporal
and the power-law exponent (Fig. 9). As such, the broad- characteristics of fractal dynamics during not only simple
band fractal activity is at least not pure noise, and likely behavioral tasks but also complex cognitive tasks in a rich
bears some functional significance. However, where the and naturalistic behavioral context. This would help
fractal dynamics comes from and how it serves brain address the functional significance of fractal dynamics.
functioning remain poorly understood and await future However, the proposed IRASA method is not intended
research (He et al. 2010; He 2014). For example, it is of to address the generative mechanism of the scale-free
potential interest to investigate inter-regional correlations broadband signal or the frequency-specific narrowband

123
Brain Topogr (2016) 29:13–26 25

signal observed in neural recordings. As a signal process- Buzsaki G, Anastassious CA, Koch C (2012) The origin of
ing method, the IRASA method by itself offers little extracellular fields and currents—EEG, ECoG, LFP and spikes.
Nat Rev Neurosci 13(6):407–420
information about the theoretical basis and functional role Ciuciu P, Varoquaux G, Abry P, Sadaghiani S, Kleinschmidt A
of the apparent scale-free or oscillatory activity in the (2012) Scale-free and multifractal time dynamics of fMRI
context of neural synchrony, coupling or computation. signals during rest and task. Front Physiol 3:186
El Boustani S, Marre O, Béhuret S, Baudot P, Yger P, Bal T, Frégnac
Y (2009) Network-state modulation of power-law frequency-
Differentiate Functional States with Fractal scaling in visual cortical neurons. PLoS Comput Biol
Characteristics 5(9):e1000519
Engel AK, Gerloff C, Hilgetag CC, Nolte G (2013) Intrinsic coupling
A large amount of literature has reported on neural oscilla- modes: multiscale interactions in ongoing brain activity. Neuron
80(4):867–886
tions in a variety of brain states (Başar et al. 2001; Buzsáki Fransson P, Metsäranta M, Blennow M, Åden U, Lagercrantz H,
and Draguhn 2004; Gray 1994; Hanslmayr et al. 2011; Vanhatalo S (2013) Early development of spatial patterns of
Lakatos et al. 2008; Rinzel and Ermentrout 1998; Schabus power-law frequency scaling in fMRI resting-state and EEG data
et al. 2011). For example, the brain’s vigilance and attention in the newborn brain. Cereb Cortex 23(3):638–646
Freeman WJ (2007) Scale-free neocortical dynamics. Scholarpedia
level has been measured mainly by the oscillation power of a 2:1357
single frequency band, e.g. the alpha band (8–13 Hz), or the Gao J, Hu J, Tung WW (2011) Facilitating joint chaos and fractal
relative amplitudes of several frequency bands, especially analysis of biosignals through nonlinear adaptive filtering. PLoS
the alpha, theta and delta bands (Olbrich et al. 2009; Wong One 6(9):e24331. doi:10.1371/[Link].0024331
González J, Gamundi A, Rial R, Nicolau MC, de Vera L, Pereda E
et al. 2013). Combined with other event-like signatures, e.g. (1999) Nonlinear, fractal, and spectral analysis of the EEG of
k-complex, spindles, oscillatory patterns have been serving lizard, Gallotia galloti. Am J Physiol Regul Integr Comp Physiol
as the main clues for sleep staging. 277(1):R86–R93
However, how the broadband fractal dynamics varies Gray CM (1994) Synchronous oscillations in neuronal systems:
mechanisms and functions. J Comput Neurosci 1(1–2):11–38
across behavioral states remains poorly understood. Fol- Hanslmayr S, Gross J, Klimesch W, Shapiro KL (2011) The role of
lowing IRASA, one can quantify the broadband fractal alpha oscillations in temporal attention. Brain Res Rev
characteristics by using the power-law exponent (He et al. 67(1):331–343
2010; He 2014), the intercept of the power-law function at He BJ (2014) Scale-free brain activity: past, present, and future.
Trends Cognit Sci 18(9):480–487
zero log-frequency, and the broadband power (Miller et al. He BJ, Zempel JM, Snyder AZ, Raichle ME (2010) The temporal
2009) as a linear function of the exponent and the intercept. structures and functional significance of scale-free brain activity.
These quantities collectively define the feature space for Neuron 66(3):353–369
potential classification of various functional states based Hwa RC, Ferree TC (2002) Scaling properties of fluctuations in the
human electroencephalogram. Phys Rev E 66(2):021901
merely on the fractal characteristics. As an initial proof of Lakatos P, Karmos G, Mehta AD, Ulbert I, Schroeder CE (2008)
concepts, the broadband power of fractal dynamics differed Entrainment of neuronal oscillations as a mechanism of atten-
notably across behavioral and arousal states: higher in tional selection. Science 320(5872):110–113
sleep than wakefulness, and higher in eyes-closed condi- Liu Z, Fukunaga M, de Zwart JA, Duyn JH (2010) Large-scale
spontaneous fluctuations and correlations in brain electrical
tion than the eyes-open condition. activity observed with magnetoencephalography. Neuroimage
51(1):102–111
Acknowledgments The research was supported in part by NIH Mandelbrot BB, Van Ness JW (1968) Fractional Brownian motions,
R01MH104402. The authors are thankful to Dr. Shao-Chin Hung for fractional noises and applications. SIAM Rev 10(4):422–437
proof reading and constructive comments, to Drs. Masaki Fukunaga Manning JR, Jacobs J, Fried I, Kahana MJ (2009) Broadband shifts in
and Jeff Duyn for assistance in collecting the MEG data, and to Dr. local field potential power spectra are correlated with single-
Naotaka Fujii for making the macaque ECoG data publicly available. neuron spiking in humans. J Neurosci 29(43):13613–13620
Miller KJ, Sorensen LB, Ojemann JG, den Nijs M (2009) Power-law
scaling in the brain surface electric potential. PLoS Comput Biol
5(12):e1000609
Nunez PL, Srinivasan R (2006) Electrical fields of the brain: the
References neurophysics of EEG. Oxford University Press, Oxford
Olbrich S, Mulert C, Karch S, Trenner M, Leicht G, Pogarell O, Hegerl
Başar E, Başar-Eroglu C, Karakaş S, Schürmann M (2001) Gamma, U (2009) EEG-vigilance and BOLD effect during simultaneous
alpha, delta, and theta oscillations govern cognitive processes. EEG/fMRI measurement. Neuroimage 45(2):319–332
Int J Psychophysiol 39(2):241–248 Peng CK, Havlin S, Stanley HE, Goldberger AL (1995) Quantifica-
Bassett GW Jr (1991) Equivariant, monotonic, 50 % breakdown tion of scaling exponents and crossover phenomena in nonsta-
estimators. Am Stat 45(2):135–137 tionary heartbeat time series. Chaos 5:82–87
Bullmore E, Sporns O (2009) Complex brain networks: graph Pereda E, Gamundi A, Rial R, Gonzalez J (1998) Non-linear behaviour
theoretical analysis of structural and functional systems. Nat Rev of human EEG: fractal exponent versus correlation dimension in
Neurosci 10(3):186–198 awake and sleep stages. Neurosci Lett 250(2):91–94
Buzsáki G, Draguhn A (2004) Neuronal oscillations in cortical Rinzel J, Ermentrout GB (1998) Analysis of neural excitability and
networks. Science 304(5679):1926–1929 oscillations. Methods Neuronal Model 2:251–292

123
26 Brain Topogr (2016) 29:13–26

Robinson PA (2003) Interpretation of scaling properties of electroen- Wong CW, Olafsson V, Tal O, Liu TT (2013) The amplitude of the
cephalographic fluctuations via spectral analysis and underlying resting-state fMRI global signal is related to EEG vigilance
physiology. Phys Rev E 67(3):032902 measures. Neuroimage 83:983–990
Schabus M, Pelikan C, Chwala-Schlegel N, Weilhart K, Roehm D, Yamamoto YOSHIHARU, Hughson RL (1991) Coarse-graining
Donis J, Klimesch W (2011) Oscillatory brain activity in spectral analysis: new method for studying heart rate variability.
vegetative and minimally conscious state during a sentence J Appl Physiol 71(3):1143–1150
comprehension task. Funct Neurol 26(1):31 Yamamoto Y, Hughson RL (1993) Extracting fractal components
Siegel M, Donner TH, Engel AK (2012) Spectral fingerprints of large- from time series. Phys D 68(2):250–264
scale neuronal interactions. Nat Rev Neurosci 13(2):121–134

123