Python分析離散心率訊號（下）
如何使用動態閾值，訊號過濾和離群值檢測來改善峰值檢測。
一些理論和背景
到目前為止，一直在研究如何分析心率訊號並從中提取最廣泛使用的時域和頻域度量。但是，使用的訊號是理想的。現在考慮這個訊號：
一個挑戰！這是遇到的訊號質量的另一個極端。老實說，當將感測器連線到手指上時（在0到4000之間），透過測量產生了該訊號。在此之後，手指中的血管需要立即適應感測器的壓縮（大約4000-5000），此後訊號變得穩定。在大約7500、9000和12000時，用力將感測器移到手指上方以產生額外的噪音。對感測器本身能很好地抑制噪聲感到驚訝。
儘管此訊號已被手動“破壞”，但實際上可能會發生部分資料包含噪音或偽影的情況。感測器可能會移動，從而產生噪聲，在測量過程中可能未正確連線或分離，參與者可能打噴嚏，移動或其任何噪聲誘發因素都可能干擾。
在幾個階段中看到如何處理此問題：
§ 開始，將該訊號傳遞給演算法的結果；
§ 注意：取樣頻率；
§ 過濾訊號以去除不想要的頻率（噪聲）；
§ 動態閾值提高檢測精度；
§ 檢測錯誤/錯過的峰；
§ 消除錯誤並將RR訊號重構為無錯誤。
開始啟動
首先，讓看看演算法如何處理該訊號。
（。。。） line 37， in detect_peaks
maximum = max（window）
ValueError： max（） arg is an empty sequence
很好，事實並非如此。怎麼了？
有人可能已經注意到了這一點，
detect_peaks（）
函式將在一種情況下中斷：當心率訊號從小於移動平均線變為變為等於而連續至少兩個資料點都沒有移動到其上方時，發生這種情況的最可能情況是訊號在一段時間內
下降到零
。然後，該函式將跳到
else
語句，並嘗試檢測ROI中的峰值，但是由於訊號從未上升到移動平均線之上，因此未標記ROI，因此
max（window）
丟擲錯誤，因為空列表沒有最大值。有時候是這樣的因為是個白痴因為在迴圈中，沒有考慮到兩個訊號的高度相等的情況，而這種情況恰好發生在訊號的早期，將感測器固定在手指上：
如果之前在程式碼中沒有注意到這一點，不用擔心，一開始也沒有。現在要更新
detect_peaks（）
函式：
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）： #Here is the update in （datapoint <= rollingmean）
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］
現在，還註釋掉
第19行
，稍後再討論。
def rolmean（dataset， hrw， fs）：
mov_avg = pd。rolling_mean（dataset。hart， window=（hrw*fs））
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
然後使用峰值檢測再次繪圖。
該模組不再丟擲錯誤並檢測峰值，但這幾乎不是想要的結果。讓從降低噪聲開始，看看是否可以改善訊號。
取樣頻率
在進行訊號過濾之前，讓確定檔案的取樣率。前兩部分的檔案為100Hz，這部分呢？實際上，實際的記錄頻率可能隨不同的裝置或系統而變化。也可能與裝置的額定值略有不同。所有計算得出的度量的準確性取決於準確瞭解取樣率，因此進行檢查非常重要。如果組合各種來源的資料，即使100Hz和101Hz的差異也可能導致明顯不同的結果。通常在每條資料行中記錄“世界時間”或“自記錄開始以來的時間”，以便隨後計算和驗證取樣率。
這樣就可以很容易地計算出取樣率：
#Simple way to get sample rate
sampletimer = ［x for x in dataset。timer］ #dataset。timer is a ms counter with start of recording at ‘0’
measures［‘fs’］ = （（len（sampletimer） / sampletimer［-1］）*1000） #Divide total length of dataset by last timer entry。 This is in ms， so multiply by 1000 to get Hz value
#If your timer is a date time string， convert to UNIX timestamp to more easily calculate with， use something like this：
unix_time = ［］
for x in dataset。datetime：
dt = datetime。datetime。strptime（Datum， “%Y-%m-%d %H：%M：%S。%f”）
unix_time。append（time。mktime（dt。timetuple（）） + （dt。microsecond / 1000000。0））
measures［‘fs’］ = （len（unix_time） / （unix_time［-1］ - unix_time［0］））
此處提供的檔案的取樣率實際上是117Hz！現在可以使用確定的取樣率自動計算度量。
請注意
，這還不是全部，除了取樣率之外，還應該檢查資料以瞭解取樣範圍。是否所有樣本均等地間隔開，例如，資料流中是否沒有“間隙”或“跳過”？如果資料包含間隙和跳躍（如果最多隻有幾個樣本），請在處理之前考慮對缺失值進行插值。如果取樣率隨時間變化很大，那麼請使用其資料記錄裝置。
現在更加準確地瞭解了取樣頻率，接下來可以進行訊號濾波。
訊號過濾
遇到偽像或嘈雜的訊號時，應該做的第一件事就是嘗試過濾訊號。為什麼？因為在任何現實的測量情況下，訊號（無論可能是什麼）都將由
訊號部分
和
誤差部分組成
。這是因為不存在理想的感測器，因此將始終拾取來自要測量的源以外的干擾。
電力線噪聲
是常見干擾的一個例子。在大多數國家/ 地區中，來自牆壁觸點的交流電源的頻率為50Hz（某些國家（例如美國）使用60Hz）。在許多原始ECG測量中也存在這種噪聲。
巴特沃斯（Butterworth）濾波器
一種常用的降低噪音的濾波器，其特點是對指定範圍內的頻率響應非常均勻。充當“頻率門”；抑制超出指定截止範圍的頻率。這個臨界點不是一條絕對線。意思是，如果將截止頻率設定為例如5Hz，則5。1Hz的訊號仍將透過濾波器，其幅度將稍微減小（或“音量”，如果更有意義）。另一方面，10Hz的訊號只會非常弱地透過（如果有的話）。這也取決於過濾器的順序。
還是很難理解？這樣考慮：正在參加一個聚會，並且同時與兩個人交談，導致無法理解其中任何一個的情況。現在在和兩個人之間放置一個過濾器。過濾器將減少人1的講話音量而不會改變人2的音量。現在可以很好地理解人2了。這是濾波器的工作，但頻率除外。
無論如何，開始編碼，看看訊號是否可以從濾波中受益。首先下載並開啟，如果還沒有這麼做過的資料集，透過定義過濾器
scipy.signal
，過濾，最後繪製訊號。假設正在使用上一部分中的程式碼，則按如下方式定義過濾器和繪圖：
from scipy。signal import butter， lfilter #Import the extra module required
#Define the filter
def butter_lowpass（cutoff， fs， order=5）：
nyq = 0。5 * fs #Nyquist frequeny is half the sampling frequency
normal_cutoff = cutoff / nyq
b， a = butter（order， normal_cutoff， btype=‘low’， analog=False）
return b， a
def butter_lowpass_filter（data， cutoff， fs， order）：
b， a = butter_lowpass（cutoff， fs， order=order）
y = lfilter（b， a， data）
return y
dataset = get_data（‘data2。csv’）
dataset = dataset［6000：12000］。reset_index（drop=True） #For visibility take a subselection of the entire signal from samples 6000 - 12000 （00：01：00 - 00：02：00）
filtered = butter_lowpass_filter（dataset。hart， 2。5， 100。0， 5）#filter the signal with a cutoff at 2。5Hz and a 5th order Butterworth filter
#Plot it
plt。subplot（211）
plt。plot（dataset。hart， color=‘Blue’， alpha=0。5， label=‘Original Signal’）
plt。legend（loc=4）
plt。subplot（212）
plt。plot（filtered， color=‘Red’， label=‘Filtered Signal’）
plt。ylim（200，800） #limit filtered signal to have same y-axis as original （filter response starts at 0 so otherwise the plot will be scaled）
plt。legend（loc=4）
plt。show（）
現在，這個訊號似乎並沒有太大改善。如果仔細觀察，線條會更平滑一些，但是一開始並沒有很多高振幅，高頻噪聲。甚至可能會爭辯說，濾波會稍微降低具有較低頻率噪聲的部分，因為在那兒也會抑制R峰。這是一個很好的例子，說明了為什麼在決定過濾
訊號
時
總是要繪製訊號
。過濾訊號會改變，這取決於決定此改變是否更好。
在另一個專案中使用過的這種嘈雜訊號就是一個Butterworth濾波器非常有用的例子：
毫無疑問，這可以改善訊號以進一步處理。
透過動態閾值提高檢測精度
減少次峰標記不正確的第一種也是最明顯的方法是提高用作閾值的移動平均值。但是提高到什麼水平？對於許多不同的訊號，這將是不同的。需要採取措施來幫助確定哪個閾值水平可能是最準確的。
一些簡單的措施可以提供幫助，可以：
§ 檢視訊號長度並計算有多少個峰與預期的峰數；
§ 確定並使用RR間隔的標準偏差（將其
稱為RRSD
）。
檢測到的峰值數量可保留有價值的資訊，但只能用來去掉明顯不正確的閾值。根據應用程式（大多數人都坐著不動），可以過濾不太可能的bpm值。例如，過濾閾值導致平均bpm <30bpm和> 130bpm。在不同情況下（體育鍛煉），閾值可以不同。
RRSD告訴有關所有RR間隔之間的差異如何分散的一些資訊。通常，如果沒有太多的噪聲，則最低RRSD
不為零
且可信的BPM值的檢測將是最合適的。RRSD不能為零，因為這意味著沒有心率變異性，這強烈表明了R峰檢測中的錯誤。
這種簡單的方法之所以有效，是因為缺少節拍會導致RR間隔大約是平均RR間隔的兩倍，而錯誤地標記一個節拍會導致RR間隔
最多
是平均RR間隔的一半左右，但通常較小。兩種情況都會導致
RRSD升高
。本質上，利用心率訊號包含一個相當穩定的重複模式這一事實。
為了說明這一點，在資料集的子選擇中繪製四個峰值檢測回合，移動平均值分別提高0％，10％，25％和35％（從上到下）。
在倒數第二個圖中，所有R峰均被正確檢測到，沒有任何東西被錯誤地標記為R峰。請注意，儘管BPM本身可能對所有四個圖都有效（都不是絕對混亂），但RRSD強烈指出了最正確的圖。說“最正確”是因為在某些情況下，無論如何設定閾值，都會殘留一些錯誤，稍後將對此進行更多介紹。還要注意，與第三個圖相比，最後一個圖中缺少單個峰如何導致RRSD大大增加。
現在如何實施呢？不能簡單地執行10。000個變化，每個變化均具有略高的移動平均線。除此之外，這將使演算法整體效能嚴重受損，也將是非常多餘的，因為許多迭代會產生相同的正確結果（而其許多迭代也會產生相同的錯誤結果！）。可能的解決方案是定期檢查，然後評估最可能的RRSD和BPM對，如下所示：
def detect_peaks（dataset， ma_perc， fs）： #Change the function to accept a moving average percentage ‘ma_perc’ argument
rolmean = ［（x+（（x/100）*ma_perc）） for x in dataset。hart_rollingmean］ #Raise moving average with passed ma_perc
window = ［］
peaklist = ［］
listpos = 0
for datapoint in dataset。hart：
rollingmean = rolmean［listpos］
if （datapoint <= rollingmean） and （len（window） <= 1）： #Here is the update in （datapoint <= rollingmean）
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］
measures［‘rolmean’］ = rolmean
calc_RR（dataset， fs）
measures［‘rrsd’］ = np。std（measures［‘RR_list’］）
def fit_peaks（dataset， fs）：
ma_perc_list = ［5， 10， 15， 20， 25， 30， 40， 50， 60， 70， 80， 90， 100， 110， 120， 150， 200， 300］ #List with moving average raise percentages， make as detailed as you like but keep an eye on speed
rrsd = ［］
valid_ma = ［］
for x in ma_perc_list： #Detect peaks with all percentages， append results to list ‘rrsd’
detect_peaks（dataset， x， fs）
bpm = （（len（measures［‘peaklist’］）/（len（dataset。hart）/fs））*60）
rrsd。append（［measures［‘rrsd’］， bpm， x］）
for x，y，z in rrsd： #Test list entries and select valid measures
if （（x > 1） and （（y > 30） and （y < 130）））：
valid_ma。append（［x， z］）
measures［‘best’］ = min（valid_ma， key = lambda t： t［0］）［1］ #Save the ma_perc for plotting purposes later on （not needed）
detect_peaks（dataset， min（valid_ma， key = lambda t： t［0］）［1］， fs） #Detect peaks with ‘ma_perc’ that goes with lowest rrsd
現在對資料集執行和情節（定時，整個演算法到目前為止，包括裝載在上i7-4790預處理約151ms，單核效能，所以仍然還是相當快的多執行緒將加快這一漲了不少。）：
這已經是好了很多。並沒有丟棄任何正確的R峰，但是仍然存在一些不正確的檢測，還有從0到大約5000的部分，其中幾乎沒有心率訊號。將回到這個嘈雜的部分，並在以後瞭解如何檢測和排除噪聲。
現在，讓從一開始就去掉嘈雜的部分，看看如何檢測和過濾離群值。
查詢錯誤檢測的峰
最後要看的是如何檢測和過濾異常峰位置。一種方法是利用心率逐漸而不是突然變化的事實。例如，bpm不會在單個拍子中從60bpm跳到120bpm，反之亦然。
同樣，這意味著返回到RR間隔。如果在檢測中跳過了一個峰，或者標記了兩個峰而不是一個峰，那麼所得的RR間隔將比平均間隔大得多或小得多。可以設定一個閾值並標記超過該閾值的間隔，類似於檢測峰的方式。對於收集的資料來說，這就足夠了。
然而，存在一種潛在的併發症。如果一次分析很長的訊號，其中心率隨時間變化很大，則可能導致錯誤的過濾。想象一下從60 bpm開始到180bpm結束時心率逐漸增加的訊號。這意味著減少RR間隔的趨勢是穩定的，這表明心率加快，而不是R峰檢測中的錯誤。僅使用基於平均RR間隔的閾值，這將導致訊號的第一部分和最後一部分被過濾。如果在資料中發生這種情況，可以先
降低RR_list的趨勢
。使用scipy。signal，這很容易：
from scipy import signal
RR_list = measures［‘RR_list’］ #First retrieve list from dictionary
RR_list_detrended = signal。detrend（RR_list， type=‘linear’）
但是，如果訊號包含一段大幅度增加的週期，然後心率又出現類似的大幅度下降，則將需要採用其方法。解決方案超出了本教程系列的範圍，但是如果遇到此問題，則可能要使用截止頻率非常低的高通濾波器。另一種方法是將訊號分成較小的部分（以便將上升趨勢和下降趨勢分開），線性下降趨勢並平均測量值。如果較小部分的長度不相等，請確保在平均之前權衡小節。
自然，
不要
使用去趨勢化的RR訊號計算任何量度，而僅將其用於檢測峰標記中的錯誤。
返回異常值過濾。對於閾值，在實踐中，發現在兩端具有250-300ms頻帶的RR差均值的閾值水平很好。使用先前的程式碼並將資料集限制為［5000：15000］（目前要排除嘈雜的開始），應像這樣實現：
RR_list = measures［‘RR_list’］ #Get measures
peaklist = measures［‘peaklist’］
ybeat = measures［‘ybeat’］
upper_threshold = （np。mean（RR_list） + 300） #Set thresholds
lower_threshold = （np。mean（RR_list） - 300）
#detect outliers
cnt = 0
removed_beats = ［］
removed_beats_y = ［］
RR2 = ［］
while cnt < len（RR_list）：
if （RR_list［cnt］ < upper_threshold） and （RR_list［cnt］ > lower_threshold）：
RR2。append（RR_list［cnt］）
cnt += 1
else：
removed_beats。append（peaklist［cnt］）
removed_beats_y。append（ybeat［cnt］）
cnt += 1
measures［‘RR_list_cor’］ = RR2 #Append corrected RR-list to dictionary
plt。subplot（211）
plt。title（‘Marked Uncertain Peaks’）
plt。plot（dataset。hart， color=‘blue’， alpha=0。6， label=‘heart rate signal’）
plt。plot（measures［‘rolmean’］， color=‘green’）
plt。scatter（measures［‘peaklist’］， measures［‘ybeat’］， color=‘green’）
plt。scatter（removed_beats， removed_beats_y， color=‘red’， label=‘Detection uncertain’）
plt。legend（framealpha=0。6， loc=4）
plt。subplot（212）
plt。title（“RR-intervals with thresholds”）
plt。plot（RR_list）
plt。axhline（y=upper_threshold， color=‘red’）
plt。axhline（y=lower_threshold， color=‘red’）
plt。show（）
結果：
似乎遇到了所有的小問題。結果是將所有正確檢測到的R峰標記為綠色的圖。不正確的標記為紅色。生成的列表
測量['RR_list_cor']
的RR列表中沒有屬於錯誤峰的那些。
以後將研究如何標記和排除噪聲段以及其一些最佳化。
彙總
整理所有程式碼，並更新一些函式以接受新變數並插入峰值抑制函式。
