EEG 8-Channel Feature Analysis Pipeline

A MATLAB pipeline for computing quantitative EEG features from multi-channel neonatal EEG recordings, designed around an 8-channel (4 monopolar + 4 bipolar) neonatal montage.

Maintained by: Saeed Montazeri, BabaCenter, 2026 (Original pipeline: BabaCenter, 2019 — individual algorithm credits are kept in the relevant files; see Credits below.)

---

What it computes

For each subject and each epoch, the pipeline extracts:

Feature	Description
PSD	Power Spectral Density, per channel and frequency band
cPSD	Cross Power Spectral Density, per channel pair and frequency band
ASI	Activation Synchrony Index (interhemispheric synchrony)
wPLI	weighted Phase Lag Index (connectivity), per channel pair and frequency band
NC	Nestedness Coefficient / Phase-Amplitude Coupling (PAC)
aEEG	amplitude-integrated EEG — mean and interquartile range
rEEG	range-EEG — mean, interquartile range, and lower 5th percentile
SC	Suppression Curve

---

Requirements

• MATLAB (built and tested on R2019a; later versions should work)
• Signal Processing Toolbox (butter, filtfilt, pwelch, cpsd, envelope, hilbert, resample)
• Statistics and Machine Learning Toolbox (kmeans, prctile, iqr, pdist2)
• No third-party downloads needed — all helper functions (including the edfread-based EDF reader and the multiWaitbar progress bar) are bundled in this repository.

---

1. Preparing your data

1.1 Folder structure

Create the following folders in the same directory as feature_analysis_pipeline.m (or update the paths at the top of that script to point elsewhere):

your_project/
├── feature_analysis_pipeline.m
├── subfunctions/              <- put all the other .m files here
├── data/                      <- your .edf recordings go here
│   ├── subject01.edf
│   ├── subject02.edf
│   └── ...
├── artifact_annotations/      <- optional, see 1.3
└── results/                   <- created automatically if missing

1.2 EDF file requirements

Each .edf file must:

• Have a sampling rate ≥ 250 Hz.

• Contain (or allow derivation of) the following 8 channels, in this exact order, as they will be assembled by get_montages.m:

  1. F3 (monopolar/referential)
  2. F4 (monopolar/referential)
  3. P3 (monopolar/referential)
  4. P4 (monopolar/referential)
  5. F3–P3 ("Left" bipolar)
  6. F4–P4 ("Right" bipolar)
  7. F3–F4 ("Frontal" bipolar)
  8. P3–P4 ("Parietal" bipolar)

  Channel labels are matched case-insensitively against common variants (e.g. F3, F3-Ref, EEGF3-REF, F3P3, F3-P3, etc.) — see get_montages.m if your file uses a label not already recognized, and add it to the matching list there.

• If a required channel is missing from a file, that channel will simply be left as zeros for that subject — always sanity-check the console output and your channel labels beforehand.

1.3 Artifact annotations (optional but recommended)

The pipeline supports excluding artifact-contaminated samples from feature calculation. By default (feature_analysis_pipeline.m, section "Construct artifact masks"), it uses a dummy all-zero mask (i.e. no artifacts are excluded). To use real artifact annotations:

1. Produce, for each subject, a binary matrix of size [time_samples x 8 channels], where 1 = artifact/missing data and 0 = clean data, in the same channel order as above and at the same sampling rate as the raw EDF (i.e. before resampling).
2. Store these matrices in a cell array named artifact_masks, one cell per subject, in the same subject order as montages (i.e. the same order get_montages.m reads the .edf files in — alphabetical by filename).
3. Replace the dummy-mask section in feature_analysis_pipeline.m with code that loads your artifact_masks cell array (e.g. from a .mat file in artifact_annotations/) instead of generating zeros.

If you don't have an artifact detector, you can leave the dummy mask in place — no epochs will be excluded — but be aware this means artifacts in your raw data will affect the computed features.

---

2. Running the pipeline

1. Open feature_analysis_pipeline.m in MATLAB.
2. Edit the "FOLDER & FILE PATHS, OTHER VARIABLES" section at the top:
   • data_path, results_path, annotation_path — adjust if your folder names differ from the structure in §1.1.
   • epoch_length_seconds / epoch_overlap_seconds — epoch window and overlap (defaults: 120 s epochs, 60 s overlap).
   • global_highpass, global_lowpass, global_gate, new_sampling_rate — global filtering/resampling settings applied to all channels before feature extraction (defaults: 0.4–30 Hz bandpass, ±1000 µV artifact gate, resampled to 250 Hz).
   • artifact_percentage_thresholds — per-feature tolerance (0–1) for how much of an epoch may be marked as artifact before that feature is skipped for that epoch (see the comments in that section for details).
3. Run the script section by section (recommended, since importing and preprocessing large EDF datasets can take a while and it's useful to inspect intermediate results):
   • Import: reads all .edf files in data_path and assembles the 8-channel montage per subject (get_montages.m).
   • Artifact masks: loads/generates the artifact mask per subject.
   • Preprocessing: applies the global bandpass filter, amplitude gate, and resampling (preprocess.m).
   • Epoching: splits the continuous data (and artifact mask) into overlapping epochs (epoch_data.m).
   • Feature extraction: computes all features per subject/epoch (get_features.m) and saves them to features.mat in the current working directory.
4. Progress bars (via multiWaitbar) show import/preprocessing/epoching progress and per-subject/per-epoch progress during feature extraction, since this last step is the most time-consuming.

---

3. Interpreting the output

Running the pipeline saves a single file, features.mat, containing one variable:

features   % 1 x N_subjects struct array (N_subjects = number of .edf files)

Each element features(i) corresponds to one subject (in the same order the .edf files were read, i.e. alphabetical by filename in data_path), and has the following fields:

Field	Size	Meaning
ASI	[8 x 8 x epochs]	ASI between every channel pair, per epoch
wPLI	[8 x 8 x 4 x epochs]	wPLI per channel pair, in 4 frequency bands (0.4–3, 3–8, 8–13, 13–22 Hz), per epoch
NC	[2 x 3 x epochs]	Nestedness/PAC: row 1 = Frontal (avg. F3/F4), row 2 = Parietal (avg. P3/P4); 3 frequency bands (3–8, 8–15, 15–30 Hz)
aEEG	struct with .mean, .iqr, each [8 x epochs]	amplitude-integrated EEG per monopolar/bipolar channel
rEEG	struct with .mean, .iqr, .li, each [8 x epochs]	range-EEG mean, IQR, and lower 5th percentile per channel
PSD	[8 x 4 x epochs]	Power spectral density per channel, in 4 bands (1–3, 3–8, 8–15, 15–30 Hz)
cPSD	[8 x 8 x 4 x epochs]	Cross power spectral density per channel pair, same 4 bands as PSD
SC	[1 x epochs]	Suppression Curve value per epoch

Notes:

• The 8 channels are always in the fixed order: [F3, F4, P3, P4, Left(F3-P3), Right(F4-P4), Frontal(F3-F4), Parietal(P3-P4)].

• NaN values indicate that an epoch (or channel/channel-pair, for cross-channel features) was skipped because its artifact percentage exceeded the threshold set in artifact_percentage_thresholds for that feature.

• If a subject's .edf file could not be read, that subject's fields will simply be NaN instead of arrays.

• To load and inspect results in a new MATLAB session:

  load('features.mat')
  features(1).PSD          % PSD for subject 1
  squeeze(features(1).ASI(1,2,:))   % ASI between F3 and F4 across all epochs

---

Credits

• Pipeline maintained by Saeed Montazeri, BabaCenter, 2026.
• Original pipeline structure: BabaCenter, 2019.
• get_artifact_prcnt.m, get_line_length.m: original algorithms by Timo Vehviläinen, BabaCenter, 2019.
• my_aeeg.m: original algorithm by Nathan Stevenson, University of Helsinki, Finland, 2017.
• getASI.m / getFFT.m / ETDF.m: Activation Synchrony Index method by Okko Räsänen & Sampsa Vanhatalo (see Räsänen & Vanhatalo, NeuroImage, 2013).
• edfreadUntilDone.m: adapted from edfread.m by Brett Shoelson, PhD, © 2009–2012 The MathWorks, Inc.
• multiWaitbar.m: by Ben Tordoff, © 2007–2014 The MathWorks, Inc.

Please retain the above attributions in any redistribution of this code.