2016年10月13日

python分析水蚤心跳數與頻率


延續上篇《用Tracker分析水蚤心跳》,繼續來談分析。
上回利用Tracker分析水蚤心跳的影片之後,得到這張圖。
水蚤心跳


接下來想問的就是,心臟搏動了幾下,心跳頻率又是怎樣?
如果要算數量,應該就是做出上面那張圖,然後一個峰一個峰算吧?算個這樣一張還好,如果今天要算的是十秒鐘內的心跳數,或是整整一分鐘的心跳數呢?有沒有簡單一點的算法?

有的,那就要寫程式去運算囉。
寫程式之前,要先教電腦判斷哪個叫做峰,這用肉眼看很明顯,但是電腦讀的是數字,是看不出這圖的。

所謂峰,就是一端高高,兩邊的比它還矮。知道這個概念之後,就是把它轉變成電腦能懂的語言:從0秒開始去掃每個時間點的亮度值,每次都去比看看,這個值是不是比前後兩個值都大,如果是,那這個點就是峰。

而要計算頻率,就是把每個峰出現的時間點紀錄下來,每一個都去減掉前一個峰的時間,這樣得出週期的時間,接著再取倒數,就可以得到頻率了。

想過這些過程之後,就是開始實做,我用python2來作,用了兩個初次使用的lib,分別是numpy和matplotlib。


程式碼如下,執行前先照著前篇文章的影片說明,把時間與亮度的對應值存成文字檔(這裡用heartBeat為例),接著就執行了。

=========================================

# -*- coding: utf-8 -*-

import numpy as ny
import matplotlib.pyplot as plt

a = ny.loadtxt('heartBeat')
dataN=a.shape[0]
HBduration=0
HBfeq=0
preHBtime=0
heartBeatNum=0
xlist=[]
ylist=[]

#for k in range(0+1,dataN-3,1):
for k in range(0+1,1784,1):
    if (a[k+1,1]<a[k,1] and a[k-1,1]<a[k,1]) or \
    (a[k+2,1]<a[k,1] and a[k-1,1]<a[k,1]) :
        HBduration=a[k,0]-preHBtime #count HeatBeat duration
        HBfeq=1/HBduration
        preHBtime=a[k,0]

        xlist.append(a[k,0])
        ylist.append(HBfeq)      
        heartBeatNum+=1
        #print(a[k,0])      
        #print(a[k,1])
     
print(heartBeatNum)
#print(xlist)
#print(ylist)


plt.plot(xlist,ylist)
#plt.xlim(0,120)
plt.xlabel('time(sec)')
plt.ylabel('Fequency(Hz)')
plt.title('Heartbeat Frequency of a Waterflea')
plt.savefig('frequency.png',dpi=300)


=========================================

執行之後,會計算出這段時間中,共計搏動了幾次心跳,同時也會存一張頻率和時間的對應圖。從結果可以發現用錄影方式來擷取水蚤心跳,很容易受到幀率影響。



frequency