Neurokit2

NeuroKit2

Overview

NeuroKit2 is a comprehensive Python toolkit for processing and analyzing physiological signals (biosignals). Use this skill to process cardiovascular, neural, autonomic, respiratory, and muscular signals for psychophysiology research, clinical applications, and human-computer interaction studies.

When to Use This Skill

Apply this skill when working with:

• Cardiac signals: ECG, PPG, heart rate variability (HRV), pulse analysis
• Brain signals: EEG frequency bands, microstates, complexity, source localization
• Autonomic signals: Electrodermal activity (EDA/GSR), skin conductance responses (SCR)
• Respiratory signals: Breathing rate, respiratory variability (RRV), volume per time
• Muscular signals: EMG amplitude, muscle activation detection
• Eye tracking: EOG, blink detection and analysis
• Multi-modal integration: Processing multiple physiological signals simultaneously
• Complexity analysis: Entropy measures, fractal dimensions, nonlinear dynamics

Core Capabilities

1. Cardiac Signal Processing (ECG/PPG)

Process electrocardiogram and photoplethysmography signals for cardiovascular analysis. See references/ecg_cardiac.md for detailed workflows.

Primary workflows:

• ECG processing pipeline: cleaning → R-peak detection → delineation → quality assessment
• HRV analysis across time, frequency, and nonlinear domains
• PPG pulse analysis and quality assessment
• ECG-derived respiration extraction

Key functions:

import neurokit2 as nk

# Complete ECG processing pipeline
signals, info = nk.ecg_process(ecg_signal, sampling_rate=1000)

# Analyze ECG data (event-related or interval-related)
analysis = nk.ecg_analyze(signals, sampling_rate=1000)

# Comprehensive HRV analysis
hrv = nk.hrv(peaks, sampling_rate=1000)  # Time, frequency, nonlinear domains

2. Heart Rate Variability Analysis

Compute comprehensive HRV metrics from cardiac signals. See references/hrv.md for all indices and domain-specific analysis.

Supported domains:

• Time domain: SDNN, RMSSD, pNN50, SDSD, and derived metrics
• Frequency domain: ULF, VLF, LF, HF, VHF power and ratios
• Nonlinear domain: Poincaré plot (SD1/SD2), entropy measures, fractal dimensions
• Specialized: Respiratory sinus arrhythmia (RSA), recurrence quantification analysis (RQA)

Key functions:

# All HRV indices at once
hrv_indices = nk.hrv(peaks, sampling_rate=1000)

# Domain-specific analysis
hrv_time = nk.hrv_time(peaks)
hrv_freq = nk.hrv_frequency(peaks, sampling_rate=1000)
hrv_nonlinear = nk.hrv_nonlinear(peaks, sampling_rate=1000)
hrv_rsa = nk.hrv_rsa(peaks, rsp_signal, sampling_rate=1000)

3. Brain Signal Analysis (EEG)

Analyze electroencephalography signals for frequency power, complexity, and microstate patterns. See references/eeg.md for detailed workflows and MNE integration.

Primary capabilities:

• Frequency band power analysis (Delta, Theta, Alpha, Beta, Gamma)
• Channel quality assessment and re-referencing
• Source localization (sLORETA, MNE)
• Microstate segmentation and transition dynamics
• Global field power and dissimilarity measures

Key functions:

# Power analysis across frequency bands
power = nk.eeg_power(eeg_data, sampling_rate=250, channels=['Fz', 'Cz', 'Pz'])

# Microstate analysis
microstates = nk.microstates_segment(eeg_data, n_microstates=4, method='kmod')
static = nk.microstates_static(microstates)
dynamic = nk.microstates_dynamic(microstates)

4. Electrodermal Activity (EDA)

Process skin conductance signals for autonomic nervous system assessment. See references/eda.md for detailed workflows.

Primary workflows:

• Signal decomposition into tonic and phasic components
• Skin conductance response (SCR) detection and analysis
• Sympathetic nervous system index calculation
• Autocorrelation and changepoint detection

Key functions:

# Complete EDA processing
signals, info = nk.eda_process(eda_signal, sampling_rate=100)

# Analyze EDA data
analysis = nk.eda_analyze(signals, sampling_rate=100)

# Sympathetic nervous system activity
sympathetic = nk.eda_sympathetic(signals, sampling_rate=100)

5. Respiratory Signal Processing (RSP)

Analyze breathing patterns and respiratory variability. See references/rsp.md for detailed workflows.

Primary capabilities:

• Respiratory rate calculation and variability analysis
• Breathing amplitude and symmetry assessment
• Respiratory volume per time (fMRI applications)
• Respiratory amplitude variability (RAV)

Key functions:

# Complete RSP processing
signals, info = nk.rsp_process(rsp_signal, sampling_rate=100)

# Respiratory rate variability
rrv = nk.rsp_rrv(signals, sampling_rate=100)

# Respiratory volume per time
rvt = nk.rsp_rvt(signals, sampling_rate=100)

6. Electromyography (EMG)

Process muscle activity signals for activation detection and amplitude analysis. See references/emg.md for workflows.

Key functions:

# Complete EMG processing
signals, info = nk.emg_process(emg_signal, sampling_rate=1000)

# Muscle activation detection
activation = nk.emg_activation(signals, sampling_rate=1000, method='threshold')

7. Electrooculography (EOG)

Analyze eye movement and blink patterns. See references/eog.md for workflows.

Key functions:

# Complete EOG processing
signals, info = nk.eog_process(eog_signal, sampling_rate=500)

# Extract blink features
features = nk.eog_features(signals, sampling_rate=500)

8. General Signal Processing

Apply filtering, decomposition, and transformation operations to any signal. See references/signal_processing.md for comprehensive utilities.

Key operations:

• Filtering (lowpass, highpass, bandpass, bandstop)
• Decomposition (EMD, SSA, wavelet)
• Peak detection and correction
• Power spectral density estimation
• Signal interpolation and resampling
• Autocorrelation and synchrony analysis

Key functions:

# Filtering
filtered = nk.signal_filter(signal, sampling_rate=1000, lowcut=0.5, highcut=40)

# Peak detection
peaks = nk.signal_findpeaks(signal)

# Power spectral density
psd = nk.signal_psd(signal, sampling_rate=1000)

9. Complexity and Entropy Analysis

Compute nonlinear dynamics, fractal dimensions, and information-theoretic measures. See references/complexity.md for all available metrics.

Available measures:

• Entropy: Shannon, approximate, sample, permutation, spectral, fuzzy, multiscale
• Fractal dimensions: Katz, Higuchi, Petrosian, Sevcik, correlation dimension
• Nonlinear dynamics: Lyapunov exponents, Lempel-Ziv complexity, recurrence quantification
• DFA: Detrended fluctuation analysis, multifractal DFA
• Information theory: Fisher information, mutual information

Key functions:

# Multiple complexity metrics at once
complexity_indices = nk.complexity(signal, sampling_rate=1000)

# Specific measures
apen = nk.entropy_approximate(signal)
dfa = nk.fractal_dfa(signal)
lyap = nk.complexity_lyapunov(signal, sampling_rate=1000)

10. Event-Related Analysis

Create epochs around stimulus events and analyze physiological responses. See references/epochs_events.md for workflows.

Primary capabilities:

• Epoch creation from event markers
• Event-related averaging and visualization
• Baseline correction options
• Grand average computation with confidence intervals

Key functions:

# Find events in signal
events = nk.events_find(trigger_signal, threshold=0.5)

# Create epochs around events
epochs = nk.epochs_create(signals, events, sampling_rate=1000,
                          epochs_start=-0.5, epochs_end=2.0)

# Average across epochs
grand_average = nk.epochs_average(epochs)

11. Multi-Signal Integration

Process multiple physiological signals simultaneously with unified output. See references/bio_module.md for integration workflows.

Key functions:

# Process multiple signals at once
bio_signals, bio_info = nk.bio_process(
    ecg=ecg_signal,
    rsp=rsp_signal,
    eda=eda_signal,
    emg=emg_signal,
    sampling_rate=1000
)

# Analyze all processed signals
bio_analysis = nk.bio_analyze(bio_signals, sampling_rate=1000)

Analysis Modes

NeuroKit2 automatically selects between two analysis modes based on data duration:

Event-related analysis (< 10 seconds):

• Analyzes stimulus-locked responses
• Epoch-based segmentation
• Suitable for experimental paradigms with discrete trials

Interval-related analysis (≥ 10 seconds):

• Characterizes physiological patterns over extended periods
• Resting state or continuous activities
• Suitable for baseline measurements and long-term monitoring

Most *_analyze() functions automatically choose the appropriate mode.

Installation

uv pip install neurokit2

For development version:

uv pip install https://github.com/neuropsychology/NeuroKit/zipball/dev

Common Workflows

Quick Start: ECG Analysis

import neurokit2 as nk

# Load example data
ecg = nk.ecg_simulate(duration=60, sampling_rate=1000)

# Process ECG
signals, info = nk.ecg_process(ecg, sampling_rate=1000)

# Analyze HRV
hrv = nk.hrv(info['ECG_R_Peaks'], sampling_rate=1000)

# Visualize
nk.ecg_plot(signals, info)

Multi-Modal Analysis

# Process multiple signals
bio_signals, bio_info = nk.bio_process(
    ecg=ecg_signal,
    rsp=rsp_signal,
    eda=eda_signal,
    sampling_rate=1000
)

# Analyze all signals
results = nk.bio_analyze(bio_signals, sampling_rate=1000)

Event-Related Potential

# Find events
events = nk.events_find(trigger_channel, threshold=0.5)

# Create epochs
epochs = nk.epochs_create(processed_signals, events,
                          sampling_rate=1000,
                          epochs_start=-0.5, epochs_end=2.0)

# Event-related analysis for each signal type
ecg_epochs = nk.ecg_eventrelated(epochs)
eda_epochs = nk.eda_eventrelated(epochs)

References

This skill includes comprehensive reference documentation organized by signal type and analysis method:

• ecg_cardiac.md: ECG/PPG processing, R-peak detection, delineation, quality assessment
• hrv.md: Heart rate variability indices across all domains
• eeg.md: EEG analysis, frequency bands, microstates, source localization
• eda.md: Electrodermal activity processing and SCR analysis
• rsp.md: Respiratory signal processing and variability
• ppg.md: Photoplethysmography signal analysis
• emg.md: Electromyography processing and activation detection
• eog.md: Electrooculography and blink analysis
• signal_processing.md: General signal utilities and transformations
• complexity.md: Entropy, fractal, and nonlinear measures
• epochs_events.md: Event-related analysis and epoch creation
• bio_module.md: Multi-signal integration workflows

Load specific reference files as needed using the Read tool to access detailed function documentation and parameters.

Additional Resources

• Official Documentation: https://neuropsychology.github.io/NeuroKit/
• GitHub Repository: https://github.com/neuropsychology/NeuroKit
• Publication: Makowski et al. (2021). NeuroKit2: A Python toolbox for neurophysiological signal processing. Behavior Research Methods. https://doi.org/10.3758/s13428-020-01516-y

Multi-Signal Integration (Bio Module)

Overview

The Bio module provides unified functions for processing and analyzing multiple physiological signals simultaneously. It acts as a wrapper that coordinates signal-specific processing functions and enables integrated multi-modal analysis.

Multi-Signal Processing

bio_process()

Process multiple physiological signals simultaneously with a single function call.

bio_signals, bio_info = nk.bio_process(ecg=None, rsp=None, eda=None, emg=None,
                                       ppg=None, eog=None, sampling_rate=1000)

Parameters:

• ecg: ECG signal array (optional)
• rsp: Respiratory signal array (optional)
• eda: EDA signal array (optional)
• emg: EMG signal array (optional)
• ppg: PPG signal array (optional)
• eog: EOG signal array (optional)
• sampling_rate: Sampling rate in Hz (must be consistent across signals or specify per signal)

Returns:

• bio_signals: Unified DataFrame containing all processed signals with columns:
• Signal-specific features (e.g., ECG_Clean, ECG_Rate, EDA_Phasic, RSP_Rate)
• All detected events/peaks
• Derived measures
• bio_info: Dictionary with signal-specific information (peak locations, parameters)

Example:

# Process ECG, respiration, and EDA simultaneously
bio_signals, bio_info = nk.bio_process(
    ecg=ecg_signal,
    rsp=rsp_signal,
    eda=eda_signal,
    sampling_rate=1000
)

# Access processed signals
ecg_clean = bio_signals['ECG_Clean']
rsp_rate = bio_signals['RSP_Rate']
eda_phasic = bio_signals['EDA_Phasic']

# Access detected peaks
ecg_peaks = bio_info['ECG']['ECG_R_Peaks']
rsp_peaks = bio_info['RSP']['RSP_Peaks']

Internal workflow: 1. Each signal is processed by its dedicated processing function:

• ecg_process() for ECG
• rsp_process() for respiration
• eda_process() for EDA
• emg_process() for EMG
• ppg_process() for PPG
• eog_process() for EOG

2. Results are merged into unified DataFrame 3. Cross-signal features computed (e.g., RSA if both ECG and RSP present)

Advantages:

• Simplified API for multi-modal recording
• Unified time base for all signals
• Automatic cross-signal feature computation
• Consistent output format

Multi-Signal Analysis

bio_analyze()

Perform comprehensive analysis on processed multi-modal signals.

bio_results = nk.bio_analyze(bio_signals, sampling_rate=1000)

Parameters:

• bio_signals: DataFrame from bio_process() or custom processed signals
• sampling_rate: Sampling rate (Hz)

Returns:

• DataFrame with analysis results for all detected signal types:
• Interval-related metrics if duration ≥ 10 seconds
• Event-related metrics if duration < 10 seconds
• Cross-signal indices (e.g., RSA if ECG + RSP available)

Computed metrics by signal:

• ECG: Heart rate statistics, HRV indices (time, frequency, nonlinear domains)
• RSP: Respiratory rate statistics, RRV, amplitude measures
• EDA: SCR count, amplitude, tonic level, sympathetic indices
• EMG: Activation count, amplitude statistics
• PPG: Similar to ECG (heart rate, HRV)
• EOG: Blink count, blink rate

Cross-signal metrics:

• RSA (Respiratory Sinus Arrhythmia): If ECG + RSP present
• Cardiorespiratory coupling: Phase synchronization indices
• Multi-modal arousal: Combined autonomic indices

Example:

# Analyze processed signals
results = nk.bio_analyze(bio_signals, sampling_rate=1000)

# Access results
heart_rate_mean = results['ECG_Rate_Mean']
hrv_rmssd = results['HRV_RMSSD']
breathing_rate = results['RSP_Rate_Mean']
scr_count = results['SCR_Peaks_N']
rsa_value = results['RSA']  # If both ECG and RSP present

Cross-Signal Features

When multiple signals are processed together, NeuroKit2 can compute integrated features:

Respiratory Sinus Arrhythmia (RSA)

Automatically computed when both ECG and respiratory signals are present.

bio_signals, bio_info = nk.bio_process(ecg=ecg, rsp=rsp, sampling_rate=1000)
results = nk.bio_analyze(bio_signals, sampling_rate=1000)

rsa = results['RSA']  # Automatically included

Computation:

• High-frequency HRV modulation by breathing
• Requires synchronized ECG R-peaks and respiratory signal
• Methods: Porges-Bohrer or peak-to-trough

Interpretation:

• Higher RSA: greater parasympathetic (vagal) influence
• Marker of cardiac-respiratory coupling
• Health indicator and emotion regulation capacity

ECG-Derived Respiration (EDR)

If respiratory signal unavailable, NeuroKit2 can estimate from ECG:

ecg_signals, ecg_info = nk.ecg_process(ecg, sampling_rate=1000)

# Extract EDR
edr = nk.ecg_rsp(ecg_signals['ECG_Clean'], sampling_rate=1000)

Use case:

• Estimate respiration when direct measurement unavailable
• Cross-validate respiratory measurements

Cardio-EDA Integration

Synchronized cardiac and electrodermal activity:

bio_signals, bio_info = nk.bio_process(ecg=ecg, eda=eda, sampling_rate=1000)

# Both signals available for integrated analysis
ecg_rate = bio_signals['ECG_Rate']
eda_phasic = bio_signals['EDA_Phasic']

# Compute correlations or coupling metrics
correlation = ecg_rate.corr(eda_phasic)

Practical Workflows

Complete Multi-Modal Recording Analysis

import neurokit2 as nk
import pandas as pd

# 1. Load multi-modal physiological data
ecg = load_ecg()        # Your data loading function
rsp = load_rsp()
eda = load_eda()
emg = load_emg()

# 2. Process all signals simultaneously
bio_signals, bio_info = nk.bio_process(
    ecg=ecg,
    rsp=rsp,
    eda=eda,
    emg=emg,
    sampling_rate=1000
)

# 3. Visualize processed signals
import matplotlib.pyplot as plt

fig, axes = plt.subplots(4, 1, figsize=(15, 12), sharex=True)

# ECG
axes[0].plot(bio_signals.index / 1000, bio_signals['ECG_Clean'])
axes[0].set_ylabel('ECG')
axes[0].set_title('Multi-Modal Physiological Recording')

# Respiration
axes[1].plot(bio_signals.index / 1000, bio_signals['RSP_Clean'])
axes[1].set_ylabel('Respiration')

# EDA
axes[2].plot(bio_signals.index / 1000, bio_signals['EDA_Phasic'])
axes[2].set_ylabel('EDA (Phasic)')

# EMG
axes[3].plot(bio_signals.index / 1000, bio_signals['EMG_Amplitude'])
axes[3].set_ylabel('EMG Amplitude')
axes[3].set_xlabel('Time (s)')

plt.tight_layout()
plt.show()

# 4. Analyze all signals
results = nk.bio_analyze(bio_signals, sampling_rate=1000)

# 5. Extract key metrics
print("Heart Rate (mean):", results['ECG_Rate_Mean'])
print("HRV (RMSSD):", results['HRV_RMSSD'])
print("Breathing Rate:", results['RSP_Rate_Mean'])
print("SCRs (count):", results['SCR_Peaks_N'])
print("RSA:", results['RSA'])

Event-Related Multi-Modal Analysis

# 1. Process signals
bio_signals, bio_info = nk.bio_process(ecg=ecg, rsp=rsp, eda=eda, sampling_rate=1000)

# 2. Detect events
events = nk.events_find(trigger_channel, threshold=0.5)

# 3. Create epochs for all signals
epochs = nk.epochs_create(bio_signals, events, sampling_rate=1000,
                          epochs_start=-1.0, epochs_end=10.0,
                          event_labels=event_labels,
                          event_conditions=event_conditions)

# 4. Signal-specific event-related analysis
ecg_eventrelated = nk.ecg_eventrelated(epochs)
rsp_eventrelated = nk.rsp_eventrelated(epochs)
eda_eventrelated = nk.eda_eventrelated(epochs)

# 5. Merge results
all_results = pd.merge(ecg_eventrelated, rsp_eventrelated,
                       left_index=True, right_index=True)
all_results = pd.merge(all_results, eda_eventrelated,
                       left_index=True, right_index=True)

# 6. Statistical comparison by condition
all_results['Condition'] = event_conditions
condition_means = all_results.groupby('Condition').mean()

Different Sampling Rates

Handle signals with different native sampling rates:

# ECG at 1000 Hz, EDA at 100 Hz
bio_signals, bio_info = nk.bio_process(
    ecg=ecg_1000hz,
    eda=eda_100hz,
    sampling_rate=1000  # Target sampling rate
)
# EDA will be automatically resampled to 1000 Hz internally

Or process separately and merge:

# Process at native sampling rates
ecg_signals, ecg_info = nk.ecg_process(ecg, sampling_rate=1000)
eda_signals, eda_info = nk.eda_process(eda, sampling_rate=100)

# Resample to common rate
eda_resampled = nk.signal_resample(eda_signals, sampling_rate=100,
                                   desired_sampling_rate=1000)

# Merge manually
bio_signals = pd.concat([ecg_signals, eda_resampled], axis=1)

Use Cases and Applications

Comprehensive Psychophysiology Research

Capture multiple dimensions of physiological arousal:

• Cardiac: Orienting, attention, emotional valence
• Respiratory: Arousal, stress, emotion regulation
• EDA: Sympathetic arousal, emotional intensity
• EMG: Muscle tension, facial expression, startle

Example: Emotional picture viewing

• ECG: Heart rate deceleration during picture viewing (attention)
• EDA: SCRs reflect emotional arousal intensity
• RSP: Breath-holding or changes reflect emotional engagement
• Facial EMG: Corrugator (frown), zygomaticus (smile) for valence

Stress and Relaxation Assessment

Multi-modal markers provide convergent evidence:

• Increased stress: ↑ HR, ↓ HRV, ↑ EDA, ↑ respiration rate, ↑ muscle tension
• Relaxation: ↓ HR, ↑ HRV, ↓ EDA, ↓ respiration rate, slow breathing, ↓ muscle tension

Intervention effectiveness:

• Compare multi-modal indices before vs. after intervention
• Identify which modalities respond to specific techniques

Clinical Assessment

Anxiety disorders:

• Heightened baseline EDA, HR
• Exaggerated responses to stressors
• Reduced HRV, respiratory variability

Depression:

• Altered autonomic balance (↓ HRV)
• Blunted EDA responses
• Irregular respiratory patterns

PTSD:

• Hyperarousal (↑ HR, ↑ EDA baseline)
• Exaggerated startle (EMG)
• Altered RSA

Human-Computer Interaction

Unobtrusive user state assessment:

• Cognitive load: ↓ HRV, ↑ EDA, suppressed blinks
• Frustration: ↑ HR, ↑ EDA, ↑ muscle tension
• Engagement: Moderate arousal, synchronized responses
• Boredom: Low arousal, irregular patterns

Athletic Performance and Recovery

Monitor training load and recovery:

• Resting HRV: Daily monitoring for overtraining
• EDA: Sympathetic activation and stress
• Respiration: Breathing patterns during exercise/recovery
• Multi-modal integration: Comprehensive recovery assessment

Advantages of Multi-Modal Recording

Convergent validity:

• Multiple indices converge on same construct (e.g., arousal)
• More robust than single measure

Discriminant validity:

• Different signals dissociate under certain conditions
• ECG reflects both sympathetic and parasympathetic
• EDA reflects primarily sympathetic

System integration:

• Understand whole-body physiological coordination
• Cross-signal coupling metrics (RSA, coherence)

Redundancy and robustness:

• If one signal quality poor, others available
• Cross-validate findings across modalities

Richer interpretation:

• HR deceleration + SCR increase = orienting with arousal
• HR acceleration + no SCR = cardiac response without sympathetic arousal

Considerations

Hardware and Synchronization

• Same device: Signals inherently synchronized
• Different devices: Requires common trigger/timestamp
• Use hardware trigger to mark simultaneous events
• Software alignment based on event markers
• Verify synchronization quality (cross-correlate redundant signals)

Signal Quality Across Modalities

• Not all signals may have equal quality
• Prioritize based on research question
• Document quality issues per signal

Computational Cost

• Processing multiple signals increases computation time
• Consider processing in batches for large datasets
• Downsample appropriately to reduce load

Analysis Complexity

• More signals = more variables = more statistical comparisons
• Risk of Type I error (false positives) without correction
• Use multivariate approaches or pre-registered analyses

Interpretation

• Avoid over-interpretation of complex multi-modal patterns
• Ground in physiological theory
• Replicate findings before strong claims

References

• Berntson, G. G., Cacioppo, J. T., & Quigley, K. S. (1993). Respiratory sinus arrhythmia: autonomic origins, physiological mechanisms, and psychophysiological implications. Psychophysiology, 30(2), 183-196.
• Cacioppo, J. T., Tassinary, L. G., & Berntson, G. (Eds.). (2017). Handbook of psychophysiology (4th ed.). Cambridge University Press.
• Kreibig, S. D. (2010). Autonomic nervous system activity in emotion: A review. Biological psychology, 84(3), 394-421.
• Laborde, S., Mosley, E., & Thayer, J. F. (2017). Heart rate variability and cardiac vagal tone in psychophysiological research–recommendations for experiment planning, data analysis, and data reporting. Frontiers in psychology, 8, 213.

Complexity and Entropy Analysis

Overview

Complexity measures quantify the irregularity, unpredictability, and multiscale structure of time series signals. NeuroKit2 provides comprehensive entropy, fractal dimension, and nonlinear dynamics measures for assessing physiological signal complexity.

Main Function

complexity()

Compute multiple complexity metrics simultaneously for exploratory analysis.

complexity_indices = nk.complexity(signal, sampling_rate=1000, show=False)

Returns:

• DataFrame with numerous complexity measures across categories:
• Entropy indices
• Fractal dimensions
• Nonlinear dynamics measures
• Information-theoretic metrics

Use case:

• Exploratory analysis to identify relevant measures
• Comprehensive signal characterization
• Comparative studies across signals

Parameter Optimization

Before computing complexity measures, optimal embedding parameters should be determined:

complexity_delay()

Determine optimal time delay (τ) for phase space reconstruction.

optimal_tau = nk.complexity_delay(signal, delay_max=100, method='fraser1986', show=False)

Methods:

• 'fraser1986': Mutual information first minimum
• 'theiler1990': Autocorrelation first zero crossing
• 'casdagli1991': Cao's method

Use for: Embedding delay in entropy, attractor reconstruction

complexity_dimension()

Determine optimal embedding dimension (m).

optimal_m = nk.complexity_dimension(signal, delay=None, dimension_max=20,
                                    method='afn', show=False)

Methods:

• 'afn': Average False Nearest Neighbors
• 'fnn': False Nearest Neighbors
• 'correlation': Correlation dimension saturation

Use for: Entropy calculations, phase space reconstruction

complexity_tolerance()

Determine optimal tolerance (r) for entropy measures.

optimal_r = nk.complexity_tolerance(signal, method='sd', show=False)

Methods:

• 'sd': Standard deviation-based (0.1-0.25 × SD typical)
• 'maxApEn': Maximize ApEn
• 'recurrence': Based on recurrence rate

Use for: Approximate entropy, sample entropy

complexity_k()

Determine optimal k parameter for Higuchi fractal dimension.

optimal_k = nk.complexity_k(signal, k_max=20, show=False)

Use for: Higuchi fractal dimension calculation

Entropy Measures

Entropy quantifies randomness, unpredictability, and information content.

entropy_shannon()

Shannon entropy - classical information-theoretic measure.

shannon_entropy = nk.entropy_shannon(signal)

Interpretation:

• Higher: more random, less predictable
• Lower: more regular, predictable
• Units: bits (information)

Use cases:

• General randomness assessment
• Information content
• Signal irregularity

entropy_approximate()

Approximate Entropy (ApEn) - regularity of patterns.

apen = nk.entropy_approximate(signal, delay=1, dimension=2, tolerance='sd')

Parameters:

• delay: Time delay (τ)
• dimension: Embedding dimension (m)
• tolerance: Similarity threshold (r)

Interpretation:

• Lower ApEn: more regular, self-similar patterns
• Higher ApEn: more complex, irregular
• Sensitive to signal length (≥100-300 points recommended)

Physiological applications:

• HRV: reduced ApEn in heart disease
• EEG: altered ApEn in neurological disorders

entropy_sample()

Sample Entropy (SampEn) - improved ApEn.

sampen = nk.entropy_sample(signal, delay=1, dimension=2, tolerance='sd')

Advantages over ApEn:

• Less dependent on signal length
• More consistent across recordings
• No self-matching bias

Interpretation:

• Same as ApEn but more reliable
• Preferred in most applications

Typical values:

• HRV: 0.5-2.5 (context-dependent)
• EEG: 0.3-1.5

entropy_multiscale()

Multiscale Entropy (MSE) - complexity across temporal scales.

mse = nk.entropy_multiscale(signal, scale=20, dimension=2, tolerance='sd',
                            method='MSEn', show=False)

Methods:

• 'MSEn': Multiscale Sample Entropy
• 'MSApEn': Multiscale Approximate Entropy
• 'CMSE': Composite Multiscale Entropy
• 'RCMSE': Refined Composite Multiscale Entropy

Interpretation:

• Entropy at different coarse-graining scales
• Healthy/complex systems: high entropy across multiple scales
• Diseased/simpler systems: reduced entropy, especially at larger scales

Use cases:

• Distinguish true complexity from randomness
• White noise: constant across scales
• Pink noise/complexity: structured variation across scales

entropy_fuzzy()

Fuzzy Entropy - uses fuzzy membership functions.

fuzzen = nk.entropy_fuzzy(signal, delay=1, dimension=2, tolerance='sd', r=0.2)

Advantages:

• More stable with noisy signals
• Fuzzy boundaries for pattern matching
• Better performance with short signals

entropy_permutation()

Permutation Entropy - based on ordinal patterns.

perment = nk.entropy_permutation(signal, delay=1, dimension=3)

Method:

• Encodes signal into ordinal patterns (permutations)
• Counts pattern frequencies
• Robust to noise and non-stationarity

Interpretation:

• Lower: more regular ordinal structure
• Higher: more random ordering

Use cases:

• EEG analysis
• Anesthesia depth monitoring
• Fast computation

entropy_spectral()

Spectral Entropy - based on power spectrum.

spec_ent = nk.entropy_spectral(signal, sampling_rate=1000, bands=None)

Method:

• Normalized Shannon entropy of power spectrum
• Quantifies frequency distribution regularity

Interpretation:

• 0: Single frequency (pure tone)
• 1: White noise (flat spectrum)

Use cases:

• EEG: spectral distribution changes with states
• Anesthesia monitoring

entropy_svd()

Singular Value Decomposition Entropy.

svd_ent = nk.entropy_svd(signal, delay=1, dimension=2)

Method:

• SVD on trajectory matrix
• Entropy of singular value distribution

Use cases:

• Attractor complexity
• Deterministic vs. stochastic dynamics

entropy_differential()

Differential Entropy - continuous analog of Shannon entropy.

diff_ent = nk.entropy_differential(signal)

Use for: Continuous probability distributions

Other Entropy Measures

Tsallis Entropy:

tsallis = nk.entropy_tsallis(signal, q=2)

• Generalized entropy with parameter q
• q=1 reduces to Shannon entropy

Rényi Entropy:

renyi = nk.entropy_renyi(signal, alpha=2)

• Generalized entropy with parameter α

Additional specialized entropies:

• entropy_attention(): Attention entropy
• entropy_grid(): Grid-based entropy
• entropy_increment(): Increment entropy
• entropy_slope(): Slope entropy
• entropy_dispersion(): Dispersion entropy
• entropy_symbolicdynamic(): Symbolic dynamics entropy
• entropy_range(): Range entropy
• entropy_phase(): Phase entropy
• entropy_quadratic(), entropy_cumulative_residual(), entropy_rate(): Specialized variants

Fractal Dimension Measures

Fractal dimensions characterize self-similarity and roughness.

fractal_katz()

Katz Fractal Dimension - waveform complexity.

kfd = nk.fractal_katz(signal)

Interpretation:

• 1: straight line
• >1: increasing roughness and complexity
• Typical range: 1.0-2.0

Advantages:

• Simple, fast computation
• No parameter tuning

fractal_higuchi()

Higuchi Fractal Dimension - self-similarity.

hfd = nk.fractal_higuchi(signal, k_max=10)

Method:

• Constructs k new time series from original
• Estimates dimension from length-scale relationship

Interpretation:

• Higher HFD: more complex, irregular
• Lower HFD: smoother, more regular

Use cases:

• EEG complexity
• HRV analysis
• Epilepsy detection

fractal_petrosian()

Petrosian Fractal Dimension - rapid estimation.

pfd = nk.fractal_petrosian(signal)

Advantages:

• Fast computation
• Direct calculation (no curve fitting)

fractal_sevcik()

Sevcik Fractal Dimension - normalized waveform complexity.

sfd = nk.fractal_sevcik(signal)

fractal_nld()

Normalized Length Density - curve length-based measure.

nld = nk.fractal_nld(signal)

fractal_psdslope()

Power Spectral Density Slope - frequency-domain fractal measure.

slope = nk.fractal_psdslope(signal, sampling_rate=1000)

Method:

• Linear fit to log-log power spectrum
• Slope β relates to fractal dimension

Interpretation:

• β ≈ 0: White noise (random)
• β ≈ -1: Pink noise (1/f, complex)
• β ≈ -2: Brown noise (Brownian motion)

fractal_hurst()

Hurst Exponent - long-range dependence.

hurst = nk.fractal_hurst(signal, show=False)

Interpretation:

• H < 0.5: Anti-persistent (mean-reverting)
• H = 0.5: Random walk (white noise)
• H > 0.5: Persistent (trending, long-memory)

Use cases:

• Assess long-term correlations
• Financial time series
• HRV analysis

fractal_correlation()

Correlation Dimension - attractor dimensionality.

corr_dim = nk.fractal_correlation(signal, delay=1, dimension=10, radius=64)

Method:

• Grassberger-Procaccia algorithm
• Estimates dimension of attractor in phase space

Interpretation:

• Low dimension: deterministic, low-dimensional chaos
• High dimension: high-dimensional chaos or noise

fractal_dfa()

Detrended Fluctuation Analysis - scaling exponent.

dfa_alpha = nk.fractal_dfa(signal, multifractal=False, q=2, show=False)

Interpretation:

• α < 0.5: Anti-correlated
• α = 0.5: Uncorrelated (white noise)
• α = 1.0: 1/f noise (pink noise, healthy complexity)
• α = 1.5: Brownian noise
• α > 1.0: Persistent long-range correlations

HRV applications:

• α1 (short-term, 4-11 beats): Reflects autonomic regulation
• α2 (long-term, >11 beats): Long-range correlations
• Reduced α1: Cardiac pathology

fractal_mfdfa()

Multifractal DFA - multiscale fractal properties.

mfdfa_results = nk.fractal_mfdfa(signal, q=None, show=False)

Method:

• Extends DFA to multiple q-orders
• Characterizes multifractal spectrum

Returns:

• Generalized Hurst exponents h(q)
• Multifractal spectrum f(α)
• Width indicates multifractality strength

Use cases:

• Detect multifractal structure
• HRV multifractality in health vs. disease
• EEG multiscale dynamics

fractal_tmf()

Multifractal Nonlinearity - deviation from monofractal.

tmf = nk.fractal_tmf(signal)

Interpretation:

• Quantifies departure from simple scaling
• Higher: more multifractal structure

fractal_density()

Density Fractal Dimension.

density_fd = nk.fractal_density(signal)

fractal_linelength()

Line Length - total variation measure.

linelength = nk.fractal_linelength(signal)

Use case:

• Simple complexity proxy
• EEG seizure detection

Nonlinear Dynamics

complexity_lyapunov()

Largest Lyapunov Exponent - chaos and divergence.

lyap = nk.complexity_lyapunov(signal, delay=None, dimension=None,
                              sampling_rate=1000, show=False)

Interpretation:

• λ < 0: Stable fixed point
• λ = 0: Periodic orbit
• λ > 0: Chaotic (nearby trajectories diverge exponentially)

Use cases:

• Detect chaos in physiological signals
• HRV: positive Lyapunov suggests nonlinear dynamics
• EEG: epilepsy detection (decreased λ before seizure)

complexity_lempelziv()

Lempel-Ziv Complexity - algorithmic complexity.

lz = nk.complexity_lempelziv(signal, symbolize='median')

Method:

• Counts number of distinct patterns
• Coarse-grained measure of randomness

Interpretation:

• Lower: repetitive, predictable patterns
• Higher: diverse, unpredictable patterns

Use cases:

• EEG: consciousness levels, anesthesia
• HRV: autonomic complexity

complexity_rqa()

Recurrence Quantification Analysis - phase space recurrences.

rqa_indices = nk.complexity_rqa(signal, delay=1, dimension=3, tolerance='sd')

Metrics:

• Recurrence Rate (RR): Percentage of recurrent states
• Determinism (DET): Percentage of recurrent points in lines
• Laminarity (LAM): Percentage in vertical structures (laminar states)
• Trapping Time (TT): Average vertical line length
• Longest diagonal/vertical: System predictability
• Entropy (ENTR): Shannon entropy of line length distribution

Interpretation:

• High DET: deterministic dynamics
• High LAM: system trapped in specific states
• Low RR: random, non-recurrent dynamics

Use cases:

• Detect transitions in system dynamics
• Physiological state changes
• Nonlinear time series analysis

complexity_hjorth()

Hjorth Parameters - time-domain complexity.

hjorth = nk.complexity_hjorth(signal)

Metrics:

• Activity: Variance of signal
• Mobility: Proportion of standard deviation of derivative to signal
• Complexity: Change in mobility with derivative

Use cases:

• EEG feature extraction
• Seizure detection
• Signal characterization

complexity_decorrelation()

Decorrelation Time - memory duration.

decorr_time = nk.complexity_decorrelation(signal, show=False)

Interpretation:

• Time lag where autocorrelation drops below threshold
• Shorter: rapid fluctuations, short memory
• Longer: slow fluctuations, long memory

complexity_relativeroughness()

Relative Roughness - smoothness measure.

roughness = nk.complexity_relativeroughness(signal)

Information Theory

fisher_information()

Fisher Information - measure of order.

fisher = nk.fisher_information(signal, delay=1, dimension=2)

Interpretation:

• High: ordered, structured
• Low: disordered, random

Use cases:

• Combine with Shannon entropy (Fisher-Shannon plane)
• Characterize system complexity

fishershannon_information()

Fisher-Shannon Information Product.

fs = nk.fishershannon_information(signal)

Method:

• Product of Fisher information and Shannon entropy
• Characterizes order-disorder balance

mutual_information()

Mutual Information - shared information between variables.

mi = nk.mutual_information(signal1, signal2, method='knn')

Methods:

• 'knn': k-nearest neighbors (nonparametric)
• 'kernel': Kernel density estimation
• 'binning': Histogram-based

Use cases:

• Coupling between signals
• Feature selection
• Nonlinear dependence

Practical Considerations

Signal Length Requirements

Parameter Selection

General guidelines:

• Use parameter optimization functions first
• Or use conventional defaults:
• Delay (τ): 1 for HRV, autocorrelation first minimum for EEG
• Dimension (m): 2-3 typical
• Tolerance (r): 0.2 × SD common

Sensitivity:

• Results can be parameter-sensitive
• Report parameters used
• Consider sensitivity analysis

Normalization and Preprocessing

Standardization:

• Many measures sensitive to signal amplitude
• Z-score normalization often recommended
• Detrending may be necessary

Stationarity:

• Some measures assume stationarity
• Check with statistical tests (e.g., ADF test)
• Segment non-stationary signals

Interpretation

Context-dependent:

• No universal "good" or "bad" complexity
• Compare within-subject or between groups
• Consider physiological context

Complexity vs. randomness:

• Maximum entropy ≠ maximum complexity
• True complexity: structured variability
• White noise: high entropy but low complexity (MSE distinguishes)

Applications

Cardiovascular:

• HRV complexity: reduced in heart disease, aging
• DFA α1: prognostic marker post-MI

Neuroscience:

• EEG complexity: consciousness, anesthesia depth
• Entropy: Alzheimer's, epilepsy, sleep stages
• Permutation entropy: anesthesia monitoring

Psychology:

• Complexity loss in depression, anxiety
• Increased regularity under stress

Aging:

• "Complexity loss" with aging across systems
• Reduced multiscale complexity

Critical transitions:

• Complexity changes before state transitions
• Early warning signals (critical slowing down)

References

• Pincus, S. M. (1991). Approximate entropy as a measure of system complexity. Proceedings of the National Academy of Sciences, 88(6), 2297-2301.
• Richman, J. S., & Moorman, J. R. (2000). Physiological time-series analysis using approximate entropy and sample entropy. American Journal of Physiology-Heart and Circulatory Physiology, 278(6), H2039-H2049.
• Peng, C. K., et al. (1995). Quantification of scaling exponents and crossover phenomena in nonstationary heartbeat time series. Chaos, 5(1), 82-87.
• Costa, M., Goldberger, A. L., & Peng, C. K. (2005). Multiscale entropy analysis of biological signals. Physical review E, 71(2), 021906.
• Grassberger, P., & Procaccia, I. (1983). Measuring the strangeness of strange attractors. Physica D: Nonlinear Phenomena, 9(1-2), 189-208.

ECG and Cardiac Signal Processing

Overview

Process electrocardiogram (ECG) and photoplethysmography (PPG) signals for cardiovascular analysis. This module provides comprehensive tools for R-peak detection, waveform delineation, quality assessment, and heart rate analysis.

Main Processing Pipeline

ecg_process()

Complete automated ECG processing pipeline that orchestrates multiple steps.

signals, info = nk.ecg_process(ecg_signal, sampling_rate=1000, method='neurokit')

Pipeline steps: 1. Signal cleaning (noise removal) 2. R-peak detection 3. Heart rate calculation 4. Quality assessment 5. QRS delineation (P, Q, S, T waves) 6. Cardiac phase determination

Returns:

• signals: DataFrame with cleaned ECG, peaks, rate, quality, cardiac phases
• info: Dictionary with R-peak locations and processing parameters

Common methods:

• 'neurokit' (default): Comprehensive NeuroKit2 pipeline
• 'biosppy': BioSPPy-based processing
• 'pantompkins1985': Pan-Tompkins algorithm
• 'hamilton2002', 'elgendi2010', 'engzeemod2012': Alternative methods

Preprocessing Functions

ecg_clean()

Remove noise from raw ECG signals using method-specific filtering.

cleaned_ecg = nk.ecg_clean(ecg_signal, sampling_rate=1000, method='neurokit')

Methods:

• 'neurokit': High-pass Butterworth filter (0.5 Hz) + powerline filtering
• 'biosppy': FIR filtering between 0.67-45 Hz
• 'pantompkins1985': Band-pass 5-15 Hz + derivative-based
• 'hamilton2002': Band-pass 8-16 Hz
• 'elgendi2010': Band-pass 8-20 Hz
• 'engzeemod2012': FIR band-pass 0.5-40 Hz

Key parameters:

• powerline: Remove 50 or 60 Hz powerline noise (default: 50)

ecg_peaks()

Detect R-peaks in ECG signals with optional artifact correction.

peaks_dict, info = nk.ecg_peaks(cleaned_ecg, sampling_rate=1000, method='neurokit', correct_artifacts=False)

Available methods (13+ algorithms):

• 'neurokit': Hybrid approach optimized for reliability
• 'pantompkins1985': Classic Pan-Tompkins algorithm
• 'hamilton2002': Hamilton's adaptive threshold
• 'christov2004': Christov's adaptive method
• 'gamboa2008': Gamboa's approach
• 'elgendi2010': Elgendi's two moving averages
• 'engzeemod2012': Modified Engelse-Zeelenberg
• 'kalidas2017': XQRS-based
• 'martinez2004', 'rodrigues2021', 'koka2022', 'promac': Advanced methods

Artifact correction: Set correct_artifacts=True to apply Lipponen & Tarvainen (2019) correction:

• Detects ectopic beats, long/short intervals, missed beats
• Uses threshold-based detection with configurable parameters

Returns:

• Dictionary with 'ECG_R_Peaks' key containing peak indices

ecg_delineate()

Identify P, Q, S, T waves and their onsets/offsets.

waves, waves_peak = nk.ecg_delineate(cleaned_ecg, rpeaks, sampling_rate=1000, method='dwt')

Methods:

• 'dwt' (default): Discrete wavelet transform-based detection
• 'peak': Simple peak detection around R-peaks
• 'cwt': Continuous wavelet transform (Martinez et al., 2004)

Detected components:

• P waves: ECG_P_Peaks, ECG_P_Onsets, ECG_P_Offsets
• Q waves: ECG_Q_Peaks
• S waves: ECG_S_Peaks
• T waves: ECG_T_Peaks, ECG_T_Onsets, ECG_T_Offsets
• QRS complex: onsets and offsets

Returns:

• waves: Dictionary with all wave indices
• waves_peak: Dictionary with peak amplitudes

ecg_quality()

Assess ECG signal integrity and quality.

quality = nk.ecg_quality(ecg_signal, rpeaks=None, sampling_rate=1000, method='averageQRS')

Methods:

• 'averageQRS' (default): Template matching correlation (Zhao & Zhang, 2018)
• Returns quality scores 0-1 for each heartbeat
• Threshold: >0.6 = good quality
• 'zhao2018': Multi-index approach using kurtosis, power spectrum distribution

Use cases:

• Identify low-quality signal segments
• Filter out noisy heartbeats before analysis
• Validate R-peak detection accuracy

Analysis Functions

ecg_analyze()

High-level analysis that automatically selects event-related or interval-related mode.

analysis = nk.ecg_analyze(signals, sampling_rate=1000, method='auto')

Mode selection:

• Duration < 10 seconds → event-related analysis
• Duration ≥ 10 seconds → interval-related analysis

Returns: DataFrame with cardiac metrics appropriate for the analysis mode.

ecg_eventrelated()

Analyze stimulus-locked ECG epochs for event-related responses.

results = nk.ecg_eventrelated(epochs)

Computed metrics:

• ECG_Rate_Baseline: Mean heart rate before stimulus
• ECG_Rate_Min/Max: Minimum/maximum heart rate during epoch
• ECG_Phase_Atrial/Ventricular: Cardiac phase at stimulus onset
• Rate dynamics across epoch time windows

Use case: Experimental paradigms with discrete trials (e.g., stimulus presentations, task events).

ecg_intervalrelated()

Analyze continuous ECG recordings for resting state or extended periods.

results = nk.ecg_intervalrelated(signals, sampling_rate=1000)

Computed metrics:

• ECG_Rate_Mean: Average heart rate over interval
• Comprehensive HRV metrics (delegates to hrv() function)
• Time domain: SDNN, RMSSD, pNN50, etc.
• Frequency domain: LF, HF, LF/HF ratio
• Nonlinear: Poincaré, entropy, fractal measures

Minimum duration:

• Basic rate: Any duration
• HRV frequency metrics: ≥60 seconds recommended, 1-5 minutes optimal

Utility Functions

ecg_rate()

Compute instantaneous heart rate from R-peak intervals.

heart_rate = nk.ecg_rate(peaks, sampling_rate=1000, desired_length=None)

Method:

• Calculates inter-beat intervals (IBIs) between consecutive R-peaks
• Converts to beats per minute (BPM): 60 / IBI
• Interpolates to match signal length if desired_length specified

Returns:

• Array of instantaneous heart rate values

ecg_phase()

Determine atrial and ventricular systole/diastole phases.

cardiac_phase = nk.ecg_phase(ecg_cleaned, rpeaks, delineate_info)

Phases computed:

• ECG_Phase_Atrial: Atrial systole (1) vs. diastole (0)
• ECG_Phase_Ventricular: Ventricular systole (1) vs. diastole (0)
• ECG_Phase_Completion_Atrial/Ventricular: Percentage of phase completion (0-1)

Use case:

• Cardiac-locked stimulus presentation
• Psychophysiology experiments timing events to cardiac cycle

ecg_segment()

Extract individual heartbeats for morphological analysis.

heartbeats = nk.ecg_segment(ecg_cleaned, rpeaks, sampling_rate=1000)

Returns:

• Dictionary of epochs, each containing one heartbeat
• Centered on R-peak with configurable pre/post windows
• Useful for beat-to-beat morphology comparison

ecg_invert()

Detect and correct inverted ECG signals automatically.

corrected_ecg, is_inverted = nk.ecg_invert(ecg_signal, sampling_rate=1000)

Method:

• Analyzes QRS complex polarity
• Flips signal if predominantly negative
• Returns corrected signal and inversion status

ecg_rsp()

Extract ECG-derived respiration (EDR) as respiratory proxy signal.

edr_signal = nk.ecg_rsp(ecg_cleaned, sampling_rate=1000, method='vangent2019')

Methods:

• 'vangent2019': Bandpass filtering 0.1-0.4 Hz
• 'charlton2016': Bandpass 0.15-0.4 Hz
• 'soni2019': Bandpass 0.08-0.5 Hz

Use case:

• Estimate respiration when direct respiratory signal unavailable
• Multi-modal physiological analysis

Simulation and Visualization

ecg_simulate()

Generate synthetic ECG signals for testing and validation.

synthetic_ecg = nk.ecg_simulate(duration=10, sampling_rate=1000, heart_rate=70, method='ecgsyn', noise=0.01)

Methods:

• 'ecgsyn': Realistic dynamical model (McSharry et al., 2003)
• Simulates P-QRS-T complex morphology
• Physiologically plausible waveforms
• 'simple': Faster wavelet-based approximation
• Gaussian-like QRS complexes
• Less realistic but computationally efficient

Key parameters:

• heart_rate: Average BPM (default: 70)
• heart_rate_std: Heart rate variability magnitude (default: 1)
• noise: Gaussian noise level (default: 0.01)
• random_state: Seed for reproducibility

ecg_plot()

Visualize processed ECG with detected R-peaks and signal quality.

nk.ecg_plot(signals, info)

Displays:

• Raw and cleaned ECG signals
• Detected R-peaks overlaid
• Heart rate trace
• Signal quality indicators

ECG-Specific Considerations

Sampling Rate Recommendations

• Minimum: 250 Hz for basic R-peak detection
• Recommended: 500-1000 Hz for waveform delineation
• High-resolution: 2000+ Hz for detailed morphology analysis

Signal Duration Requirements

• R-peak detection: Any duration (≥2 beats minimum)
• Basic heart rate: ≥10 seconds
• HRV time domain: ≥60 seconds
• HRV frequency domain: 1-5 minutes (optimal)
• Ultra-low frequency HRV: ≥24 hours

Common Issues and Solutions

Poor R-peak detection:

• Try different methods: method='pantompkins1985' often robust
• Ensure adequate sampling rate (≥250 Hz)
• Check for inverted ECG: use ecg_invert()
• Apply artifact correction: correct_artifacts=True

Noisy signal:

• Use appropriate cleaning method for noise type
• Adjust powerline frequency if outside US/Europe
• Consider signal quality assessment before analysis

Missing waveform components:

• Increase sampling rate (≥500 Hz for delineation)
• Try different delineation methods ('dwt', 'peak', 'cwt')
• Verify signal quality with ecg_quality()

Integration with Other Signals

ECG + RSP: Respiratory Sinus Arrhythmia

# Process both signals
ecg_signals, ecg_info = nk.ecg_process(ecg, sampling_rate=1000)
rsp_signals, rsp_info = nk.rsp_process(rsp, sampling_rate=1000)

# Compute RSA
rsa = nk.hrv_rsa(ecg_info['ECG_R_Peaks'], rsp_signals['RSP_Clean'], sampling_rate=1000)

Multi-modal Integration

# Process multiple signals at once
bio_signals, bio_info = nk.bio_process(
    ecg=ecg_signal,
    rsp=rsp_signal,
    eda=eda_signal,
    sampling_rate=1000
)

References

• Pan, J., & Tompkins, W. J. (1985). A real-time QRS detection algorithm. IEEE transactions on biomedical engineering, 32(3), 230-236.
• Hamilton, P. (2002). Open source ECG analysis. Computers in cardiology, 101-104.
• Martinez, J. P., Almeida, R., Olmos, S., Rocha, A. P., & Laguna, P. (2004). A wavelet-based ECG delineator: evaluation on standard databases. IEEE Transactions on biomedical engineering, 51(4), 570-581.
• Lipponen, J. A., & Tarvainen, M. P. (2019). A robust algorithm for heart rate variability time series artefact correction using novel beat classification. Journal of medical engineering & technology, 43(3), 173-181.

Electrodermal Activity (EDA) Analysis

Overview

Electrodermal Activity (EDA), also known as Galvanic Skin Response (GSR) or Skin Conductance (SC), measures the electrical conductance of the skin, reflecting sympathetic nervous system arousal and sweat gland activity. EDA is widely used in psychophysiology, affective computing, and lie detection.

Main Processing Pipeline

eda_process()

Automated processing of raw EDA signals returning tonic/phasic decomposition and SCR features.

signals, info = nk.eda_process(eda_signal, sampling_rate=100, method='neurokit')

Pipeline steps: 1. Signal cleaning (low-pass filtering) 2. Tonic-phasic decomposition 3. Skin conductance response (SCR) detection 4. SCR feature extraction (onset, peak, amplitude, rise/recovery times)

Returns:

• signals: DataFrame with:
• EDA_Clean: Filtered signal
• EDA_Tonic: Slow-varying baseline
• EDA_Phasic: Fast-varying responses
• SCR_Onsets, SCR_Peaks, SCR_Height: Response markers
• SCR_Amplitude, SCR_RiseTime, SCR_RecoveryTime: Response features
• info: Dictionary with processing parameters

Methods:

• 'neurokit': cvxEDA decomposition + neurokit peak detection
• 'biosppy': Median smoothing + biosppy approach

Preprocessing Functions

eda_clean()

Remove noise through low-pass filtering.

cleaned_eda = nk.eda_clean(eda_signal, sampling_rate=100, method='neurokit')

Methods:

• 'neurokit': Low-pass Butterworth filter (3 Hz cutoff)
• 'biosppy': Low-pass Butterworth filter (5 Hz cutoff)

Automatic skipping:

• If sampling rate < 7 Hz, cleaning is skipped (already low-pass)

Rationale:

• EDA frequency content typically 0-3 Hz
• Remove high-frequency noise and motion artifacts
• Preserve slow SCRs (typical rise time 1-3 seconds)

eda_phasic()

Decompose EDA into tonic (slow baseline) and phasic (rapid responses) components.

tonic, phasic = nk.eda_phasic(eda_cleaned, sampling_rate=100, method='cvxeda')

Methods:

1. cvxEDA (default, recommended):

tonic, phasic = nk.eda_phasic(eda_cleaned, sampling_rate=100, method='cvxeda')

• Convex optimization approach (Greco et al., 2016)
• Sparse phasic driver model
• Most physiologically accurate
• Computationally intensive but superior decomposition

2. Median smoothing:

tonic, phasic = nk.eda_phasic(eda_cleaned, sampling_rate=100, method='smoothmedian')

• Median filter with configurable window
• Fast, simple
• Less accurate than cvxEDA

3. High-pass filtering (Biopac's Acqknowledge):

tonic, phasic = nk.eda_phasic(eda_cleaned, sampling_rate=100, method='highpass')

• High-pass filter (0.05 Hz) extracts phasic
• Fast computation
• Tonic derived by subtraction

4. SparsEDA:

tonic, phasic = nk.eda_phasic(eda_cleaned, sampling_rate=100, method='sparseda')

• Sparse deconvolution approach
• Alternative optimization method

Returns:

• tonic: Slow-varying skin conductance level (SCL)
• phasic: Fast skin conductance responses (SCRs)

Physiological interpretation:

• Tonic (SCL): Baseline arousal, general activation, hydration
• Phasic (SCR): Event-related responses, orienting, emotional reactions

eda_peaks()

Detect Skin Conductance Responses (SCRs) in phasic component.

peaks, info = nk.eda_peaks(eda_phasic, sampling_rate=100, method='neurokit',
                           amplitude_min=0.1)

Methods:

• 'neurokit': Optimized for reliability, configurable thresholds
• 'gamboa2008': Gamboa's algorithm
• 'kim2004': Kim's approach
• 'vanhalem2020': Van Halem's method
• 'nabian2018': Nabian's algorithm

Key parameters:

• amplitude_min: Minimum SCR amplitude (default: 0.1 µS)
• Too low: false positives from noise
• Too high: miss small but valid responses
• rise_time_max: Maximum rise time (default: 2 seconds)
• rise_time_min: Minimum rise time (default: 0.01 seconds)

Returns:

• Dictionary with:
• SCR_Onsets: Indices where SCR begins
• SCR_Peaks: Indices of peak amplitude
• SCR_Height: Peak height above baseline
• SCR_Amplitude: Onset-to-peak amplitude
• SCR_RiseTime: Onset-to-peak duration
• SCR_RecoveryTime: Peak-to-recovery duration (50% decay)

SCR timing conventions:

• Latency: 1-3 seconds after stimulus (typical)
• Rise time: 0.5-3 seconds
• Recovery time: 2-10 seconds (to 50% recovery)
• Minimum amplitude: 0.01-0.05 µS (detection threshold)

eda_fixpeaks()

Correct detected SCR peaks (currently placeholder for EDA).

corrected_peaks = nk.eda_fixpeaks(peaks)

Note: Less critical for EDA than cardiac signals due to slower dynamics.

Analysis Functions

eda_analyze()

Automatically select appropriate analysis type based on data duration.

analysis = nk.eda_analyze(signals, sampling_rate=100)

Mode selection:

• Duration < 10 seconds → eda_eventrelated()
• Duration ≥ 10 seconds → eda_intervalrelated()

Returns:

• DataFrame with EDA metrics appropriate for analysis mode

eda_eventrelated()

Analyze stimulus-locked EDA epochs for event-related responses.

results = nk.eda_eventrelated(epochs)

Computed metrics (per epoch):

• EDA_SCR: Presence of SCR (binary: 0 or 1)
• SCR_Amplitude: Maximum SCR amplitude during epoch
• SCR_Magnitude: Mean phasic activity
• SCR_Peak_Amplitude: Onset-to-peak amplitude
• SCR_RiseTime: Time to peak from onset
• SCR_RecoveryTime: Time to 50% recovery
• SCR_Latency: Delay from stimulus to SCR onset
• EDA_Tonic: Mean tonic level during epoch

Typical parameters:

• Epoch duration: 0-10 seconds post-stimulus
• Baseline: -1 to 0 seconds pre-stimulus
• Expected SCR latency: 1-3 seconds

Use cases:

• Emotional stimulus processing (images, sounds)
• Cognitive load assessment (mental arithmetic)
• Anticipation and prediction error
• Orienting responses

eda_intervalrelated()

Analyze extended EDA recordings for overall arousal and activation patterns.

results = nk.eda_intervalrelated(signals, sampling_rate=100)

Computed metrics:

• SCR_Peaks_N: Number of SCRs detected
• SCR_Peaks_Amplitude_Mean: Average SCR amplitude
• EDA_Tonic_Mean, EDA_Tonic_SD: Tonic level statistics
• EDA_Sympathetic: Sympathetic nervous system index
• EDA_SympatheticN: Normalized sympathetic index
• EDA_Autocorrelation: Temporal structure (lag 4 seconds)
• EDA_Phasic_*: Mean, SD, min, max of phasic component

Recording duration:

• Minimum: 10 seconds
• Recommended: 60+ seconds for stable SCR rate
• Sympathetic index: ≥64 seconds required

Use cases:

• Resting state arousal assessment
• Stress level monitoring
• Baseline sympathetic activity
• Long-term affective state

Specialized Analysis Functions

eda_sympathetic()

Derive sympathetic nervous system activity from frequency band (0.045-0.25 Hz).

sympathetic = nk.eda_sympathetic(signals, sampling_rate=100, method='posada',
                                  show=False)

Methods:

• 'posada': Posada-Quintero method (2016)
• Spectral power in 0.045-0.25 Hz band
• Validated against other autonomic measures
• 'ghiasi': Ghiasi method (2018)
• Alternative frequency-based approach

Requirements:

• Minimum duration: 64 seconds
• Sufficient for frequency resolution in target band

Returns:

• EDA_Sympathetic: Sympathetic index (absolute)
• EDA_SympatheticN: Normalized sympathetic index (0-1)

Interpretation:

• Higher values: increased sympathetic arousal
• Reflects tonic sympathetic activity, not phasic responses
• Complements SCR analysis

Use cases:

• Stress assessment
• Arousal monitoring over time
• Cognitive load measurement
• Complementary to HRV for autonomic balance

eda_autocor()

Compute autocorrelation to assess temporal structure of EDA signal.

autocorr = nk.eda_autocor(eda_phasic, sampling_rate=100, lag=4)

Parameters:

• lag: Time lag in seconds (default: 4 seconds)

Interpretation:

• High autocorrelation: persistent, slowly-varying signal
• Low autocorrelation: rapid, uncorrelated fluctuations
• Reflects temporal regularity of SCRs

Use case:

• Assess signal quality
• Characterize response patterns
• Distinguish sustained vs. transient arousal

eda_changepoints()

Detect abrupt shifts in mean and variance of EDA signal.

changepoints = nk.eda_changepoints(eda_phasic, penalty=10000, show=False)

Method:

• Penalty-based segmentation
• Identifies transitions between states

Parameters:

• penalty: Controls sensitivity (default: 10,000)
• Higher penalty: fewer, more robust changepoints
• Lower penalty: more sensitive to small changes

Returns:

• Indices of detected changepoints
• Optional visualization of segments

Use cases:

• Identify state transitions in continuous monitoring
• Segment data by arousal level
• Detect phase changes in experiments
• Automated epoch definition

Visualization

eda_plot()

Create static or interactive visualizations of processed EDA.

nk.eda_plot(signals, info, static=True)

Displays:

• Raw and cleaned EDA signal
• Tonic and phasic components
• Detected SCR onsets, peaks, and recovery
• Sympathetic index time course (if computed)

Interactive mode (`static=False`):

• Plotly-based interactive exploration
• Zoom, pan, hover for details
• Export to image formats

Simulation and Testing

eda_simulate()

Generate synthetic EDA signals with configurable parameters.

synthetic_eda = nk.eda_simulate(duration=10, sampling_rate=100, scr_number=3,
                                noise=0.01, drift=0.01)

Parameters:

• duration: Signal length in seconds
• sampling_rate: Sampling frequency (Hz)
• scr_number: Number of SCRs to include
• noise: Gaussian noise level
• drift: Slow baseline drift magnitude
• random_state: Seed for reproducibility

Returns:

• Synthetic EDA signal with realistic SCR morphology

Use cases:

• Algorithm testing and validation
• Educational demonstrations
• Method comparison

Practical Considerations

Sampling Rate Recommendations

• Minimum: 10 Hz (adequate for slow SCRs)
• Standard: 20-100 Hz (most commercial systems)
• High-resolution: 1000 Hz (research-grade, oversampled)

Recording Duration

• SCR detection: ≥10 seconds (depends on stimulus)
• Event-related: Typically 10-20 seconds per trial
• Interval-related: ≥60 seconds for stable estimates
• Sympathetic index: ≥64 seconds (frequency resolution)

Electrode Placement

• Standard sites:
• Palmar: distal/middle phalanges (fingers)
• Plantar: sole of foot
• High density: Thenar/hypothenar eminence
• Avoid: Hairy skin, low sweat gland density areas
• Bilateral: Left vs. right hand (typically similar)

Signal Quality Issues

Flat signal (no variation):

• Check electrode contact and gel
• Verify proper placement on sweat gland-rich areas
• Allow 5-10 minute adaptation period

Excessive noise:

• Movement artifacts: minimize participant motion
• Electrical interference: check grounding, shielding
• Thermal effects: control room temperature

Baseline drift:

• Normal: slow changes over minutes
• Excessive: electrode polarization, poor contact
• Solution: use eda_phasic() to separate tonic drift

Non-responders:

• ~5-10% of population have minimal EDA
• Genetic/physiological variation
• Not indicative of equipment failure

Best Practices

Preprocessing workflow:

# 1. Clean signal
cleaned = nk.eda_clean(eda_raw, sampling_rate=100, method='neurokit')

# 2. Decompose tonic/phasic
tonic, phasic = nk.eda_phasic(cleaned, sampling_rate=100, method='cvxeda')

# 3. Detect SCRs
signals, info = nk.eda_peaks(phasic, sampling_rate=100, amplitude_min=0.05)

# 4. Analyze
analysis = nk.eda_analyze(signals, sampling_rate=100)

Event-related workflow:

# 1. Process signal
signals, info = nk.eda_process(eda_raw, sampling_rate=100)

# 2. Find events
events = nk.events_find(trigger_channel, threshold=0.5)

# 3. Create epochs (-1 to 10 seconds around stimulus)
epochs = nk.epochs_create(signals, events, sampling_rate=100,
                          epochs_start=-1, epochs_end=10)

# 4. Event-related analysis
results = nk.eda_eventrelated(epochs)

# 5. Statistical analysis
# Compare SCR amplitude across conditions

Clinical and Research Applications

Emotion and affective science:

• Arousal dimension of emotion (not valence)
• Emotional picture viewing
• Music-induced emotion
• Fear conditioning

Cognitive processes:

• Mental workload and effort
• Attention and vigilance
• Decision-making and uncertainty
• Error processing

Clinical populations:

• Anxiety disorders: heightened baseline, exaggerated responses
• PTSD: fear conditioning, extinction deficits
• Autism: atypical arousal patterns
• Psychopathy: reduced fear responses

Applied settings:

• Lie detection (polygraph)
• User experience research
• Driver monitoring
• Stress assessment in real-world settings

Neuroimaging integration:

• fMRI: EDA correlates with amygdala, insula activity
• Concurrent recording during brain imaging
• Autonomic-brain coupling

Interpretation Guidelines

SCR amplitude:

• 0.01-0.05 µS: Small but detectable
• 0.05-0.2 µS: Moderate response
• >0.2 µS: Large response
• Context-dependent: Normalize within-subject

SCR frequency:

• Resting: 1-3 SCRs per minute (typical)
• Stressed: >5 SCRs per minute
• Non-specific SCRs: Spontaneous (no identifiable stimulus)

Tonic SCL:

• Range: 2-20 µS (highly variable across individuals)
• Within-subject changes more interpretable than absolute levels
• Increases: arousal, stress, cognitive load
• Decreases: relaxation, habituation

References

• Boucsein, W. (2012). Electrodermal activity (2nd ed.). Springer Science & Business Media.
• Greco, A., Valenza, G., & Scilingo, E. P. (2016). cvxEDA: A convex optimization approach to electrodermal activity processing. IEEE Transactions on Biomedical Engineering, 63(4), 797-804.
• Posada-Quintero, H. F., Florian, J. P., Orjuela-Cañón, A. D., Aljama-Corrales, T., Charleston-Villalobos, S., & Chon, K. H. (2016). Power spectral density analysis of electrodermal activity for sympathetic function assessment. Annals of biomedical engineering, 44(10), 3124-3135.
• Dawson, M. E., Schell, A. M., & Filion, D. L. (2017). The electrodermal system. In Handbook of psychophysiology (pp. 217-243). Cambridge University Press.

EEG Analysis and Microstates

Overview

Analyze electroencephalography (EEG) signals for frequency band power, channel quality assessment, source localization, and microstate identification. NeuroKit2 integrates with MNE-Python for comprehensive EEG processing workflows.

Core EEG Functions

eeg_power()

Compute power across standard frequency bands for specified channels.

power = nk.eeg_power(eeg_data, sampling_rate=250, channels=['Fz', 'Cz', 'Pz'],
                     frequency_bands={'Delta': (0.5, 4),
                                     'Theta': (4, 8),
                                     'Alpha': (8, 13),
                                     'Beta': (13, 30),
                                     'Gamma': (30, 45)})

Standard frequency bands:

• Delta (0.5-4 Hz): Deep sleep, unconscious processes
• Theta (4-8 Hz): Drowsiness, meditation, memory encoding
• Alpha (8-13 Hz): Relaxed wakefulness, eyes closed
• Beta (13-30 Hz): Active thinking, focus, anxiety
• Gamma (30-45 Hz): Cognitive processing, binding

Returns:

• DataFrame with power values for each channel × frequency band combination
• Columns: Channel_Band (e.g., 'Fz_Alpha', 'Cz_Beta')

Use cases:

• Resting state analysis
• Cognitive state classification
• Sleep staging
• Meditation or neurofeedback monitoring

eeg_badchannels()

Identify problematic channels using statistical outlier detection.

bad_channels = nk.eeg_badchannels(eeg_data, sampling_rate=250, bad_threshold=2)

Detection methods:

• Standard deviation outliers across channels
• Correlation with other channels
• Flat or dead channels
• Channels with excessive noise

Parameters:

• bad_threshold: Z-score threshold for outlier detection (default: 2)

Returns:

• List of channel names identified as problematic

Use case:

• Quality control before analysis
• Automatic bad channel rejection
• Interpolation or exclusion decisions

eeg_rereference()

Re-express voltage measurements relative to different reference points.

rereferenced = nk.eeg_rereference(eeg_data, reference='average', robust=False)

Reference types:

• 'average': Average reference (mean of all electrodes)
• 'REST': Reference Electrode Standardization Technique
• 'bipolar': Differential recording between electrode pairs
• Specific channel name: Use single electrode as reference

Common references:

• Average reference: Most common for high-density EEG
• Linked mastoids: Traditional clinical EEG
• Vertex (Cz): Sometimes used in ERP research
• REST: Approximates infinity reference

Returns:

• Re-referenced EEG data

eeg_gfp()

Compute Global Field Power - the standard deviation of all electrodes at each time point.

gfp = nk.eeg_gfp(eeg_data)

Interpretation:

• High GFP: Strong, synchronized brain activity across regions
• Low GFP: Weak or desynchronized activity
• GFP peaks: Points of stable topography, used for microstate detection

Use cases:

• Identify periods of stable topographic patterns
• Select time points for microstate analysis
• Event-related potential (ERP) visualization

eeg_diss()

Measure topographic dissimilarity between electric field configurations.

dissimilarity = nk.eeg_diss(eeg_data1, eeg_data2, method='gfp')

Methods:

• GFP-based: Normalized difference
• Spatial correlation
• Cosine distance

Use case:

• Compare topographies between conditions
• Microstate transition analysis
• Template matching

Source Localization

eeg_source()

Perform source reconstruction to estimate brain-level activity from scalp recordings.

sources = nk.eeg_source(eeg_data, method='sLORETA')

Methods:

• 'sLORETA': Standardized Low-Resolution Electromagnetic Tomography
• Zero localization error for point sources
• Good spatial resolution
• 'MNE': Minimum Norm Estimate
• Fast, well-established
• Bias toward superficial sources
• 'dSPM': Dynamic Statistical Parametric Mapping
• Normalized MNE
• 'eLORETA': Exact LORETA
• Improved localization accuracy

Requirements:

• Forward model (lead field matrix)
• Co-registered electrode positions
• Head model (boundary element or spherical)

Returns:

• Source space activity estimates

eeg_source_extract()

Extract activity from specific anatomical brain regions.

regional_activity = nk.eeg_source_extract(sources, regions=['PFC', 'MTL', 'Parietal'])

Region options:

• Standard atlases: Desikan-Killiany, Destrieux, AAL
• Custom ROIs
• Brodmann areas

Returns:

• Time series for each region
• Averaged or principal component across voxels

Use cases:

• Region-of-interest analysis
• Functional connectivity
• Source-level statistics

Microstate Analysis

Microstates are brief (80-120 ms) periods of stable brain topography, representing coordinated neural networks. Typically 4-7 microstate classes (often labeled A, B, C, D) with distinct functions.

microstates_segment()

Identify and extract microstates using clustering algorithms.

microstates = nk.microstates_segment(eeg_data, n_microstates=4, sampling_rate=250,
                                      method='kmod', normalize=True)

Methods:

• 'kmod' (default): Modified k-means optimized for EEG topographies
• Polarity-invariant clustering
• Most common in microstate literature
• 'kmeans': Standard k-means clustering
• 'kmedoids': K-medoids (more robust to outliers)
• 'pca': Principal component analysis
• 'ica': Independent component analysis
• 'aahc': Atomize and agglomerate hierarchical clustering

Parameters:

• n_microstates: Number of microstate classes (typically 4-7)
• normalize: Normalize topographies (recommended: True)
• n_inits: Number of random initializations (increase for stability)

Returns:

• Dictionary with:
• 'maps': Microstate template topographies
• 'labels': Microstate label at each time point
• 'gfp': Global field power
• 'gev': Global explained variance

microstates_findnumber()

Estimate the optimal number of microstates.

optimal_k = nk.microstates_findnumber(eeg_data, show=True)

Criteria:

• Global Explained Variance (GEV): Percentage of variance explained
• Elbow method: find "knee" in GEV curve
• Typically 70-80% GEV achieved
• Krzanowski-Lai (KL) Criterion: Statistical measure balancing fit and parsimony
• Maximum KL indicates optimal k

Typical range: 4-7 microstates

• 4 microstates: Classic A, B, C, D states
• 5-7 microstates: Finer-grained decomposition

microstates_classify()

Reorder microstates based on anterior-posterior and left-right channel values.

classified = nk.microstates_classify(microstates)

Purpose:

• Standardize microstate labels across subjects
• Match conventional A, B, C, D topographies:
• A: Left-right orientation, parieto-occipital
• B: Right-left orientation, fronto-temporal
• C: Anterior-posterior orientation, frontal-central
• D: Fronto-central, anterior-posterior (inverse of C)

Returns:

• Reordered microstate maps and labels

microstates_clean()

Preprocess EEG data for microstate extraction.

cleaned_eeg = nk.microstates_clean(eeg_data, sampling_rate=250)

Preprocessing steps:

• Bandpass filtering (typically 2-20 Hz)
• Artifact rejection
• Bad channel interpolation
• Re-referencing to average

Rationale:

• Microstates reflect large-scale network activity
• High-frequency and low-frequency artifacts can distort topographies

microstates_peaks()

Identify GFP peaks for microstate analysis.

peak_indices = nk.microstates_peaks(eeg_data, sampling_rate=250)

Purpose:

• Microstates typically analyzed at GFP peaks
• Peaks represent moments of maximal, stable topographic activity
• Reduces computational load and noise sensitivity

Returns:

• Indices of GFP local maxima

microstates_static()

Compute temporal properties of individual microstates.

static_metrics = nk.microstates_static(microstates)

Metrics:

• Duration (ms): Mean time spent in each microstate
• Typical: 80-120 ms
• Reflects stability and persistence
• Occurrence (per second): Frequency of microstate appearances
• How often each state is entered
• Coverage (%): Percentage of total time in each microstate
• Relative dominance
• Global Explained Variance (GEV): Variance explained by each class
• Quality of template fit

Returns:

• DataFrame with metrics for each microstate class

Interpretation:

• Changes in duration: altered network stability
• Changes in occurrence: shifting state dynamics
• Changes in coverage: dominance of specific networks

microstates_dynamic()

Analyze transition patterns between microstates.

dynamic_metrics = nk.microstates_dynamic(microstates)

Metrics:

• Transition matrix: Probability of transitioning from state i to state j
• Reveals preferential sequences
• Transition rate: Overall transition frequency
• Higher rate: more rapid switching
• Entropy: Randomness of transitions
• High entropy: unpredictable switching
• Low entropy: stereotyped sequences
• Markov test: Are transitions history-dependent?

Returns:

• Dictionary with transition statistics

Use cases:

• Identify abnormal microstate sequences in clinical populations
• Network dynamics and flexibility
• State-dependent information processing

microstates_plot()

Visualize microstate topographies and time course.

nk.microstates_plot(microstates, eeg_data)

Displays:

• Topographic maps for each microstate class
• GFP trace with microstate labels
• Transition plot showing state sequences
• Statistical summary

MNE Integration Utilities

mne_data()

Access sample datasets from MNE-Python.

raw = nk.mne_data(dataset='sample', directory=None)

Available datasets:

• 'sample': Multi-modal (MEG/EEG) example
• 'ssvep': Steady-state visual evoked potentials
• 'eegbci': Motor imagery BCI dataset

mne_to_df() / mne_to_dict()

Convert MNE objects to NeuroKit-compatible formats.

df = nk.mne_to_df(raw)
data_dict = nk.mne_to_dict(epochs)

Use case:

• Work with MNE-processed data in NeuroKit2
• Convert between formats for analysis

mne_channel_add() / mne_channel_extract()

Manage individual channels in MNE objects.

# Extract specific channels
subset = nk.mne_channel_extract(raw, ['Fz', 'Cz', 'Pz'])

# Add derived channels
raw_with_eog = nk.mne_channel_add(raw, new_channel_data, ch_name='EOG')

mne_crop()

Trim recordings by time or samples.

cropped = nk.mne_crop(raw, tmin=10, tmax=100)

mne_templateMRI()

Provide template anatomy for source localization.

subjects_dir = nk.mne_templateMRI()

Use case:

• Source analysis without individual MRI
• Group-level source localization
• fsaverage template brain

eeg_simulate()

Generate synthetic EEG signals for testing.

synthetic_eeg = nk.eeg_simulate(duration=60, sampling_rate=250, n_channels=32)

Practical Considerations

Sampling Rate Recommendations

• Minimum: 100 Hz for basic power analysis
• Standard: 250-500 Hz for most applications
• High-resolution: 1000+ Hz for detailed temporal dynamics

Recording Duration

• Power analysis: ≥2 minutes for stable estimates
• Microstates: ≥2-5 minutes, longer preferred
• Resting state: 3-10 minutes typical
• Event-related: Depends on trial count (≥30 trials per condition)

Artifact Management

• Eye blinks: Remove with ICA or regression
• Muscle artifacts: High-pass filter (≥1 Hz) or manual rejection
• Bad channels: Detect and interpolate before analysis
• Line noise: Notch filter at 50/60 Hz

Best Practices

Power analysis:

# 1. Clean data
cleaned = nk.signal_filter(eeg_data, sampling_rate=250, lowcut=0.5, highcut=45)

# 2. Identify and interpolate bad channels
bad = nk.eeg_badchannels(cleaned, sampling_rate=250)
# Interpolate bad channels using MNE

# 3. Re-reference
rereferenced = nk.eeg_rereference(cleaned, reference='average')

# 4. Compute power
power = nk.eeg_power(rereferenced, sampling_rate=250, channels=channel_list)

Microstate workflow:

# 1. Preprocess
cleaned = nk.microstates_clean(eeg_data, sampling_rate=250)

# 2. Determine optimal number of states
optimal_k = nk.microstates_findnumber(cleaned, show=True)

# 3. Segment microstates
microstates = nk.microstates_segment(cleaned, n_microstates=optimal_k,
                                     sampling_rate=250, method='kmod')

# 4. Classify to standard labels
microstates = nk.microstates_classify(microstates)

# 5. Compute temporal metrics
static = nk.microstates_static(microstates)
dynamic = nk.microstates_dynamic(microstates)

# 6. Visualize
nk.microstates_plot(microstates, cleaned)

Clinical and Research Applications

Cognitive neuroscience:

• Attention, working memory, executive function
• Language processing
• Sensory perception

Clinical populations:

• Epilepsy: seizure detection, localization
• Alzheimer's disease: slowing of EEG, microstate alterations
• Schizophrenia: altered microstates, especially state C
• ADHD: increased theta/beta ratio
• Depression: frontal alpha asymmetry

Consciousness research:

• Anesthesia monitoring
• Disorders of consciousness
• Sleep staging

Neurofeedback:

• Real-time frequency band training
• Alpha enhancement for relaxation
• Beta enhancement for focus

References

• Michel, C. M., & Koenig, T. (2018). EEG microstates as a tool for studying the temporal dynamics of whole-brain neuronal networks: A review. Neuroimage, 180, 577-593.
• Pascual-Marqui, R. D., Michel, C. M., & Lehmann, D. (1995). Segmentation of brain electrical activity into microstates: model estimation and validation. IEEE Transactions on Biomedical Engineering, 42(7), 658-665.
• Gramfort, A., Luessi, M., Larson, E., Engemann, D. A., Strohmeier, D., Brodbeck, C., ... & Hämäläinen, M. (2013). MEG and EEG data analysis with MNE-Python. Frontiers in neuroscience, 7, 267.

Electromyography (EMG) Analysis

Overview

Electromyography (EMG) measures electrical activity produced by skeletal muscles during contraction. EMG analysis in NeuroKit2 focuses on amplitude estimation, muscle activation detection, and temporal dynamics for psychophysiology and motor control research.

Main Processing Pipeline

emg_process()

Automated EMG signal processing pipeline.

signals, info = nk.emg_process(emg_signal, sampling_rate=1000)

Pipeline steps: 1. Signal cleaning (high-pass filtering, detrending) 2. Amplitude envelope extraction 3. Muscle activation detection 4. Onset and offset identification

Returns:

• signals: DataFrame with:
• EMG_Clean: Filtered EMG signal
• EMG_Amplitude: Linear envelope (smoothed rectified signal)
• EMG_Activity: Binary activation indicator (0/1)
• EMG_Onsets: Activation onset markers
• EMG_Offsets: Activation offset markers
• info: Dictionary with activation parameters

Typical workflow:

• Process raw EMG → Extract amplitude → Detect activations → Analyze features

Preprocessing Functions

emg_clean()

Apply filtering to remove noise and prepare for amplitude extraction.

cleaned_emg = nk.emg_clean(emg_signal, sampling_rate=1000)

Filtering approach (BioSPPy method):

• Fourth-order Butterworth high-pass filter (100 Hz)
• Removes low-frequency movement artifacts and baseline drift
• Removes DC offset
• Signal detrending

Rationale:

• EMG frequency content: 20-500 Hz (dominant: 50-150 Hz)
• High-pass at 100 Hz isolates muscle activity
• Removes ECG contamination (especially in trunk muscles)
• Removes motion artifacts (<20 Hz)

EMG signal characteristics:

• Random, zero-mean oscillations during contraction
• Higher amplitude = stronger contraction
• Raw EMG: both positive and negative deflections

Feature Extraction

emg_amplitude()

Compute linear envelope representing muscle contraction intensity.

amplitude = nk.emg_amplitude(cleaned_emg, sampling_rate=1000)

Method: 1. Full-wave rectification (absolute value) 2. Low-pass filtering (smooth envelope) 3. Downsampling (optional)

Linear envelope:

• Smooth curve following EMG amplitude modulation
• Represents muscle force/activation level
• Suitable for further analysis (activation detection, integration)

Typical smoothing:

• Low-pass filter: 10-20 Hz cutoff
• Moving average: 50-200 ms window
• Balance: responsiveness vs. smoothness

Activation Detection

emg_activation()

Detect periods of muscle activation (onsets and offsets).

activity, info = nk.emg_activation(emg_amplitude, sampling_rate=1000, method='threshold',
                                   threshold='auto', duration_min=0.05)

Methods:

1. Threshold-based (default):

activity = nk.emg_activation(amplitude, method='threshold', threshold='auto')

• Compares amplitude to threshold
• threshold='auto': Automatic based on signal statistics (e.g., mean + 1 SD)
• threshold=0.1: Manual absolute threshold
• Simple, fast, widely used

2. Gaussian Mixture Model (GMM):

activity = nk.emg_activation(amplitude, method='mixture', n_clusters=2)

• Unsupervised clustering: active vs. rest
• Adaptive to signal characteristics
• More robust to varying baseline

3. Changepoint detection:

activity = nk.emg_activation(amplitude, method='changepoint')

• Detects abrupt transitions in signal properties
• Identifies activation/deactivation points
• Useful for complex temporal patterns

4. Bimodality (Silva et al., 2013):

activity = nk.emg_activation(amplitude, method='bimodal')

• Tests for bimodal distribution (active vs. rest)
• Determines optimal separation threshold
• Statistically principled

Key parameters:

• duration_min: Minimum activation duration (seconds)
• Filters brief spurious activations
• Typical: 50-100 ms
• threshold: Activation threshold (method-dependent)

Returns:

• activity: Binary array (0 = rest, 1 = active)
• info: Dictionary with onset/offset indices

Activation metrics:

• Onset: Transition from rest to activity
• Offset: Transition from activity to rest
• Duration: Time between onset and offset
• Burst: Single period of continuous activation

Analysis Functions

emg_analyze()

Automatically select event-related or interval-related analysis.

analysis = nk.emg_analyze(signals, sampling_rate=1000)

Mode selection:

• Duration < 10 seconds → event-related
• Duration ≥ 10 seconds → interval-related

emg_eventrelated()

Analyze EMG responses to discrete events/stimuli.

results = nk.emg_eventrelated(epochs)

Computed metrics (per epoch):

• EMG_Activation: Presence of activation (binary)
• EMG_Amplitude_Mean: Average amplitude during epoch
• EMG_Amplitude_Max: Peak amplitude
• EMG_Bursts: Number of activation bursts
• EMG_Onset_Latency: Time from event to first activation (if applicable)

Use cases:

• Startle response (orbicularis oculi EMG)
• Facial EMG during emotional stimuli (corrugator, zygomaticus)
• Motor response latencies
• Muscle reactivity paradigms

emg_intervalrelated()

Analyze extended EMG recordings.

results = nk.emg_intervalrelated(signals, sampling_rate=1000)

Computed metrics:

• EMG_Bursts_N: Total number of activation bursts
• EMG_Amplitude_Mean: Mean amplitude across entire interval
• EMG_Activation_Duration: Total time in active state
• EMG_Rest_Duration: Total time in rest state

Use cases:

• Resting muscle tension assessment
• Chronic pain or stress-related muscle activity
• Fatigue monitoring during sustained tasks
• Postural muscle assessment

Simulation and Visualization

emg_simulate()

Generate synthetic EMG signals for testing.

synthetic_emg = nk.emg_simulate(duration=10, sampling_rate=1000, burst_number=3,
                                noise=0.1, random_state=42)

Parameters:

• burst_number: Number of activation bursts to include
• noise: Background noise level
• random_state: Reproducibility seed

Generated features:

• Random EMG-like oscillations during bursts
• Realistic frequency content
• Variable burst timing and amplitude

Use cases:

• Algorithm validation
• Detection parameter tuning
• Educational demonstrations

emg_plot()

Visualize processed EMG signal.

nk.emg_plot(signals, info, static=True)

Displays:

• Raw and cleaned EMG signal
• Amplitude envelope
• Detected activation periods
• Onset/offset markers

Interactive mode: Set static=False for Plotly visualization

Practical Considerations

Sampling Rate Recommendations

• Minimum: 500 Hz (Nyquist for 250 Hz upper frequency)
• Standard: 1000 Hz (most research applications)
• High-resolution: 2000-4000 Hz (detailed motor unit studies)
• Surface EMG: 1000-2000 Hz typical
• Intramuscular EMG: 10,000+ Hz for single motor units

Recording Duration

• Event-related: Depends on paradigm (e.g., 2-5 seconds per trial)
• Sustained contraction: Seconds to minutes
• Fatigue studies: Minutes to hours
• Chronic monitoring: Days (wearable EMG)

Electrode Placement

Surface EMG (most common):

• Bipolar configuration (two electrodes over muscle belly)
• Reference/ground electrode over electrically neutral site (bone)
• Skin preparation: clean, abrade, reduce impedance
• Inter-electrode distance: 10-20 mm (SENIAM standards)

Muscle-specific guidelines:

• Follow SENIAM (Surface EMG for Non-Invasive Assessment of Muscles) recommendations
• Palpate muscle during contraction to locate belly
• Align electrodes with muscle fiber direction

Common muscles in psychophysiology:

• Corrugator supercilii: Frowning, negative affect (above eyebrow)
• Zygomaticus major: Smiling, positive affect (cheek)
• Orbicularis oculi: Startle response, fear (around eye)
• Masseter: Jaw clenching, stress (jaw muscle)
• Trapezius: Shoulder tension, stress (upper back)
• Frontalis: Forehead tension, surprise

Signal Quality Issues

ECG contamination:

• Common in trunk and proximal muscles
• High-pass filtering (>100 Hz) usually sufficient
• If persistent: template subtraction, ICA

Motion artifacts:

• Low-frequency disturbances
• Electrode cable movement
• Secure electrodes, minimize cable motion

Electrode issues:

• Poor contact: high impedance, low amplitude
• Sweat: gradual amplitude increase, instability
• Hair: clean or shave area

Cross-talk:

• Adjacent muscle activity bleeding into recording
• Careful electrode placement
• Small inter-electrode distance

Best Practices

Standard workflow:

# 1. Clean signal (high-pass filter, detrend)
cleaned = nk.emg_clean(emg_raw, sampling_rate=1000)

# 2. Extract amplitude envelope
amplitude = nk.emg_amplitude(cleaned, sampling_rate=1000)

# 3. Detect activation periods
activity, info = nk.emg_activation(amplitude, sampling_rate=1000,
                                   method='threshold', threshold='auto')

# 4. Comprehensive processing (alternative)
signals, info = nk.emg_process(emg_raw, sampling_rate=1000)

# 5. Analyze
analysis = nk.emg_analyze(signals, sampling_rate=1000)

Normalization:

# Maximum voluntary contraction (MVC) normalization
mvc_amplitude = np.max(mvc_emg_amplitude)  # From separate MVC trial
normalized_emg = (amplitude / mvc_amplitude) * 100  # Express as % MVC

# Common in ergonomics, exercise physiology
# Allows comparison across individuals and sessions

Clinical and Research Applications

Psychophysiology:

• Facial EMG: Emotional valence (smile vs. frown)
• Startle response: Fear, surprise, defensive reactivity
• Stress: Chronic muscle tension (trapezius, masseter)

Motor control and rehabilitation:

• Gait analysis
• Movement disorders (tremor, dystonia)
• Stroke rehabilitation (muscle re-activation)
• Prosthetic control (myoelectric)

Ergonomics and occupational health:

• Work-related musculoskeletal disorders
• Postural assessment
• Repetitive strain injury risk

Sports science:

• Muscle activation patterns during exercise
• Fatigue assessment (median frequency shift)
• Training optimization

Biofeedback:

• Relaxation training (reduce muscle tension)
• Neuromuscular re-education
• Chronic pain management

Sleep medicine:

• Chin EMG for REM sleep atonia
• Periodic limb movements
• Bruxism (teeth grinding)

Advanced EMG Analysis (Beyond NeuroKit2 Basic Functions)

Frequency domain:

• Median frequency shift during fatigue
• Power spectrum analysis
• Requires longer segments (≥1 second per analysis window)

Motor unit identification:

• Intramuscular EMG
• Spike detection and sorting
• Firing rate analysis
• Requires high sampling rates (10+ kHz)

Muscle coordination:

• Co-contraction indices
• Synergy analysis
• Multi-muscle integration

Interpretation Guidelines

Amplitude (linear envelope):

• Higher amplitude ≈ stronger contraction (not perfectly linear)
• Relationship to force: sigmoid, influenced by many factors
• Within-subject comparisons most reliable

Activation threshold:

• Automatic thresholds: convenient but verify visually
• Manual thresholds: may be needed for non-standard muscles
• Resting baseline: should be near zero (if not, check electrodes)

Burst characteristics:

• Phasic: Brief bursts (startle, rapid movements)
• Tonic: Sustained activation (postural, sustained grip)
• Rhythmic: Repeated bursts (tremor, walking)

References

• Fridlund, A. J., & Cacioppo, J. T. (1986). Guidelines for human electromyographic research. Psychophysiology, 23(5), 567-589.
• Hermens, H. J., Freriks, B., Disselhorst-Klug, C., & Rau, G. (2000). Development of recommendations for SEMG sensors and sensor placement procedures. Journal of electromyography and Kinesiology, 10(5), 361-374.
• Silva, H., Scherer, R., Sousa, J., & Londral, A. (2013). Towards improving the ssability of electromyographic interfaces. Journal of Oral Rehabilitation, 40(6), 456-465.
• Tassinary, L. G., Cacioppo, J. T., & Vanman, E. J. (2017). The skeletomotor system: Surface electromyography. In Handbook of psychophysiology (pp. 267-299). Cambridge University Press.

Electrooculography (EOG) Analysis

Overview

Electrooculography (EOG) measures eye movements and blinks by detecting electrical potential differences generated by eye position changes. EOG is used in sleep studies, attention research, reading analysis, and artifact correction for EEG.

Main Processing Pipeline

eog_process()

Automated EOG signal processing pipeline.

signals, info = nk.eog_process(eog_signal, sampling_rate=500, method='neurokit')

Pipeline steps: 1. Signal cleaning (filtering) 2. Blink detection 3. Blink rate calculation

Returns:

• signals: DataFrame with:
• EOG_Clean: Filtered EOG signal
• EOG_Blinks: Binary blink markers (0/1)
• EOG_Rate: Instantaneous blink rate (blinks/min)
• info: Dictionary with blink indices and parameters

Methods:

• 'neurokit': NeuroKit2 optimized approach (default)
• 'agarwal2019': Agarwal et al. (2019) algorithm
• 'mne': MNE-Python method
• 'brainstorm': Brainstorm toolbox approach
• 'kong1998': Kong et al. (1998) method

Preprocessing Functions

eog_clean()

Prepare raw EOG signal for blink detection.

cleaned_eog = nk.eog_clean(eog_signal, sampling_rate=500, method='neurokit')

Methods:

• 'neurokit': Butterworth filtering optimized for EOG
• 'agarwal2019': Alternative filtering
• 'mne': MNE-Python preprocessing
• 'brainstorm': Brainstorm approach
• 'kong1998': Kong's method

Typical filtering:

• Low-pass: 10-20 Hz (remove high-frequency noise)
• High-pass: 0.1-1 Hz (remove DC drift)
• Preserves blink waveform (typical duration 100-400 ms)

EOG signal characteristics:

• Blinks: Large amplitude, stereotyped waveform (200-400 ms)
• Saccades: Rapid step-like deflections (20-80 ms)
• Smooth pursuit: Slow ramp-like changes
• Baseline: Stable when eyes fixated

Blink Detection

eog_peaks()

Detect eye blinks in EOG signal.

blinks, info = nk.eog_peaks(cleaned_eog, sampling_rate=500, method='neurokit',
                            threshold=0.33)

Methods:

• 'neurokit': Amplitude and duration criteria (default)
• 'mne': MNE-Python blink detection
• 'brainstorm': Brainstorm approach
• 'blinker': BLINKER algorithm (Kleifges et al., 2017)

Key parameters:

• threshold: Amplitude threshold (fraction of max amplitude)
• Typical: 0.2-0.5
• Lower: more sensitive (may include false positives)
• Higher: more conservative (may miss small blinks)

Returns:

• Dictionary with 'EOG_Blinks' key containing blink peak indices

Blink characteristics:

• Frequency: 15-20 blinks/min (resting, comfortable)
• Duration: 100-400 ms (mean ~200 ms)
• Amplitude: Varies with electrode placement and individual factors
• Waveform: Biphasic or triphasic

eog_findpeaks()

Low-level blink detection with multiple algorithms.

blinks_dict = nk.eog_findpeaks(cleaned_eog, sampling_rate=500, method='neurokit')

Use cases:

• Custom parameter tuning
• Algorithm comparison
• Research method development

Feature Extraction

eog_features()

Extract characteristics of individual blinks.

features = nk.eog_features(signals, sampling_rate=500)

Computed features:

• Amplitude velocity ratio (AVR): Peak velocity / amplitude
• Discriminates blinks from artifacts
• Blink-amplitude ratio: Consistency of blink amplitudes
• Duration metrics: Blink duration statistics (mean, SD)
• Peak amplitude: Maximum deflection
• Peak velocity: Maximum rate of change

Use cases:

• Blink quality assessment
• Drowsiness detection (blink duration increases when sleepy)
• Neurological assessment (altered blink dynamics in disease)

eog_rate()

Compute blink frequency (blinks per minute).

blink_rate = nk.eog_rate(blinks, sampling_rate=500, desired_length=None)

Method:

• Calculate inter-blink intervals
• Convert to blinks per minute
• Interpolate to match signal length

Typical blink rates:

• Resting: 15-20 blinks/min
• Reading/visual tasks: 5-10 blinks/min (suppressed)
• Conversation: 20-30 blinks/min
• Stress/dry eyes: >30 blinks/min
• Drowsiness: Variable, longer blinks

Analysis Functions

eog_analyze()

Automatically select event-related or interval-related analysis.

analysis = nk.eog_analyze(signals, sampling_rate=500)

Mode selection:

• Duration < 10 seconds → event-related
• Duration ≥ 10 seconds → interval-related

eog_eventrelated()

Analyze blink patterns relative to specific events.

results = nk.eog_eventrelated(epochs)

Computed metrics (per epoch):

• EOG_Blinks_N: Number of blinks during epoch
• EOG_Rate_Mean: Average blink rate
• EOG_Blink_Presence: Binary (any blinks occurred)
• Temporal distribution of blinks across epoch

Use cases:

• Blink-locked ERP contamination assessment
• Attention and engagement during stimuli
• Visual task difficulty (suppressed blinks during demanding tasks)
• Spontaneous blinks after stimulus offset

eog_intervalrelated()

Analyze blink patterns over extended periods.

results = nk.eog_intervalrelated(signals, sampling_rate=500)

Computed metrics:

• EOG_Blinks_N: Total number of blinks
• EOG_Rate_Mean: Average blink rate (blinks/min)
• EOG_Rate_SD: Blink rate variability
• EOG_Duration_Mean: Average blink duration (if available)
• EOG_Amplitude_Mean: Average blink amplitude (if available)

Use cases:

• Resting state blink patterns
• Drowsiness or fatigue monitoring (increased duration)
• Sustained attention tasks (suppressed rate)
• Dry eye assessment (increased rate, incomplete blinks)

Simulation and Visualization

eog_plot()

Visualize processed EOG signal and detected blinks.

nk.eog_plot(signals, info)

Displays:

• Raw and cleaned EOG signal
• Detected blink markers
• Blink rate time course

Practical Considerations

Sampling Rate Recommendations

• Minimum: 100 Hz (basic blink detection)
• Standard: 250-500 Hz (research applications)
• High-resolution: 1000 Hz (detailed waveform analysis, saccades)
• Sleep studies: 200-250 Hz typical

Recording Duration

• Blink detection: Any duration (≥1 blink)
• Blink rate estimation: ≥60 seconds for stable estimate
• Event-related: Depends on paradigm (seconds per trial)
• Sleep EOG: Hours (full night)

Electrode Placement

Standard configurations:

Horizontal EOG (HEOG):

• Two electrodes: lateral canthi (outer corners) of left and right eyes
• Measures horizontal eye movements (saccades, smooth pursuit)
• Bipolar recording (left - right)

Vertical EOG (VEOG):

• Two electrodes: above and below one eye (typically right)
• Measures vertical eye movements and blinks
• Bipolar recording (above - below)

Sleep EOG:

• Often uses slightly different placement (temple area)
• E1: 1 cm lateral and 1 cm below outer canthus of left eye
• E2: 1 cm lateral and 1 cm above outer canthus of right eye
• Captures both horizontal and vertical movements

EEG contamination removal:

• Frontal electrodes (Fp1, Fp2) can serve as EOG proxies
• ICA-based EOG artifact removal common in EEG preprocessing

Common Issues and Solutions

Electrode issues:

• Poor contact: low amplitude, noise
• Skin preparation: clean, light abrasion
• Conductive gel: ensure good contact

Artifacts:

• Muscle activity (especially frontalis): high-frequency noise
• Movement: cable artifacts, head motion
• Electrical noise: 50/60 Hz hum (ground properly)

Saturation:

• Large saccades may saturate amplifier
• Adjust gain or voltage range
• More common with low-resolution systems

Best Practices

Standard workflow:

# 1. Clean signal
cleaned = nk.eog_clean(eog_raw, sampling_rate=500, method='neurokit')

# 2. Detect blinks
blinks, info = nk.eog_peaks(cleaned, sampling_rate=500, method='neurokit')

# 3. Extract features
features = nk.eog_features(signals, sampling_rate=500)

# 4. Comprehensive processing (alternative)
signals, info = nk.eog_process(eog_raw, sampling_rate=500)

# 5. Analyze
analysis = nk.eog_analyze(signals, sampling_rate=500)

EEG artifact correction workflow:

# Option 1: Regression-based removal
# Identify EOG components from cleaned EOG signal
# Regress out EOG from EEG channels

# Option 2: ICA-based removal (preferred)
# 1. Run ICA on EEG data including EOG channels
# 2. Identify ICA components correlated with EOG
# 3. Remove EOG components from EEG data
# NeuroKit2 integrates with MNE for this workflow

Clinical and Research Applications

EEG artifact correction:

• Blinks contaminate frontal EEG channels
• ICA or regression methods remove EOG artifacts
• Essential for ERP studies

Sleep staging:

• Rapid eye movements (REMs) during REM sleep
• Slow rolling eye movements during drowsiness
• Sleep onset and stage transitions

Attention and cognitive load:

• Blink rate suppressed during demanding tasks
• Blinks cluster at task boundaries (natural breakpoints)
• Spontaneous blink as indicator of attention shifts

Fatigue and drowsiness monitoring:

• Increased blink duration when sleepy
• Slower eyelid closures
• Partial or incomplete blinks
• Driver monitoring applications

Reading and visual processing:

• Blinks suppressed during reading
• Eye movements during saccades (line changes)
• Fatigue effects on reading efficiency

Neurological disorders:

• Parkinson's disease: Reduced spontaneous blink rate
• Schizophrenia: Increased blink rate
• Tourette syndrome: Excessive blinking (tics)
• Dry eye syndrome: Increased, incomplete blinks

Affective and social cognition:

• Blink synchrony in social interaction
• Emotional modulation of blink rate
• Blink-related potentials in ERPs

Human-computer interaction:

• Gaze tracking preprocessing
• Attention monitoring
• User engagement assessment

Eye Movement Types Detectable with EOG

Blinks:

• Large amplitude, brief duration (100-400 ms)
• NeuroKit2 primary focus
• Vertical EOG most sensitive

Saccades:

• Rapid, ballistic eye movements (20-80 ms)
• Step-like voltage deflections
• Horizontal or vertical
• Require higher sampling rates for detailed analysis

Smooth pursuit:

• Slow tracking of moving objects
• Ramp-like voltage changes
• Lower amplitude than saccades

Fixations:

• Stable gaze
• Baseline EOG with small oscillations
• Duration varies (200-600 ms typical in reading)

Note: Detailed saccade/fixation analysis typically requires eye tracking (infrared, video-based). EOG useful for blinks and gross eye movements.

Interpretation Guidelines

Blink rate:

• Normal resting: 15-20 blinks/min
• <10 blinks/min: Visual task engagement, concentration
• >30 blinks/min: Stress, dry eyes, fatigue
• Context-dependent: Task demands, lighting, screen use

Blink duration:

• Normal: 100-400 ms (mean ~200 ms)
• Prolonged: Drowsiness, fatigue (>500 ms)
• Short: Normal alertness

Blink amplitude:

• Varies with electrode placement and individuals
• Within-subject comparisons most reliable
• Incomplete blinks: reduced amplitude (dry eye, fatigue)

Temporal patterns:

• Clustered blinks: Transitions between tasks or cognitive states
• Suppressed blinks: Active visual processing, sustained attention
• Post-stimulus blinks: After completing visual processing

References

• Kleifges, K., Bigdely-Shamlo, N., Kerick, S. E., & Robbins, K. A. (2017). BLINKER: Automated extraction of ocular indices from EEG enabling large-scale analysis. Frontiers in Neuroscience, 11, 12.
• Agarwal, M., & Sivakumar, R. (2019). Blink: A fully automated unsupervised algorithm for eye-blink detection in EEG signals. In 2019 57th Annual Allerton Conference on Communication, Control, and Computing (pp. 1113-1121). IEEE.
• Kong, X., & Wilson, G. F. (1998). A new EOG-based eyeblink detection algorithm. Behavior Research Methods, Instruments, & Computers, 30(4), 713-719.
• Schleicher, R., Galley, N., Briest, S., & Galley, L. (2008). Blinks and saccades as indicators of fatigue in sleepiness warnings: Looking tired? Ergonomics, 51(7), 982-1010.

Epochs and Event-Related Analysis

Overview

Event-related analysis examines physiological responses time-locked to specific stimuli or events. NeuroKit2 provides tools for event detection, epoch creation, averaging, and event-related feature extraction across all signal types.

Event Detection

events_find()

Automatically detect events/triggers in a signal based on threshold crossings or changes.

events = nk.events_find(event_channel, threshold=0.5, threshold_keep='above',
                        duration_min=1, inter_min=0)

Parameters:

• threshold: Detection threshold value
• threshold_keep: 'above' or 'below' threshold
• duration_min: Minimum event duration (samples) to keep
• inter_min: Minimum interval between events (samples)

Returns:

• Dictionary with:
• 'onset': Event onset indices
• 'offset': Event offset indices (if applicable)
• 'duration': Event durations
• 'label': Event labels (if multiple event types)

Common use cases:

TTL triggers from experiments:

# Trigger channel: 0V baseline, 5V pulses during events
events = nk.events_find(trigger_channel, threshold=2.5, threshold_keep='above')

Button presses:

# Detect when button signal goes high
button_events = nk.events_find(button_signal, threshold=0.5, threshold_keep='above',
                               duration_min=10)  # Debounce

State changes:

# Detect periods above/below threshold
high_arousal = nk.events_find(eda_signal, threshold='auto', duration_min=100)

events_plot()

Visualize event timing relative to signals.

nk.events_plot(events, signal)

Displays:

• Signal trace
• Event markers (vertical lines or shaded regions)
• Event labels

Use case:

• Verify event detection accuracy
• Inspect temporal distribution of events
• Quality control before epoching

Epoch Creation

epochs_create()

Create epochs (segments) of data around events for event-related analysis.

epochs = nk.epochs_create(data, events, sampling_rate=1000,
                          epochs_start=-0.5, epochs_end=2.0,
                          event_labels=None, event_conditions=None,
                          baseline_correction=False)

Parameters:

• data: DataFrame with signals or single signal
• events: Event indices or dictionary from events_find()
• sampling_rate: Signal sampling rate (Hz)
• epochs_start: Start time relative to event (seconds, negative = before)
• epochs_end: End time relative to event (seconds, positive = after)
• event_labels: List of labels for each event (optional)
• event_conditions: List of condition names for each event (optional)
• baseline_correction: If True, subtract baseline mean from each epoch

Returns:

• Dictionary of DataFrames, one per epoch
• Each DataFrame contains signal data with time relative to event (Index=0 at event onset)
• Includes 'Label' and 'Condition' columns if provided

Typical epoch windows:

• Visual ERP: -0.2 to 1.0 seconds (200 ms baseline, 1 s post-stimulus)
• Cardiac orienting: -1.0 to 10 seconds (capture anticipation and response)
• EMG startle: -0.1 to 0.5 seconds (brief response)
• EDA SCR: -1.0 to 10 seconds (1-3 s latency, slow recovery)

Event Labels and Conditions

Organize events by type and experimental conditions:

# Example: Emotional picture experiment
event_times = [1000, 2500, 4200, 5800]  # Event onsets in samples
event_labels = ['trial1', 'trial2', 'trial3', 'trial4']
event_conditions = ['positive', 'negative', 'positive', 'neutral']

epochs = nk.epochs_create(signals, events=event_times, sampling_rate=1000,
                          epochs_start=-1, epochs_end=5,
                          event_labels=event_labels,
                          event_conditions=event_conditions)

Access epochs:

# Epoch by number
epoch_1 = epochs['1']

# Filter by condition
positive_epochs = {k: v for k, v in epochs.items() if v['Condition'][0] == 'positive'}

Baseline Correction

Remove pre-stimulus baseline from epochs to isolate event-related changes:

Automatic (during epoch creation):

epochs = nk.epochs_create(data, events, sampling_rate=1000,
                          epochs_start=-0.5, epochs_end=2.0,
                          baseline_correction=True)  # Subtracts mean of entire baseline

Manual (after epoch creation):

# Subtract baseline period mean
baseline_start = -0.5  # seconds
baseline_end = 0.0     # seconds

for key, epoch in epochs.items():
    baseline_mask = (epoch.index >= baseline_start) & (epoch.index < baseline_end)
    baseline_mean = epoch[baseline_mask].mean()
    epochs[key] = epoch - baseline_mean

When to baseline correct:

• ERPs: Always (isolates event-related change)
• Cardiac/EDA: Usually (removes inter-individual baseline differences)
• Absolute measures: Sometimes not desired (e.g., analyzing absolute amplitude)

Epoch Analysis and Visualization

epochs_plot()

Visualize individual or averaged epochs.

nk.epochs_plot(epochs, column='ECG_Rate', condition=None, show=True)

Parameters:

• column: Which signal column to plot
• condition: Plot only specific condition (optional)

Displays:

• Individual epoch traces (semi-transparent)
• Average across epochs (bold line)
• Optional: Shaded error (SEM or SD)

Use cases:

• Visualize event-related responses
• Compare conditions
• Identify outlier epochs

epochs_average()

Compute grand average across epochs with statistics.

average_epochs = nk.epochs_average(epochs, output='dict')

Parameters:

• output: 'dict' (default) or 'df' (DataFrame)

Returns:

• Dictionary or DataFrame with:
• 'Mean': Average across epochs at each time point
• 'SD': Standard deviation
• 'SE': Standard error of mean
• 'CI_lower', 'CI_upper': 95% confidence intervals

Use case:

• Compute event-related potentials (ERPs)
• Grand average cardiac/EDA/EMG responses
• Group-level analysis

Condition-specific averaging:

# Separate averages by condition
positive_epochs = {k: v for k, v in epochs.items() if v['Condition'][0] == 'positive'}
negative_epochs = {k: v for k, v in epochs.items() if v['Condition'][0] == 'negative'}

avg_positive = nk.epochs_average(positive_epochs)
avg_negative = nk.epochs_average(negative_epochs)

epochs_to_df()

Convert epochs dictionary to unified DataFrame.

epochs_df = nk.epochs_to_df(epochs)

Returns:

• Single DataFrame with all epochs stacked
• Includes 'Epoch', 'Time', 'Label', 'Condition' columns
• Facilitates statistical analysis and plotting with pandas/seaborn

Use case:

• Prepare data for mixed-effects models
• Plotting with seaborn/plotly
• Export to R or statistical software

epochs_to_array()

Convert epochs to 3D NumPy array.

epochs_array = nk.epochs_to_array(epochs, column='ECG_Rate')

Returns:

• 3D array: (n_epochs, n_timepoints, n_columns)

Use case:

• Machine learning input (epoched features)
• Custom array-based analysis
• Statistical tests on array data

Signal-Specific Event-Related Analysis

NeuroKit2 provides specialized event-related analysis for each signal type:

ECG Event-Related

ecg_epochs = nk.epochs_create(ecg_signals, events, sampling_rate=1000,
                              epochs_start=-1, epochs_end=10)
ecg_results = nk.ecg_eventrelated(ecg_epochs)

Computed metrics:

• ECG_Rate_Baseline: Heart rate before event
• ECG_Rate_Min/Max: Minimum/maximum rate during epoch
• ECG_Phase_*: Cardiac phase at event onset
• Rate dynamics across time windows

EDA Event-Related

eda_epochs = nk.epochs_create(eda_signals, events, sampling_rate=100,
                              epochs_start=-1, epochs_end=10)
eda_results = nk.eda_eventrelated(eda_epochs)

Computed metrics:

• EDA_SCR: Presence of SCR (binary)
• SCR_Amplitude: Maximum SCR amplitude
• SCR_Latency: Time to SCR onset
• SCR_RiseTime, SCR_RecoveryTime
• EDA_Tonic: Mean tonic level

RSP Event-Related

rsp_epochs = nk.epochs_create(rsp_signals, events, sampling_rate=100,
                              epochs_start=-0.5, epochs_end=5)
rsp_results = nk.rsp_eventrelated(rsp_epochs)

Computed metrics:

• RSP_Rate_Mean: Average breathing rate
• RSP_Amplitude_Mean: Average breath depth
• RSP_Phase: Respiratory phase at event
• Rate/amplitude dynamics

EMG Event-Related

emg_epochs = nk.epochs_create(emg_signals, events, sampling_rate=1000,
                              epochs_start=-0.1, epochs_end=1.0)
emg_results = nk.emg_eventrelated(emg_epochs)

Computed metrics:

• EMG_Activation: Presence of activation
• EMG_Amplitude_Mean/Max: Amplitude statistics
• EMG_Onset_Latency: Time to activation onset
• EMG_Bursts: Number of activation bursts

EOG Event-Related

eog_epochs = nk.epochs_create(eog_signals, events, sampling_rate=500,
                              epochs_start=-0.5, epochs_end=2.0)
eog_results = nk.eog_eventrelated(eog_epochs)

Computed metrics:

• EOG_Blinks_N: Number of blinks during epoch
• EOG_Rate_Mean: Blink rate
• Temporal blink distribution

PPG Event-Related

ppg_epochs = nk.epochs_create(ppg_signals, events, sampling_rate=100,
                              epochs_start=-1, epochs_end=10)
ppg_results = nk.ppg_eventrelated(ppg_epochs)

Computed metrics:

• Similar to ECG: rate dynamics, phase information

Practical Workflows

Complete Event-Related Analysis Pipeline

import neurokit2 as nk

# 1. Process physiological signals
ecg_signals, ecg_info = nk.ecg_process(ecg, sampling_rate=1000)
eda_signals, eda_info = nk.eda_process(eda, sampling_rate=100)

# 2. Align sampling rates if needed
eda_signals_resampled = nk.signal_resample(eda_signals, sampling_rate=100,
                                           desired_sampling_rate=1000)

# 3. Merge signals into single DataFrame
signals = pd.concat([ecg_signals, eda_signals_resampled], axis=1)

# 4. Detect events
events = nk.events_find(trigger_channel, threshold=0.5)

# 5. Add event labels and conditions
event_labels = ['trial1', 'trial2', 'trial3', ...]
event_conditions = ['condition_A', 'condition_B', 'condition_A', ...]

# 6. Create epochs
epochs = nk.epochs_create(signals, events, sampling_rate=1000,
                          epochs_start=-1.0, epochs_end=5.0,
                          event_labels=event_labels,
                          event_conditions=event_conditions,
                          baseline_correction=True)

# 7. Signal-specific event-related analysis
ecg_results = nk.ecg_eventrelated(epochs)
eda_results = nk.eda_eventrelated(epochs)

# 8. Merge results
results = pd.merge(ecg_results, eda_results, left_index=True, right_index=True)

# 9. Statistical analysis by condition
results['Condition'] = event_conditions
condition_comparison = results.groupby('Condition').mean()

Handling Multiple Event Types

# Different event types with different markers
event_type1 = nk.events_find(trigger_ch1, threshold=0.5)
event_type2 = nk.events_find(trigger_ch2, threshold=0.5)

# Combine events with labels
all_events = np.concatenate([event_type1['onset'], event_type2['onset']])
event_labels = ['type1'] * len(event_type1['onset']) + ['type2'] * len(event_type2['onset'])

# Sort by time
sort_idx = np.argsort(all_events)
all_events = all_events[sort_idx]
event_labels = [event_labels[i] for i in sort_idx]

# Create epochs
epochs = nk.epochs_create(signals, all_events, sampling_rate=1000,
                          epochs_start=-0.5, epochs_end=3.0,
                          event_labels=event_labels)

# Separate by type
type1_epochs = {k: v for k, v in epochs.items() if v['Label'][0] == 'type1'}
type2_epochs = {k: v for k, v in epochs.items() if v['Label'][0] == 'type2'}

Quality Control and Artifact Rejection

# Remove epochs with excessive noise or artifacts
clean_epochs = {}
for key, epoch in epochs.items():
    # Example: reject if EDA amplitude too high (movement artifact)
    if epoch['EDA_Phasic'].abs().max() < 5.0:  # Threshold
        # Example: reject if heart rate change too large (invalid)
        if epoch['ECG_Rate'].max() - epoch['ECG_Rate'].min() < 50:
            clean_epochs[key] = epoch

print(f"Kept {len(clean_epochs)}/{len(epochs)} epochs")

# Analyze clean epochs
results = nk.ecg_eventrelated(clean_epochs)

Statistical Considerations

Sample Size

• ERP/averaging: 20-30+ trials per condition minimum
• Individual trial analysis: Mixed-effects models handle variable trial counts
• Group comparisons: Pilot data for power analysis

Time Window Selection

• A priori hypotheses: Pre-register time windows based on literature
• Exploratory: Use full epoch, correct for multiple comparisons
• Avoid: Selecting windows based on observed data (circular)

Baseline Period

• Should be free of anticipatory effects
• Sufficient duration for stable estimate (500-1000 ms typical)
• Shorter for fast dynamics (e.g., startle: 100 ms sufficient)

Condition Comparison

• Repeated measures ANOVA for within-subject designs
• Mixed-effects models for unbalanced data
• Permutation tests for non-parametric comparisons
• Correct for multiple comparisons (time points/signals)

Common Applications

Cognitive psychology:

• P300 ERP analysis
• Error-related negativity (ERN)
• Attentional blink
• Working memory load effects

Affective neuroscience:

• Emotional picture viewing (EDA, HR, facial EMG)
• Fear conditioning (HR deceleration, SCR)
• Valence/arousal dimensions

Clinical research:

• Startle response (orbicularis oculi EMG)
• Orienting response (HR deceleration)
• Anticipation and prediction error

Psychophysiology:

• Cardiac defense response
• Orienting vs. defensive reflexes
• Respiratory changes during emotion

Human-computer interaction:

• User engagement during events
• Surprise/violation of expectation
• Cognitive load during task events

References

• Luck, S. J. (2014). An introduction to the event-related potential technique (2nd ed.). MIT press.
• Bradley, M. M., & Lang, P. J. (2000). Measuring emotion: Behavior, feeling, and physiology. In R. D. Lane & L. Nadel (Eds.), Cognitive neuroscience of emotion (pp. 242-276). Oxford University Press.
• Boucsein, W. (2012). Electrodermal activity (2nd ed.). Springer.
• Gratton, G., Coles, M. G., & Donchin, E. (1983). A new method for off-line removal of ocular artifact. Electroencephalography and clinical neurophysiology, 55(4), 468-484.