From Neurons to Networks: A Holistic Review of Electroencephalography (EEG) from Neurophysiological Foundations to Advanced Decoding

1. Introduction

Electroencephalography (EEG) has substantially evolved throughout the years. The incorporation of mathematical tools and computational models has enabled researchers to gain a profound insight into brain function and its applications now encompass a wider array of fields. The study of biomedical signals draws upon multiple disciplines such as biology, biochemistry, neuroscience, medicine, engineering, mathematics and computer science, rendering it as one of the most fascinating areas of research. From a biological and neuroscience point of view, EEG signals elucidate neuronal dynamics that can be proven pivotal in diagnosis and treatment of many neuropsychiatric conditions. From a mathematical and signal processing point of view, they represent a highly complex and information-rich class of signals that can be challenging to interpret. This is because they are intricate by nature, for they are the result of million neurons’ activity which create a non-stationary, nonlinear and temporally dynamic signal [1]. In addition, the useful signal is often submerged in a highly noisy background originating from various sources. For this reason, sophisticated pre-processing techniques are considered imperative to isolate brain activity [2]. Finally, the relationship between the observed electrical potentials in the head and the underlying neurophysiological processes is not direct or simple, as the signal has a complex spatial topography and it is a multi-faceted outcome derived from the interactions of the brain sources [3].

These properties make clinical analysis vulnerable to subjectivity and require advanced computing tools for the reliable and quantitative extraction of diagnostically important information. For these reasons, in recent years, many researchers have been working to advance the methods for processing biomedical signals [4,5,6,7,8]. Basic descriptive statistics provide essential baselines for artifact detection and signal quality assessment, while spectral and time-frequency analysis methods allow for the precise characterisation of oscillatory power dynamics across different brain states [9]. Furthermore, time series analysis techniques [10] and advanced statistical methods [11], are utilised to model the temporal evolution of the signal for forecasting and feature extraction. To map these temporal dynamics to anatomical substrates, spatial analysis and source modeling techniques [12] are used. Connectivity and network analysis employs graph theoretical metrics to quantify the functional integration and information flow between these distributed regions [13]. Additionally, recognising the brain’s complex dynamical nature, nonlinear and chaotic analysis methods are increasingly applied to detect pathological changes in signal complexity [14]. However, the current state of the field is witnessing a paradigmatic shift from these traditional methods toward “Large EEG Models” and Generative AI [15]; while foundation models like LaBraM [16] and Gram [17] leverage massive pre-training on thousands of hours of data to achieve universal cross-subject generalisation and diffusion models [18] enable the direct reconstruction of visual stimuli from neural signals. Consequently, this review aims to describe the fundamental methods and techniques used in EEG analysis, while enabling the researcher to develop a comprehensive understanding of the whole field.

2. Biophysical Principles and Neurophysiological Basis

EEG records the temporal changes in the brain’s electric field, which arise from the sum of extracellular currents associated with the postsynaptic potentials of neurons. The recorded signal is a result not only of endogenous neuronal activity, but also of the electrical properties of the tissues (brain, bone, skin) that function as a potential conductor, as well as the orientation of the sources with respect to the surface of the head. The human brain consists of tens of billions of neurons and glial cells. The basis of electrical activity is the membrane potential $({\mathrm{V}}_{\mathrm{m}}\approx -70 \mathrm{mV})$ [19], which is due to the difference in ion concentration $\left({\mathrm{K}}^{+},{\mathrm{Na}}^{+},{\mathrm{Cl}}^{-},{\mathrm{Ca}}^{2+}\right)$ inside and outside the membrane and it is regulated by active pumps $\left(\mathrm{e}.\mathrm{g}., {\mathrm{Na}}^{+}/ {\mathrm{K}}^{+}-\mathrm{ATPase}\right)$ and passive ion channels. The action potential is a rapid change in this membrane voltage that transmits information along the neuronal terminals.

The Hodgkin-Huxley (HH) model [20,21,22] is a Nobel Prize winning mathematical model that faithfully describes the mechanism of generation and propagation of action potentials in nerve cells. This model was not just a description, but a predictive theory that laid the foundation for modern computational neuroscience. Hodgkin and Huxley’s basic idea was this: The action potential arises from changes in the permeability of the neuron’s membrane to specific ions $\left({\mathrm{Na}}^{+} \mathrm{and} {\mathrm{K}}^{+}\right)$, which are controlled by voltage-gated ion channels. Imagine the neuron’s membrane as a battery with many small switches (the channels). At rest $(\approx -70 \mathrm{mV})$: ${\mathrm{K}}^{+}$ channels are closed, ${\mathrm{Na}}^{+}$channels are tightly closed (but ready to open). When the potential exceeds a threshold $(\approx -55\mathrm{mV})$, the ${\mathrm{Na}}^{+}$ channels open rapidly (depolarisation). Positive ${\mathrm{Na}}^{+}$ ions flow into the cell, further uncoupling the potential $(\approx +40\mathrm{mV})$. This is the anodic phase of the spike. Then, the ${\mathrm{Na}}^{+}$ channels begin to close after a short period of time and the slowly activated ${\mathrm{K}}^{+}$ channels open fully, driving the potential below the normal resting potential, before returning to its original value (refractory period) [23]. This choreography of opening and closing of the channels is what produces the action potential. For an intuitive approach to the Hodgkin-Huxley model, refer to the Appendix A.

The transmission of electrical signals in the brain is explained by the principle of volume conduction [24,25]. According to this, electrical charges interact following the well-known laws of electrostatics: opposite charges attract and like repel. Volume conduction is the process by which a group of ions repels neighboring ions with the same charge. These, in turn, repel subsequent ions, creating a “wave” of charge that moves through the extracellular space. Since the brain is not homogeneous, variations in tissue density can either facilitate or hinder the flow of ions. Furthermore, a signal from a large source (large dipole) can travel a much greater distance than a signal from a small source. To record such a signal from an electrode on the scalp, it must first traverse several layers: from brain tissue, dura mater, skull bone, skin, and hair, until it finally reaches the electrode. At the boundaries between different tissues, where volume conduction ceases, the signal is transmitted through the phenomenon of capacitive coupling. At these points, energy advantages are created, where charges accumulate on one side of an insulating barrier (e.g., bone) and, through electrostatic repulsion, cause the movement of like charges on the opposite side. The succession of these layers (brain, meninges, skull, skin, etc.) creates a series of such capacitors [26]. Figure 1 illustrates multi-scale electrical conduction.

A great review [27] attempts a unified and deep understanding of the biophysical mechanisms that generate extracellular electric fields and currents (e.g., EEG, LFP, ECoG). According to this work, extracellular fields do not originate exclusively from synaptic currents but are the result of the superposition of all active transmembrane currents, including action potentials $\left({\mathrm{Na}}^{+}\mathrm{a}\mathrm{n}\mathrm{d}{\mathrm{Ca}}^{2+}\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{k}\mathrm{e}\mathrm{s}\right)$, intrinsic membrane currents (ionic currents that flow through voltage-gated channels of the neuronal membrane and are not directly induced by synaptic activity), and afterhyperpolarisations (AHPs) that immediately follow the termination of an action potential. In addition to action potentials, calcium spikes, intrinsic membrane oscillations, and neuron-glia interactions play a role in shaping the EEG signal. The overall contribution is critically determined by two factors: (a) the cellular-synaptic architectural design of the tissue (e.g., parallel colonies of pyramidal neurons generate strong “open” fields), and (b) the simultaneous synchronisation of these currents.

The EEG is not directly associated with the action potential, but postsynaptic potentials (EPSPs, IPSPs) that are evoked by neurotransmitters (e.g., glutamate, GABA) in the postsynaptic neurons. EPSPs cause depolarisation from positive ion influx, making the outside relatively more positive, while IPSPs cause hyperpolarisation from negative ion influx or positive ion efflux, making the outside relatively more negative [28]. To detect the signal in the EEG, the spatial and temporal synchronisation of a large population activity, over 10,000–50,000 pyramidal cells [29], is required with a minimum threshold of six square centimeters of synchronised cortex [30]. The main source of the signal comes from the pyramidal neurons of the neocortex, which, due to the parallel arrangement and the large size of their diverging dendrites, create equivalent electric dipoles. The geometry of these dipoles (tangential or radial to the skull surface) determines the amplitude and polarity of the measured signal. A dipole is characterised by a current “sink” (negative charge, usually at the base of the dendrites) and a current “source” (positive charge, close to the particle). Recording critically depends on the orientation of the dipole with respect to the electrode: radial dipoles are easier to detect than tangential ones [31,32].

The extracellular currents resulting from these synaptic events obey the physics of the conductor’s volume. In the quasi-static field approximation, which is valid for EEG frequencies (<100 Hz), we ignore magnetic induction effects and consider the field to be changing so slowly that at any given time it is as if it were static. Under this approach, no significant charge accumulates in the tissues; this is expressed mathematically by Equation (1), where the divergence measures how much a vector field spreads or concentrates at a point and $j$ represents the current density. We also know from Ohm’s law that the current $j$ is proportional to the electric field $E$ (2) where $\sigma$ is the tissue conductivity. Electric field $E$ can be written as the negative slope of the electric potential $\Phi$ (3) where the slope shows the direction and magnitude of the maximum increase in potential and substituting all of this into (1) we get Equation (4). This implies the electric field $E$ is curl-free and current is divergence-free except at active sources.

$$\nabla j=0,$$

(1)

$$J=\sigma E,$$

(2)

$$E=-\nabla \Phi ,$$

(3)

$$\nabla (\sigma \nabla \Phi )=0,$$

(4)

Equation (4) would hold if there were no sources. However, synapses are the sources of current. The density of these “primary” or “active” currents is denoted as $J^p\left(\mathrm{r}\right)$. The deviation of this current, $\nabla J^p\left(\mathrm{r}\right)$, tells us exactly where the sources (positive deviation) and sinks (negative deviation) are. In this way, the final equation that describes how the sources create the potential in a conductor with conductance $\sigma$is Poisson’s Equation (5) [33].

$$\nabla \left(\sigma \left(r\right)\nabla \Phi \left(r\right)\right)=-\nabla J^p\left(r\right),$$

(5)

Instead of modeling each neuron individually, an equivalent current dipole can represent the synchronised activity of a cortical region. This is a mathematical abstraction: a vector $\mathrm{p}$ that has magnitude (proportional to the sum of the currents) and direction (from source to sink). In a simple, infinite, and homogeneous conductor of conductance$\sigma$, the potential $\Phi$ at a distance $\mathrm{r}$ from such a dipole is given by the classical formula (6). Of course, the head is neither infinite nor homogeneous. For more realistic models (e.g., four-sphere head model), the formula becomes more complex [34]. Importantly, only the net (summed) dipolar moment of the active cortical patch matters at the scalp: isolated single-neuron fields cancel, but synchronised populations form a macroscopic dipole that is measurable.

$$\Phi \left(r\right)={\displaystyle \frac{1}{4\mathsf{\pi }\mathsf{\sigma }}}{\displaystyle \frac{\mathrm{p}·\mathrm{r}}{{‖\mathrm{r}‖}^3}},$$

(6)

The Forward Problem [35]: “If the location and strength of the sources (dipoles) inside the brain are known, what will be the electric potential in the scalp?” is a well-defined physical problem. Its solution involves solving Poisson’s equation for a given $J^p$. In computational terms, this is summarised in a linear matrix, the Lead Field Matrix $\left(\mathrm{L}\right)$, that “connects” the sources to the electrodes by Equation (7), where $V$ is the vector of potentials at the electrodes and $I^p$ is the vector of source strengths. The Inverse Problem [36]: “Given the measured potentials $V$ at the electrodes, what are the sources $I^p$ inside the brain that caused them?”. This problem is ill-posed, because infinitely many different source distributions can produce the same pattern in the scalp. To find a unique solution, constraints must be introduced (e.g., the solution must have the minimum possible energy) through a process called regularisation.

$$V=L{I}^p,$$

(7)

3. Recording Technique

EEG technique involves recording the electrical activity of the brain through electrodes placed on the scalp, based on the principle of differential amplification, where the voltages between pairs of electrodes are amplified and displayed as waves on a digital screen. The basic components of an EEG system include the electrodes, the orientation circuits, the analog-to-digital converters (ADCs), and the processing and feedback systems. The process begins with scalp preparation, including cleansing and sometimes light skin peeling to reduce the scalp impedance below 5–10 kΩ for traditional electrode systems to improve conductivity and minimise interference [37]. The placement of electrodes follows the international 10–20 system, which ensures standardised and reproducible placement of the electrodes over brain regions (e.g., frontal, temporal, parietal). The traditional 10–20 electrode system (Figure 2) involves 19 electrode sites and two earlobe-mounted electrodes (A1/A2) that are associated with specific anatomical regions, such that 10–20% of the distance between them is used as the interval for electrode placement. The locations are named with two characters where the first character indicates the brain region (Fp = frontopolar, F = frontal, C = central, P = parietal, O = occipital, T = temporal and A = auricular), and the second character is a number (even = right, odd = left) or “z” for the central electrodes (e.g., Cz) [37,38,39,40,41,42].

EEG technology is characterised by a wide variety of design parameters and applications, which can be classified into multiple dimensions. EEG systems can be categorised based on four main axes (Figure 3): acquisition approach and interface, architecture and electronics, data and computations and context and integration [43]. In terms of signal acquisition, systems range from completely non-invasive scalp surface EEG [44], around-ear EEG for improved portability [45], and subgaleal/epicranial placement [46], to invasive techniques such as electrocorticography (ECoG) with subcutaneous arrays [47,48], stereotactic EEG (sEEG) with approved (depth) electrodes [49], and intracortical microelectrodes for research and neuroprosthetic applications [50]. The device architecture can be wired for high-reliability laboratory applications, wireless (with Bluetooth, Wi-Fi) for mobile and ambiguous studies [51], hybrid [52], wearable (consumer headsets) [53], implantable with telemetry [54], or even modular with structured electrode arrays [55]. Electrode technology covers a wide range, from the contact method: liquid/gel [56], semi-dry/hydrogel [57], dry (e.g., needle, finger) [58,59,60], microneedle electrodes [61], to construction materials such as metal films [62] and polyimide [63], knitted/woven electrodes [64], conducting polymers/graphene [65], with a distinction between active (with built-in preamplifier) and passive electrodes, as well as redeemable versus reusable. The number of channels and montage varies from low density (1–8) for basic brain-computer interfaces (BCI), medium density (16–64) for research, high density (64–256) for source detection, to extremely high density (>256) for special applications, with layouts such as the standard 10–20, 10–10, or custom grids [66,67,68].

Amplifier and signal preprocessing characteristics are critical and include input noise (e.g., <0.5 μV RMS), input impedance (e.g., >1 GΩ), common mode rejection ratio (CMRR), bandwidth and filters, dynamic 1.4, dynamic 2, dynamic bandwidth, filters, and active isolation (driven right leg) [69,70]. Temporal and spectral capabilities are determined by the sampling rate (low: <250 Hz, medium: 250–1000 Hz and high: >1000 Hz serving as approximate guidelines). Real-time latency and synchronisation methods (TTL, PTP, GPS) [71,72] are critical for closed-loop applications where the system must process data and provide feedback within milliseconds. Data management and connectivity involve local storage, data streaming, communication protocols (BLE, Wi-Fi), security, and synchronisation. The software and analysis stack includes real-time preprocessing, AI inference on the device or in the cloud, closed-loop control, and standards compliance (e.g., BIDS) [73]. Multi-modal integration is common, such as EEG + fNIRS, EEG + motion sensors (IMU), EEG + eye tracking/ECG/EGD, and synchronisation with TMS/fMRI/MEG [74,75]. Regulatory compliance and robustness determine whether a system is for research only, a medical device (Class II/III, CE, FDA), MRI compatible, or deployed for industrial use. Ergonomics are tailored to the user (infant, pediatric, adult, animal) and the duration of use. Finally, the power system can be mains-connected, battery-powered, or even energy harvesting systems, with low-power modes for extended studies, while the ecosystem can be closed (vendor-dependent) or open source/hardware [76,77,78].

4. Applications

Contemporary research is witnessing an unprecedented expansion in the application range of EEG (Figure 4) which now sufficiently provide valuable information about a wide spectrum of brain functions employed in a plethora of fields, depending on the purpose [79]. One of the most advanced areas is brain-computer interfaces (BCI), where EEG signals are used to control machines, robots, drones and prosthetic limbs, aiding people with motor disabilities [80]. Moreover, its technology is leveraged to create video games providing both entertaining and therapeutic effects. Beyond rehabilitation, it is used in neuropsychiatric diagnosis (e.g., epilepsy, Alzheimer’s disease) and in neuroscience to measure cognitive load, attention, stress and emotional state [81]. Recent research extends applications to neuroadaptive environments [82] that adjust lighting, sound or content based on the user’s brain activity, digital psychiatry [83] for detecting anxiety and depression through AI-EEG combinations [84] and neuroergonomics [85], where EEG is used in real-world or simulated environments to assess cognitive states like mental fatigue, workload and attention. This can be utilised in a new emerging research field called neuromarketing [86], where the goal is to create more effective marketing campaigns and desirable products. In the field of education, EEG is used in the evaluation of learning effort and concentration, leading to the development of intelligent and adaptive teaching systems [87]. In sports and wellness, it is utilised for neurofeedback and emotion regulation [88]. Importantly, EEG-biometric systems offer secure identification based on brain patterns [89].

5. Basic Characteristics

EEG signals are classified by their frequency into distinct bands, each of which is associated with specific brain functions and states (Figure 5). Delta waves (0.5–4 Hz) are the predominant activity during deep sleep (stage N3) and are normal in infants, while their presence in awake adults may indicate pathology [90,91]. Theta waves (4–8 Hz) are observed during states of deep relaxation and meditation and are normal in children [92,93]. Alpha waves (8–13 Hz) predominate in the occipital region when a person is awake and calm with eyes closed and are associated with a state of calm and inner alertness [94,95]. Beta waves (13–30 Hz) are associated with active, alert thinking, concentration, and problem solving [96,97]. Gamma waves (30–100 Hz) are thought to be involved in higher cognitive functions, such as committing information to memory, perception, and integrating sensory information [98,99]. Finally, EEG frequencies above 100 Hz are not typically categorised into a standard brainwave band, but are considered part of high-frequency activity, often associated with the upper end of the gamma band and beyond. This range can be difficult to measure reliably due to low amplitude, and it can be influenced by factors like anxiety, neurostimulation, and certain neurological conditions [100,101].

In addition, EEGs show a variety of morphological patterns (Figure 5) that reflect different physiological functions or pathological conditions. The Mu (μ) rhythm occurs in the alpha frequency range (8–13 Hz) but has a characteristic arcuate shape and is recorded over the central (Rolandic) areas. In contrast to the alpha rhythm, which is suppressed by opening the eyes, the μ rhythm reacts and is suppressed during actual or even thought of movement, as it is associated with the motor cortex [102]. During sleep, several distinct patterns appear. Sleep spindles are rhythmic activity patterns in the 11–15 Hz range and are one of the defining characteristics of stage N2 (light) sleep, believed to reflect the synchronised activity of thalamo-cortical neuronal networks [103]. The K-complex is a large, high-amplitude biphasic wave that also characterises stage N2 sleep. It is a sudden delta wave, often triggered by external stimuli, and is thought to play a role in maintaining sleep by suppressing arousal [104]. Vertex waves (or V-waves) are short, sharp waves that occur in the central-frontal regions during stage N1 sleep and signal the onset of sleep [105]. In the awake state, lambda waves are transient, triangular waves that occur in the occipital region and are associated with visual exploration or saccadic eye movements, often activated during reading or watching a scene [106]. Finally, spike waves are an important EEG finding that is often associated with epileptic activity [107].

EEG is a complex, spatiotemporal signal that is often treated as a stochastic process. Although the underlying biophysical mechanisms may not be truly random, their high level of intricacy makes statistical description the usual and practical choice for quantitative analysis. EEG signals can be considered quasi-stationary over short temporal windows (on the order of a few to tens of seconds) [108]. This property is fundamental for EEG analysis, as it allows the application of classical signal processing techniques to short segments of the signal, over which it can be considered stationary and permits the estimation of instantaneous measures such as mean, variance, skewness, and kurtosis. These summary statistics and their corresponding distribution tables describe the probability of amplitude values and provide snapshots of the energy distribution; however, their estimates are strongly affected by the possible non-independence of successive samples and by the non-stationarity of the signal, so tests for normality (e.g., Kolmogorov–Smirnov) and corrections for correlation of neighboring samples must be applied with caution [109,110]. Digitisation of EEG requires sampling and smoothing/quantisation. The choice of sampling rate follows the Nyquist theorem stating that to accurately reconstruct an analog signal from its discrete sampling, the sampling frequency (sampling rate) must be at least twice the highest frequency contained in the signal. EEG analysis must be accompanied by appropriate pre-filtering to prevent aliasing, while the quantisation width (e.g., 9–11 bits) affects the signal-to-noise ratio and the fidelity of the representation [111]. Temporal correlation between samples is described by autocorrelation/coherence functions and spectral properties are projected via the power spectral density, which is the Fourier transform of the autocorrelation. Techniques such as periodogram/FFT and smoothing or ensemble-averaging modes are used for reliable spectrum estimates, while spectral moments and parameters such as Hjorth indices summarise morphological features of the spectrum [112]. When the distribution is non-Gaussian, higher-order analysis (e.g., bispectrum/bicoherence) reveals nonlinear correlations and spectral phase coupling between harmonic components [113,114]. Alternatively, interval or period analysis based on zero-crossing machines, measurement of intervals between peaks or half-waves, and interval-amplitude diagrams offers a simple but powerful temporal description of wave structure and is mathematically related to spectral moments and zero-crossing rates (N0, N1, N2), although it is sensitive to high-frequency noise (hysteresis/dead-band is commonly introduced to avoid spurious crossings) and tends to overestimate fast components while underestimating rare long periods [115]. In summary, the comprehensive description of EEG for research combines statistical snapshots (distribution shapes and moments), temporal relationships (auto-/intercorrelation), and spectral representations (PSD, high-order spectra), complementing -where necessary- parametric approaches (e.g., AR models) or pattern detection methods for feature extraction [116].

6. Methods for EEG Processing

6.1. Basic Statistical & Descriptive Analysis

Basic statistical and descriptive methods are the cornerstone of quantitative analysis of EEG signals, as they allow an initial understanding of their shape, variability and statistical behavior. EEG, as a stochastic and nonlinear signal, exhibits fluctuations that reflect the dynamic interaction of large populations of neurons. As mentioned previously, descriptive parameters, such as mean, variance, standard deviation, skewness and kurtosis, are handy tools to get key insights into the statistical properties of the signal [117,118]. At the same time, the investigation of the probability distributions of the EEG provides information on the nature and nonlinearity of the underlying processes and it is the deciding factor in the selection of appropriate statistical tests or models.

The correct choice of sampling frequency, according to the Nyquist theorem, and the use of appropriate anti-aliasing filters are vital for the prevention of distortions as well as the loss of information, while the analysis in bits governs the accuracy of the recording. The autocorrelation and cross-correlation functions can be used to estimate temporal dependencies and delays between channels, providing a first takeaway of synchronisation or functional correlation [119,120]. Finally, methods such as interval or period analysis pave the way for the study of wave periodicity, utilising techniques such as calculating zero-crossings or half-wave intervals to extract simple but insightful temporal features [121,122]. Overall, basic statistical and descriptive methods constitute a necessary first level of analysis, which not only provides an initial overview of the data, but also guides the selection of more complex spectral, parametric or nonlinear approaches.

6.2. Spectral & Time-Frequency Analysis

There are a lot of transformations that reveal the energy structure of the EEG in terms of frequency. The Continuous Time Fourier Transform (CTFT) transforms a time signal into a frequency signal, representing it as a sum of sine and cosine functions [123]. For a continuous signal $x\left(t\right)$, it is defined by Equation (8) where $f$ is the frequency spanning from 0 Hz to half the Nyquist frequency and $j$ the imaginary unit.

$$F\left(f\right)=\underset{-\infty }{\overset{\infty }{\int }}x\left(t\right)e^{-j2\pi ft}dt,$$

(8)

According to Fourier analysis, any physical signal can be decomposed into discrete frequencies or a continuous frequency spectrum. The statistical average of the signal is called EEG spectrum and allows the analysis of the signal in the frequency domain, showing the contribution of each frequency to the signal. EEG signals are discrete due to digitisation since the sampling process converts the original analog signal into a sequence of discrete time instants and quantisation converts the actual voltages into discrete numerical values. Since the input is a sequence of numbers, the mathematically appropriate tool for its analysis in the frequency domain is the Discrete Fourier Transform (DFT) (9) which acts on a discrete-time signal consisting of $N$ samples. The index of the frequency component calculated $\left(k\right)$ corresponds to a specific frequency $f_k=k\times {\displaystyle \frac{F_s}{N}}$, where $F_s$ is the sampling frequency [124].

$$F\left(f\right)=\sum _{n=0}^{N-1}x\left(n\right)e^{{\displaystyle \frac{-j2\pi kn}{N}}},$$

(9)

DFT is implemented via the Fast Fourier Transform (FFT) algorithm which takes advantage of the symmetries and periodicities of the exponential term to break the calculation into smaller, repeated steps, reducing the complexity from $O\left(N^2\right)$ to $O\left(N{\log }_2\left(N\right)\right)$ [125]. The Short-Time Fourier Transform (SFFT) calculation (10) procedure is as follows: 1) The signal is divided into segments and a sliding window is applied ($m$ is the window shift index in time). 2) Each segment is multiplied by a window function $w\left[n\right]$[126]. 3) The DFT (via FFT) is computed on each windowed segment, giving the spectral coefficients for each position. Small windows offer good time but poor frequency resolution, while large windows do the opposite, a fact that expresses the Heisenberg uncertainty limit. Although simple and flexible, the STFT does not automatically adapt to changes in the frequency structure of the signal, and for this reason it is often combined with other techniques (such as multitaper or wavelets) for improved statistical accuracy and spectral leakage reduction.

$$\mathrm{STFT}\left\{x\left(t\right)\right\}\left(m, \omega \right)=\sum _{n=-\infty }^{\infty }x\left[n\right]w\left[n-m\right]e^{-j\omega n},$$

(10)

This transformation allows us to calculate the power spectrum and quantify the energy in the classic brain waves (delta, theta, alpha, beta, gamma), providing fundamental information about brain function. By applying the Fourier transform we can obtain the power spectral density (PSD), which is directly related to functional activity and cognitive states [127]. It has been used as a basic tool for extracting spectral features by detecting frequencies associated with physiological states such as sleep [128], mental load [129] or pathological states such as epilepsy [130]. Nevertheless, the FFT assumes that the signal is stationary within the analysis window. This results in weaknesses in detecting transient phenomena. For example, it cannot detect short signal changes well (such as the characteristic “spikes” in epileptic seizures). It also requires a relatively long signal length to provide reliable spectral estimation, as very short segments have low frequency resolution and cause spectral leakage due to the window boundaries. The FFT is also sensitive to noise.

In contrast to the FFT, which can be used to obtain the amplitude and phase, the Periodogram is a simpler method for a statistical estimation of the spectrum. It is utilised to obtain information only about how the power (energy) of the signal is distributed in the frequency spectrum. It is the basic non-parametric estimation of the PSD and it is calculated from the square of the magnitude of the DFT of the signal, normalised with the sample length (11) [131]. This facilitates the comparison of the power spectrum of different signals, or of the same signal under different conditions, the identification of periodicities as well as the frequency with the maximum contribution to the total power of the signal. The periodogram is statistically inconsistent since the variance of the estimate does not decrease with increasing number of samples, making it sensitive to noise [132]. Because of this variability, smoothing techniques were developed to achieve a more stable estimate, a development that led to the Welch method. The Welch method takes many “statistical samples” (from overlapping segments) and combines them to give a much more reliable picture. The original signal $x\left[n\right]$ of length $N$ is divided into $K$ overlapping segments. Each segment $i$ has length $L$. The $i$-th segment is: $x_i\left[m\right]=x\left[m+i\times D\right], m=0, 1,\dots ,L-1$, where $D$ is the number of offset samples between two consecutive segments (if the overlap is 50%, then $D=L/2$). The total number of segments is $K=\left[{\displaystyle \frac{N-L}{D}}\right]+1$. A window function $w\left[n\right]$is applied to each segment, so the window section is $u_i\left[m\right]=x_i\left[m\right]\times w\left[m\right]$. The periodogram is calculated for each windowed section (12) and the final PSD estimate from the average of the individual periodograms (13) [133].

$$P\left(f\right)={\displaystyle \frac{1}{N}}{\left|\sum _{n=0}^{N-1}x\left(n\right)e^{{\displaystyle \frac{-j2\pi kn}{N}}}\right|}^2,$$

(11)

$$P_{xx}^i\left(f\right)={\displaystyle \frac{1}{L\times U}}{\left|\sum _{n=0}^{L-1}{u_i\left[m\right]\times e}^{{\displaystyle \frac{-j2\pi fm}{L}}}\right|}^2,$$

(12)

$$P_{xx}^{Welch}\left(f\right)={\displaystyle \frac{1}{K}}\sum _{i=0}^{K-1}\left[{\displaystyle \frac{1}{L\times U}}{\left|\sum _{m=0}^{N-1}{u_i\left[m\right]\times e}^{{\displaystyle \frac{-j2\pi fm}{L}}}\right|}^2\right],$$

(13)

Unlike the Welch method, which relies on dividing the signal into overlapping segments and applying a single window to control spectral leakage, the Multitaper method applies multiple optimal mathematical windows (Slepian Tapers), which are designed to maximise the concentration of spectral information within a predefined frequency bandwidth, while minimising leakage outside it [134]. Thus, eliminating the need for subjective window selection and segmentation. This fundamental difference offers a statistically superior estimate of the power spectrum, with minimised variance and minimal leakage at the same time, making it particularly robust for the analysis of short and dynamically changing EEG, where the Welch method can introduce significant bias or miss critical temporal details.

The Wavelet Transform is a fundamental tool in the study of non-stationary signals that overcomes the limitations of the classical Fourier. The mathematical core of the method is based on the concept of the mother wavelet, a function of short duration and zero mean, which resembles a small wave or pulse. The central idea is the following: instead of decomposing the signal into sinusoidal functions of infinite duration, the wavelet transform uses the mother wavelet, which shifts and scales in time, to investigate the signal at different time scales. The significant difference from the STFT is that instead of using a fixed-size window, the window size is adjustable. For high frequencies, the window is narrow, providing high temporal resolution, while for low frequencies, the window is widened, providing high frequency resolution. The Continuous Wavelet Transform (CWT) of a signal $x\left(t\right)$ is defined according to Equation (14) as the inner product of the signal with a shifted and scaled version of the parent wavelet $\psi \left(t\right)$ (e.g., Morlet, Mexican Hat, Daubechies), where $a$ is the scale factor which controls the “stretching” of the wavelet, $b$the translation factor which handles the temporal position of the wavelet along the signal and ${\psi }^{*}$ denotes the complex conjugate.

$$W\left(a,b\right)={\displaystyle \frac{1}{\sqrt{\mid a\mid }}}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}x\left(t\right)\hspace{0.17em}{\psi }^{*}\left({\displaystyle \frac{t-b}{a}}\right)\hspace{0.17em}dt,$$

(14)

CWT calculates the coefficients $W\left(a,b\right)$ for a continuous range of scales and time shifts, creating a two-dimensional power map (scalogram) that illustrates at which times and frequencies EEG activity is increased. The Discrete Wavelet Transform (DWT) is based on the sampling of the scale and shift parameters at discrete values, usually $a=2^j$ and $b=k\times 2^j$, where $j,k\in \mathbb{Z}$. The DWT analyses the signal through binary multiresolution analysis (MRA), separating the signal into two main parts at each level: low frequencies (approximation coefficients), representing the general form of the signal and high frequencies (detail coefficients), representing the details and transient features. The process is implemented through low-pass and high-pass filters, which are applied iteratively to successive decomposition levels [135,136,137]. The inverse process (inverse DWT) allows the reconstruction of the signal without loss of information. Simply put, the CWT acts like a “magnifying lens” that scans each moment in time at all frequencies, while the DWT “breaks” the signal into structural levels, revealing the content of the EEG in an economical but complete way. Wavelet Packet Decomposition (WPD) is a generalisation of the DWT. While in DWT for all subsequent decomposition levels, only the low-frequency approximation component is further split into new approximation and detail components, in WPD, both the low-frequency approximation and the high-frequency detail components are recursively decomposed in subsequent levels. In this way, WPD offers a more uniform and flexible analysis of the entire frequency spectrum, making it extremely useful for EEG signals that contain information at multiple, non-hierarchically organised frequencies [138]. Wavelet Transform has been used in various ways, including feature extraction [139], EEG signal denoising [140], seizure [141] and Event Related Potentials (ERPs) detection [142].

The Synchrosqueezing Transform (SST) is an improved and re-distributed version of the Wavelet Transform (WT), designed to provide better spectral clarity and detection of time-varying frequencies. Initially, CWT is calculated (14) and then the instantaneous frequency is computed (15). Then, the wavelet coefficients are “squeezed”, redistributing the energy from the region $(a,b)$ to the frequency $\omega$, through the kernel (16), where $\delta$ is Dirac delta function [143].

$$\omega \left(a, b\right)=-i{\displaystyle \frac{{\partial }_{t}W\left(a,b\right)}{W\left(a,b\right)}},$$

(15)

$$T\left(\omega ,b\right)=\int W\left(a,b\right)\delta \left(\omega -\omega \left(a,b\right)\right)da,$$

(16)

SST refines the wavelet scalogram moving the energy from areas of ambiguity to the correct frequencies, achieving highly accurate time-frequency representations [144,145]. Stockwell Transform (S-transform) (17) uses a frequency-changing window, like wavelet transform, but in addition, preserves the phase. It uses the Gaussian window function, which depends on the frequency and the time shift parameter $\tau$, which determines the center of the Gaussian window in time [146].

$$S\left(\tau ,f\right)={\int }_{-\infty }^{+\infty }x\left(t\right){\displaystyle \frac{\left|f\right|}{\sqrt{2\pi }}}{e}^{-{\displaystyle \frac{{\left(\tau -t\right)}^2{f}^2}{2}}}{e}^{-i2\pi ft}dt,$$

(17)

Stockwell Transform (S-transform) analyses the signal as “an orchestra of oscillations” arising from the EEG itself, without imposing an external basis. It consists of two stages: (a) Empirical Mode Decomposition (EMD) decomposes the signal $x\left(t\right)$ into a finite number of Intrinsic Mode Functions (IMFs) (18), where each $c_i\left(t\right)$ is an IMF and $r_n\left(t\right)$ is the residual (trend). (b) Then, each IMF is subjected to a Hilbert Transform (19) ($P.V.$ is the Cauchy Principal Value which ensures that the integral converges) and the analytical representation is defined (20) where $a_i\left(t\right)$is the instantaneous amplitude and ${\phi }_i\left(t\right)$ is the instantaneous phase [147,148].

$$x\left(t\right)=\sum _{i=1}^n{c}_i\left(t\right)+r_n\left(t\right),$$

(18)

$${\tilde{c}}_i\left(t\right)={\displaystyle \frac{1}{\pi }}P.V.{\int }_{-\infty }^{+\infty }{\displaystyle \frac{c_i\left(\tau \right)}{t-\tau }}{d}_{\tau },$$

(19)

$$z_i\left(t\right)=c_i\left(t\right)+j{\tilde{c}}_i\left(t\right)=a_i\left(t\right)e^{j{\phi }_i\left(t\right)},$$

(20)

6.3. Time-Series Analysis

EEG signal is a time-series as it captures measurements of the brain’s electrical activity over a consistent interval of time. Thus, time-series methods can be applied. The Prony method is an ancestor of parametric autoregressive models, using exponential fitting to the signal (21), where $A_i$ is the amplitude and $s_i={\alpha }_i+j{\omega }_i$ is the complex frequency [149,150].

$$x_t=\sum _{i=1}^N{A}_i{e}^{s_{i}t},$$

(21)

The method fits a sum of damped exponential functions and is useful for the estimation of dominant frequency components, analysis of damped oscillations and modeling of transient episodes. In modern applications, it has been extended to Matrix Prony and Subspace Prony [151], which utilise algebraic and eigenspace (subspace) techniques to estimate exponential components in multichannel EEG signals, allowing more stable and accurate identification of dominant oscillatory modes and temporally consistent spectral analysis even in highly noisy data.

A signal $x_t$ can be described by an Autoregressive (AR) model as a linear combination of its $p$ previous values (22), where $a_i$ is the autoregressive coefficient and ${\epsilon }_t$ is white noise. Hence the AR model predicts the current value based on the memory of the signal. In a Moving Average (MA) model the signal depends on the current and past errors (23), where $b_i$ corresponds to a noise filter. The MA smooths out the stochastic variability, capturing transient or random fluctuations. Autoregressive Moving Average (ARMA) combines the previous two (24) and offers a balance between memory and random component, which is ideal for EEG that exhibits both internal dynamics and noisy effects. Parameter estimation is executed with methods such as Yule-Walker, Levinson-Durbin, or Burg’s algorithm [152,153,154]. Extending ARMA by introducing an integration term for non-stationarity, we obtain the Autoregressive Integrated Moving Average (ARIMA) (25), where $B$ is the lag operator and $d$ is the degree of integration. Thus, ARIMA(p,d,q) can describe EEG series with long-term variations or trends, while preserving local dynamics [155]. In Autoregressive Moving Average with Exogenous Inputs (ARMAX) models, external variables $u_t$ are introduced that do not belong to the system itself but affect the behavior of the time series, such as sensory stimuli or experimental conditions that affect the EEG during recording (26). In practice, it allows the modeling of EEG with stimuli (e.g., sensory input, task events), combining temporal memory, stochasticity and external influence [156]. Fractional ARIMA (FARIMA or ARFIMA) generalises ARIMA by allowing fractional integration $d\in \mathbb{R}$ (27). This facilitates the capture o long-range dependence or memory effects in EEG. This means that current samples are influenced not only by recent, but also by incredibly old events [157]. Vector Autoregressive (VAR) models extend AR to multivariate time series (28), where ${\mathrm{x}}_t$ is a vector of EEG channels, ${\mathrm{A}}_i$ is the coefficient matrix, and ${\epsilon }_t$ is an error vector. The basic idea is that each channel can influence and be influenced by the others, thus allowing the analysis of interdependencies and causality between brain regions [158]. Multivariate VAR (MVAR) extensions allow the calculation of the Granger causality which estimates whether the activity of one channel can predict the future behavior of another, the Directed Transfer Function (DTF) and the Partial Directed Coherence (PDC) which are based on the spectral representation of the MVAR model and support the estimation of frequency-dependent directionality, while their time-varying extensions capture how these causal connections change dynamically during cognitive processes or states [159,160,161]. Volatility, how the noise or energy of the signal changes over time, can be modeled with the Generalised Autoregressive Conditional Heteroskedasticity (GARCH) model. The central idea of GARCH(p, q) is that the variance of the noise depends on previous values of the variance and the errors. For a time-series $x_t$, ${\sigma }_t^2$is the conditional variance, which evolves according to Equation (29) [162]. Choosing the optimal p, q order of a parametric model is critical to the accuracy and stability of the model. The most common criteria are the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC) [163]. All these methods have been used for instance in dealing with artifacts [164], the prediction of epileptic seizures [165], the diagnosis of autism [166], and the real-time monitoring of the effects of drugs (such as anesthetics) [167].

$$x_t=\sum _{i=1}^p{a}_i{x}_{t-i}+{\epsilon }_t,$$

(22)

$$x_t=\sum _{i=0}^q{b}_i{\epsilon }_{t-i},$$

(23)

$$x_t=\sum _{i=1}^p{a}_i{x}_{t-i}+\sum _{j=0}^q{b}_j{\epsilon }_{t-j},$$

(24)

$${\nabla }^d{x}_t=(1-B{)}^d{x}_t,$$

(25)

$$x_t=\sum _{i=1}^p{a}_i{x}_{t-i}+\sum _{j=0}^q{b}_j{\epsilon }_{t-j}+\sum _{k=1}^r{c}_k{u}_{t-k},$$

(26)

$$(1-B{)}^d{x}_t={\epsilon }_t, \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}(1-B{)}^d={\sum }_{k=0}^{\mathrm{\infty }}{\displaystyle \frac{\mathsf{\Gamma }\left(k-d\right)}{\mathsf{\Gamma }\left(-d\right)\mathsf{\Gamma }\left(k+1\right)}}{B}^k,$$

(27)

$${\mathrm{x}}_t=\sum _{i=1}^p{\mathrm{A}}_i{\mathrm{x}}_{t-i}+{\epsilon }_t,$$

(28)

$${\sigma }_t^2={\alpha }_0+\sum _{i=1}^q{\alpha }_i{\epsilon }_{t-i}^2+\sum _{j=1}^p{\beta }_j{\sigma }_{t-j}^2,$$

(29)

6.4. Advanced Statistical Analysis

The goal of advanced statistical and probabilistic methods is not only to describe the signal, but also to estimate the possible generating processes behind it. Bayesian approach is based on Bayes’ Theorem (30) [168].

$$P\left(\theta ∣D\right)={\displaystyle \frac{P\left(D∣\theta \right)P\left(\theta \right)}{P\left(D\right)}},$$

(30)

The logic is that we combine prior knowledge (priors) with observed data to get updated estimates (posterior estimates). Bayesian approach is instrumental to Bayesian source localisation, which estimates the most likely distribution of brain sources that produce the observed signal, hierarchical bayesian models, which estimates connectivity by integrating information between subjects and uncertainty quantification, which calculates probabilities for the reliability of each parameter [169]. Another useful tool is the Expectation-Maximisation (EM) algorithm which is an iterative optimisation approach used when data contains latent (unobservable) variables. The objective is to maximise the likelihood $P(X\mid \theta )$ based on the available data $X$ and the algorithm alternates between two stages: E-step calculates the expected value of the latent variables (Expectation) and M-step updates the parameters by maximising the expected likelihood (Maximisation) [170,171]. EM is exploited in gaussian mixtures for clustering EEG states as well as in Hidden Markov Models (HMM) and Dynamic Causal Modelling (DCM) for parametric estimation. HMMs assume that there is a latent sequence of states $S_t$ that produces the observed data $X_t$ and the probability of the entire sequence is calculated by Equation (31) [172].

$$P\left(X,S\right)=P\left(S_1\right)\prod _{t}P\left(S_t∣S_{t-1}\right)P\left(X_t∣S_t\right),$$

(31)

This way, the EEG is considered as an alternation of latent “brain states”, which are not directly observed, but leave an imprint on the EEG waves. The estimation of the parameters is usually done with the EM algorithm and these models have been used in the detection of brain states (e.g., sleep, attention, judgment), brain-state decoding and dynamic functional connectivity and event segmentation in cognitive experiments [173]. DCM is a Bayesian state-space model that describes how brain regions interact causally. The system is described by differential Equation (32), where $x$ is the internal brain states, $u$ is the external stimuli and $\theta$ is the connection strengths. The parameters are estimated via variational Bayes (VB), which approximates the posterior distribution of connections [174]. DCM does not simply describe correlations, but how one region causes the response of another. Thus, it is a causal generative model. It is useful in the analysis of directional connectivity (effective connectivity) [175].

$$\dot{x}=f\left(x,u,\theta \right),$$

(32)

6.5. Spatial Analysis & Source Modeling

Spatial analysis and source modeling includes methods for estimating EEG brain sources (inverse problem). Independent Component Analysis (ICA) is based on the model (33), where the observed EEG signal x is a linear mixture of independent sources $s$, through a mixing matrix $A$.

$$x=A\times s,$$

(33)

The method aims to estimate $A^{-1}$, maximising the statistical independence of the components (e.g., through maximising non-Gaussianity or minimising mutual information). This technique is used in EEG denoising (EOG/EMG subtraction), detection of independent rhythms, and functional source analysis (EEG-fMRI integration) [176]. Joint ICA (EEG-fMRI) is applied to unified datasets $X=[X_{EEG},X_{fMRI}]$, where the aim is to find common independent factors that explain both modalities simultaneously. Instead of separating only EEG components, it finds common spatiotemporal factors that explain both the dynamics of the EEG and the spatial distribution of the fMRI. It has been used in combined spatiotemporal interpretation of functional brain activity (neurovascular coupling) [177]. Independent Vector Analysis (IVA) extends ICA to multiple datasets (34) where each $s_k$ is independent of the others.

$$x_k=A_k{s}_k, k=1,\dots ,K,$$

(34)

The goal is to estimate $A_k$ that preserve internal interdependence and external independence. It has been utilized in group-level EEG analysis, multimodal analysis (EEG-MEG) and comparative EEG studies in populations [178,179]. Principal Component Analysis (PCA) finds orthogonal components that maximise the variance of the data (35).

$$X=U\mathsf{\Sigma }{V}^T,$$

(35)

The principal components are the columns of $V$, and the eigenvalues of $X^{T}X$ give the contribution of each component. This method reduces the dimension of the EEG to a few axes that describe the largest variation (energy) of the signal. It has been used in pre-processing for ICA, data compression, EEG power factor analysis and pattern classification [180].

Non-Negative Factorisation (NMF) assumes that the data $X$ can be approximated as a sum of positive bases $W$ that are time activated by $H$ (36). Practically, this means we look for repeating spatial/temporal patterns that add up to form the signal. This is useful for isolating alpha/β patterns and artefacts [181,182]. In Sparse Coding or Dictionary Learning we assume that each instance of the signal is written as a linear combination of some columns of a dictionary $D$, with sparse coefficients $\alpha$ (37) [183,184]. In other words, we have a “dictionary of waveforms” and at any given time a few of them are activated. EMD mentioned before decomposes the signal into intrinsic oscillations (IMFs) plus a residual (38). It peels the signal to isolate the natural, local time-varying frequencies, suitable for transient/nonstationary features. In the inverse source localisation problem, the linear model $x = Lj + n$ is used [185]. Beamformers construct spatial filters $w$ that maximise the gain at the location of interest and minimise the power elsewhere, the formulation and closed-form solution are given in (39) [186,187]. The Minimum-Norm Estimate (MNE) favours the solution with the lowest total energy and is written as regularised least-squares (40), so when several scenarios explain the data, we choose the most parsimonious explanation [188]. Variants such as LORETA/sLORETA/eLORETA introduce spatial smoothness terms (matrix $R$) in the normalised problem for smoother and more reliable source allocations (41) [189]. Dipole fitting is expressed as a nonlinear optimisation for fitting a few dipoles $(\theta ,q)$ that explain the signal and formulated in (42) [190,191]. These methodologies are chosen or combined depending on the objective (e.g., beamformer for targeted spatial isolation, MNE/LORETA for distributed evoked maps, dipole fitting for focal sources), while decomposition techniques (NMF, dictionary learning, EMD) are often used as pre-processing or to extract features.

$$\underset{W,H\ge 0}{\min }\parallel X-WH{\parallel }_F^2,$$

(36)

$$\underset{D,A}{\min }\parallel X-\mathit{DA}{\parallel }_F^2+\lambda \parallel A{\parallel }_1,\hspace{0.25em} \mathrm{s}.\mathrm{t}. \parallel d_k{\parallel }_2 \le 1\forall k,$$

(37)

$$x\left(t\right)=\sum _{k=1}^K{\mathrm{I}\mathrm{M}\mathrm{F}}_k\left(t\right) +\hspace{0.25em} r\left(t\right),$$

(38)

$$\underset{w}{\min }{w}^{\mathrm{\top }}Rw, \mathrm{s}.\mathrm{t}.\hspace{0.25em} w^{\mathrm{\top }}l=1\Rightarrow w={\displaystyle \frac{R^{-1}l}{l^{\mathrm{\top }}{R}^{-1}l}},$$

(39)

$$\widehat{j}=\mathrm{a}\mathrm{r}\mathrm{g}\underset{j}{\min }\parallel x-Lj{\parallel }_2^2+\lambda \parallel j{\parallel }_2^2\Rightarrow \widehat{j}=(L^{\mathrm{\top }}L+\lambda I{)}^{-1}{L}^{\mathrm{\top }}x,$$

(40)

$$\widehat{j}=(L^{\mathrm{\top }}L+\lambda R{)}^{-1}{L}^{\mathrm{\top }}x,$$

(41)

$$\underset{\theta ,q}{\min }\parallel x-L\left(\theta \right)q{\parallel }_2^2,$$

(42)

6.6. Connectivity & Network Analysis

Another domain of EEG analysis is the study of functional and causal relationships between brain regions based on temporal, spectral and probabilistic correlations. Functional connectivity measures statistical correlation or coherence between pairs of signals: Pearson satisfies the need for linear correlation between two time series (43), while Coherence expresses the frequency-dependent linear relationship through spectral densities (44) [192]. In practice, high $r$ or coherence values in a specific band indicate possible involvement in common cognitive processes and are used in studies of hemispheric cooperation, sleep, and pathologies such as epilepsy [193]. Causal connectivity is based on MVAR models that predict the present from the past; in Granger’s formulation, if the predictive power of one signal is improved by adding the past terms of another, we consider that Granger causality exists (45). Spectral extensions of this framework give the frequency transfer function $H\left(f\right)$ (46) and directivity measures such as the DTF (47) and PDC (48) mentioned before, which represent the flow of information per frequency; these are useful for identifying information flow in sensory networks and pathological conditions [194].

$$r_{xy}={\displaystyle \frac{\mathbb{E}\left[\left(x-{\mu }_x\right)\left(y-{\mu }_y\right)\right]}{{\sigma }_x{\sigma }_y}},$$

(43)

$$C_{xy}\left(f\right)={\displaystyle \frac{\mid S_{xy}\left(f\right){\mid }^2}{S_{xx}\left(f\right)\hspace{0.17em}{S}_{yy}\left(f\right)}},$$

(44)

$$x_t=\sum _{k=1}^p{A}_{xx}\left(k\right)\hspace{0.17em}{x}_{t-k}+\sum _{k=1}^p{A}_{xy}\left(k\right)\hspace{0.17em}{y}_{t-k}+e_t,$$

(45)

$$H(f)=(I-\sum _{k=1}^{p}A\left(k\right)e^{-i2\pi fk}{)}^{-1},$$

(46)

$${\mathrm{D}\mathrm{T}\mathrm{F}}_{ij}\left(f\right)={\displaystyle \frac{\mid H_{ij}\left(f\right){\mid }^2}{\sum _k\mid H_{ik}\left(f\right){\mid }^2}},$$

(47)

$${\mathrm{P}\mathrm{D}\mathrm{C}}_{ij}\left(f\right)={\displaystyle \frac{A_{ij}\left(f\right)}{\sqrt{\sum _k\mid A_{kj}\left(f\right){\mid }^2}}},$$

(48)

Nonlinear connectivity analysis includes metrics such as Transfer Entropy (TE) [195], which measures the directional nonlinear information from Y to X (49), and Mutual Information (MI) [196], which measures the general mutual dependence between two variables (50); the intuition here is that these metrics detect relationships that are not constrained by linearity or Gaussian assumptions and are therefore applicable to complex cognitive processes, judgments, and consciousness studies.

$${\mathrm{T}\mathrm{E}}_{Y\to X}=\sum p\left(x_{t+1},x_t,y_t\right)\log {\displaystyle \frac{p\left(x_{t+1}∣x_t,y_t\right)}{p\left(x_{t+1}∣x_t\right)}},$$

(49)

$$\mathrm{M}\mathrm{I}\left(X,Y\right)=\sum _{x,y}p\left(x,y\right)\log {\displaystyle \frac{p\left(x,y\right)}{p\left(x\right)p\left(y\right)}},$$

(50)

Phase methods focus on phase synchronisation: Phase-Locking Value (PLV) measures the stability of the phase difference (51), Phase-Lag Index (PLI) excludes zero phase lags to reduce the effect of volume conduction (52), while weighted-PLI (53) and Amplitude Envelope Correlation (AEC) (54) provide variants that emphasise reliable, non-zero phase contributions or amplitude envelope correlation; these measures are widely used to study phase locking in attention, perception, neurofeedback, and epilepsy. Cross-Frequency Coupling (CFC), and especially Phase-Amplitude Coupling (PAC), examines whether the phase of a low frequency modulates the amplitude of a higher one; this can be measured simply as phase-amplitude correlation or with the Modulation Index (Tort) based on Kullback-Leibler divergence (55), and serves to highlight hierarchical neural interactions (e.g., memory, attention) [197,198,199].

$$\mathrm{P}\mathrm{L}\mathrm{V}=\left|{\displaystyle \frac{1}{N}}\sum _{n=1}^N{e}^{i\left({\varphi }_x\left(n\right)-{\varphi }_y\left(n\right)\right)}\right|,$$

(51)

$$\mathrm{P}\mathrm{L}\mathrm{I}=\left|\mathrm{sign}\left(\sin \left({\varphi }_x-{\varphi }_y\right)\right)\right|,$$

(52)

$$\mathrm{w}\mathrm{P}\mathrm{L}\mathrm{I}={\displaystyle \frac{\left|\mathbb{E}\left[\mathrm{Im}\left\{S_{xy}\right\}\right]\right|}{\mathbb{E}\left[\left|\mathrm{Im}\left\{S_{xy}\right\}\right|\right]}},$$

(53)

$$\mathrm{A}\mathrm{E}\mathrm{C}=\mathrm{corr}\left({\mathrm{e}\mathrm{n}\mathrm{v}}_x\left(t\right),{\mathrm{e}\mathrm{n}\mathrm{v}}_y\left(t\right)\right),$$

(54)

$$\mathrm{P}\mathrm{A}\mathrm{C}\approx \mathrm{corr}\left({\mathrm{p}\mathrm{h}\mathrm{a}\mathrm{s}\mathrm{e}}_{\mathrm{low}}\left(t\right),{\mathrm{a}\mathrm{m}\mathrm{p}}_{\mathrm{high}}\left(t\right)\right),\mathrm{M}\mathrm{I}={\displaystyle \frac{D_{KL}\left(P\parallel U\right)}{\log N}},$$

(55)

Symbolising connectivity as a graph G(V,E) allows the calculation of indicator graphs such as modularity $Q$ (56), global efficiency (57), clustering coefficient (58), and centrality metrics, which quantify the organisation, efficiency, and roles of nodes in the network. Such indices are applied in resting-state studies and in diseases such as Alzheimer’s [200]. Since connectivity varies over time, time-varying approaches are applied: sliding windows for time-local estimates or Kalman-based/recursive MVAR models [201] where the parameters $A_t$ change over time and are estimated using Bayesian or recursive filters (59); these are crucial for studying dynamic FC, microstates, and task-related changes. DCM model mentioned previously approaches causality with biophysically parameterised state-space models (60) and Bayesian inversion to estimate the effect of experimental factors on connections; DCM allows comparison of causal structure hypotheses and is widely used in SPM for EEG/MEG [202]. Finally, the Phase Synchronisation Index (PSI) generally measures the stability of the phase difference (61) and is used to evaluate synchronisation in sensory, motor, and cognitive processes [203]. In all of the above approaches, the choice of the appropriate metric depends on the background of the hypothesis (linear vs. nonlinear, statistical vs. causal), sensitivity to volume conduction, desired spatial/temporal resolution, and the semantics of the study; practical factors such as spectral density estimation, noise model, normalisation, and correct estimation of hyperparameters are critical for reliable conclusions.

$$Q={\displaystyle \frac{1}{2m}}\sum _{i,j}\left(A_{ij}-{\displaystyle \frac{k_i{k}_j}{2m}}\right)\delta \left(c_i,c_j\right),$$

(56)

$$E_{\mathrm{g}\mathrm{l}\mathrm{o}\mathrm{b}}={\displaystyle \frac{1}{N\left(N-1\right)}}\sum _{i\ne j}{\displaystyle \frac{1}{d_{ij}}},$$

(57)

$$C_i={\displaystyle \frac{2T_i}{k_i\left(k_i-1\right)}},$$

(58)

$$x_t=A_t{x}_{t-1}+w_t,$$

(59)

$$\dot{x}=Ax+Bu+C,$$

(60)

$$\mathrm{P}\mathrm{S}\mathrm{I}={\displaystyle \frac{1}{N}}\sum _{n=1}^N{e}^{i\left({\varphi }_x\left(n\right)-{\varphi }_y\left(n\right)\right)},$$

(61)

6.7. Nonlinear & Chaotic Analysis

EEG activity is an extremely complex and nonlinear signal that reflects the collective action of large neural populations. Brain dynamics are often characterised by nonlinear interactions, transient states, and chaotic behavior, which makes the use of purely linear models inadequate for their description [204]. For this reason, a wide range of nonlinear and chaotic analysis methods have been developed, which allow the investigation of the stability, complexity, and predictability of EEG time-series, offering a deeper understanding of the underlying neural dynamics. One of the fundamental measures in chaos theory is the Lyapunov exponents, which describe the average rate of divergence of two neighboring trajectories in phase space. If two initial states differ by $dx\left(0\right)$, then the Lyapunov exponent is defined as the logarithm of the ratio of the final to the initial distance, divided by time (62):

$$\lambda =\underset{t\to \infty }{\lim }{\displaystyle \frac{1}{t}}\ln {\displaystyle \frac{\mid \delta x\left(t\right)\mid }{\mid \delta x\left(0\right)\mid }},$$

(62)

Positive values of λ indicate sensitivity to initial conditions and therefore chaotic behavior. In EEG analysis, the presence of positive Lyapunov exponents has been linked to unstable and chaotic neural states, such as in epileptic seizures, while negative or zero values correspond to stable or periodic states, such as sleep or anesthesia [205,206]. Fractal dimension (D) is another fundamental measure of complexity, as it expresses how the detail of a signal changes with the scale of observation. Using the box-counting approach, this dimension is defined as the limit of the ratio of the logarithms of the number of squares $N\left(ϵ\right)$ required to cover the signal trajectory and their inverse size (63).

$$D=\underset{ϵ\to 0}{\lim }{\displaystyle \frac{\log N\left(ϵ\right)}{\log \left(\frac{1}{ϵ}\right)}},$$

(63)

High $D$values correspond to more complex EEG signals, which have been observed during cognitive processes or in pathological conditions such as epilepsy. To investigate the dynamic structure of the EEG, reconstruction of the attractor in phase space is often used, according to Takens’ theorem. Each point of the signal can be represented as a state vector $\mathrm{x}\left(t\right)$ through the delay τ and the embedding dimension m, as shown in Equation (64). The geometry of the reconstructed tractor can be described by the correlation dimension $(D₂)$, defined by the correlation function $C\left(r\right)$ in Equation (65).

$$\mathrm{x}\left(t\right)=\left[x\left(t\right),x\left(t+\tau \right),\dots ,x\left(t+\left(m-1\right)\tau \right)\right],$$

(64)

$$C\left(r\right)={\displaystyle \frac{2}{N\left(N-1\right)}}\sum _{i<j}H\left(r-‖x_i-x_j‖\right),D_2={\displaystyle \frac{\mathrm{dln}C\left(r\right)}{\mathrm{dln}r}},$$

(65)

The existence of a low-dimensional attractor $(D₂ < 5)$ is indicative of chaotic but limited dynamics, which has been observed in epileptic seizures and cognitive transitions [207,208,209]. Approximate entropy (ApEn) introduces a statistical measure of signal regularity. It is defined as the difference between the probability that similar patterns of length $m$ remain similar when extended by one dimension, as given in Equation (66). Low ApEn values indicate predictable signals, while high values indicate greater complexity. Sample Entropy (SampEn) improves ApEn by excluding self-correlations and is defined as the negative logarithm of the ratio of probabilities A and B of two similar patterns of length $m+1$ and$m$ (67). SampEn is less biased and more reliable in short EEG sequences, which is why it is used in studies of sleep, consciousness, and cognitive load. Multiscale Sample Entropy (MSE) extends SampEn by calculating entropy at different time levels through “coarse-graining” of the signal. Each new sequence is obtained as the average of the values within windows of length $\tau$ according to Equation (68). In this way, MSE provides a profile of complexity at multiple time scales, while the Refined Composite MSE (RCMSE) version offers improved stability for short EEG time-series [210]. Kolmogorov entropy (K) links chaos theory with information theory, expressing the rate at which new information is generated in a dynamic system. It is defined as the sum of all positive Lyapunov exponents, as shown in Equation (69). High $K$values correspond to increased chaos and rapid loss of predictability, characteristics that have been observed in EEG states of increased neural excitation. Permutation entropy (PE) offers a quick and efficient assessment of EEG complexity. It converts the time series into permutation sequences and calculates the Shannon entropy of their probabilities, as given in Equation (70). This measure is particularly robust to noise and has been used extensively in applications such as seizure detection and anesthesia monitoring. TE (49) mentioned before measures the nonlinear, directed flow of information between two time series and is defined as the conditional information difference between probability distributions. Unlike Granger’s Linear Causality, TE does not assume linearity or normality and has been widely used to study the effective connectivity of the brain [211,212].

$$ApEn\left(m,r,N\right)={\varphi }_m\left(r\right)-{\varphi }_{m+1}\left(r\right),$$

(66)

$$SampEn\left(m,r,N\right)=-\ln {\displaystyle \frac{A}{B}},$$

(67)

$$y_j^{\left(\tau \right)}={\displaystyle \frac{1}{\tau }}\sum _{i=\left(j-1\right)\tau +1}^{j\tau }{x}_i,$$

(68)

$$K=\sum _{{\lambda }_i>0}{\lambda }_i,$$

(69)

$$PE=-\sum _{i=1}^{n!}{p}_i\ln p_i,$$

(70)

Another fundamental measure is the Hurst exponent (H), which quantifies the long-term dependence and self-similarity of the signal. The relationship between the $R/S$ ratio and the sample size $N$ follows Equation (71). $H$values greater than 0.5 indicate persistent behavior, while values less than 0.5 correspond to mean-reverting processes. The measurement of H has been applied in the analysis of fatigue, aging, and cognitive levels [213,214]. Finally, Recurrence Quantification Analysis (RQA) is based on the creation of a recurrence plot, i.e., a table that records when the system’s trajectory returns to previous states. This table is defined as (72) where $\Theta$ is the Heaviside function and $ϵ$ is the proximity radius. This table yields indices such as the repetition rate, determinism, and stationarity, which reveal transitions between different brain states [215,216,217]. Overall, the above nonlinear and chaotic methods provide a multidimensional approach to the study of brain activity.

$${\displaystyle \frac{R}{S}}\propto N^H,$$

(71)

$$R_{ij}=\mathsf{\Theta }\left(ϵ-‖x_i-x_j‖\right),$$

(72)

6.8. Machine Learning & Deep Learning

The analysis of EEG signals has been reshaped in the last half-decade. Τhe manual feature engineering that defined the late 20th century paved the way for sophisticated data-driven approaches, more suitable for the dynamic human brain signals [218]. A plethora of ML algorithms have been used, and the choice of classifier depends heavily on the dataset size and the requirement for interpretability. Support Vector Machines (SVM) are the most robust classifier for “small data” regimes, which characterises many BCI pilot studies (where N < 20 subjects). SVMs map input vectors into a high-dimensional space using a kernel function (e.g., Radial Basis Function) to find the optimal hyperplane separating two classes (e.g., “Left Hand Movement” vs. “Right Hand Movement”). Linear Discriminant Analysis (LDA) is the standard for Brain-Computer Interfaces (BCI), particularly in P300 speller paradigms. It assumes that each class follows a Gaussian distribution with a shared covariance matrix. It is less powerful than nonlinear classifiers but its simplicity allows for adaptive implementations, where the classifier updates in real-time as the user’s brain signals drift during a session. For offline analysis where computational speed is less critical, ensemble methods like Random Forest and Gradient Boosting (XGBoost, LightGBM) have seen increased use [219,220,221]. These methods aggregate the decisions of hundreds of decision trees, handling nonlinear relationships and interacting features better than SVMs. They also provide future ranking, offering a degree of intrinsic interpretability that facilitates the identification of the frequency bands or channels that are driving the classification.

While conventional ML operates in Euclidean space, a significant theoretical advancement in the 2020–2025 period has been the widespread adoption of Riemannian Geometry [222]. This approach represents a “middle ground” between classical ML and Deep Learning (DL) as it is mathematically rigorous, highly accurate, but does not require the massive parameters of a neural network. The core insight of this framework is that the spatial covariance matrix of the EEG signal contains the essential information regarding brain state. Covariance matrices are Symmetric Positive Definite (SPD). Geometrically, SPD matrices do not lie on a flat plane (Euclidean space) but form a curved manifold (a convex cone). Applying Euclidean operations (like standard averaging) to these matrices introduces distortions, known as the “swelling effect”. Manifold-Based Classification Riemannian methods utilise the Affine Invariant Riemannian Metric (AIRM) to calculate distances on this manifold. The geodesic distance (the shortest path along the curve) between two covariance matrices supply with a much more robust measure of similarity than the Euclidean distance [223,224]. The most successful implementation of this is the Tangent Space Mapping (TSM) approach. Covariance matrices on the manifold are projected onto a tangent space which is a flat Euclidean vector space tangent to the geometric mean of the dataset. Once in the tangent space, the matrices are treated as vectors and then standard Euclidean classifiers (like SVM or LDA) are applied to these tangent vectors. This “Riemannian-Geometry-based” pipeline currently holds state-of-the-art accuracy for Motor Imagery BCI tasks on small datasets, outperforming complex Convolutional Neural Networks when training data is limited. Techniques like Riemannian Procrustes Analysis (RPA) allow researchers to geometrically align the data manifolds of different subjects [225,226,227]. By translating, rotating, and scaling the covariance matrices of a new user to match the “average user” manifold, models can be trained on a database of existing subjects and applied to a new subject with minimal calibration. This capability is critical for the development of “plug-and-play” BCI systems.

The fundamental advantage of DL is the end-to-end learning as the model learns the feature extraction steps itself, without the need for manual selection of frequency bands or entropy measures. Convolutional Neural Networks (CNNs) are the most reliable tool of modern EEG analysis, utilised in over 90% of DL studies until 2023 [218]. However, the architecture used for EEG differs significantly from the 2D CNNs used in computer vision because EEG data has two dimensions, time, and channels. Since these dimensions have different physical units and correlation structures, standard square filters are suboptimal. EEGNet [228] is a lightweight architecture that uses depthwise separable convolutions. It first applies a temporal convolution (acting as a frequency filter) and then a spatial convolution (acting as a spatial filter) separately. DeepConvNet [229] and ShallowConvNet [230] use varying depths to capture different levels of abstraction. ShallowConvNet is designed to extract logarithmic band power features, while DeepConvNet can learn more complex, hierarchical representations. Multi-Branch CNNs [231] process the raw signal in parallel streams with different kernel sizes (e.g., small kernels for gamma waves, large kernels for delta waves). The branches are then fused, allowing the model to capture multi-scale temporal dynamics simultaneously. This approach has shown success in emotion recognition and sleep staging. While CNNs excel at extracting local features, they struggle with long-term dependencies (e.g., a sleep cycle lasting 90 minutes). To address this, Long Short-Term Memory (LSTM) networks and Gated Recurrent Units (GRU) are employed. The standard architecture is the C-RNN (Convolutional-Recurrent Neural Network) [232,233,234,235]. In this hybrid model a CNN front-end extracts spatial-spectral features from short EEG chunks (e.g., 1 second). These feature vectors are fed as a sequence into an LSTM back-end. The LSTM captures the temporal evolution of the brain over longer periods (e.g., 30 seconds). The introduction of the Attention Mechanism has refined EEG models. Attention allows the network to weigh the importance of different input parts dynamically. For instance, channel attention can be utilised so the model learns which electrodes are most relevant for the current task and suppresses noise from irrelevant channels. Temporal attention makes the model focus on the specific time segments containing the event of interest while ignoring the background EEG activity. Now, attention blocks are routinely inserted into CNN architectures (e.g., SE-Net blocks) to improve performance and interpretability [236].

One of the most theoretically grounded advancements in the last years is the application of Graph Neural Networks (GNNs) [237]. GNNs treat the electrodes as nodes in a graph, connected by edges representing their relationship, taking advantage of the actual 3D geometry of the head. In static graphs, edges are defined by physical Euclidean distance (e.g., Fp1 is connected to Fp2 but not O1). In functional graphs, edges are defined by functional connectivity measures (e.g., Phase Locking Value or Coherence). If two regions oscillate in sync, they are connected, regardless of physical distance. The most advanced GNNs employ Dynamic Graph Learning. Since functional connectivity in the brain is not static but evolves on a millisecond timescale, these models learn the adjacency matrix from the data itself at each time step. The network dynamically strengthens or weakens connections between nodes, effectively learning the evolving functional network of the brain during a cognitive task. GNNs have demonstrated superior performance in Emotion Recognition, a task heavily dependent on the functional interplay between frontal and temporal lobes. By explicitly modeling these interactions, GNNs capture the global network patterns of emotion that local CNN filters miss [238,239,240,241].

The success of Large Language Models (LLMs) like GPT-4 has inspired a wave of “EEG Foundation Models”. These models leverage the Transformer architecture [242,243] to learn universal representations of brain activity from massive, diverse datasets. The continuous signal is sliced into patches (e.g., 200 ms windows). These patches are projected into embedding vectors. The model calculates the relationship between every pair of tokens. This allows it to capture global dependencies across the entire recording duration, solving the “receptive field” limitation of CNNs. LaBraM [16] addresses the heterogeneity of EEG data (different datasets use different electrode caps). It segments signals into “EEG channel patches” and uses a Neural Tokenizer trained via Vector-Quantised Neural Spectrum Prediction. It is pre-trained on approximately 2,500 hours of data using a Masked Modeling objective (like BERT). Random patches of the EEG are masked, and the model must reconstruct the missing neural codes. This forces the model to learn the grammar of neural dynamics. Similarly, BrainBERT [244] and NeuroBERT [245] utilise masked auto-encoding on intracranial and scalp EEG, respectively. These foundation models aim to create a “universal embedding” space for brain signals. A key capability of these models is Cross-Subject Generalisation. Because they are trained on thousands of individuals, they learn features that are invariant to the anatomical differences between subjects, addressing the “calibration problem” that has plagued BCI for decades.

One rapidly evolving area is the application of Generative Diffusion Models to EEG [246]. These models, which power image generators, are being adapted for two distinct purposes in neuroscience: Data Augmentation and Visual Decoding. The scarcity of high-quality labeled EEG data (especially for rare pathologies) is a major bottleneck. Denoising Diffusion Probabilistic Models (DDPMs) are now used to generate synthetic EEG signals. The model learns the data distribution by reversing a gradual noise addition process. Starting from pure Gaussian noise, it iteratively “denoises” the signal to produce a realistic EEG waveform. Recent studies confirm that this synthetic data preserves the spectral and temporal characteristics of real EEG (e.g., Visual Evoked Potentials) and can be used to train classifiers, improving accuracy [247,248,249]. DreamDiffusion [250] is a framework that uses a Stable Diffusion backbone but replaces the text encoder with a pre-trained EEG encoder. It generates images directly from EEG signals recorded while a user looks at a picture. NeuroDM [251] uses an EEG-Visual-Transformer (EV-Transformer) to extract visual-semantic features from the EEG. These features then guide a diffusion model (EEG-Guided Diffusion) to reconstruct the image. NeuroDM has achieved state-of-the-art results, generating images that not only match the category (e.g., “dog”) but also the structural composition of the stimulus seen by the user. Beyond generating images, diffusion is used for feature learning. EEGDM [252] employs a Structured State-Space Model (SSM) within a diffusion framework. By training the model to denoise EEG signals, the SSM captures the complex temporal dynamics of the brain. The latent representations learned during this process are then fused and used for highly accurate classification of seizures, outperforming standard Transformer baselines.

Self-Supervised and Contrastive Learning are emerging paradigms used to further reduce reliance on labeled data [253,254]. Contrastive learning trains a model to recognise that two augmented versions of the same EEG segment are “similar” (positive pair) while distinguishing them from other segments (negative pairs). The challenge in EEG is defining “augmentation” because unlike images, it cannot be rotated or cropped without destroying information. Signal-Transformation Contrastive Learning, applies physiologically plausible augmentations like permutation of symmetric channels, time-warping, or adding band-limited noise. Frameworks like MoCo (Momentum Contrast) [255] allow for massive queues of negative samples, enabling the model to learn highly discriminative features without a single label. Contrastive Predictive Coding (CPC) trains the model to predict the future latent representation of the signal based on the past. This forces the model to encode information that is slow-varying and predictive (like the underlying brain state) while discarding fast-varying noise [256].

For clinical adoption, a doctor must know for example why the AI diagnosed epilepsy. Explainable AI (XAI) has thus become a critical research area. Saliency and Relevance Mapping Methods like SHAP (SHapley Additive exPlanations) and LRP (Layer-wise Relevance Propagation) are now routinely applied to EEG models [257,258]. They generate heatmaps showing which time points and electrodes contribute to the prediction. For example, in a depression detection model, XAI might reveal that the model is focusing on Alpha asymmetry in the frontal electrodes, confirming a known biomarker. SHERPA (SHAP-based ERP Analysis) [259] uses XAI not just for validation but for discovery. By calculating the “importance score” of every point in an ERP, SHERPA has identified subtle cognitive processes (like negative selection mechanisms) that standard statistical analysis missed. This marks a shift where AI is used to advance neuroscience theory, not just automate diagnosis. Figure 6 summarises the diverse methods applied to EEG data.

7. Future Trajectory and Evolution

EEG analysis is fast advancing towards the incorporation of integrated approaches that combine multimodal data and real-time edge computing, improve the interpretability of algorithms, and address ethical issues. These efforts aim to bridge the gap between EEG’s excellent temporal resolution and its limited spatial information and move neuroscientific discoveries from the lab into practical, everyday use. The goal is to guarantee responsible, equitable, and safe use of neurotechnologies. Integrating other imaging modalities into EEG analysis combines the millisecond neuronal fluctuations recorded by EEG with the high spatial precision of fMRI and the complementary sensitivity of MEG. Successful fusion requires careful safety and artifact management. The integration of other biomarkers (fNIRS, EMG, ECG, GSR, eye-tracking, accelerometers) in polymorphic frameworks provides complementary information for neurophysiological and psychophysiological states [260,261,262,263]. The transition to real-time applications at the edge faces latency, privacy, and energy consumption challenges. Classic algorithms with a small footprint (LDA, linear SVM) remain useful, while architectures such as EEGNet show that deep models can be sufficiently compressed (pruning, quantisation) to run on ARM processors via TensorFlow Lite Micro or ONNX Runtime Mobile [264]. Design practices for the edge include stream-oriented buffer management with sliding windows, incremental feature updates to avoid re-computation, event-driven subroutine activation to save energy, and dynamic online model adaptation with trust mechanisms to prevent catastrophic forgetting. Local processing limits the exposure of raw neural data, while encrypted channels and limited transmission of fused features ensure privacy. Interpretability is a key prerequisite for clinical trust and responsible implementation. Strategies include integrating attention mechanisms into networks to visualise temporal/spatial “zones” of importance, using Grad-CAM on spectrograms, Layer-wise Relevance Propagation (LRP) to decompose contributions down to the millisecond/channel level, and self-attention analysis in transformer models to uncover long-term dependencies [265]. Also, methods such as SHAP and fidelity tests (retraining on top-features) assess the reliability of explanations. Finally, the ethical and social consequences require careful policy: strong privacy measures (local preprocessing, encryption, retention limitation), dynamic and informed consent, algorithm evaluation and mitigating bias due to disparities in signal quality. In applications of job control or educational assessment, rules are needed that prevent coercive use or stigmatisation. It is evident that data ownership and sharing agreements need to be carefully controlled to strike a balance between the benefit of open science and the rights of participants [266]. The field of EEG analysis is poised to transition from task-specific deep learning architectures to Universal Brain Foundation Models. Future methodologies will likely center on calibration-free systems, were massive pre-training on thousands of hours of diverse neural data, exemplified by Large EEG Models (LEMs), effectively eliminating the historic bottleneck of subject-specific training sessions. Additionally, the integration of Generative Diffusion Models is expected to evolve from simple data augmentation to sophisticated Neural Decoding, enabling the direct, high-fidelity reconstruction of visual and semantic content from brain activity. This paradigm shift will be anchored by advanced XAI frameworks, which are critical for decoding the “black box” of these foundation models to discover novel neurophysiological biomarkers, thereby bridging the gap between computational power and clinical trust.

Overall, the future trajectory of EEG analysis will shift towards more advanced, all-in-one solutions that are faster, easier to understand, and ethically appropriate. Combined with better sensors and processing devices, this will enable a responsible transition from the laboratory to a wide range of applications.