一些理論和背景
心率包含許多有關信息。如果擁有心率傳感器和一些數據，那么當然可以購買分析包或嘗試一些可用的開源產品，但是並非所有產品都可以滿足需求。也是這種情況。那么，為什么不嘗試自己做一個人呢？如果正在閱讀本文，那么可能想嘗試一下。本文適合那些有興趣了解更多關於心率分析以及如何使用一些基本模塊用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內標記數據點，找到最高點的位置。
后兩種方法在計算上便宜得多，但精度也較低。這些方法的高精度比曲線擬合方法更加依賴於高采樣率。畢竟，曲線的實際最大值可能會位於兩個數據點之間，而不是位於實際數據點上，如果采樣率較高，則誤差容限將減小。還將使用曲線擬合方法來更精確地近似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值和計算出的列表放在一個可以調用的字典中，並可以附加將在第2部分中計算出的度量。還讓編寫一個包裝函數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（）函數與包裝器分開。這樣，模塊還可以接受已經存儲在內存中的數據幀對象，例如在由另一個程序生成這些數據幀對象時很有用。這使模塊保持靈活。
