Python分析离散心率信号（上）
Python分析離散心率信號（上）
一些理論和背景
心率包含許多有關信息。如果擁有心率傳感器和一些數據，那么當然可以購買分析包或嘗試一些可用的開源產品，但是并非所有產品都可以滿足需求。也是這種情況。那么，為什么不嘗試自己做一個人呢？如果正在閱讀本文，那么可能想嘗試一下。本文適合那些有興趣了解更多關于心率分析以及如何使用一些基本模塊用Python編寫簡單但有效的分析算法的人。在談論心率信號的形狀時，每個搏動都具有QRS復數的特征，如a。，I-III（QRS）所示。還有一個P波（IV）和一個T波（V）。R成分是大峰（II）。這是感興趣的一個each beat is characterized by a QRS-complex as shown below in
a.), I-III (Q-R-S). There is also a P-wave (IV), and a T-wave (V). The R
component is the large peak (II). This is the one of interest to us:：
在光體積描記圖信號（PPG）中，信號略有不同。PPG信號顯示在圖像的b。）中。這里的組成部分是舒張壓峰值（I）和舒張壓峰值（III），舒張壓峰值是最高的血壓點。這兩個波被所謂的重癥陷波（II）分開。在這種情況下，收縮峰（I）用于心率提取。
為了開發峰值檢測算法，這是一個很好的巧合：對于兩種信號，都需要標記每個復合物中的最高峰。
首先，首先讓下載數據集并繪制信號圖，以便對數據有所了解并開始尋找有意義地分析數據的方法。將熊貓用于大多數數據任務，并將matplotlib用于大多數繪圖需求。
import pandas as pd import matplotlib.pyplot as plt
dataset = pd.read_csv(“data.csv”) #Read data
from CSV datafile
plt.title(“Heart Rate Signal”) #The title of our plot
plt.plot(dataset.hart) #Draw the plot object
plt.show() #Display the plot
復制
信號看起來很干凈。請記住，即使使用良好的傳感器，也并非總是如此，尤其是在實驗室外進行測量時！以后會分析如何處理噪聲和自動確定信號質量。
檢測第一個峰
第一步是找到所有R峰的位置。為此，需要確定感興趣區域（ROI），即信號中的每個R峰。得到這些之后，需要確定最大值。有幾種方法可以做到這一點：
§ 在ROI數據點上擬合曲線，求解最大值的x位置；
§ 確定ROI中每組點之間的斜率，找到斜率反轉的位置；
§ 在ROI內標記數據點，找到最高點的位置。
第一個提供了精確的數學解決方案（更精確：根據需要精確），但是也是計算上最昂貴的。對于示例數據集，差異并不重要，但是在開發此算法時，不得不快速遍歷近十億個數據點，而未來的數據集似乎只會變得更大。大部分代碼還將被重寫為C以便在嵌入式ARM芯片上使用，因此需要保持其簡單高效。
后兩種方法在計算上便宜得多，但精度也較低。這些方法的高精度比曲線擬合方法更加依賴于高采樣率。畢竟，曲線的實際最大值可能會位于兩個數據點之間，而不是位于實際數據點上，如果采樣率較高，則誤差容限將減小。還將使用曲線擬合方法來更精確地近似R峰，因為現在，將ROI中最高點的位置用作拍子的位置。
現在開始工作：首先將不同的峰彼此分離。為此，繪制一個移動平均線，標記心率信號位于移動平均線上方的ROI，并最終在每個ROI中找到最高點，如下所示：
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import math
dataset = pd.read_csv(“data.csv”)
#Calculate moving average with 0.75s
in both directions, then append do dataset
hrw = 0.75 #One-sided
window size, as proportion of the sampling frequency
fs = 100 #The example
dataset was recorded at 100Hz
mov_avg = dataset[‘hart’].rolling(int(hrw*fs)).mean() #Calculate
moving average
#Impute where moving average function
returns NaN, which is the beginning of the signal where x hrw
avg_hr = (np.mean(dataset.hart))
mov_avg = [avg_hr if math.isnan(x) else x for x in mov_avg]
mov_avg = [x*1.2 for x in mov_avg] #For now we
raise the average by 20% to prevent the secondary heart contraction from
interfering, in part 2 we will do this dynamically
dataset[‘hart_rollingmean’] = mov_avg #Append the
moving average to the dataframe
#Mark regions of interest
window = []
peaklist = []
listpos = 0 #We use a
counter to move over the different data columns
for datapoint in dataset.hart:
rollingmean = dataset.hart_rollingmean[listpos] #Get local
mean
if (datapoint < rollingmean) and (len(window) < 1): #If no
detectable R-complex activity -> do nothing
listpos += 1
elif (datapoint > rollingmean): #If signal
comes above local mean, mark ROI
window.append(datapoint)
listpos += 1
else: #If signal
drops below local mean -> determine highest point
maximum = max(window)
beatposition = listpos - len(window) + (window.index(max(window))) #Notate the
position of the point on the X-axis
peaklist.append(beatposition) #Add detected
peak to list
window = [] #Clear marked
ROI
listpos += 1
ybeat = [dataset.hart[x] for x in peaklist] #Get the
y-value of all peaks for plotting purposes
plt.title(“Detected
peaks in signal”)
plt.xlim(0,2500)
plt.plot(dataset.hart, alpha=0.5, color=‘blue’) #Plot
semi-transparent HR
plt.plot(mov_avg, color =‘green’) #Plot moving
average
plt.scatter(peaklist, ybeat, color=‘red’) #Plot detected
peaks
plt.show()
復制
已經在信號中標記了每個R復合體的最高點，還不錯！但是，這是一個理想的信號。將了解如何處理噪聲，測量誤差，檢測誤差以及如何防止模塊在質量差的信號部分拋出誤差。
有人可能會問：為什么要移動平均線，為什么不僅僅使用700左右的水平線作為閾值？這是一個有效的問題，可以很好地處理此信號。為什么使用移動平均線與理想化程度較低的信號有關。R峰值的幅度會隨時間變化，尤其是當傳感器稍微移動時。較小的次要峰的振幅也可以獨立于R峰的振幅而變化，有時具有幾乎相同的振幅。減少錯誤檢測到的峰的一種方法是通過擬合不同高度處的移動平均值并確定哪種擬合效果最佳。稍后再詳細介紹。
計算心率
知道每個峰值在時間上的位置，因此計算此信號的平均“每分鐘心跳數”（BPM）度量很簡單。只需計算峰之間的距離，取平均值并轉換為每分鐘的值，如下所示：
RR_list = []
cnt = 0
while (cnt < (len(peaklist)-1)):
RR_interval = (peaklist[cnt+1] - peaklist[cnt]) #Calculate distance
between beats in # of samples
ms_dist = ((RR_interval / fs) * 1000.0) #Convert
sample distances to ms distances
RR_list.append(ms_dist) #Append to
list
cnt += 1
bpm = 60000 / np.mean(RR_list) #60000 ms (1
minute) / average R-R interval of signal
print “Average
Heart Beat is: %.01f” %bpm #Round off to
1 decimal and print
復制
?
還要更新繪圖方法以在圖例中顯示BPM：
plt.title(“Detected
peaks in signal”)
plt.xlim(0,2500)
plt.plot(dataset.hart, alpha=0.5, color=‘blue’, label=“raw
signal”) #Plot semi-transparent HR
plt.plot(mov_avg, color =‘green’, label=“moving
average”) #Plot moving average
plt.scatter(peaklist, ybeat, color=‘red’, label=“average:
%.1f BPM” %bpm) #Plot detected
peaks
plt.legend(loc=4, framealpha=0.6)
plt.show()
復制
已邁出分析心率信號的第一步。每分鐘的心跳數是一種非常有用的量度，在科學研究中經常使用，并且還有許多非科學用途，但是信號包含的信息要多得多。將處理從心率信號中提取更復雜的信息。
四舍五入
最后，整理一下代碼并將其放入可調用函數中。這將使生活在下一部分變得更加輕松，并且代碼將更加有條理和可重用。請注意，可能要做的是整潔的事情，使函數成為類的一部分，但為了使本教程也可供那些不太熟悉Python（并且可能對類不熟悉或不熟悉）的人訪問，選擇從此處省略本教程系列中的所有代碼。
讓將BPM值和計算出的列表放在一個可以調用的字典中，并可以附加計算出的度量。還讓編寫一個包裝函數process（），以便可以使用盡可能少的代碼來調用分析：
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import math
measures = {}
def get_data(filename):
dataset = pd.read_csv(filename)
return dataset
def rolmean(dataset, hrw, fs):
mov_avg = dataset[‘hart’].rolling(int(hrw*fs)).mean()
avg_hr= (np.mean(dataset.hart))
mov_avg = [avg_hr if math.isnan(x) else x for x in mov_avg]
mov_avg = [x*1.2 for x in mov_avg]
dataset[‘hart_rollingmean’] = mov_avg
def detect_peaks(dataset):
window = []
peaklist = []
listpos = 0
for datapoint in dataset.hart:
rollingmean = dataset.hart_rollingmean[listpos]
if (datapoint < rollingmean) and (len(window) < 1):
listpos += 1
elif (datapoint > rollingmean):
window.append(datapoint)
listpos += 1
else:
maximum = max(window)
beatposition = listpos - len(window) + (window.index(max(window)))
peaklist.append(beatposition)
window = []
listpos += 1
measures[‘peaklist’] = peaklist
measures[‘ybeat’] = [dataset.hart[x] for x in peaklist]
def calc_RR(dataset, fs):
RR_list = []
peaklist = measures[‘peaklist’]
cnt = 0
while (cnt < (len(peaklist)-1)):
RR_interval = (peaklist[cnt+1] - peaklist[cnt])
ms_dist = ((RR_interval / fs) * 1000.0)
RR_list.append(ms_dist)
cnt += 1
measures[‘RR_list’] = RR_list
def calc_bpm():
RR_list = measures[‘RR_list’]
measures[‘bpm’] = 60000 / np.mean(RR_list)
def plotter(dataset, title):
peaklist = measures[‘peaklist’]
ybeat = measures[‘ybeat’]
plt.title(title)
plt.plot(dataset.hart, alpha=0.5, color=‘blue’, label=“raw signal”)
plt.plot(dataset.hart_rollingmean, color =‘green’, label=“moving average”)
plt.scatter(peaklist, ybeat, color=‘red’, label=“average: %.1f BPM” %measures[‘bpm’])
plt.legend(loc=4, framealpha=0.6)
plt.show()
def process(dataset, hrw, fs): #Remember; hrw was the one-sided window size (we used 0.75) and fs was the sample rate (file is recorded at 100Hz)
rolmean(dataset, hrw, fs)
detect_peaks(dataset)
calc_RR(dataset, fs)
calc_bpm()
plotter(dataset, “My Heartbeat Plot”)
復制
這樣調用：
import heartbeat as hb #Assuming we named the file ‘heartbeat.py’
dataset = hb.get_data(“data.csv”)
hb.process(dataset, 0.75, 100)
#We have imported our Python module as
an object called ‘hb’
#This object contains the dictionary
‘measures’ with all values in it
#Now we can also retrieve the BPM
value (and later other values) like this:
bpm = hb.measures[‘bpm’]
#To view all objects in the
dictionary, use “keys()” like so:
print hb.measures.keys()
復制
請注意，將get_data（）函數與包裝器分開。這樣，模塊還可以接受已經存儲在內存中的數據幀對象，例如在由另一個程序生成這些數據幀對象時很有用。這使模塊保持靈活。
總結
以上是生活随笔為你收集整理的Python分析离散心率信号（上）的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 适用于Windows和Linux的Yol
- 下一篇： Python分析离散心率信号（中）