neurokit2

NeuroKit2

概述

NeuroKit2 提供用于生理信号处理的统一高级 API。每种信号类型（ECG、EEG、EMG、EDA、PPG、RSP）都遵循相同的 nk.{signal}_process() → nk.{signal}_analyze() 工作流：输入原始信号，输出清理后的信号和特征。该库会自动处理去趋势、滤波、峰值检测、伪影校正和特征提取，参数针对生物信号的特性进行了调整。结果以 pandas DataFrame 的形式返回，便于进行后续统计分析。NeuroKit2 还提供 nk.{signal}_simulate() 函数，用于生成合成测试信号。

适用场景

• 从 ECG 记录中提取 HRV（心率变异性）特征，用于压力或自主神经系统分析
• 检测 ECG 中的 R 波峰，并计算 RR 间期、SDNN、RMSSD、pNN50、LF/HF 比值
• 处理 EDA/GSR 信号，将其分解为紧张性（SCL）和相位性（SCR）成分，用于心理生理学研究
• 清理并分段 EMG 信号，以便在生物力学中检测肌肉激活的起始和结束
• 处理可穿戴传感器采集的 PPG 信号，将其作为 ECG 的替代信号来测量心率和 SpO2
• 生成合成生理信号，用于算法验证和单元测试
• 处理多通道 EEG/MEG 并进行源定位时，请改用 MNE；处理原始底层 DSP 时，请使用 scipy.signal

前置条件

• Python 软件包：neurokit2、numpy、pandas、matplotlib
• 数据要求：生理信号的一维 NumPy 数组或 pandas Series；已知采样率（Hz）
• 典型采样率：ECG 250–1000 Hz，EEG 256–2048 Hz，EDA 4–64 Hz，EMG 1000–2000 Hz

pip install neurokit2 numpy pandas matplotlib scipy

快速开始

import neurokit2 as nk
import matplotlib.pyplot as plt

# Generate and process synthetic ECG (10 seconds at 500 Hz)
ecg_signal = nk.ecg_simulate(duration=10, sampling_rate=500, heart_rate=70)
signals, info = nk.ecg_process(ecg_signal, sampling_rate=500)

# Plot processed ECG
nk.ecg_plot(signals, info)
plt.savefig("ecg_processed.pdf", bbox_inches="tight")
print(f"R-peaks detected: {len(info['ECG_R_Peaks'])}")
print(signals[["ECG_Clean", "ECG_Rate", "ECG_Quality"]].describe())

核心 API

模块 1：ECG 处理

从原始 ECG 到清理后的信号、R 波峰和瞬时心率的完整处理流程。

import neurokit2 as nk
import numpy as np

# Load real data (example: CSV with one ECG column at 250 Hz)
# ecg_raw = pd.read_csv("ecg_recording.csv")["ecg"].values
# Or use synthetic:
ecg_raw = nk.ecg_simulate(duration=60, sampling_rate=250, heart_rate=72, noise=0.05)

# Full processing pipeline
signals, info = nk.ecg_process(ecg_raw, sampling_rate=250)

# signals: DataFrame with columns ECG_Raw, ECG_Clean, ECG_Rate, ECG_R_Peaks, ...
# info: dict with R-peak indices, P/Q/S/T wave locations

print(f"R-peaks: {len(info['ECG_R_Peaks'])} detected")
print(f"Mean heart rate: {signals['ECG_Rate'].mean():.1f} bpm")
print(f"ECG quality (0-1): {signals['ECG_Quality'].mean():.2f}")

# Get delineated waves (P, Q, S, T)
_, waves_dict = nk.ecg_delineate(signals["ECG_Clean"], info["ECG_R_Peaks"],
                                   sampling_rate=250, method="dwt")
print(f"P-wave peaks found: {np.sum(~np.isnan(waves_dict['ECG_P_Peaks']))}")

# HRV analysis from ECG
hrv_time = nk.hrv_time(info["ECG_R_Peaks"], sampling_rate=250)
hrv_freq = nk.hrv_frequency(info["ECG_R_Peaks"], sampling_rate=250)
hrv_nonlinear = nk.hrv_nonlinear(info["ECG_R_Peaks"], sampling_rate=250)

print("Time-domain HRV:")
print(f"  SDNN  : {hrv_time['HRV_SDNN'].values[0]:.2f} ms")
print(f"  RMSSD : {hrv_time['HRV_RMSSD'].values[0]:.2f} ms")
print(f"  pNN50 : {hrv_time['HRV_pNN50'].values[0]:.2f}%")

print("\nFrequency-domain HRV:")
print(f"  LF power  : {hrv_freq['HRV_LF'].values[0]:.4f} ms²")
print(f"  HF power  : {hrv_freq['HRV_HF'].values[0]:.4f} ms²")
print(f"  LF/HF     : {hrv_freq['HRV_LFHF'].values[0]:.3f}")

模块 2：EDA / GSR 处理

将皮肤电活动（EDA）/电皮肤反应（GSR）分解为紧张性和相位性成分。

import neurokit2 as nk

# Simulate EDA signal (10 min at 4 Hz with 3 events)
eda_raw = nk.eda_simulate(duration=600, sampling_rate=4, scr_number=10, noise=0.01)

# Process: detrend, filter, decompose into SCL (tonic) + SCR (phasic)
signals, info = nk.eda_process(eda_raw, sampling_rate=4)

print(f"SCR peaks detected: {len(info['SCR_Peaks'])}")
print(f"SCR recovery times (mean): {signals['SCR_RecoveryTime'].mean():.2f} s")
print(f"SCL (tonic) mean: {signals['EDA_Tonic'].mean():.4f} μS")
print(f"SCR (phasic) amplitude mean: {signals['EDA_Phasic'].mean():.4f} μS")

# Analyze epochs around events
events = nk.events_create(event_onsets=[60, 120, 240], event_durations=3,
                           desired_length=len(eda_raw))
epoch_df = nk.epochs_create(signals, events, sampling_rate=4,
                              epochs_start=-5, epochs_end=10)

模块 3：EMG 处理

从表面 EMG 中检测肌肉激活。

import neurokit2 as nk

# Simulate EMG with 3 activations at 1000 Hz
emg_raw = nk.emg_simulate(duration=10, sampling_rate=1000, burst_number=3)

# Process: rectify, envelope, detect activation periods
signals, info = nk.emg_process(emg_raw, sampling_rate=1000)

print(f"EMG activation onsets: {len(info['EMG_Onsets'])}")
print(f"EMG activation offsets: {len(info['EMG_Offsets'])}")

# Activation periods
for onset, offset in zip(info["EMG_Onsets"], info["EMG_Offsets"]):
    duration_ms = (offset - onset) / 1000 * 1000  # samples → ms
    print(f"  Activation: {onset/1000:.2f}s – {offset/1000:.2f}s ({duration_ms:.0f} ms)")

# Plot
nk.emg_plot(signals)
import matplotlib.pyplot as plt
plt.savefig("emg_processed.pdf", bbox_inches="tight")

模块 4：EEG / 复杂度分析

适用于 EEG 和其他生物信号的信号复杂度测量。

import neurokit2 as nk
import numpy as np

# Generate EEG-like signal
eeg = nk.eeg_simulate(duration=30, sampling_rate=256, brain_information=0.7)

# Band power (delta, theta, alpha, beta, gamma)
bands = nk.eeg_power(eeg, sampling_rate=256,
                      frequency_band=["Delta", "Theta", "Alpha", "Beta", "Gamma"])
print("EEG Band Power:")
for band, power in bands.items():
    print(f"  {band}: {power:.4f}")

# Complexity measures (for any 1D biosignal)
signal = np.random.randn(1000)  # test signal

entropy_sample = nk.entropy_sample(signal, order=2, r=0.2 * np.std(signal))
entropy_fuzzy = nk.entropy_fuzzy(signal, order=2)
dfa = nk.fractal_dfa(signal)

print(f"Sample entropy: {entropy_sample[0]:.4f}")
print(f"Fuzzy entropy: {entropy_fuzzy[0]:.4f}")
print(f"DFA exponent: {dfa[0]:.4f}")

模块 5：PPG 处理

处理来自可穿戴传感器的光电容积描记（PPG）信号。

import neurokit2 as nk

# Simulate PPG at 100 Hz (common wearable sampling rate)
ppg_raw = nk.ppg_simulate(duration=30, sampling_rate=100, heart_rate=65, noise=0.01)

# Process: clean, detect peaks
signals, info = nk.ppg_process(ppg_raw, sampling_rate=100)

print(f"PPG peaks detected: {len(info['PPG_Peaks'])}")
print(f"Mean heart rate: {signals['PPG_Rate'].mean():.1f} bpm")

# Compute HRV from PPG peaks (as ECG surrogate)
hrv = nk.hrv(info["PPG_Peaks"], sampling_rate=100)
print(f"RMSSD from PPG: {hrv['HRV_RMSSD'].values[0]:.2f} ms")

模块 6：信号模拟

生成合成生理信号，用于测试和算法验证。

import neurokit2 as nk
import matplotlib.pyplot as plt

fig, axes = plt.subplots(5, 1, figsize=(12, 10))

signals_dict = {
    "ECG":  nk.ecg_simulate(duration=5, sampling_rate=500, heart_rate=70),
    "EDA":  nk.eda_simulate(duration=5, sampling_rate=64, scr_number=2),
    "EMG":  nk.emg_simulate(duration=5, sampling_rate=1000, burst_number=2),
    "PPG":  nk.ppg_simulate(duration=5, sampling_rate=100, heart_rate=70),
    "RSP":  nk.rsp_simulate(duration=5, sampling_rate=100, respiratory_rate=15),
}

for ax, (name, signal) in zip(axes, signals_dict.items()):
    sr = {"ECG": 500, "EDA": 64, "EMG": 1000, "PPG": 100, "RSP": 100}[name]
    import numpy as np
    t = np.arange(len(signal)) / sr
    ax.plot(t, signal, lw=0.8)
    ax.set_ylabel(name)

axes[-1].set_xlabel("Time (s)")
plt.tight_layout()
plt.savefig("physiological_signals.pdf", bbox_inches="tight")
print("Synthetic signals plotted")

核心概念

事件相关分析（Epochs）

NeuroKit2 使用 nk.events_create() 标记刺激开始时间，然后使用 nk.epochs_create() 将连续信号分割为围绕各个事件的固定长度窗口。这使得能够根据刺激时间锁定 ERP（事件相关电位）和 SCR（皮肤电反应）分析。

信号质量指数（SQI）

ECG 和 PPG 处理会为每个采样点计算质量指数（0–1）。低于 0.5 的质量值表示该片段存在噪声或伪影。计算 HRV 特征之前，可使用 signals["ECG_Quality"] > 0.5 屏蔽不可靠区域。

常见工作流

工作流 1：多信号生理记录分析

import neurokit2 as nk
import pandas as pd
import numpy as np

# Simulate 5-minute multi-modal recording
sr_ecg, sr_eda, sr_ppg = 250, 4, 100
duration = 300  # seconds

ecg = nk.ecg_simulate(duration=duration, sampling_rate=sr_ecg, heart_rate=72)
eda = nk.eda_simulate(duration=duration, sampling_rate=sr_eda, scr_number=15)
ppg = nk.ppg_simulate(duration=duration, sampling_rate=sr_ppg, heart_rate=72)

# Process each modality
ecg_signals, ecg_info = nk.ecg_process(ecg, sampling_rate=sr_ecg)
eda_signals, eda_info = nk.eda_process(eda, sampling_rate=sr_eda)
ppg_signals, ppg_info = nk.ppg_process(ppg, sampling_rate=sr_ppg)

# Extract HRV from ECG
hrv = nk.hrv(ecg_info["ECG_R_Peaks"], sampling_rate=sr_ecg, show=False)

# Summary table
summary = pd.Series({
    "HR_mean_bpm":       ecg_signals["ECG_Rate"].mean(),
    "HRV_RMSSD_ms":     hrv["HRV_RMSSD"].values[0],
    "HRV_LF_HF":        hrv["HRV_LFHF"].values[0],
    "SCL_mean_uS":      eda_signals["EDA_Tonic"].mean(),
    "SCR_count":        len(eda_info["SCR_Peaks"]),
    "PPG_HR_mean_bpm":  ppg_signals["PPG_Rate"].mean(),
})
print(summary.to_string())

工作流 2：压力检测特征提取

import neurokit2 as nk
import pandas as pd
import numpy as np

def extract_stress_features(ecg_array, eda_array, sr_ecg=250, sr_eda=4):
    """Extract HRV + EDA features for stress classification."""
    # ECG features
    ecg_signals, ecg_info = nk.ecg_process(ecg_array, sampling_rate=sr_ecg)
    hrv_time = nk.hrv_time(ecg_info["ECG_R_Peaks"], sampling_rate=sr_ecg)
    hrv_freq = nk.hrv_frequency(ecg_info["ECG_R_Peaks"], sampling_rate=sr_ecg)

    # EDA features
    eda_signals, eda_info = nk.eda_process(eda_array, sampling_rate=sr_eda)

    return {
        "HR_mean":    ecg_signals["ECG_Rate"].mean(),
        "RMSSD":      hrv_time["HRV_RMSSD"].values[0],
        "SDNN":       hrv_time["HRV_SDNN"].values[0],
        "LF_HF":      hrv_freq["HRV_LFHF"].values[0],
        "SCL_mean":   eda_signals["EDA_Tonic"].mean(),
        "SCR_count":  len(eda_info["SCR_Peaks"]),
        "SCR_amp_mean": eda_signals["EDA_Phasic"].mean(),
    }

# Test with simulated data
ecg = nk.ecg_simulate(duration=120, sampling_rate=250, heart_rate=85)  # elevated HR
eda = nk.eda_simulate(duration=120, sampling_rate=4, scr_number=8)

features = extract_stress_features(ecg, eda)
df = pd.DataFrame([features])
print(df.T.to_string(header=False))

关键参数

参数	模块/函数	默认值	范围 / 选项	作用
sampling_rate	所有 *_process	必填	64–2048 Hz（取决于模态）	必须与实际记录的采样率一致，才能确保计时准确
method	ecg_peaks	"neurokit"	"neurokit"、"pantompkins1985"、"hamilton2002"、"elgendi2010"	R 峰检测算法；"neurokit" 最稳健
heart_rate	ecg_simulate	70	30–200 bpm	模拟的 ECG 心率
noise	ecg_simulate	0.01	0–0.5	添加到模拟信号中的高斯噪声幅度
scr_number	eda_simulate	1	0–50	模拟 EDA 中 SCR 事件的数量
order	entropy_sample	2	1–5	计算熵时使用的嵌入维度
r	entropy_sample	0.2×SD	0.1–0.5×SD	样本熵的容差
epochs_start	epochs_create	-1	负浮点数（s）	分段开始时间，即事件发生前的秒数
epochs_end	epochs_create	3	正浮点数（s）	分段结束时间，即事件发生后的秒数

常用方案

方案：批量处理多个 ECG 文件

import neurokit2 as nk
import pandas as pd
from pathlib import Path
import numpy as np

results = []
for ecg_file in Path("ecg_data").glob("*.csv"):
    try:
        ecg = pd.read_csv(ecg_file)["ecg"].values
        _, info = nk.ecg_process(ecg, sampling_rate=250)
        hrv = nk.hrv_time(info["ECG_R_Peaks"], sampling_rate=250)
        results.append({
            "file": ecg_file.stem,
            "RMSSD": hrv["HRV_RMSSD"].values[0],
            "SDNN":  hrv["HRV_SDNN"].values[0],
            "pNN50": hrv["HRV_pNN50"].values[0],
        })
    except Exception as e:
        print(f"FAILED {ecg_file.name}: {e}")

df = pd.DataFrame(results)
df.to_csv("hrv_results.csv", index=False)
print(df.describe())

方案：按质量指数剔除伪迹

import neurokit2 as nk
import numpy as np

ecg = nk.ecg_simulate(duration=60, sampling_rate=250, noise=0.1)
signals, info = nk.ecg_process(ecg, sampling_rate=250)

# Identify low-quality R-peaks
r_peaks = info["ECG_R_Peaks"]
quality_at_peaks = signals["ECG_Quality"].iloc[r_peaks].values

good_peaks = r_peaks[quality_at_peaks > 0.5]
print(f"R-peaks before cleaning: {len(r_peaks)}")
print(f"High-quality R-peaks  : {len(good_peaks)}")

# Recompute HRV on clean peaks only
hrv_clean = nk.hrv_time(good_peaks, sampling_rate=250)
print(f"RMSSD (cleaned): {hrv_clean['HRV_RMSSD'].values[0]:.2f} ms")

故障排除

问题	原因	解决方案
R 峰检测位置错误	sampling_rate 不正确	仔细检查记录元数据；错误的采样率会导致所有峰值时间偏移
HRV RMSSD 极高（>200 ms）	漏检心搏或额外检测	更换检测方法：尝试 method="pantompkins1985"；检查原始信号
ecg_delineate 对 P/Q/S/T 返回的全是 NaN	记录时长过短（<5 秒）或噪声过大	使用时长 ≥10 s 的记录；调用 delineate 前先降低噪声
真实数据中的 EDA 未检测到 SCR 峰	低通滤波器截止频率设置过高，或 EDA 采样率错误	验证 EDA 采样率；在 eda_phasic() 中设置 phasic_method="cvx"
熵函数返回 NaN	信号过短或为恒定值（方差为零）	至少使用 200–500 个样本；确保信号不是恒定值
EMG 激活起始检测漏掉激活	对低幅值信号自动设置的阈值过高	传入自定义阈值：nk.emg_activation(signal, threshold=0.5)
长时间记录出现内存错误	DataFrame 包含数百万行	按 60 秒窗口处理；拼接 HRV 分段结果

相关技能

• matplotlib-scientific-plotting — 绘制处理后的 neurokit2 信号和 HRV Poincaré 图
• statsmodels-statistical-modeling — 对提取的 HRV/EDA 特征进行 ANOVA 和混合效应模型分析
• scikit-learn-machine-learning — 使用提取的特征进行压力/情绪分类

参考资料

• NeuroKit2 文档 — API 参考、教程和算法对比
• NeuroKit2 论文：Makowski 等（2021），Behavior Research Methods — 方法与验证
• HRV 标准：Task Force（1996），Circulation — HRV 时域和频域指标的权威定义
• NeuroKit2 GitHub — 源代码、示例和社区
• EDA 综述：Boucsein（2012），Electrodermal Activity — EDA 信号解读参考书