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 信号をクリーンアップしてセグメント化する
• 心拍数や SpO2 の ECG 代替指標として、ウェアラブルセンサーの PPG 信号を処理する
• アルゴリズムの検証や単体テストのために、生理信号の合成データを生成する
• マルチチャンネル EEG/MEG や電流源推定を扱う場合は代わりに MNE を使用します。低レベルの生の DSP には scipy.signal を使用します

前提条件

• Python パッケージ：neurokit2、numpy、pandas、matplotlib
• データ要件：生理信号の 1D 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")

重要な概念

事象関連解析（エポック）

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	負の浮動小数点数（秒）	イベント開始前からエポック開始までの秒数です
epochs_end	epochs_create	3	正の浮動小数点数（秒）	イベント開始後からエポック終了までの秒数です

よく使うレシピ

レシピ：複数の 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 が正しくありません	記録のメタデータを再確認してください。SR が誤っていると、すべてのピーク時刻がずれます
HRV RMSSD が極端に高い（>200 ms）	拍動の見落とし、または余分な検出	検出方法を切り替えて、method="pantompkins1985" を試してください。生信号も確認してください
ecg_delineate が P/Q/S/T すべてに NaN を返す	記録時間が短すぎる（<5 秒）、またはノイズが非常に多い	10 秒以上の記録を使用し、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 のポアンカレプロットを作成します
• statsmodels-statistical-modeling — 抽出した HRV／EDA 特徴量に対する ANOVA および混合効果モデル
• scikit-learn-machine-learning — 抽出した特徴量を使用したストレス／感情の分類

参考資料

• NeuroKit2 ドキュメント — API リファレンス、チュートリアル、アルゴリズムの比較
• NeuroKit2 論文：Makowski et al. (2021), Behavior Research Methods — 方法論と検証
• HRV 標準：Task Force (1996), Circulation — HRV の時間領域および周波数領域指標の正式な定義
• NeuroKit2 GitHub — ソース、例、コミュニティ
• EDA レビュー：Boucsein (2012), Electrodermal Activity — EDA 信号の解釈に関する参考書