基于心音香農能量包絡的實時心率檢測的研究和實現*

范能勝， 張懷岺

(廣東醫科大學生物醫學工程學院，東莞 523808)

1 引 言

心臟是人體的重要器官，其節律性跳動是由竇房結產生的電信號引起的，稱之為心電信號。心臟的一些疾病會引起心電波形的改變，但也有些疾病只會引起心音的改變，如瓣膜類的疾病。心臟在工作時，由于瓣膜的開啟和關閉，血液對心臟內壁、瓣膜和血管的沖擊會發出節律性的聲音，稱之為“心音”。一般認為有四種心音[1]，分別是第一心音(S1)、第二心音(S2)、第三心音(S3)和第四心音(S4)。S1代表心室收縮期的開始，其振幅大，持續時間較長；S2代表心室舒張期的開始，時間較短，振幅低；S3和S4不易聽到。在分析心音信號時，S1和S2是最重要的信號來源。

聽診器是常用的心音探測設備，其結構簡單，攜帶方便，是醫生的首選。但通過聽診診斷心臟疾病主觀性較強，醫生需要經過長期的實踐，經驗越豐富的醫生診斷陽性率越高；隨著科技的發展，用特定的傳感器將心音拾取出來，通過放大和模數轉換后變成數字信號，利用分析算法從信號中提取有用信息，可以避免醫生主觀因素的影響。一些典型病例的心音信號在數字化后，有利于保存和傳播，可以用來培訓醫生；此外，在日常的診療活動中，有時需要知道患者的心率變化情況，雖然通過心電信號計算心率最為準確，但是需要連接心電圖機，操作不便，醫生一般通過計數患者頸部或腕部動脈在單位時間內的搏動次數來計算心率，方法簡單，但耗時長。

本研究介紹了一種利用心音信號來計算心率的方法，通過心音傳感器拾取心音信號并將其數字化，按照一定算法提取心音包絡信號，根據包絡信號計算心率，理論上只要測量兩個以上的心動周期就可以將心率計算出來，最終在聽診的同時完成心率的測量。

2 材料和設備

2.1 心音傳感器

正常人的心音信號頻率一般在5～600 Hz[2]，測試系統選用了安徽合肥華科電子技術研究所生產的HKY-06B+型心音傳感器作為心音拾取設備，內含微音傳感元件，可以采集心臟和其它體表動脈搏動信號，經過內部集成化信號處理電路處理，輸出模擬電信號，幅值為0.5～1.5 V，頻率響應為10～600 Hz，使用單電源供電，電壓范圍3～5 V，傳感器實物見圖1。

2.2 實驗硬件平臺

由于在實現算法時需要做浮點數運算，同時還需要顯示波形和其它信息，故本研究選用型號STM32F103VET6的單片機作為核心處理器，該單片機具有32位數據處理寬度，內部資源豐富，處理速度快，并自帶靜態存儲控制器(flexible static memory controller,FSMC)接口功能，方便連接LCD液晶顯示；采用2.8英寸并帶觸摸控制功能的液晶顯示器進行操作控制和顯示波形，核心實驗電路見圖2。

圖1 心音傳感器Fig.1 Heart sound sensor

由于單片機自帶的12位A/D轉換器，可以滿足實驗要求，不必再外擴A/D轉換電路，心音傳感器輸出的模擬信號經單片機的PC1引腳輸入，前輟為LCD的網絡標號連接液晶顯示器，前輟為TP和SPI的網絡標號連接觸摸屏。

2.3 實驗數據來源

根據心音傳感器頻率響應和奈圭斯特(Nyquist)準則，采樣頻率必須是信號最高頻率的兩倍以上才可以還原出原始信號，因此設置采樣頻率為1 250 Hz。將心音傳感器固定在體表合適位置，采樣8 s，獲得10 000個心音采樣數據，每個采樣數據占用兩個字節，將采樣數據保存在一個數組中，并在采樣結束后導入matlab軟件，作為算法分析驗證的實驗數據。

3 算法

在對心音信號進行采集時，還可能引入其它的噪聲信號，如呼吸音信號、傳感器和胸壁摩擦產生的噪聲信號等。處理算法的目的是衰減各種噪聲信號的同時，將心音信號的包絡提取出來作為分析處理使用。由于本研究的目的是利用心音信號來計算心率，所以需要從心音信號包絡中準確識別心音峰值。根據相鄰峰值之間的距離以及采樣參數，計算瞬時心率值。有關心音信號處理的流程和算法，前人做了許多研究，包括心音信號的預處理、求取香農能量包絡、尋找心音峰值等。結合這些研究成果和實際要求，我們采用如下方法對原始心音信號進行處理。

3.1 信號歸一化預處理

為便于后續處理數據，需對原始輸入信號進行歸一化處理[4]，按式(1)進行:

(1)

圖2 實驗電路Fig.2 Main circuit

式中，xk為原始采樣數據，max (|xi|)為輸入序列數據絕對值的最大值，xnorm為歸一化計算后的輸出信號。

由于傳感器的輸出電壓在0.5～1.5 V，所以輸入信號振幅的峰-峰值為1.0 V，為了實時處理的需要，將式(1)中的max(|xi|)固定設為峰-峰值的一半，即0.5 V。在實時處理時，先對數字化的信號減去一個固定數值(基線下移)，處理過程見圖3，然后再按照式(1)進行歸一化運算。

3.2 計算香農能量包絡

使用式(1)對原始信號進行歸一化處理后，可以采用不同的方法來計算心音信號的包絡，這些方法[4-8]包括：

絕對值：E=|x|

平方能量：E=x2

香農熵：E=-|x|log |x|

香農能量：E=-x2log (x2)

圖3 原始信號處理過程(a).原始信號;(b).基線下移Fig.3 Process of original signal (a).original signal;(b).baseline down

圖4顯示了不同計算方法的情況下歸一化原始信號幅值和包絡輸出值之間的關系曲線，由于對稱性，負值部分被忽略了。

圖4 不同包絡偵測方法的比較Fig.4 The comparision of different envelope detection methods

由圖4可知，絕對值的方法對輸入信號的所有組成部分采用了相同的權重，難以將需要的信號凸顯出來；平方能量對高強度信號部分給予了較大的權重，不利于低強度信號的識別；香農熵方法對低強度信號給予了較大的權重，衰減了高強度信號；香農能量的方法對中等強度信號給予了更多權重，對低強度信號比高強度信號衰減得更多，這樣增強了中等偏高強度的信號。

對于心率測量只需關注第一心音S1即可，其它信號可以盡可能的衰減，根據文獻[9]的研究,采用式(2)求取香農能量的平均值Eavg：

(2)

式中，N為平均濾波的數據點數，x(i)為輸入信號，n為冪次。冪次n越高，計算所得的包絡波形越平滑。在本研究中，經過驗證，當n=4時可以取得較好的抑噪效果。

3.3 香農能量平均濾波的實現

用式(2)進行平均濾波時，N取值越大，得到的波形越平滑[10]，圖5是N=201時的平均濾波效果，由圖中可以看出，大部分的小幅波動被有效抑制，突出了第一心音信號。

圖5 信號處理過程(a).歸一化處理后的信號;(b).E=-x4log (x4);(c).Fig.5 Signal processing(a).normalization;(b).E=-x4log (x4);(c).

在實時處理時，可以設置一個環型數組用來保存圖4(b)中計算出來的香農能量數據，數組的長度為平均濾波的點數，每次有新的數據進入時，讓數組下標加一再保存。在此以7點平均濾波為例說明移動平均濾波的實現過程。

定義一個長度為7個元素的浮點數組Array，輸入的數據按先后順序依次保存在數組Array[0]～Array[6]，當第8個數據到來時覆蓋最早保存的數據，將其存在Array[0]單元，以后的數據保存方法以此類推，環型數組見圖6。

采用迭代遞歸的方式進行運算，環型數組每次存入的數值由式(3)給出：

(3)

此處的out(i)是本次運算的輸出結果，out(i-1)為前一次運算的輸出結果，m為平均濾波長度，in(i+p)表示在(i+p)刻的輸入值，可理解為最新的輸入值，in(i-q)表示在(i-q)時刻的輸入值，in(i+p)-in(i-q)在環型數組中表示(i+p)刻的輸入值覆蓋掉(i-q)刻的輸入值，當數組長度m=7時，由式(3)計算得到：

圖6 環型數組Fig.6 Circular array

(4)

式(5)-式(8)按順序演示了環型數組保存時的迭代過程，如果括弧內的值為負，表示該項的值為0。

out(0)=out(-1)+in(3)-in(-4)=in(3)

(5)

out(1)=out(0)+in(4)-in(-3)=in(3)+in(4)

(6)

out(2)=out(1)+in(5)-in(-2)=in(3)+in(4)+in(5)

(7)

…

out(7)=out(6)+in(10)-in(0)

(8)

in(10)表示新輸入的第8個數據，in(10)-in(0)表示第8個輸入的數據覆蓋掉下標為0的數組元素，平均濾波的結果為out(i)/m。由于整個運算過程都是簡單的四則運算，所以運算速度很快，完全滿足信號實時處理的要求。

3.4 尋找峰值

要計算實時心率，必須尋找心音包絡信號的峰值并確定它在數列中的坐標[11-12]，由于采樣率已知，則只需知道相鄰兩個峰之間距離(數據點數)就可計算出瞬時心率，計算公式如下：

(9)

式中，HR表示心率，N為相鄰兩個心音峰值之間的數據點數，等于兩個峰值的坐標之差，SR為采樣率。

峰值指單位時間內的最大值，此“單位時間”稱之為時間窗口。窗口的寬度需要根據采樣率和心率來確定，本研究中將窗口寬度設置為0.32 S。同時考慮噪聲的影響，還需設置一個閾值，在窗口范圍內，只有大于這個閾值的最大值才能認為是心音峰值，見圖7。

圖7 尋找峰值(方框為時間窗口)Fig.7 Finding peak(red rectangle present time window)

在程序實現時，將新的輸入值與前一次的值進行比較，前者大于后者，則進行最大值替換，保存數據坐標，并清零窗口計數器，反之，則不替換，并啟動窗口計數器；窗口計數器達到窗口計數寬度，則認為當前的最大值為心音峰值，尋找峰值的算法如下：

//WinCnt: 窗口計數器；WINWIDTH: 窗口寬度；cnt：輸入值計數器(數值坐標)

//CurrVal: 當前值; MaxVal: 最大值；MaxPos: 最大值坐標；ThresholdVal: 閾值If(WinCnt

{

if(CurrVal >MaxVal AND CurrVal>ThresholdVal)

{

MaxVal = CurrVal;

MaxPos = cnt;

WinCnt = 0;

}

else

{

WinCnt++;

}

}

else

{

return MaxPos;//返回峰值的坐標

}

根據以上算法，每找到一個心音峰值，返回其坐標，同時復位相關變量，然后根據式(9)即可計算出心率。實際應用時，可以每三次瞬時心率取平均值作為心率，既可滿足實時顯示的要求，又可以平滑計算引入的誤差。

4 結果和討論

心音信號因受到測量位置、呼吸音、環境等因素的影響，信號噪聲較大，容易受干擾，理論上通過心音信號計算出來的心率比心電信號計算出來的心率誤差大。為了驗證算法的準確性，使用一臺心電圖機(型號ECG-101，深圳邦健公司)作為對照設備，實驗設備和對照設備見圖8。

圖8 對照設備和實驗設備(a).對照設備;(b)實驗設備;(c).結果對照Fig.8 Reference device and experimental device(a).reference device;(b)experimental device;(c).results contrast

選擇了六位測量對象，使用兩種設備同時對每個人的心率進行測量，實驗時測量的結果見表1。

由表1可以看出，與心電圖機測量的心率對照，實驗設備測量值的準確率達到97%以上。由于條件的限制，本次測試的對象均為健康成年人，未對有心臟疾病的患者進行測試，假如對這類患者進行測試，需要進一步完善算法，以提高測量的準確性。