How we build clinical-grade EEG processing from the ground up, and why it matters for early neurological detection.
Most EEG research happens in MATLAB with EEGLAB, or in Python with MNE. These are fine tools for offline analysis in a lab setting. They are not suitable for building clinical products.
EEGLAB is fundamentally a MATLAB GUI application. Its processing functions are designed for interactive use: load a recording, click through menus, visually inspect results, adjust parameters, re-run. This workflow makes sense when a researcher is exploring data. It falls apart when you need deterministic, automated processing of continuous EEG streams in a clinical environment. EEGLAB also requires a MATLAB license, which adds cost and deployment constraints that are unnecessary for a production system.
MNE-Python is more programmatic, but it carries its own baggage. The internal pipeline has implicit assumptions about data structure and processing order that are poorly documented. Memory management is delegated to NumPy and Python's garbage collector, which is acceptable for batch processing but problematic for real-time applications where latency spikes from GC pauses are unacceptable. The abstractions are also heavy: MNE's Raw and Epochs objects carry significant overhead per operation.
We write our signal processing in C and C++. Not because it is fashionable, but because the clinical application demands it. When you are processing 64 channels at 1024 Hz with a latency budget of a few milliseconds, you need direct control over memory layout, SIMD vectorization, and thread scheduling. You cannot get that through a Python wrapper around a Fortran library.
Building from the ground up also means we own every numerical decision. We know exactly how our bandpass filter handles edge effects. We know the precision characteristics of our FFT implementation at every stage. When a clinician asks why the system flagged a particular segment, we can trace the decision back through every transform and threshold, which is a regulatory requirement that black-box toolkits cannot easily satisfy.
The brain operates across multiple frequency scales simultaneously. Slow oscillations in the theta range (4-8 Hz) coordinate large neuronal populations across distributed cortical areas. Fast gamma oscillations (30-80 Hz) reflect local computation within small neural assemblies. These are not independent processes. The phase of the slow oscillation modulates the amplitude of the fast oscillation, a phenomenon called phase-amplitude coupling, or PAC. This is how the brain coordinates "what" is being processed (gamma) with "when" and "where" it should be routed (theta).
Quantifying PAC requires extracting instantaneous phase and amplitude from narrowband-filtered signals. We apply the Hilbert transform to obtain the analytic signal for each band, then compute the coupling strength using the Modulation Index as described by Tort and colleagues. The method works as follows: for each time point, record the instantaneous phase of the low-frequency signal and the instantaneous amplitude of the high-frequency signal. Bin the phase values into N equal-width bins spanning -pi to pi (typically 18 bins). For each bin, compute the mean high-frequency amplitude. If there is no coupling, the amplitude distribution across phase bins will be uniform. If there is coupling, certain phase bins will have systematically higher amplitudes.
The Modulation Index quantifies the deviation from uniformity using Kullback-Leibler divergence. A value of 0 means no coupling; values approaching 1 indicate strong, phase-locked amplitude modulation. Statistical significance is established by comparing against surrogate distributions generated by time-shifting the phase and amplitude series relative to each other, which preserves the spectral properties of each signal while destroying the temporal relationship between them.
Why this matters clinically: Patients with mild cognitive impairment (MCI) who go on to develop Alzheimer's show significantly reduced theta-gamma PAC compared to patients whose MCI remains stable. This reduction correlates with cognitive test scores (MMSE, MoCA), but critically, it does not correlate with CSF biomarkers like amyloid-beta or phosphorylated tau.
This dissociation is important. It means CFC is detecting functional network degradation that precedes the structural neurodegeneration visible through molecular biomarkers. The brain's coordination machinery is failing before the tissue itself shows classical Alzheimer's pathology. That window is where early intervention has the most potential, and CFC gives us a non-invasive way to see into it.
Spectral analysis is the foundation of any EEG processing pipeline. The question is not whether to estimate the power spectrum, but how to do it well. The standard EEG frequency bands are well-established:
| Band | Range | Associated Activity |
|---|---|---|
| Delta | 0.5 - 4 Hz | Deep sleep, pathological slowing |
| Theta | 4 - 8 Hz | Memory encoding, drowsiness |
| Alpha | 8 - 13 Hz | Relaxed wakefulness, posterior dominant rhythm |
| Beta | 13 - 30 Hz | Active thinking, motor planning |
| Gamma | >30 Hz | Local cortical processing, binding |
The naive approach is the periodogram: take an FFT of the entire signal, compute the squared magnitude. This gives an unbiased estimate of the power spectral density, but with terrible variance. Every spectral estimate is noisy because you have exactly one realization of the signal. In a clinical recording with 50/60 Hz line noise, muscle artifacts, and non-stationary brain dynamics, a raw periodogram is essentially unusable for quantitative analysis.
Welch's method improves on this by dividing the signal into overlapping segments (typically 50-75% overlap), applying a window function (Hanning or Hamming) to each segment, computing the periodogram of each, and averaging the results. Averaging reduces variance at the cost of frequency resolution: shorter segments mean wider frequency bins. This is the bias-variance tradeoff at its most explicit. For most clinical EEG applications, Welch's method with 2-4 second windows and 50% overlap provides adequate resolution and stability.
Multitaper spectral estimation (Thomson's method) takes a different approach. Instead of averaging over time segments, it applies multiple orthogonal window functions (Discrete Prolate Spheroidal Sequences, also called Slepian tapers) to the same data segment and averages the resulting spectra. This provides near-optimal bias-variance tradeoff for a given data length and bandwidth. The number of tapers controls the smoothing: more tapers reduce variance but widen the effective bandwidth. For moderate-length EEG recordings (10-60 seconds), multitaper estimation with 3-5 tapers often outperforms Welch's method, particularly for resolving closely spaced spectral peaks.
In practice, we use Welch's method for continuous monitoring where computational efficiency matters, and multitaper estimation for diagnostic segments where spectral accuracy is paramount. The choice is application-driven, not ideological.
The intuitive approach to artifact removal is bandpass filtering. Eye blinks produce large, slow deflections, so filter them out with a high-pass at 1 Hz. Muscle contamination is high-frequency, so low-pass at 40 Hz. This logic is wrong in a subtle but important way.
Bandpass filtering only works when the artifact and the signal of interest occupy non-overlapping frequency bands. In practice, this is almost never the case. Eye blinks generate potentials that extend well into the theta and alpha bands, which are exactly where the diagnostically relevant brain activity lives. Muscle EMG is broadband, covering 0 Hz to 200+ Hz, overlapping with every EEG band. A bandpass filter that removes an eye blink also destroys the frontal theta signal underneath it. You cannot separate two signals that share the same frequency band using frequency-domain methods alone.
Independent Component Analysis (ICA) solves this by separating signals based on statistical independence, not frequency. ICA assumes the recorded EEG at each electrode is a linear mixture of independent source signals: brain activity from different regions, eye movements, heartbeat, muscle tension. It finds the unmixing matrix that maximizes the statistical independence of the estimated sources. Once you have the components, you can identify and remove the artifact components (by their spatial topography and temporal characteristics) and reconstruct the EEG from the remaining neural components. An eye blink component can be removed without affecting the theta-band brain signal, because they are separated in the source domain, not the frequency domain.
Two ICA variants dominate in EEG processing. FastICA uses a fixed-point iteration to maximize non-Gaussianity (via negentropy approximation) and converges quickly, but assumes the sources are super-Gaussian. Extended Infomax handles both super- and sub-Gaussian source distributions by adapting its nonlinearity, making it more robust for EEG where source distributions vary. We use Extended Infomax for offline analysis and FastICA where speed is the priority.
For real-time applications, ICA is impractical because it requires a stationary data segment of sufficient length (typically 30+ seconds) to estimate the unmixing matrix reliably. This is where Artifact Subspace Reconstruction (ASR) becomes essential. ASR maintains a running model of "clean" EEG statistics from calibration data. During online processing, it uses a sliding-window PCA decomposition to identify principal components whose variance exceeds a threshold relative to the clean baseline. Those components are interpolated from neighboring channels rather than simply removed. The result is continuous, low-latency artifact correction without the stationarity requirements of ICA.
Writing correct algorithms is necessary but not sufficient. For real-time clinical use, the implementation must meet hard latency constraints. A 64-channel system sampled at 1024 Hz produces roughly 250 KB of data per second. That is not a lot of data, but the processing must complete within the sample interval, and the processing is not simple: filtering, artifact detection, spectral estimation, and feature extraction all need to happen before the next batch arrives.
SIMD vectorization is the primary optimization lever. Modern processors (ARM NEON on embedded systems, Intel SSE/AVX on workstations) can process 4 to 8 single-precision floats in a single instruction. For an FIR filter, this means computing 4-8 output samples per cycle instead of one. For FFT butterfly operations, the speedup is similar. In practice, we see a 4x throughput improvement for FFT and approximately 60% CPU reduction for FIR filtering when moving from scalar to SIMD code.
The key to effective SIMD is memory layout. Data must be contiguous and aligned. We use 16-byte alignment (for NEON) or 32-byte alignment (for AVX) with block processing: instead of processing one sample at a time through the entire pipeline, we process blocks of 256-1024 samples through each stage. This keeps the working set in L1 cache and avoids the latency penalty of main memory access, which on modern hardware can be 100x slower than cache.
For FPGA and microcontroller deployment, we use fixed-point arithmetic. EEG signals have limited dynamic range (typically 12-16 bits from the ADC), so 16-bit or 32-bit fixed-point representation loses negligible precision while eliminating the need for floating-point hardware. This is critical for low-power wearable devices where every milliwatt counts.
These optimizations are not academic exercises. They are what make it possible to run the full processing pipeline, from raw acquisition through spectral decomposition and CFC analysis, on a portable device with a power budget measured in hundreds of milliwatts. That is the difference between a laboratory instrument and a clinical tool that can be deployed at the point of care.