Python分析离散心率信号（中）
Python分析離散心率信號（中）
一些理論和背景
心率信號不僅包含有關心臟的信息，還包含有關呼吸，短期血壓調節(jié)，體溫調節(jié)和荷爾蒙血壓調節(jié)（長期）的信息。也（盡管不總是始終如一）與精神努力相關聯(lián)，這并不奇怪，因為大腦是一個非常饑餓的器官，因此消耗了總葡萄糖的25％和氧氣消耗的20％。如果活動增加，心臟需要更加努力地工作以保持其供應。
感興趣的是這些措施可以被分為時間序列數(shù)據(jù)連接頻域數(shù)據(jù)。如果熟悉傅立葉變換，則頻率部分會很有意義。如果不是，請參閱維基百科頁面具有很好的解釋，并且對過程也非常直觀。基本思想是，要獲取隨時間重復的信號（例如心率信號），并確定由哪些頻率構成該信號。將信號從時域轉換到頻域。這是特別喜歡的另一種可視化效果，清楚地顯示了如何將重復信號近似為隨時間重復的不同正弦波之和。
時間序列數(shù)據(jù)
對于心率信號的時間序列部分，主要關注心跳之間的間隔及其隨時間的變化。想要所有R復合物（R1，R2，… Rn）的位置，之間的間隔（RR1，RR2，… RRn，定義為）以及相鄰間隔之間的差異（RRdiff-1，…RRdiff-n，定義為）。
可視化：
在科學文獻中經(jīng)常發(fā)現(xiàn)的時間序列度量是：
§ BPM，即每分鐘的心跳量，在上一部分中進行了計算；
§ IBI（心跳間隔），即心跳之間的平均距離，在前一部分中隱式地將其計算為BPM計算的一部分；
§ SDNN，心跳間隔的標準偏差：
§ SDSD，相鄰RR間隔之間的連續(xù)差的標準偏差：
§ RMSSD，相鄰RR間隔之間的連續(xù)差的均方根：
§ pNN50 / pNN20，比例差異大于50ms / 20ms。
IBI，SDNN，SDSD，RMSSD和pNNx（以及頻域量度）通常歸為“心率變異性”（HRV）量度，因為可提供有關心率隨時間變化的信息。
頻域數(shù)據(jù)
在心率信號的頻率側，最常發(fā)現(xiàn)的量度稱為HF（高頻），MF（中頻）和LF（低頻）頻段，這是對創(chuàng)新性命名水平的永恒證明。MF en HF頻段通常合并在一起，并被標記為“ HF”。LF和HF大致對應于LF頻段的0.04-0.15Hz和HF頻段的0.16-0.5Hz。LF頻段似乎與短期血壓變化有關，HF頻段與呼吸頻率有關。
通過在RR間隔數(shù)據(jù)序列上執(zhí)行快速傅立葉變換來計算頻譜。顧名思義，該方法與離散傅立葉變換方法相比是快速的。數(shù)據(jù)集越大，方法之間的速度差異就越大。
首先對信號進行重新采樣，以計算頻譜，然后將重新采樣的信號轉換到頻域，然后以給定的間隔積分曲線下的面積，從而計算出HF和LF的量度。
時域度量–入門
查看上述時間序列度量，需要一些成分才能輕松計算所有這些度量：
§ 所有R峰的位置列表；
§ 所有后續(xù)RR對之間的間隔的列表（RR1，RR2，… RR-n）；
§ RR對之間所有后續(xù)間隔之間的差異的列表（RRdiff1，…RRdiffn）；
§ RR對之間所有后續(xù)差之間的平方差的列表。
已經(jīng)從第一部分的detect_peaks（）函數(shù)獲得了所有R峰位置的列表，該列表包含在dict [‘peaklist’]中。還有來自calc_RR（）函數(shù)的RR對差異列表，位于dict [‘RR_list’]中。大！無需編寫其代碼，已經(jīng)完成了50％的編寫。
為了獲得最后兩個成分，使用dict [‘RR_list’]并計算差值和相鄰值之間的平方差：
RR_diff = []
RR_sqdiff = []
RR_list = measures[‘RR_list’]
cnt = 1 #Use counter to iterate over RR_list
while (cnt < (len(RR_list)-1)): #Keep going as long as there are R-R intervals
RR_diff.append(abs(RR_list[cnt] - RR_list[cnt+1])) #Calculate
absolute difference between successive R-R interval
RR_sqdiff.append(math.pow(RR_list[cnt] - RR_list[cnt+1], 2)) #Calculate squared differencecnt += 1
print(RR_diff, RR_sqdiff)
復制
計算時域度量值現(xiàn)在有了所有成分，可以輕松計算所有度量值：
ibi = np.mean(RR_list) #Take the mean of RR_list to get the mean Inter Beat Interval
print(“IBI:”, ibi)
sdnn = np.std(RR_list) #Take standard deviation of all R-R intervals
print(“SDNN:”, sdnn)
sdsd = np.std(RR_diff) #Take standard deviation of the differences between all
subsequent R-R intervals
print(“SDSD:”, sdsd)
rmssd = np.sqrt(np.mean(RR_sqdiff)) #Take root of the mean of the list of squared differences
print(“RMSSD:”, rmssd)
nn20 = [x for x in RR_diff if (x>20)] #First create a list of all values over 20, 50
nn50 = [x for x in RR_diff if (x>50)]
pnn20 = float(len(NN20)) / float(len(RR_diff)) #Calculate the proportion of NN20, NN50 intervals to all
intervals
pnn50 = float(len(NN50)) / float(len(RR_diff)) #Note the use of float(), because we don’t want Python to
think we want an int() and round the proportion to 0 or 1
print(“pNN20, pNN50:”, pnn20, pnn50)
復制
時間序列度量就是這樣。讓將包裝在可調用函數(shù)中。擴展calc_RR（）函數(shù)以計算額外成分，并將其附加到字典對象中，還將calc_bpm（）與其時間序列測量值合并到一個新函數(shù)calc_ts_measures（）中，并將附加到字典中：
def calc_RR(dataset, fs):
peaklist = measures['peaklist']RR_list = []cnt = 0while (cnt < (len(peaklist)-1)):RR_interval = (peaklist[cnt+1] - peaklist[cnt])ms_dist = ((RR_interval / fs) * 1000.0)RR_list.append(ms_dist)cnt += 1RR_diff = []RR_sqdiff = []cnt = 0while (cnt < (len(RR_list)-1)):RR_diff.append(abs(RR_list[cnt] - RR_list[cnt+1]))RR_sqdiff.append(math.pow(RR_list[cnt] - RR_list[cnt+1], 2))cnt += 1measures['RR_list'] = RR_listmeasures['RR_diff'] = RR_diffmeasures['RR_sqdiff'] = RR_sqdiff
def calc_ts_measures():
RR_list = measures['RR_list']RR_diff = measures['RR_diff']RR_sqdiff = measures['RR_sqdiff']
measures[‘bpm’] = 60000 / np.mean(RR_list)
measures['ibi'] = np.mean(RR_list)measures['sdnn'] = np.std(RR_list)measures['sdsd'] = np.std(RR_diff)measures['rmssd'] = np.sqrt(np.mean(RR_sqdiff))NN20 = [x for x in RR_diff if (x>20)]NN50 = [x for x in RR_diff if (x>50)]measures['nn20'] = NN20measures['nn50'] = NN50measures['pnn20'] = float(len(NN20)) / float(len(RR_diff))measures['pnn50'] = float(len(NN50)) / float(len(RR_diff))
#Don’t forget to update our process() wrapper to include the new function
def process(dataset, hrw, fs):
rolmean(dataset, hrw, fs)detect_peaks(dataset)calc_RR(dataset, fs)calc_ts_measures()
復制
這樣稱呼：
import heartbeat as hb #Assuming we
named the file ‘heartbeat.py’
dataset = hb.get_data(‘data.csv’)
hb.process(dataset, 0.75, 100)
#The module dict now contains all the variables computed over our signal:
hb.measures[‘bpm’]
hb.measures[‘ibi’]
hb.measures[‘sdnn’]
#etcetera
#Remember that you can get a list of all dictionary entries with “keys()”:
print(hb.measures.keys())
復制
在計算機上，分析信號并計算單個線程中的所有度量大約需要25毫秒。
頻域度量–入門頻域度量
計算比較棘手。主要原因是不想要將心率信號轉換到頻域（這樣做只會返回等于BPM / 60（以Hz表示的心跳）的強頻率）。相反，希望將RR間隔轉換到頻域。很難理解？這樣思考：隨著身體需求的變化，心率隨著時間的變化而變化。這種變化表現(xiàn)為心跳之間的距離隨時間變化（之前計算的RR間隔）。RR峰之間的距離隨其自身頻率隨時間變化。為了可視化，繪制之前計算的RR間隔：
從圖中可以清楚地看到，RR間隔不會隨心跳而劇烈變化，而是隨時間呈正弦波狀變化（更準確地說是不同正弦波的組合）。想找到組成該模式的頻率。
但是，任何傅立葉變換方法都依賴于均勻間隔的數(shù)據(jù)，并且RR間隔在時間上肯定不是均勻間隔的。這是因為間隔的時間位置取決于長度，每個間隔都不同。希望這是有道理的。
要找到措施，需要：
§ 創(chuàng)建一個帶有RR間隔的均勻間隔的時間線；
§ 內插信號，既可以創(chuàng)建均勻間隔的時間序列，又可以提高分辨率；
§ 在某些研究中，此插值步驟也稱為重新采樣。
§ 將信號轉換到頻域；
§ 積分頻譜的LF和HF部分下方的面積。
計算頻域量度
首先，為RR間隔創(chuàng)建均勻間隔的時間線。為此，獲取所有R峰的樣本位置，這些樣本位置位于在第1部分中計算的列表dict [‘peaklist’]中。然后，將y值分配給列表dict [‘中的這些樣本位置RR_list’]，其中包含所有RR間隔的持續(xù)時間。最后，對信號進行插值。
from scipy.interpolate import interp1d #Import the
interpolate function from SciPy
peaklist = measures[‘peaklist’] #First retrieve the lists we need
RR_list = measures[‘RR_list’]
RR_x = peaklist[1:] #Remove the first entry, because first interval is
assigned to the second beat.
RR_y = RR_list #Y-values are equal to interval lengths
RR_x_new = np.linspace(RR_x[0],RR_x[-1],RR_x[-1]) #Create evenly spaced timeline starting at the second peak, its endpoint and length equal to position of last peak
f = interp1d(RR_x, RR_y, kind=‘cubic’) #Interpolate the signal with cubic spline interpolation
復制
請注意，時間序列不是從第一個峰值開始，而是從第二個R峰值的采樣位置開始。因為使用間隔，所以第一個間隔在第二個峰值可用。
現(xiàn)在，可以使用創(chuàng)建的函數(shù)f（）查找信號范圍內任何x位置的y值：
print f(250)
#Returns 997.619845418, the Y value at x=250
復制
同樣，可以將整個時間序列RR_x_new傳遞給該函數(shù)并進行繪制：
plt.title(“Original and Interpolated Signal”)
plt.plot(RR_x, RR_y, label=“Original”, color=‘blue’)
plt.plot(RR_x_new, f(RR_x_new), label=“Interpolated”, color=‘red’)
plt.legend()
plt.show()
復制
現(xiàn)在，要查找組成插值信號的頻率，請使用numpy的快速傅立葉變換np.fft.fft（）方法，計算采樣間隔，將采樣倉轉換為Hz并作圖：
#Set variables
n = len(dataset.hart) #Length of the
signal
frq = np.fft.fftfreq(len(dataset.hart), d=((1/fs))) #divide the bins into frequency categories
frq = frq[range(n/2)] #Get single side of the frequency range
#Do FFT
Y = np.fft.fft(f(RR_x_new))/n #Calculate FFT
Y = Y[range(n/2)] #Return one side of the FFT
#Plot
plt.title(“Frequency Spectrum of Heart Rate Variability”)
plt.xlim(0,0.6) #Limit X axis to frequencies of interest (0-0.6Hz for
visibility, we are interested in 0.04-0.5)
plt.ylim(0, 50) #Limit Y axis for visibility
plt.plot(frq, abs(Y)) #Plot it
plt.xlabel(“Frequencies in Hz”)
plt.show()
復制
可以清楚地看到信號中的LF和HF頻率峰值。
剩下的最后一件事是對LF（0.04 – 0.15Hz）和HF（0.16 – 0.5Hz）頻帶下的曲線下面積進行積分。需要找到與感興趣的頻率范圍相對應的數(shù)據(jù)點。在FFT期間，計算了單側頻率范圍frq，因此可以在其中搜索所需的數(shù)據(jù)點位置。
lf = np.trapz(abs(Y[(frq>=0.04) & (frq<=0.15)])) #Slice frequency spectrum where x is between 0.04 and
0.15Hz (LF), and use NumPy’s trapezoidal integration function to find the area
print(“LF:”, lf)
hf = np.trapz(abs(Y[(frq>=0.16) & (frq<=0.5)])) #Do the same for 0.16-0.5Hz (HF)
print(“HF:”, hf)
復制
返回：
LF: 38.8900414093
HF: 47.8933830871
復制
這些是感興趣頻譜處頻譜下的區(qū)域。請記住，從理論上講，HF與呼吸有關，而LF與短期血壓調節(jié)有關。這些措施也與增加精神活動有關。
全面包裝
已經(jīng)走了很長一段路，現(xiàn)在可以從心率信號中提取許多有意義的指標。但是，如果輸入其心率數(shù)據(jù)，則該模塊很可能會失敗，因為可能比理想數(shù)據(jù)更嘈雜，或者可能包含偽像。將處理信號過濾，錯誤檢測和離群值剔除。
總結
以上是生活随笔為你收集整理的Python分析离散心率信号（中）的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： Python分析离散心率信号（上）
- 下一篇： Python分析离散心率信号（下）