超详细分析如何用python处理ECG&PPG数据

这篇博文主要和大家分享一下如何使用python处理ECG和PPG的数据，从而使用PPG和ECG的数据进行血压的推测。
首先普及一下ECG和PPG，首先ECG 心电图(electrocardiogram)心脏在每个心动周期中,由起搏点、心房、心室相继兴奋,伴随着生物电的变化,通过心电描记器从体表引出多种形式的电位变化的图形。PPG(photoplethysmograph)是利用光电容积描记(PPG)技术进行人体运动心率的检测是红外无损检测技术在生物医学中的一个应用.它利用光电传感器,检测经过人体血液和组织吸收后的反射光强度的不同,描记出血管容积在心动周期内的变化,从得到的脉搏波形中计算出心率。
我们这里数据是一组采样频率为250Hz的采样点数据，我们首先看ECG的数据，我们首先读取ECG的数据观察一下：

# ECG 数据处理
ecg_data = []
with open('../Desktop/ECG.txt') as ecg_file:
    lines = ecg_file.readlines()
for line in lines:
    ecg_data.append(int(line))    
sample_rate = 250              #采样频率,每秒250个样本
x = ecg_data
sample_count = 6136        #采样数
t = np.linspace(0,sample_count/sample_rate,sample_count)
# 作图
plt.figure(figsize=(20,10))
ax0 = plt.subplot(211)
#画时域信号
#ax0.set_xlim([0, 5])
ECG_line = ax0.plot(t, x, label='ECG')
ax0.set_xlabel("Time(s)")
ax0.set_ylabel("Amp(μV)")
plt.show()

可以得到ECG的图像如下所示：

可以看到原始的ECG数据比较杂乱，也很不稳定，上下波动也很大，首先我们想到的是将这个振动的波形拉到水平线上，所以我们需要找到这个原始数据波形的包络线，然后将原始数据减去这个包络线，就可以将这个波形拉成水平。我们这里使用低通滤波器将原始数据进行滤波就可以得到包络线，滤波器的参数需要自行调整：

b, a = signal.butter(8, 0.01, 'lowpass')   #配置滤波器 8 表示滤波器的阶数
ecg_outline_data = signal.filtfilt(b, a, ecg_data)  #data为要过滤的信号

得到的原始数据包络线如下所示：

然后我们就可以将原始数据拉平和均值归零，代码如下：

ecg_norm_data = ecg_data - ecg_outline_data
avg = sum(ecg_data)/len(ecg_data)
for i in range(len(ecg_data)):
    ecg_data[i] = ecg_data[i] - avg

然后我们就可以得到拉平之后的数据：

可以从图中看出，整个波形含有大量的高频噪声，我们需要过滤掉高频信号，这里我们使用低通滤波器进行去噪处理：

b, a = signal.butter(8, 0.18, 'lowpass')   #配置滤波器 8 表示滤波器的阶数
ecg_filt_data = signal.filtfilt(b, a, ecg_norm_data)#为要过滤的信号

然后我们可以得到滤波之后的图像：

我们如果想计算ECG的心率，首先我们将滤波之后的时域数据进行傅里叶变换得到频域的数据：

xFFT = np.abs(np.fft.rfft(x)/sample_count)  #快速傅里叶变换
xFreqs = np.linspace(0, sample_rate//2, sample_count//2+1) * 60

然后我们看看ECG的频域图像：

从上图频域的图像可以看到，ECG的波形含有大量的低频噪声成分，很难看到主要的频率成分。所以我们再来看看PPG的数据：

同样的，我们将两组滤波之后的数据放在一起，绿色的是PPG的处理之后的数据，蓝色的是ECG的处理之后的数据，可以从图中发现前面和后面两端的数据并不是很稳定，那么我们就取中间稳定的部分，15-20秒的波形是比较稳定的，我们放大来看看：

我们再来看看PPG做傅里叶变换之后的频域图像：

这张图就很明显地显示PPG的频率主要的成分在75Hz靠右边一点，那么，我们可以计算得到这组波形的频率：

print(int(xFreqs[np.argmax(xFFT, axis=0)]))   

答案是78，那么就可以认为这组波形的心率F是78Hz，周期T就是1/78s。接下来是求解相位，我们来看一下标准的PPG数据图像，我在网上找到的：

图中可以看出PPG数据有两个波谷，分别叫Systolic Peak和Diastolic Peak，再来看我们的PPG数据，对应到我们的PPG数据中：

所谓的相位就是ECG数据的波峰和紧跟着的PPG的Systolic Peak之间的时间间隔，在图中可以表示成：

图中两条虚线所对应的时间差就是相位，我们如何去求这个相位呢？首先波峰就是一个周期内最高的点，波谷是一个周期内最低的点，但是因为波形有上下波动，我们还需要在求解之后进行一个校验的处理：

frequency = 78
T = 60 / frequency
sample_num = int(T * sample_rate)
print(sample_num)
length = 10*sample//sample_num
ecg_cut_data = ecg_filter_data[10*sample_rate:20*sample_rate]
ecg_peak = []
for i in range(length):
    peak = np.argmax(ecg_cut_data[i*sample_num:(i+1)*sample_num], axis=0) / sample_rate + i*sample_num/sample_rate + 10
    ecg_peak.append(peak)
# 检验
for i in range(length-1):
    if ecg_peak[i+1] - ecg_peak[i] < T * 0.7:
        ecg_peak[i+1] = ecg_peak[i] + T
print(ecg_peak)

ppg_cut_data = ppg_filter_data[10*sample_rate:20*sample_rate]
ppg_peak = []
for i in range(length):
    peak = np.argmax(ppg_cut_data[i*sample_num:(i+1)*sample_num], axis=0) / sample_rate+ i*sample_num/sample_rate + 10
    ppg_peak.append(peak)
# 检验
for i in range(length-1):
    if ppg_peak[i+1] - ppg_peak[i] < T * 0.7:
        ppg_peak[i+1] = ppg_peak[i] + T
print(ppg_peak)

得到的结果：

ecg_peak = [10.244, 10.972, 11.748, 12.532, 13.308, 14.052, 14.84, 15.632000000000001, 16.387999999999998, 17.144, 17.92, 18.704, 19.448]
ppg_peak = [10.356, 11.136, 11.868, 12.64, 13.436, 14.187999999999999, 15.2, 15.96923076923077, 16.532, 17.292, 18.024, 18.808, 19.588]

我们来计算十个稳定周期的相位求平均，就是最终的相位值。因为是ECG的波先来，所以我们要做一个判断：

total_diff = 0
count = 0
if ecg_peak[0] <= ppg_peak[0]:
    for i in range(min(len(ecg_peak), len(ppg_peak))):
        total_diff += ppg_peak[i] - ecg_peak[i]
        count += 1
elif ecg_peak[0] > ppg_peak[0]:
    for i in range(min(len(ecg_peak), len(ppg_peak)-1)):
        total_diff += ppg_peak[i+1] - ecg_peak[i]
        count += 1
avg_phase = total_diff/count
print('相位: ', avg_phase)

这里得到的相位差是：0.1619408284023671s。然后我们可以通过专业医疗器械进行血压测量，我们就能够调整频率和相位的权重，可以多测量一些然后求平均值。这篇博文主要是像大家分享一下PPG/ECG数据处理的流程，希望能够给大家一些帮助，谢谢。