神经生理信号处理的Python工具箱_neurokit2

NeuroKit2是一个用户友好的包，提供了先进的生物信号处理程序的方便访问。没有丰富的编程知识或生物医学信号处理知识的研究人员和临床医生只需要两行代码就可以分析生理数据。

举个栗子：

import neurokit2 as nk
 
# Download example data
data = nk.data("bio_eventrelated_100hz")
 
# Preprocess the data (filter, find peaks, etc.)
processed_data, info = nk.bio_process(ecg=data["ECG"], rsp=data["RSP"], eda=data["EDA"], sampling_rate=100)
 
# Compute relevant features
results = nk.bio_analyze(processed_data, sampling_rate=100)

然后就分析完信号了。

安装：

你可以从PyPI安装NeuroKit2

pip install neurokit2

或conda-forge

conda install -c conda-forge neurokit2

文件：

通用

①10分钟内熟悉Python
②记录高质量信号
③多好的生理信号处理软件啊
④安装Python和NeuroKit
⑤包括数据集
⑥额外的资源

例子

①模拟人工生理信号
②定制处理管道
③与事件相关的分析
④Interval-related分析
⑤分析皮肤电活动(EDA)
⑥分析呼吸率变异性(RRV)
⑦提取和可视化个人心跳
⑧定位心电图P、Q、S、T波
⑨生理信号的复杂性分析
⑩分析眼电图EOG数据
⑪将函数与信号匹配

不知道哪个适合你的情况?参照下面的流程图:

生理信号预处理

模拟生理信号

import numpy as np
import pandas as pd
import neurokit2 as nk
 
# Generate synthetic signals
ecg = nk.ecg_simulate(duration=10, heart_rate=70)
ppg = nk.ppg_simulate(duration=10, heart_rate=70)
rsp = nk.rsp_simulate(duration=10, respiratory_rate=15)
eda = nk.eda_simulate(duration=10, scr_number=3)
emg = nk.emg_simulate(duration=10, burst_number=2)
 
# Visualise biosignals
data = pd.DataFrame({"ECG": ecg,
                     "PPG": ppg,
                     "RSP": rsp,
                     "EDA": eda,
                     "EMG": emg})
nk.signal_plot(data, subplots=True)

皮肤电活动(EDA / GSR，Electrodermal Activity)

# Generate 10 seconds of EDA signal (recorded at 250 samples / second) with 2 SCR peaks
eda = nk.eda_simulate(duration=10, sampling_rate=250, scr_number=2, drift=0.01)
 
# Process it
signals, info = nk.eda_process(eda, sampling_rate=250)
 
# Visualise the processing
nk.eda_plot(signals, sampling_rate=250)

 心脏活动(ECG，Cardiac activity)

# Generate 15 seconds of ECG signal (recorded at 250 samples / second)
ecg = nk.ecg_simulate(duration=15, sampling_rate=250, heart_rate=70)
 
# Process it
signals, info = nk.ecg_process(ecg, sampling_rate=250)
 
# Visualise the processing
nk.ecg_plot(signals, sampling_rate=250)

 呼吸(RSP，Respiration )

# Generate one minute of respiratory (RSP) signal (recorded at 250 samples / second)
rsp = nk.rsp_simulate(duration=60, sampling_rate=250, respiratory_rate=15)
 
# Process it
signals, info = nk.rsp_process(rsp, sampling_rate=250)
 
# Visualise the processing
nk.rsp_plot(signals, sampling_rate=250)

肌电图(EMG，Electromyography )

# Generate 10 seconds of EMG signal (recorded at 250 samples / second)
emg = nk.emg_simulate(duration=10, sampling_rate=250, burst_number=3)
 
# Process it
signal, info = nk.emg_process(emg, sampling_rate=250)
 
# Visualise the processing
nk.emg_plot(signals, sampling_rate=250)

光电容积脉搏波 (PPG/BVP，Photoplethysmography)

# Generate 15 seconds of PPG signal (recorded at 250 samples / second)
ppg = nk.ppg_simulate(duration=15, sampling_rate=250, heart_rate=70)
 
# Process it
signals, info = nk.ppg_process(ppg, sampling_rate=250)
 
# Visualize the processing
nk.ppg_plot(signals, sampling_rate=250)

眼动电图(EOG， Electrooculography )

# Import EOG data
eog_signal = nk.data("eog_100hz")
 
# Process it
signals, info = nk.eog_process(eog_signal, sampling_rate=100)
 
# Plot
plot = nk.eog_plot(signals, info, sampling_rate=100)

 生理数据分析

生理数据的分析通常分为两种类型，事件相关或间隔相关。

 事件相关

这种类型的分析指的是在对某一事件作出反应时立即发生的生理变化。例如，上图中虚线所示的刺激(例如，情绪刺激)呈现后的生理变化。在这种情况下，分析是基于时代的。epoch是生理信号的一小块(通常小于10秒)，锁定于特定的刺激，因此感兴趣的生理信号被相应的时间分段。这由上图中的橙色框表示。在本例中，使用bio_analyze() will。

间隔相关

这种类型的分析指的是在较长的时间内(从几秒钟到几天的活动)出现的生理特征和特征。典型的用例是休息状态的时期，在此期间活动被记录了几分钟，而参与者正在休息，或者在不同的条件下，没有特定的时间锁定事件(例如，看电影，听音乐，从事体育活动，等等)。例如，当人们想要比较不同强度的物理运动下的生理活动时，就会使用这种类型的分析。

心率变异性

# Download data
data = nk.data("bio_resting_8min_100hz")
 
# Find peaks
peaks, info = nk.ecg_peaks(data["ECG"], sampling_rate=100)
 
# Compute HRV indices
nk.hrv(peaks, sampling_rate=100, show=True)
>>>    HRV_RMSSD  HRV_MeanNN   HRV_SDNN  ...   HRV_CVI  HRV_CSI_Modified  HRV_SampEn
>>> 0  69.697983  696.395349  62.135891  ...  4.829101        592.095372    1.259931

 其他

心电图描述

描绘心电信号(ECG)的QRS波，包括p峰、t峰以及它们的起止和偏移。

# Download data
ecg_signal = nk.data(dataset="ecg_3000hz")['ECG']
 
# Extract R-peaks locations
_, rpeaks = nk.ecg_peaks(ecg_signal, sampling_rate=3000)
 
# Delineate
signal, waves = nk.ecg_delineate(ecg_signal, rpeaks, sampling_rate=3000, method="dwt", show=True, show_type='all')

 信号处理

信号处理功能：
①过滤:使用不同的方法。
②去趋势:消除基线漂移或趋势。
③失真:添加噪音和工件。

# Generate original signal
original = nk.signal_simulate(duration=6, frequency=1)
 
# Distort the signal (add noise, linear trend, artifacts etc.)
distorted = nk.signal_distort(original,
                              noise_amplitude=0.1,
                              noise_frequency=[5, 10, 20],
                              powerline_amplitude=0.05,
                              artifacts_amplitude=0.3,
                              artifacts_number=3,
                              linear_drift=0.5)
 
# Clean (filter and detrend)
cleaned = nk.signal_detrend(distorted)
cleaned = nk.signal_filter(cleaned, lowcut=0.5, highcut=1.5)
 
# Compare the 3 signals
plot = nk.signal_plot([original, distorted, cleaned])

 复杂性(熵，分形维数，…)

优化复杂度参数(延迟tau，尺寸m，公差r)

# Generate signal
signal = nk.signal_simulate(frequency=[1, 3], noise=0.01, sampling_rate=100)
 
# Find optimal time delay, embedding dimension and r
parameters = nk.complexity_optimize(signal, show=True)

 计算复杂性特征

①熵:样本熵(SampEn)、近似熵(ApEn)、模糊熵(FuzzEn)、多尺度熵(MSE)、Shannon熵(ShEn)

②分形维数:相关维数D2，…

③消除趋势波动分析法

信号分解

# Create complex signal
signal = nk.signal_simulate(duration=10, frequency=1)  # High freq
signal += 3 * nk.signal_simulate(duration=10, frequency=3)  # Higher freq
signal += 3 * np.linspace(0, 2, len(signal))  # Add baseline and linear trend
signal += 2 * nk.signal_simulate(duration=10, frequency=0.1, noise=0)  # Non-linear trend
signal += np.random.normal(0, 0.02, len(signal))  # Add noise
 
# Decompose signal using Empirical Mode Decomposition (EMD)
components = nk.signal_decompose(signal, method='emd')
nk.signal_plot(components)  # Visualize components
 
# Recompose merging correlated components
recomposed = nk.signal_recompose(components, threshold=0.99)
nk.signal_plot(recomposed)  # Visualize components

信号功率谱密度（PSD）

# Generate complex signal
signal = nk.signal_simulate(duration=20, frequency=[0.5, 5, 10, 15], amplitude=[2, 1.5, 0.5, 0.3], noise=0.025)
 
# Get the PSD using different methods
welch = nk.signal_psd(signal, method="welch", min_frequency=1, max_frequency=20, show=True)
multitaper = nk.signal_psd(signal, method="multitapers", max_frequency=20, show=True)
lomb = nk.signal_psd(signal, method="lomb", min_frequency=1, max_frequency=20, show=True)
burg = nk.signal_psd(signal, method="burg", min_frequency=1, max_frequency=20, order=10, show=True)

 统计数字

最高密度间隔（HDI）

x = np.random.normal(loc=0, scale=1, size=100000)
 
ci_min, ci_max = nk.hdi(x, ci=0.95, show=True)