Combined EEG-fNIRS Pipeline: Multimodal Neuroimaging from Recording to Analysis

Combined EEG-fNIRS Pipeline

EEG measures the scalp potential generated by synchronised postsynaptic currents in cortical pyramidal neurons. It has millisecond temporal resolution but poor spatial resolution, and the signal is smeared by the skull and scalp. fNIRS measures changes in cortical haemoglobin concentration (HbO and HbR) via near-infrared light that penetrates the skull; it has centimetre-scale spatial resolution and a haemodynamic response that lags neural activity by ~5 seconds. Combining the two modalities gives complementary information: EEG reveals the when and what frequency of neural events, fNIRS reveals the where and the underlying metabolic demand.

The central challenge in a joint pipeline is not just preprocessing each modality correctly, but ensuring they are aligned in time, share the same event structure, and are fused in a way that exploits rather than discards the complementary information.

Step 1: Hardware Synchronisation

Both modalities must be time-locked to the same experimental events. The most common approaches are:

Shared trigger channel. A TTL pulse sent from the stimulus computer is recorded simultaneously by both the EEG amplifier and the fNIRS device. The trigger timestamps can then be used to epoch both signals to the same events.

Hardware synchroniser. Devices such as the NIRx Aurora or Artinis OxyMon have a dedicated synchronisation port that accepts the EEG trigger signal and stamps it in the fNIRS data stream.

Integrated cap. Several commercial systems (NIRx NIRScout, Artinis Brite) use a cap that houses both EEG electrodes and fNIRS optodes in fixed spatial positions, reducing co-registration error and ensuring no cable-induced motion artifacts are modality-specific.

Regardless of method, the key requirement is that the same event code appears in both data files at a time offset no larger than one sample of the slower modality (typically fNIRS at ~10–50 Hz).

Step 2: EEG Preprocessing

EEG is sampled at 250–2000 Hz and contains a wide range of artifacts — eye movements, muscle activity, heartbeat, line noise, and electrode pops. The standard pipeline in MNE-Python:

Filtering

import mne

raw_eeg = mne.io.read_raw_fif("sub-01_eeg.fif", preload=True)

# band-pass: remove slow drift and high-frequency muscle noise
raw_eeg.filter(l_freq=0.5, h_freq=50.0, method="fir")

# notch: suppress power-line interference
raw_eeg.notch_filter(freqs=50.0)

Bad channel rejection and interpolation

raw_eeg.plot(block=True)          # mark bads interactively
raw_eeg.interpolate_bads(reset_bads=True)  # spherical spline interpolation

Re-referencing

raw_eeg.set_eeg_reference("average", projection=True)
raw_eeg.apply_proj()

ICA for ocular and cardiac artifacts

from mne.preprocessing import ICA

ica = ICA(n_components=25, method="fastica", random_state=42)
ica.fit(raw_eeg, picks="eeg")

# find EOG components automatically
eog_indices, eog_scores = ica.find_bads_eog(raw_eeg)
ica.exclude = eog_indices

ica.apply(raw_eeg)

Epoching

events, event_id = mne.events_from_annotations(raw_eeg)

epochs_eeg = mne.Epochs(
    raw_eeg, events, event_id,
    tmin=-0.2, tmax=1.5,
    baseline=(-0.2, 0.0),
    preload=True,
    reject={"eeg": 100e-6},   # reject epochs with peak-to-peak > 100 µV
)

Step 3: fNIRS Preprocessing

fNIRS is sampled at 10–50 Hz and contains slow physiological fluctuations (heartbeat ~1 Hz, respiration ~0.3 Hz, Mayer waves ~0.1 Hz) that dwarf the haemodynamic response of interest. The pipeline must separate task-related HbO/HbR changes from these nuisance signals.

Load and convert to optical density

raw_fnirs = mne.io.read_raw_snirf("sub-01_fnirs.snirf", preload=True)

# convert raw light intensity to optical density changes
raw_od = mne.preprocessing.nirs.optical_density(raw_fnirs)

Bad channel detection

from mne.preprocessing.nirs import scalp_coupling_index

sci = scalp_coupling_index(raw_od)
raw_od.info["bads"] = [
    ch for ch, s in zip(raw_od.ch_names, sci) if s < 0.5
]

Motion artifact correction

Temporal Derivative Distribution Repair (TDDR) is robust to spiky motion artifacts without requiring manual threshold setting:

from mne.preprocessing.nirs import temporal_derivative_distribution_repair

raw_od_clean = temporal_derivative_distribution_repair(raw_od)

For larger slow-drift artifacts, a wavelet-based or spline interpolation approach is preferred. Short-channel regression removes global scalp hemodynamics when short-separation (<8 mm) channels are available:

from mne.preprocessing.nirs import short_channel_regression

raw_od_clean = short_channel_regression(raw_od_clean)

Convert to haemoglobin concentration (Modified Beer-Lambert Law)

The MBLL relates optical density change ΔOD to haemoglobin concentration change ΔC:

where ε is the molar extinction coefficient and L is the differential path length factor (DPF), which accounts for the scattering of light in tissue. Solving for two wavelengths gives ΔHbO and ΔHbR independently.

raw_haemo = mne.preprocessing.nirs.beer_lambert_law(raw_od_clean, ppf=0.1)

Band-pass filtering and epoching

raw_haemo.filter(l_freq=0.01, h_freq=0.1, method="fir")

events_fnirs, event_id_fnirs = mne.events_from_annotations(raw_haemo)

epochs_fnirs = mne.Epochs(
    raw_haemo, events_fnirs, event_id_fnirs,
    tmin=-5.0, tmax=20.0,       # fNIRS needs a long window for the HRF
    baseline=(-5.0, 0.0),
    preload=True,
)

Step 4: Temporal Alignment and Joint Epoching

Once both modalities are preprocessed, the shared trigger channel is the reference point for alignment. Because EEG and fNIRS are typically stored in separate files and may have been recorded on devices with slightly different clocks, a clock-drift correction is sometimes needed:

# compute the sample offset between modalities using shared triggers
eeg_trigger_times  = events_eeg[:, 0] / raw_eeg.info["sfreq"]
fnirs_trigger_times = events_fnirs[:, 0] / raw_fnirs.info["sfreq"]

drift_per_second = (
    (fnirs_trigger_times[-1] - eeg_trigger_times[-1]) -
    (fnirs_trigger_times[0]  - eeg_trigger_times[0])
) / (eeg_trigger_times[-1] - eeg_trigger_times[0])

After alignment, epochs from both modalities are indexed by the same event list, so trial N in epochs_eeg and trial N in epochs_fnirs correspond to the same stimulus presentation.

Step 5: Joint Analysis

GLM-based mass-univariate analysis

The canonical haemodynamic response function (HRF) is convolved with the stimulus boxcar to form a design matrix, then fit channel-by-channel:

from nilearn.glm.first_level import make_first_level_design_matrix
from nilearn.glm.first_level import run_glm

frame_times = epochs_fnirs.times
X = make_first_level_design_matrix(frame_times, events=events_df, hrf_model="spm")
labels, estimates = run_glm(epochs_fnirs.get_data().mean(axis=0).T, X.values)

EEG event-related potentials (ERPs) are computed in parallel from the EEG epochs:

evoked = epochs_eeg.average()
evoked.plot_topomap(times=[0.1, 0.2, 0.3])

The two maps — fNIRS beta weights and EEG topographies — can then be compared spatially if electrode and optode positions are co-registered to the same head model.

EEG band-power envelope correlated with fNIRS

A common analysis computes the analytic amplitude of an EEG frequency band (e.g. alpha 8–12 Hz) and correlates it with the HbO timecourse across trials or subjects:

import numpy as np
from scipy.signal import hilbert

# alpha envelope from EEG
alpha = epochs_eeg.copy().filter(8, 12).get_data()  # (n_trials, n_ch, n_times)
alpha_env = np.abs(hilbert(alpha, axis=-1))

# HbO from fNIRS (resample to EEG sampling rate for direct comparison)
hbo = epochs_fnirs.copy().pick("hbo").resample(epochs_eeg.info["sfreq"]).get_data()

# correlation across time for each EEG-fNIRS channel pair
corr = np.corrcoef(alpha_env.mean(axis=(0, 1)), hbo.mean(axis=(0, 1)))

Data-driven fusion

Step 6: Source-level integration

Both modalities can be projected to the cortical surface for a direct spatial comparison. EEG uses a minimum-norm estimate or LCMV beamformer (see the MEG beamforming post); fNIRS uses a diffuse optical tomography (DOT) reconstruction that inverts the photon transport model.

With MNE-Python and a shared FreeSurfer head model the workflow is:

# EEG source estimate
inv = mne.minimum_norm.make_inverse_operator(epochs_eeg.info, fwd, noise_cov)
stc_eeg = mne.minimum_norm.apply_inverse(evoked, inv, lambda2=1.0/9.0)

# fNIRS source estimate (requires MNE-NIRS or Cedalion)
# optodes projected to cortex via sensitivity matrix

Once both are on the same cortical surface, you can compute the spatial correlation of activity maps or constrain the EEG inverse solution with the fNIRS spatial prior.

Practical considerations

Physiological noise. The biggest challenge in fNIRS is separating task-related HbO changes (~1 µM) from cardiac (~50 µM) and respiratory (~10 µM) fluctuations. The 0.01–0.1 Hz band-pass helps, but is not sufficient on its own. Short-channel regression is the gold standard when the hardware supports it.

Scalp artifacts in EEG. When optodes and electrodes share the scalp, the fNIRS light sources can introduce broadband noise into EEG channels. Physically separating light sources from EEG electrodes by at least 1–2 cm and using shielded cables reduces cross-talk.

Head model co-registration. Optodes and electrodes must be digitised in the same coordinate system — either via a 3D digitiser (Polhemus), photogrammetry, or probabilistic atlases (for template-based analysis). Errors here propagate into any source-level fusion.

Haemodynamic lag. The fNIRS HRF peaks at ~5–8 seconds post-stimulus. Correlating EEG and fNIRS in a single epoch window requires the fNIRS epoch to extend far beyond the EEG window of interest (hence tmax=20.0 above vs tmax=1.5 for EEG).

Neurovascular coupling assumptions. The rationale for combining EEG and fNIRS rests on neurovascular coupling — the assumption that neural activity drives local blood flow. This coupling is well-established in healthy adults at rest, but can break down in pathology (stroke, dementia) or under pharmacological intervention.

Tooling summary

Further reading

• Zander et al. (2024). Emerging Neuroimaging Approach of Hybrid EEG-fNIRS Recordings: Data Collection and Analysis Challenges. Tandfonline
• Ortega-Martinez et al. (2025). Multimodal fNIRS-EEG sensor fusion: Review of data-driven methods and perspective for naturalistic brain imaging. Imaging Neuroscience. MIT Press
• Shin et al. (2023). Simultaneous multimodal fNIRS-EEG recordings reveal new insights in neural activity during motor execution, observation, and imagery. Scientific Reports. Nature
• MNE-Python fNIRS preprocessing tutorial
• NIRSTORM Brainstorm tutorial