EEG analysis turns many correlated time series (channels × time × trials × subjects) into a few interpretable measures and statistical conclusions. This tutorial gives an overview of the domains and measures of EEG analysis and of the statistical problems specific to EEG data: correlated samples, the multiple comparisons problem and non-stationarity.
Keep analyses simple: prefer straightforward, interpretable methods; complexity should serve clarity, not obscure it.
Data integrity
Stay close to the data: avoid transformations and assumptions that distance you from the raw observations; always look at the raw signals.
Record clean signals: check electrode impedances before recording (and between blocks in long sessions); no processing can fully recover a bad recording.
Understanding and transparency
Understand your analyses: be able to explain every step of your methodology and its rationale.
Write a clear Methods section: report all parameters (filters, epochs, rejection criteria, software versions) so that others can replicate your work.
Avoiding pitfalls
Prevent over-fitting: make sure that models generalize beyond the data used to fit them (e.g., cross-validation).
Validate the pipeline on simulated data with a known effect, and perform sanity checks (visualization, consistency tests).
Plan the analysis (time windows, channels, frequency bands, statistical tests) before looking at the data, ideally in a pre-registration.
Reproducibility and validation
Design the study so that others can reproduce your findings; share data and code when possible.
Split-half replication: analyze one half of the data and replicate the result with the other half.
Collect extra data: plan for at least 10% more participants or trials than needed to compensate for artifacts and exclusions.
Exploring parameters
Test whether the results are robust to reasonable parameter changes (e.g., a high-pass cutoff of 0.1 Hz vs 0.5 Hz, an alpha band of 8–12 Hz vs 8–13 Hz).
Pilot testing
Run pilot studies to test stimuli, equipment, timing and the participants’ understanding of the instructions.
Informed decision-making
Base methods and analyses on evidence, best practices and the specific goals of your research.
Domains of Signal Analysis
A signal can be analyzed in different domains, each showing a different aspect of the data:
Time domain: the signal as a function of time, \(x(t)\) (continuous-time) or \(x[n]\) (discrete-time); shows how the signal evolves over time (amplitudes, latencies, waveforms).
Frequency domain: the signal as a function of frequency, \(X(f)\); shows which frequencies are present and with what amplitude and phase.
Time-frequency domain: the frequency content as a function of time, \(X(t, f)\); shows how oscillations change over time.
Spatial domain: the signal as a function of position on the scalp (topography) or in the brain (sources).
Transforms move a signal between domains:
Fourier transform: time domain → frequency domain; the main tool for analyzing the frequency content of a signal. Its inverse transforms the spectrum back to the time domain.
Short-time Fourier transform, wavelet transform and Hilbert transform (of band-pass filtered signals): time domain → time-frequency domain.
Laplace transform (continuous-time) and Z-transform (discrete-time): generalizations of the Fourier transform used to describe systems such as filters (transfer functions, stability, transient responses).
time domain x[n] ──── Fourier transform ────► frequency domain X(f)
◄─── inverse transform ─────
EEG Measures
Amplitude-Based Measures
Voltage: the amplitude of the signal, in microvolts (μV), e.g., the amplitude and latency of an ERP component.
Band amplitude: the amplitude of activity in a frequency band. NeuroAnalyzer’s band definitions (band_frq()):
Band
Frequency range
Delta
0.5–4 Hz
Theta
4–8 Hz
Alpha
8–13 Hz
Beta
13–30 Hz
Gamma
> 30 Hz
Band limits are conventions that differ between studies (e.g., alpha 8–12 Hz); report the ones you use. Sub-bands are also available (e.g., :alpha_lower, :beta_higher, :gamma_1); band_frq() clamps the upper limit to the Nyquist frequency:
bands = [:delta, :theta, :alpha, :beta, :gamma]for b in bandsprintln(rpad(b, 6), band_frq(eeg; band = b), " Hz")end
Pattern amplitude: the amplitude of characteristic EEG patterns (graphoelements), such as:
spike-and-wave complexes,
sharp waves,
K-complexes,
sleep spindles (sigma activity, 11–16 Hz).
Spectral Measures
After transforming the signal to the frequency domain (e.g., with the FFT or Welch’s method, see frequency analysis), the activity can be described by:
Power spectral density (PSD): power per frequency bin (μV²/Hz).
Band power (μV²): the PSD integrated over a frequency band. Power is the squared amplitude, so it emphasizes the largest components (e.g., the alpha peak) and has a skewed distribution; log-power or decibels (dB) are often used for statistics.
Amplitude spectrum (μV): the square root of power; it compresses the differences between large and small components, making the small ones more visible.
Peak frequency: the frequency of maximum power within a band, e.g., the individual alpha frequency (IAF), see band_mpower().
Relative power: the power of a band divided by the total (broadband) power, usually expressed in %:
Relative power removes differences in overall signal strength between subjects and recordings (e.g., due to skull thickness or electrode impedance), but a change in one band changes the relative power of all the other bands too.
Absolute and relative band power at O1, relative to the total power in 0.5–30 Hz:
# broadband power: the reference for relative powerbp_total =band_power(eeg; ch ="O1", flim = (0.5, 30))[1]for b in bands[1:4] bp =band_power(eeg; ch ="O1", flim =band_frq(eeg; band = b))[1]println(rpad(b, 6),round(bp; digits =1), " μV² (",round(100* bp / bp_total; digits =1), "%)")end
The shares add up to slightly more than 100%, as the neighbouring bands share their boundary frequency bins (e.g., 4 Hz belongs to both delta and theta).
Coherence
Coherence measures the linear relationship between two signals as a function of frequency. Magnitude-squared coherence is the cross-spectrum normalized by the power spectra of both signals:
It ranges from 0 (no linear relationship at \(f\)) to 1 (a constant phase difference and amplitude ratio at \(f\) across all segments or trials). Coherence between nearby channels is inflated by volume conduction, see coherence.
Domains of EEG Analysis
Time-domain analysis: maps the temporal sequence of brain processes with millisecond precision.
Event-related potentials (ERPs) / fields (ERFs), see Introduction to ERPs/ERFs: averaging many trials time-locked to an event removes activity that is not time- and phase-locked to it. Limitations of averaging:
it assumes that the response is identical in every trial; latency jitter smears and attenuates the components,
activity that is not phase-locked to the event (induced oscillations) cancels out,
trial-to-trial variability, needed e.g. for single-trial connectivity measures, is lost,
the signal-to-noise ratio of single trials is low, so many trials are needed.
Microstates: quasi-stable scalp topographies lasting about 60–120 ms, which are thought to reflect successive states of large-scale brain networks.
Complexity measures: fractal dimension, entropy and other measures of signal irregularity and self-similarity.
Frequency-domain analysis: the frequency content of the signal (power spectrum, band power). It assumes that the signal is stationary within the analysis window and loses all information about timing within it.
Time-frequency analysis: power and phase as a function of time and frequency (short-time Fourier transform, wavelet convolution, filter-Hilbert), see time-frequency analysis. It captures both evoked (phase-locked) and induced (non-phase-locked) activity and is close to the neurophysiology of oscillations, but is limited by the time-frequency trade-off: better frequency resolution requires longer windows and thus worse temporal resolution.
Spatial-domain analysis: the distribution of activity across the scalp or the brain.
Topography: maps of the scalp distribution of a measure, see topographical plots.
Spatial filters: re-referencing, surface Laplacian / current source density (CSD).
Source separation: independent component analysis (ICA), principal component analysis (PCA).
Multiway analysis: analyzes data across several dimensions at once (e.g., channel × time × frequency), e.g., with tensor decompositions.
Major Types of EEG Analysis
Core analysis techniques
Spectral analysis: power in frequency bands (delta, theta, alpha, beta, gamma).
Event-related potentials (ERPs): brain responses to specific events, such as the P300 or the error-related negativity (ERN).
EEG source imaging (ESI): estimation of the location, orientation and distribution of the EEG generators by solving the inverse problem.
Additional analysis types
Steady-state responses: responses to periodic stimulation at the stimulation frequency and its harmonics, e.g., steady-state visual evoked potentials (SSVEP) and auditory steady-state responses (ASSR).
Hemispheric asymmetry: differences between homologous left and right channels, e.g., frontal alpha asymmetry \(\ln P_{\alpha}(F4) - \ln P_{\alpha}(F3)\), see lateralization and symmetry.
Coherence and other connectivity measures: synchrony between brain regions.
Power spectral density (PSD): power at each frequency.
Entropy-based measures:
Entropy: irregularity or unpredictability of the signal (e.g., Shannon, approximate, sample entropy).
Multiscale entropy (MSE): sample entropy computed at several time scales (of the signal progressively coarse-grained by averaging consecutive samples).
Differential entropy (DE): the entropy of a continuous random variable; for a Gaussian signal of variance \(\sigma^2\) (e.g., a band-pass filtered EEG), \(DE = \frac{1}{2} \ln(2 \pi e \sigma^2)\), i.e., the logarithm of band power up to a constant.
Differential asymmetry (DASM): the difference of DE between symmetric left and right channels, \(DE_{left} - DE_{right}\).
Rational asymmetry (RASM): the ratio of DE between symmetric left and right channels, \(DE_{left} / DE_{right}\).
Asymmetry (ASM): DASM and RASM combined into one feature vector.
Differential caudality (DCAU): the difference of DE between frontal and posterior channels, \(DE_{frontal} - DE_{posterior}\).
Subject-level (first-level) analysis: analysis of the trials or epochs within each subject.
Descriptive statistics: mean, variance or standard deviation across trials (e.g., the ERP is the mean of the trials).
Test statistics: e.g., t-values of the difference between conditions within the subject.
Model parameters: parameters of a model fitted to the trials (e.g., regression coefficients of a general linear model).
Contrasts: comparisons between:
time periods within a condition (e.g., post-stimulus activity vs baseline),
conditions (e.g., task vs control).
Group-level (second-level) analysis: statistics on the subject-level measures; the subjects are the observations.
Between-subject effects: differences between groups (e.g., patients vs healthy controls).
Within-subject effects: differences between conditions within the same subjects (e.g., pre-treatment vs post-treatment).
Interactions between between-subject and within-subject factors (e.g., a treatment effect that differs between groups).
The two-level approach gives each subject the same weight, regardless of the number of trials. Mixed-effects (hierarchical) models analyze all trials at once, with random effects for subjects (and items), and use the trial-level information more efficiently, e.g., when subjects have different numbers of trials.
Example Data
The statistical methods below are demonstrated on simulated group data: 20 subjects, one channel, epochs from −200 to 800 ms. For each subject, the data are the difference between two conditions (condition B − condition A), as used in a within-subject (paired) design. Condition B has an additional component of 1.5 μV peaking at 400 ms; the background EEG is simulated as temporally smooth noise (SD 2 μV), since neighbouring EEG samples are correlated.
Random.seed!(42)fs =250t =collect(-0.2:(1/ fs):0.8)nt =length(t)n_subj =20# condition effect: Gaussian-shaped component, 1.5 μV at 400 mseffect =1.5.*exp.(-(t .-0.4) .^2./ (2*0.05^2))# temporally smooth Gaussian noise: white noise smoothed by# a Gaussian kernel (σ = 5 samples = 20 ms), standard deviation sdfunctionsmooth_noise(n; sd =2.0, σ =5) k =exp.(-(-3σ:3σ) .^2./ (2* σ^2)) k ./=sqrt(sum(abs2, k)) x =randn(n +6σ)return sd .* [sum(x[i:(i +6σ)] .* k) for i in1:n]end# subject-level differences (B - A), shape (subjects, time points)d =zeros(n_subj, nt)for s in1:n_subj d[s, :] = effect .+smooth_noise(nt)end
fig =Figure(size = fig_size)ax =Axis(fig[1, 1]; ax_style..., title ="Condition difference (B - A)", xlabel ="Time [s]", ylabel ="Amplitude [μV]")for s in1:n_subjlines!(ax, t, d[s, :]; color = (:gray, 0.4))endlines!(ax, t, vec(mean(d; dims =1)); color =:black, linewidth =3)lines!(ax, t, effect; color =:red, linestyle =:dash)fig
Gray: single subjects; black: group mean; red: the true effect.
Permutation Testing
Overview
Permutation testing is a non-parametric method to assess the significance of an effect without assumptions about the distribution of the data. The data are repeatedly reshuffled under the null hypothesis to build the null distribution of the test statistic.
The data may only be shuffled in ways that do not change their distribution under the null hypothesis (exchangeability):
Between-subject design (e.g., patients vs controls): the group labels are randomly permuted across subjects.
Within-subject design, subject level: the condition labels are randomly permuted across trials.
Within-subject design, group level: the condition labels are swapped within subjects, i.e., the sign of each subject’s difference (B − A) is randomly flipped.
How it works
Calculate the test statistic (e.g., mean difference, t-value) for the original data.
Shuffle the data (permute the labels or flip the signs) and recalculate the test statistic; repeat many times (e.g., 1000–10000).
The test statistics of all permutations form the null distribution: the range of values expected if there is no difference between the conditions.
The p-value is the proportion of permutations with a test statistic at least as extreme as the observed one. Counting the original data as one of the permutations, \(p = (k + 1) / (n_{perm} + 1)\), where \(k\) is the number of permutations with a statistic at least as extreme; this p-value is never 0.
If the p-value is below the significance level (e.g., 0.05), the null hypothesis is rejected.
Permutation test of the difference at 400 ms (sign flipping) compared with the parametric one-sample t-test:
# one-sample t-valuetval(x) =mean(x) / (std(x) /sqrt(length(x)))i_400 =argmin(abs.(t .-0.4))x = d[:, i_400]t_obs =tval(x)n_perm =5000# null distribution: random sign flips of the subject differencest_null = [tval(x .*rand([-1, 1], n_subj)) for _ in1:n_perm]k =count(abs.(t_null) .>=abs(t_obs))p_perm = (k +1) / (n_perm +1)println("t = ", round(t_obs; digits =2))println("permutation p = ", round(p_perm; digits =4))println("t-test p = ", round(pvalue(OneSampleTTest(x)); sigdigits =2))
t = 7.12
permutation p = 0.0002
t-test p = 9.1e-7
None of the 5000 permutations reached the observed t-value, so the permutation p-value is at its minimum \(1 / 5001\); a smaller p-value requires more permutations.
NeuroAnalyzer provides perm_cmp() for permutation-based comparison of 3-dimensional arrays (e.g., spectrograms of two conditions) and cmp_stat() for the proportion of a statistic distribution above or below a value.
The Multiple Comparisons Problem
Mass univariate analysis
Mass univariate analysis treats each sample (channel × time point × frequency) as a separate variable and applies a statistical test to all of them (or to a preselected subset), to find the samples with an effect. This means thousands to millions of tests, e.g., 64 channels × 250 time points × 30 frequencies = 480 000 tests.
The problem
At a significance level \(\alpha\), each test has a probability \(\alpha\) of a false positive (type I error) when there is no effect. Across \(m\) independent tests, the probability of at least one false positive, the family-wise error rate (FWER), is:
\[
FWER = 1 - (1 - \alpha)^m
\]
for m in [1, 10, 100, 1000]println(rpad("m = $m:", 10),"FWER = ", round(1-0.95^m; digits =3))end
m = 1: FWER = 0.05
m = 10: FWER = 0.401
m = 100: FWER = 0.994
m = 1000: FWER = 1.0
With enough tests, “significant” samples are found somewhere in the data even when there is no effect anywhere. To keep the FWER at \(\alpha\), the threshold of the individual tests must be adjusted.
EEG samples are not independent: neighbouring time points, frequencies and channels are correlated, so the actual FWER is lower than the formula predicts, but still far above \(\alpha\). 1000 simulated experiments with no effect at all (pure noise, 251 time points):
n_sim =1000false_pos =0for _ in1:n_sim d0 =reduce(vcat, [smooth_noise(nt)' for _ in 1:n_subj]) p0 = [pvalue(OneSampleTTest(d0[:, i])) for i in1:nt]global false_pos +=any(p0 .<0.05)end
In most of the experiments without any effect at least one time point was “significant”.
Reducing the number of comparisons
Region of interest (ROI) approach: test only predefined channels, sources, time windows or frequency bands (e.g., the P300 at Pz in 300–500 ms); this dramatically reduces the number of tests.
Collapsing the data: instead of testing each sample, summarize the data (e.g., the mean amplitude in a time window, the peak latency of a component, the band power) and test the summary measures.
Avoiding circular inference (“double dipping”)
The ROI must be chosen independently of the effect being tested. Selecting the channels, time windows or frequencies where the difference between conditions is largest and then testing that difference in the same data biases the test towards significance: even pure noise has a “largest difference” somewhere.
Valid ways to define an ROI:
from previous studies or a pre-registration,
from an independent data set (e.g., a pilot study or the other half of the data),
from data that are orthogonal to the contrast, e.g., the average of all conditions (the ERP averaged over conditions A and B) to find the latency of a component, then compare A vs B in that window.
Solutions to the Multiple Comparisons Problem
Bonferroni correction: each test uses the threshold \(\alpha / m\). It controls the FWER for any dependence between tests, but it ignores the correlations between EEG samples, so it is very conservative (low statistical power, many type II errors) when the tests are correlated.
False discovery rate (FDR): controls the expected proportion of false positives among the significant results instead of the probability of any false positive. It is less stringent than the Bonferroni correction (fewer type II errors), at the cost of some false positives among the significant results. Benjamini-Hochberg procedure:
sort the \(m\) p-values in ascending order: \(p_{(1)} \le p_{(2)} \le \dots \le p_{(m)}\),
find the largest \(k\) such that \(p_{(k)} \le \frac{k}{m} q\), where \(q\) is the desired FDR (e.g., 0.05),
reject the null hypotheses of the tests with the \(k\) smallest p-values.
Maximum statistic: each test statistic is compared with the null distribution of the maximum test statistic across all tests. If the maximum of the null data exceeds the threshold in only 5% of cases, the probability of any false positive is 5%, so this controls the FWER while taking the correlations into account (fewer type II errors than Bonferroni). For parametric tests the null distribution can be derived from random field theory; for EEG it is usually obtained by permutation (the maximum statistic of each permutation).
Cluster-based correction: only effects that form contiguous clusters (in time, frequency or space) are considered; cluster sizes or masses are compared with the null distribution of the maximum cluster statistic (see below).
Cluster-Based Permutation Test
Overview
Cluster-based analysis uses the fact that real effects in EEG extend over neighbouring samples (time points, frequencies, channels), while false positives are scattered. The test is performed on clusters of contiguous samples instead of on individual samples, which greatly increases sensitivity while controlling the FWER (Maris and Oostenveld 2007).
Key steps
Sample-level statistics: calculate a test statistic (e.g., a t-value) for each sample, as in the mass univariate analysis.
Cluster formation: samples whose statistic exceeds a cluster-forming threshold (e.g., the t-value corresponding to p < 0.05) are grouped into clusters of adjacent samples with the same sign.
Cluster statistic: each cluster is quantified, usually by the sum of its test statistics (cluster mass) or by the number of samples (cluster size).
Null distribution: the data are permuted (labels permuted or signs flipped, as above), steps 1–3 are repeated, and the maximum cluster statistic of each permutation is recorded.
Hypothesis testing: the p-value of each observed cluster is the proportion of permutations whose maximum cluster statistic is at least as large.
Interpretation
A significant cluster means that the conditions differ; the test controls the FWER at the level of clusters, not of individual samples.
It does not establish that the effect is present at each sample of the cluster, nor the exact onset, offset or spatial extent of the effect: these depend on the cluster-forming threshold and the noise (Sassenhagen and Draschkow 2019).
The cluster-forming threshold is arbitrary but must be fixed before the analysis; it affects which effects are found (high thresholds favour strong, short effects; low thresholds weak, long-lasting ones), but not the validity of the test.
Comparison of the Methods
All methods applied to the example data. First, the t-value and p-value at each time point:
# one-sample t-values of all columns (time points)tvals(x) =vec(mean(x; dims =1) ./ (std(x; dims =1) ./sqrt(size(x, 1))))tv =tvals(d)pv = [pvalue(OneSampleTTest(d[:, i])) for i in1:nt]
Benjamini-Hochberg FDR:
functionfdr_bh(p; q =0.05) m =length(p) o =sortperm(p)# largest k for which p(k) ≤ k / m × q k =findlast(p[o] .<= q .* (1:m) ./ m) sig =falses(m)isnothing(k) || (sig[o[1:k]] .=true)return sigend
Maximum statistic (permutation with sign flipping):
max_t = [maximum(abs, tvals(d .*rand([-1, 1], n_subj))) for _ in1:n_perm]# FWER-corrected threshold: 95th percentile of the maximum |t|thr_max =quantile(max_t, 0.95)println("max-t threshold: |t| > ", round(thr_max; digits =2))
max-t threshold: |t| > 3.89
Cluster-based permutation test:
# clusters: runs of adjacent time points with |t| above# the threshold and of the same signfunctionclusters(tv, thr) cl =UnitRange{Int64}[] i =1while i <=length(tv)ifabs(tv[i]) > thr j = iwhile j <length(tv) &&abs(tv[j +1]) > thr &&sign(tv[j +1]) ==sign(tv[i]) j +=1endpush!(cl, i:j) i = j +1else i +=1endendreturn clend# cluster mass: sum of the t-values in the clustercl_mass(tv, cl) = [abs(sum(tv[c])) for c in cl]# cluster-forming threshold: t-value for p < 0.05 (two-tailed)thr_cl =quantile(TDist(n_subj -1), 0.975)cl =clusters(tv, thr_cl)m_obs =cl_mass(tv, cl)# null distribution of the maximum cluster massmax_m =zeros(n_perm)for i in1:n_perm tp =tvals(d .*rand([-1, 1], n_subj)) m =cl_mass(tp, clusters(tp, thr_cl)) max_m[i] =isempty(m) ? 0.0:maximum(m)endsig_cl =falses(nt)for (c, m) inzip(cl, m_obs) p = (count(max_m .>= m) +1) / (n_perm +1) p <0.05&& (sig_cl[c] .=true)println("cluster ", round(t[c[1]]; digits =3), "–",round(t[c[end]]; digits =3), " s: mass = ",round(m; digits =1), ", p = ", round(p; digits =4))end
uncorrected 22 time points
Bonferroni 11 time points
FDR 15 time points
max-t 13 time points
cluster 22 time points
The dashed lines mark the cluster-forming threshold (uncorrected p < 0.05). The uncorrected tests find the most time points, but give no control of false positives. The Bonferroni correction is the most conservative: the 251 tests are strongly correlated, so dividing α by 251 over-corrects. The max-t and FDR corrections are less conservative, and the cluster-based test detects the whole cluster around 400 ms, at the cost of no sample-level inference.
None of the methods recovers the full extent of the true effect (about 0.3–0.5 s, red line in the first figure): where significance starts and ends depends on the noise, the number of subjects and the method, so the borders of a significant time window are not estimates of the onset and offset of the effect.
How conservative the Bonferroni correction is here can be checked by simulation: the FWER of the Bonferroni and FDR corrections in experiments without any effect:
fp_bonf =0fp_fdr =0for _ in1:n_sim d0 =reduce(vcat, [smooth_noise(nt)' for _ in 1:n_subj]) p0 = [pvalue(OneSampleTTest(d0[:, i])) for i in1:nt]global fp_bonf +=any(p0 .<0.05/ nt)global fp_fdr +=any(fdr_bh(p0))end
The Bonferroni correction keeps the FWER well below the nominal 0.05. When there is no effect at all, every significant result is a false one, so FDR control is equivalent to FWER control and the FWER of the FDR procedure is also below 0.05.
Stationarity
Definition
A signal is stationary if its statistical properties (mean, variance, autocovariance and thus the power spectrum) do not change over time. In an ensemble of signals (e.g., many trials), the distributions of amplitudes are then the same at every time point.
Strict stationarity: the whole distribution is constant over time.
Weak (wide-sense) stationarity: the mean and autocovariance are constant over time; a signal may be stationary in the mean but not in the variance.
Ergodic process
A stationary process is ergodic if the statistics of a single, long realization (time averages) equal the statistics of the ensemble (averages across realizations); then one long recording is representative of the process.
Non-stationarity in EEG
EEG signals are non-stationary: the states of neuronal assemblies change continuously with cognitive and physiological states (e.g., attention, drowsiness, sleep stages, task events).
The non-stationarities are themselves informative, as they reflect neural events and processes.
Consequences
Many methods, such as the Fourier transform, assume stationarity. The spectrum of a non-stationary signal is the average over the analysis window: it does not show when the frequency content changed, and abrupt changes (e.g., transients, artifacts) spread power over a broad range of frequencies.
Addressing non-stationarity
EEG is approximately stationary over short periods, so it is analyzed in short segments:
Adaptive segmentation: segment boundaries placed where the statistical properties of the signal change.
If the signal is not stationary in the mean (e.g., slow drifts):
Detrending: remove linear or polynomial trends.
Pre-whitening: e.g., the first derivative (difference) of the signal, which attenuates slow dynamics.
Time-varying mean subtraction: subtract a moving average from the signal.
Filtering: high-pass, low-pass or band-pass filters to keep the frequency components of interest.
Stationarity can be assessed with NeuroAnalyzer’s stationarity() function, see stationarity.
Univariate and Multivariate Analysis
Univariate analysis
Analyzes one channel at a time: each channel is evaluated independently and the results are summarized afterwards, e.g., “80% of channels showed a spectral peak at 23.4 Hz”.
Multivariate analysis
Analyzes multiple channels simultaneously:
It uses the relationships between the measurements (e.g., the covariance of all EEG channels).
It extracts features from the entire dataset, capturing interactions among channels (e.g., spatial filters, ICA, PCA, decoding).
It requires multichannel data.
Multivariate data with a given covariance (correlation) structure can be simulated using the Cholesky decomposition\(\Sigma = U^T U\): multiplying uncorrelated white noise by \(U\) gives data with covariance \(\Sigma\):
# channels 1 and 2 correlated at 0.8, channel 3 uncorrelatedΣ = [1.00.80.0;0.81.00.0;0.00.01.0]U =cholesky(Σ).U# 10 000 time points × 3 channelsx_mv =randn(10_000, 3) * U# correlation matrix of the simulated channels (+ 0.0 avoids -0.0)round.(cor(x_mv); digits =2) .+0.0
Signal decomposition separates a signal into distinct components, revealing patterns, frequencies or modes that are not apparent in the raw signal.
Common methods
Fourier transform (FT): decomposes a signal into sinusoids of different frequencies, revealing its spectral content.
Wavelet transform (WT): a time-frequency representation, showing both when and at what frequency a feature occurs.
Empirical mode decomposition (EMD): adaptively decomposes a signal into intrinsic mode functions (IMFs), oscillatory modes embedded in the data.
Singular value decomposition (SVD): factorizes a data matrix (e.g., channels × time) into orthogonal components, \(X = U S V^T\); used for dimensionality reduction and denoising.
Principal component analysis (PCA): transforms correlated variables into uncorrelated components ordered by explained variance; equivalent to the SVD of the mean-centered data.
Independent component analysis (ICA): separates the signal into statistically independent components, e.g., brain activity and artifacts (eye blinks, muscle activity, heartbeat).
Proper orthogonal decomposition (POD): extracts the spatial or temporal modes that capture most of the energy of the signal; mathematically equivalent to PCA.
Dynamic mode decomposition (DMD): approximates the dynamics by a linear system and identifies spatial modes, each with its own frequency and growth or decay rate.
Key concepts
Mode: a distinct component of the signal, often representing a specific pattern or oscillation.
Multivariate empirical mode decomposition (MEMD):
Extends EMD to multivariate signals, using projections of the signal in many directions to define multidimensional envelopes of the local extrema.
Simultaneous decomposition: all channels are decomposed together, so the scales common to several channels appear in IMFs with the same index (mode alignment).
References
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society: Series B. 1995;57(1):289–300.
Cohen MX. Analyzing Neural Time Series Data: Theory and Practice. MIT Press; 2014.
Costa M, Goldberger AL, Peng CK. Multiscale entropy analysis of complex physiologic time series. Physical Review Letters. 2002;89(6):068102.
Groppe DM, Urbach TP, Kutas M. Mass univariate analysis of event-related brain potentials/fields I: a critical tutorial review. Psychophysiology. 2011;48(12):1711–1725.
Kriegeskorte N, Simmons WK, Bellgowan PSF, Baker CI. Circular analysis in systems neuroscience: the dangers of double dipping. Nature Neuroscience. 2009;12(5):535–540.
Maris E, Oostenveld R. Nonparametric statistical testing of EEG- and MEG-data. Journal of Neuroscience Methods. 2007;164(1):177–190.
Nichols TE, Holmes AP. Nonparametric permutation tests for functional neuroimaging: a primer with examples. Human Brain Mapping. 2002;15(1):1–25.
Sassenhagen J, Draschkow D. Cluster-based permutation tests of MEG/EEG data do not establish significance of effect latency or location. Psychophysiology. 2019;56(6):e13335.
