延續上篇《用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)
=========================================
執行之後,會計算出這段時間中,共計搏動了幾次心跳,同時也會存一張頻率和時間的對應圖。從結果可以發現用錄影方式來擷取水蚤心跳,很容易受到幀率影響。