脈搏波信號特征包含人體心血管的生理病理信息, 準確識別其特征點對分析人體心血管的健康狀況有重要意義。脈搏波特征點因人而異, 且易受干擾, 因而獲取高信噪比的完整信號是準確識別特征點的前提。本文基于數學形態學, 設計了一種組合濾波器對信號作預處理, 可有效抑制基線漂移并去除高頻噪聲。利用特征點與小波變換過零點的位置對應關系, 識別預處理后信號的特征點, 并提出用微分法來校正特征點的位置偏移。采用三個參數可調的高斯函數重構出四種典型的脈搏波進行對照實驗, 證明了本文方法識別脈搏波特征點的準確性。
引用本文： 孫薇, 唐寧, 江貴平. 脈搏波信號特征點識別與預處理方法研究. 生物醫學工程學雜志, 2015, 32(1): 197-201. doi: 10.7507/1001-5515.20150036 復制
引言
脈搏是心臟和血管狀態等的重要信息的外在反映, 人體任何一個系統的狀態變化都會影響到脈搏系統。可在體表多處測量脈搏, 如手腕、頸部、足腕、左胸等，根據脈搏壓力的變化，采集到的脈搏信號會呈現出一定的波形[1]，其中手腕上采集的橈動脈信號是最普遍的。正常橈動脈脈搏波如圖 1所示。
圖1
正常橈動脈脈搏波圖形
Figure1.
Normal radial artery pulse diagram
人體脈搏波一般被認為有6個特征點，如圖 1所示，點b是主動脈瓣開放點，主要反映收縮期末血管內的壓力和容積；點c是收縮期最高壓力點，即主波峰值；點d是主動脈擴張降壓點，是左心室射血沖擊主動脈發生彈性振動造成的；點e是左心射血停止點，重搏前波，主要與外周阻力、血管彈性及降支下降速度等變化速度有關；點f是反潮波起點，重搏波波谷；點g是重搏波，可以反映主動脈瓣的功能狀況、血管彈性和血流流動狀態。這些特征點分別反映出心血管的不同狀態[2]。因此脈搏波特征點的準確提取及有效的預處理，對于分析生理病理信息會起到重要的作用。
1 小波變換過零點提取脈搏波特征點
若存在ε＞0，f(x)∈L2(R)在區間[x0-ε, x0]嚴格正(負)，而在區間[x0, x0+ε]嚴格負(正)時，稱函數f(x)在x0點有一個過零點。若小波函數是一個光滑函數θ(t)的二階導數，小波變換的過零點就表示了信號突變點的位置[3]。
設θ(t)是某一起平滑作用的低通平滑函數(濾波器)，且滿足條件：
| $ \int\limits_{-\infty }^{\infty }{\theta \left(t \right)\text{d}t=1\ \ \ \ \underset{\left| t \right|\to \infty }{\mathop{\lim }}\, }\theta \left(t \right)=0 $ |
一般地, 選高斯函數為θ(t)，即
| $ \theta \left(t \right)=\frac{1}{\sqrt{2\pi \sigma }}{{\text{e}}^{1\frac{{{t}^{2}}}{2{{\sigma }^{2}}}}} $ |
θ(t)是二次可導的，其一、二階導數分別為：
| $ {{\Psi }^{\left(1 \right)}}\left(t \right)=\frac{d\theta \left(t \right)}{\text{d}t}\ \ \ {{\Psi }^{\left(2 \right)}}\left(t \right)=\frac{{{d}^{2}}\theta \left(t \right)}{\text{d}t} $ |
則有
| $ \begin{align} & \int\limits_{-\infty }^{\infty }{{{\Psi }^{\left(1 \right)}}\left(t \right)\text{d}t=\int\limits_{-\infty }^{\infty }{\frac{d\theta \left(t \right)}{\text{d}t}=0}} \\ & \int\limits_{-\infty }^{\infty }{{{\Psi }^{\left(2 \right)}}\left(t \right)\text{d}t=\int\limits_{-\infty }^{\infty }{\frac{{{d}^{2}}\theta \left(t \right)}{\text{d}{{t}^{2}}}\text{d}t=0}} \\ \end{align} $ |
所以可以將Ψ(1)(t)和Ψ(2)(t)作為小波函數。
以它們為小波母函數，函數f(t)在尺度為a時間為t的卷積型小波變換可定義為：
| $ \begin{align} & W_{a}^{\left(1 \right)}f\left(t \right)=f\left(t \right)*{{\Psi }^{\left(1 \right)}}a\left(t \right)=\frac{1}{a}\int\limits_{-\infty }^{\infty }{f\left(\tau \right){{\Psi }^{\left(1 \right)}}\left(\frac{t-\tau }{a} \right)}\text{d}\tau \\ & \ \ \ \ \ \ \ \ \ \ \ \ =a\frac{d}{\text{d}t}\left(f*{{\theta }_{a}} \right)\left(t \right)\\ \end{align} $ |
| $ \begin{align} & W_{a}^{\left(2 \right)}f\left(t \right)=f\left(t \right)*{{\Psi }^{\left(2 \right)}}a\left(t \right)=\frac{1}{a}\int\limits_{-\infty }^{\infty }{f\left(\tau \right){{\Psi }^{\left(2 \right)}}\left(\frac{t-\tau }{a} \right)}\text{d}\tau \\ & \ \ \ \ \ \ \ \ \ \ \ ={{a}^{2}}\frac{{{d}^{2}}}{\text{d}{{t}^{2}}}\left(f*{{\theta }_{a}} \right)\left(t \right)\\ \end{align} $ |
由式(6)可看出，當小波函數可看作某一平滑函數的二階導數時，信號小波變換模的過零點對應于信號的突變點[4-5]。本文選擇平滑小波Marr(墨西哥草帽小波)為小波基函數，因為，Marr小波基是高斯函數的二階導數，滿足過零點檢測條件。第二，該小波不是緊支撐的，也不是正交的，但它是對稱的，可用于連續小波變換，得到平滑而連續的小波系數[6]。用Marr小波對脈搏信號進行各尺度分解，并計算第5層小波系數的所有過零點，根據上升支與下降支來區分提取到的點位置對應原信號的波谷和波峰點。并且判斷局部極小值與極大值之間的上升支的過零點位置對應起始點b，找到起始點后，按照順序依次識別c、d、e、f和g點。本文之所以選取第5層小波系數來計算，主要是由于硬件系統采樣頻率為400 Hz，根據Nyquist采樣定理，則信號的頻率范圍為0~200 Hz。那么小波變換對脈搏信號第4層、第5層分解的低頻頻帶分別為0~12.5 Hz和0~6.25 Hz，根據下文可知脈搏波頻譜主要分布在0~10 Hz，本文分別提取4層和5層小波系數作特征點提取，發現由于5層的幅度大一些，且沒有特征點缺失，檢測的過零點位置更準確，因此采用第5層的小波系數的過零點進行特征點的識別。
經過Marr小波分解的脈搏信號，在較小尺度上信號的小波變換幅度很小，很難對特征點進行準確定位，而且，當信號中存在噪聲時，在這些尺度上受噪聲的影響比較嚴重；而在大的尺度上，特征參量的位置會發生比較大的偏移[7-8]。
2 差分法校正小波變換位置偏移
通過時域差分法來校正小波變換后的偏移。為了確保求得的偏移量的準確性，采用起始點和峰值點兩個點進行校正。首先，提取每個周期中起始點b和峰值點c的位置B=(b1, b2, …, bm)，C=(c1, c2, …, cm)，然后與小波提取的對應的b、c點位置相減，得到差值，取差值的平均值作為偏移量。在周期識別過程中，逐個點尋找極大值和極小值點很費時，由于脈搏波周期多為700~l 200 ms，而特征點b與c的時間間隔多為70~120 ms，即整個周期的十分之一左右[9]，因此可將信號一階差分之后選取一個合適的L=λT，其中λ為常數，T為b與c點之間的時間間隔，在適當的L時間間隔內尋找極大值和極小值便可識別出b和c點。同時，識別出的峰值點也有可能是噪聲點，那么通過設定兩點之間時間差的范圍可以去除噪聲點。
3 脈搏波信號的預處理
脈搏波信號主要受基線漂移、工頻干擾以及高頻干擾的影響。針對上述干擾本文采用數學形態學的方法，因為同其他濾波器相比，形態學方法直觀且易于使用，速度快，運算量小，準確率高，不涉及較難理解的數學變換及公式，也不需要定義不同的變換域和有效帶寬問題。下面給出形態學中最基本的膨脹、腐蝕、開運算和閉運算的定義:
設一維離散化脈搏波信號為f(t)(t=0, 1, …, T-1)，結構元素g(m)(t=0, 1, …, M-1)的一維離散函數，且T≥M。f(t)關于g(m)的形態學運算子定義為：
膨脹運算：
腐蝕運算：
開運算:
閉運算:
為了同時去除信號中的“尖峰”和“低谷”，以及避免由互補的形態運算產生相反的偏差，即開運算的收縮性使形態開-閉濾波器輸出偏小，閉運算的擴張性使形態閉-開濾波器輸出偏大[10]，本文采用形態開-閉和形態閉-開平均組合濾波器。
抑制基線漂移的計算公式如下：
| $ {{p}_{1}}\left(t \right)f\left(t \right)-\frac{\left[f\left(t \right)\text{ }\!\!{}^\circ\!\!\text{ }{{g}_{1}}\cdot {{g}_{1}}+f\left(t \right)\cdot {{g}_{1}}\text{ }\!\!{}^\circ\!\!\text{ }{{g}_{1}} \right]}{2} $ |
去除高頻噪聲的計算公式如下：
| $ {{p}_{2}}\left(t \right)=\frac{f\left(t \right)\text{ }\!\!{}^\circ\!\!\text{ }{{g}_{2}}\cdot {{g}_{2}}+f\left(t \right)\cdot {{g}_{2}}\text{ }\!\!{}^\circ\!\!\text{ }{{g}_{2}}}{2} $ |
本文對不同形狀的結構元素作了比較，發現其對信號處理的結果影響不大，實際應用中采用了計算量較小的直線形結構元素。結構元素的寬度由波形寬度和采樣頻率來決定，對于采樣頻率400 Hz的脈搏波信號，其整波寬度在400左右，局部特征最小寬度在30以上，因此，抑制基線漂移的結構元素g1長度要大于脈搏波信號最大特征分量，取450；去除高頻噪聲的結構元素g2長度要小于脈搏波信號中的最大特征分量而大于噪聲的長度，取25。結構元素高度對濾波效果影響不大[11]，取為0。如圖 2、3所示，數學形態學濾波可以有效地抑制基線漂移和去除高頻噪聲，同時保留信號的有用信息。
圖2
抑制基線漂移
Figure2.
Inhibition of baseline drift
圖3
去除高頻噪聲
Figure3.
Removal of high-frequency noise
4 仿真實驗
由于個體差異，脈搏波的特征點也有很大的不同，因此本研究通過仿真幾種典型的脈搏波形來驗證本文的方法。脈搏波在形式上可以看做主波、重搏波和重搏前波等三者疊加而成。錢偉立等[12]提出將脈搏波的每個特征波形都用一個高斯函數來近似，脈搏波的特征變化由三個高斯函數的參數變化構成，仿真出不同脈搏波。本文選取幾種在臨床上既典型又有重要意義的仿真脈搏信號來驗證方法，典型仿真脈搏波形可參見文獻[13-14]。
數學解析模型表示為：
| $ f\left(t \right)=\sum\limits_{i=1}^{3}{{{h}_{i}}\cdot {{\text{e}}^{-\frac{{{\left(t-{{\mu }_{i}} \right)}^{2}}}{\delta _{i}^{2}}}}} $ |
其中h1、h2、h3分別表示主峰波、重搏前波和重搏波的高度，μ1、μ2、μ3分別表示主峰波、重搏前波和重搏波的波形偏移量，δ1、δ2、δ3分別表示主峰波、重搏前波和重搏波的波形寬度。本文仿真的幾種典型脈搏波的波形參數見表 1。
表1
幾種典型脈搏波的仿真波形參數
Table1.
Simulation parameters of several typical pulse waves
為了驗證仿真脈搏波信號的有效性，本文將仿真的平脈(健康人脈搏波)與實際采集的脈搏波分別作傅里葉變換，得到對應的頻譜，如圖 4所示。脈搏波的頻率集中在0.7~10 Hz，其中5~10 Hz是其細節特征，從圖中可以看出仿真得到的信號符合實際脈搏波的頻譜分布，并且細節特征更明顯。
圖4
脈搏波頻譜對比圖
Figure4.
Comparison chart of pulse wave spectrum
5 結果與討論
實驗數據來源于仿真和自行研發的脈搏信號采集儀。采集儀的電路系統中選擇截止頻率為20 Hz的四階巴特沃斯低通濾波器，主要是由于在15 Hz以下的脈搏波能量衰減達到-40 dB，在8 Hz以下時能量衰減可以達到-35 dB，并且其95%的能量都集中在0~6 Hz，因此既保證了脈搏波信號的重要生理信息，又濾掉了高頻信號和50 Hz的工頻干擾。對采集儀采集到的脈搏信號做特征點識別之前，采用上文提到的基于數學形態學的方法，設計組合濾波器對信號做預處理。
提取預處理后的實際脈搏波和仿真的幾種典型脈搏波的特征點，這里以典型平脈為例。按照本文提出的方法對其作小波變換5層分解，提取第5層小波系數，如圖 5上圖所示；查找并提取所有過零點的位置，并對應到原信號上，如圖 5中圖所示；過零點位置與對應的脈搏波的特征點位置存在偏移，采用本文提出的時域差分法求出偏移量，校正偏移，可以有效地識別b點和c點的位置，如圖 5下圖所示，經過校正后，可以有效地識別脈搏波的特征點。
圖5
小波變換后的偏移校正
Figure5.
Offset correction after wavelet transform
如圖 6所示，采用本文的方法分別對實際采集到的脈搏波信號、仿真的數脈信號以及滑脈信號做特征點識別。其中計算的偏移量是每5個脈搏周期的平均值，由于硬件系統的顯示模塊，每次顯示信號的長度大約為5個脈搏周期，因此為了簡化算法并減小計算量，本文設定每5個脈搏周期進行一次偏移校正。從圖中可以看出對特征點不明顯的數脈可以有效地識別6個特征點，對于特征點被淹沒的滑脈等可以有效識別信號顯現出的特征點。由于硬件系統設計還有缺陷，實際采集到的脈搏信號f點和g點被噪聲所覆蓋，因此沒有識別出這兩個特征點，這也是本文方法在今后需要改進的地方。從醫學上來說，識別起始點b、峰值點c、降中峽d和重搏波e的位置可以得到人體心血管等系統的很多生理信息，為病例分析提供了很好的依據。
圖6
偏移校正后的結果圖
Figure6.
Results of offset correction
引言
脈搏是心臟和血管狀態等的重要信息的外在反映, 人體任何一個系統的狀態變化都會影響到脈搏系統。可在體表多處測量脈搏, 如手腕、頸部、足腕、左胸等，根據脈搏壓力的變化，采集到的脈搏信號會呈現出一定的波形[1]，其中手腕上采集的橈動脈信號是最普遍的。正常橈動脈脈搏波如圖 1所示。
圖1
正常橈動脈脈搏波圖形
Figure1.
Normal radial artery pulse diagram
人體脈搏波一般被認為有6個特征點，如圖 1所示，點b是主動脈瓣開放點，主要反映收縮期末血管內的壓力和容積；點c是收縮期最高壓力點，即主波峰值；點d是主動脈擴張降壓點，是左心室射血沖擊主動脈發生彈性振動造成的；點e是左心射血停止點，重搏前波，主要與外周阻力、血管彈性及降支下降速度等變化速度有關；點f是反潮波起點，重搏波波谷；點g是重搏波，可以反映主動脈瓣的功能狀況、血管彈性和血流流動狀態。這些特征點分別反映出心血管的不同狀態[2]。因此脈搏波特征點的準確提取及有效的預處理，對于分析生理病理信息會起到重要的作用。
1 小波變換過零點提取脈搏波特征點
若存在ε＞0，f(x)∈L2(R)在區間[x0-ε, x0]嚴格正(負)，而在區間[x0, x0+ε]嚴格負(正)時，稱函數f(x)在x0點有一個過零點。若小波函數是一個光滑函數θ(t)的二階導數，小波變換的過零點就表示了信號突變點的位置[3]。
設θ(t)是某一起平滑作用的低通平滑函數(濾波器)，且滿足條件：
| $ \int\limits_{-\infty }^{\infty }{\theta \left(t \right)\text{d}t=1\ \ \ \ \underset{\left| t \right|\to \infty }{\mathop{\lim }}\, }\theta \left(t \right)=0 $ |
一般地, 選高斯函數為θ(t)，即
| $ \theta \left(t \right)=\frac{1}{\sqrt{2\pi \sigma }}{{\text{e}}^{1\frac{{{t}^{2}}}{2{{\sigma }^{2}}}}} $ |
θ(t)是二次可導的，其一、二階導數分別為：
| $ {{\Psi }^{\left(1 \right)}}\left(t \right)=\frac{d\theta \left(t \right)}{\text{d}t}\ \ \ {{\Psi }^{\left(2 \right)}}\left(t \right)=\frac{{{d}^{2}}\theta \left(t \right)}{\text{d}t} $ |
則有
| $ \begin{align} & \int\limits_{-\infty }^{\infty }{{{\Psi }^{\left(1 \right)}}\left(t \right)\text{d}t=\int\limits_{-\infty }^{\infty }{\frac{d\theta \left(t \right)}{\text{d}t}=0}} \\ & \int\limits_{-\infty }^{\infty }{{{\Psi }^{\left(2 \right)}}\left(t \right)\text{d}t=\int\limits_{-\infty }^{\infty }{\frac{{{d}^{2}}\theta \left(t \right)}{\text{d}{{t}^{2}}}\text{d}t=0}} \\ \end{align} $ |
所以可以將Ψ(1)(t)和Ψ(2)(t)作為小波函數。
以它們為小波母函數，函數f(t)在尺度為a時間為t的卷積型小波變換可定義為：
| $ \begin{align} & W_{a}^{\left(1 \right)}f\left(t \right)=f\left(t \right)*{{\Psi }^{\left(1 \right)}}a\left(t \right)=\frac{1}{a}\int\limits_{-\infty }^{\infty }{f\left(\tau \right){{\Psi }^{\left(1 \right)}}\left(\frac{t-\tau }{a} \right)}\text{d}\tau \\ & \ \ \ \ \ \ \ \ \ \ \ \ =a\frac{d}{\text{d}t}\left(f*{{\theta }_{a}} \right)\left(t \right)\\ \end{align} $ |
| $ \begin{align} & W_{a}^{\left(2 \right)}f\left(t \right)=f\left(t \right)*{{\Psi }^{\left(2 \right)}}a\left(t \right)=\frac{1}{a}\int\limits_{-\infty }^{\infty }{f\left(\tau \right){{\Psi }^{\left(2 \right)}}\left(\frac{t-\tau }{a} \right)}\text{d}\tau \\ & \ \ \ \ \ \ \ \ \ \ \ ={{a}^{2}}\frac{{{d}^{2}}}{\text{d}{{t}^{2}}}\left(f*{{\theta }_{a}} \right)\left(t \right)\\ \end{align} $ |
由式(6)可看出，當小波函數可看作某一平滑函數的二階導數時，信號小波變換模的過零點對應于信號的突變點[4-5]。本文選擇平滑小波Marr(墨西哥草帽小波)為小波基函數，因為，Marr小波基是高斯函數的二階導數，滿足過零點檢測條件。第二，該小波不是緊支撐的，也不是正交的，但它是對稱的，可用于連續小波變換，得到平滑而連續的小波系數[6]。用Marr小波對脈搏信號進行各尺度分解，并計算第5層小波系數的所有過零點，根據上升支與下降支來區分提取到的點位置對應原信號的波谷和波峰點。并且判斷局部極小值與極大值之間的上升支的過零點位置對應起始點b，找到起始點后，按照順序依次識別c、d、e、f和g點。本文之所以選取第5層小波系數來計算，主要是由于硬件系統采樣頻率為400 Hz，根據Nyquist采樣定理，則信號的頻率范圍為0~200 Hz。那么小波變換對脈搏信號第4層、第5層分解的低頻頻帶分別為0~12.5 Hz和0~6.25 Hz，根據下文可知脈搏波頻譜主要分布在0~10 Hz，本文分別提取4層和5層小波系數作特征點提取，發現由于5層的幅度大一些，且沒有特征點缺失，檢測的過零點位置更準確，因此采用第5層的小波系數的過零點進行特征點的識別。
經過Marr小波分解的脈搏信號，在較小尺度上信號的小波變換幅度很小，很難對特征點進行準確定位，而且，當信號中存在噪聲時，在這些尺度上受噪聲的影響比較嚴重；而在大的尺度上，特征參量的位置會發生比較大的偏移[7-8]。
2 差分法校正小波變換位置偏移
通過時域差分法來校正小波變換后的偏移。為了確保求得的偏移量的準確性，采用起始點和峰值點兩個點進行校正。首先，提取每個周期中起始點b和峰值點c的位置B=(b1, b2, …, bm)，C=(c1, c2, …, cm)，然后與小波提取的對應的b、c點位置相減，得到差值，取差值的平均值作為偏移量。在周期識別過程中，逐個點尋找極大值和極小值點很費時，由于脈搏波周期多為700~l 200 ms，而特征點b與c的時間間隔多為70~120 ms，即整個周期的十分之一左右[9]，因此可將信號一階差分之后選取一個合適的L=λT，其中λ為常數，T為b與c點之間的時間間隔，在適當的L時間間隔內尋找極大值和極小值便可識別出b和c點。同時，識別出的峰值點也有可能是噪聲點，那么通過設定兩點之間時間差的范圍可以去除噪聲點。
3 脈搏波信號的預處理
脈搏波信號主要受基線漂移、工頻干擾以及高頻干擾的影響。針對上述干擾本文采用數學形態學的方法，因為同其他濾波器相比，形態學方法直觀且易于使用，速度快，運算量小，準確率高，不涉及較難理解的數學變換及公式，也不需要定義不同的變換域和有效帶寬問題。下面給出形態學中最基本的膨脹、腐蝕、開運算和閉運算的定義:
設一維離散化脈搏波信號為f(t)(t=0, 1, …, T-1)，結構元素g(m)(t=0, 1, …, M-1)的一維離散函數，且T≥M。f(t)關于g(m)的形態學運算子定義為：
膨脹運算：
腐蝕運算：
開運算:
閉運算:
為了同時去除信號中的“尖峰”和“低谷”，以及避免由互補的形態運算產生相反的偏差，即開運算的收縮性使形態開-閉濾波器輸出偏小，閉運算的擴張性使形態閉-開濾波器輸出偏大[10]，本文采用形態開-閉和形態閉-開平均組合濾波器。
抑制基線漂移的計算公式如下：
| $ {{p}_{1}}\left(t \right)f\left(t \right)-\frac{\left[f\left(t \right)\text{ }\!\!{}^\circ\!\!\text{ }{{g}_{1}}\cdot {{g}_{1}}+f\left(t \right)\cdot {{g}_{1}}\text{ }\!\!{}^\circ\!\!\text{ }{{g}_{1}} \right]}{2} $ |
去除高頻噪聲的計算公式如下：
| $ {{p}_{2}}\left(t \right)=\frac{f\left(t \right)\text{ }\!\!{}^\circ\!\!\text{ }{{g}_{2}}\cdot {{g}_{2}}+f\left(t \right)\cdot {{g}_{2}}\text{ }\!\!{}^\circ\!\!\text{ }{{g}_{2}}}{2} $ |
本文對不同形狀的結構元素作了比較，發現其對信號處理的結果影響不大，實際應用中采用了計算量較小的直線形結構元素。結構元素的寬度由波形寬度和采樣頻率來決定，對于采樣頻率400 Hz的脈搏波信號，其整波寬度在400左右，局部特征最小寬度在30以上，因此，抑制基線漂移的結構元素g1長度要大于脈搏波信號最大特征分量，取450；去除高頻噪聲的結構元素g2長度要小于脈搏波信號中的最大特征分量而大于噪聲的長度，取25。結構元素高度對濾波效果影響不大[11]，取為0。如圖 2、3所示，數學形態學濾波可以有效地抑制基線漂移和去除高頻噪聲，同時保留信號的有用信息。
圖2
抑制基線漂移
Figure2.
Inhibition of baseline drift
圖3
去除高頻噪聲
Figure3.
Removal of high-frequency noise
4 仿真實驗
由于個體差異，脈搏波的特征點也有很大的不同，因此本研究通過仿真幾種典型的脈搏波形來驗證本文的方法。脈搏波在形式上可以看做主波、重搏波和重搏前波等三者疊加而成。錢偉立等[12]提出將脈搏波的每個特征波形都用一個高斯函數來近似，脈搏波的特征變化由三個高斯函數的參數變化構成，仿真出不同脈搏波。本文選取幾種在臨床上既典型又有重要意義的仿真脈搏信號來驗證方法，典型仿真脈搏波形可參見文獻[13-14]。
數學解析模型表示為：
| $ f\left(t \right)=\sum\limits_{i=1}^{3}{{{h}_{i}}\cdot {{\text{e}}^{-\frac{{{\left(t-{{\mu }_{i}} \right)}^{2}}}{\delta _{i}^{2}}}}} $ |
其中h1、h2、h3分別表示主峰波、重搏前波和重搏波的高度，μ1、μ2、μ3分別表示主峰波、重搏前波和重搏波的波形偏移量，δ1、δ2、δ3分別表示主峰波、重搏前波和重搏波的波形寬度。本文仿真的幾種典型脈搏波的波形參數見表 1。
表1
幾種典型脈搏波的仿真波形參數
Table1.
Simulation parameters of several typical pulse waves
為了驗證仿真脈搏波信號的有效性，本文將仿真的平脈(健康人脈搏波)與實際采集的脈搏波分別作傅里葉變換，得到對應的頻譜，如圖 4所示。脈搏波的頻率集中在0.7~10 Hz，其中5~10 Hz是其細節特征，從圖中可以看出仿真得到的信號符合實際脈搏波的頻譜分布，并且細節特征更明顯。
圖4
脈搏波頻譜對比圖
Figure4.
Comparison chart of pulse wave spectrum
5 結果與討論
實驗數據來源于仿真和自行研發的脈搏信號采集儀。采集儀的電路系統中選擇截止頻率為20 Hz的四階巴特沃斯低通濾波器，主要是由于在15 Hz以下的脈搏波能量衰減達到-40 dB，在8 Hz以下時能量衰減可以達到-35 dB，并且其95%的能量都集中在0~6 Hz，因此既保證了脈搏波信號的重要生理信息，又濾掉了高頻信號和50 Hz的工頻干擾。對采集儀采集到的脈搏信號做特征點識別之前，采用上文提到的基于數學形態學的方法，設計組合濾波器對信號做預處理。
提取預處理后的實際脈搏波和仿真的幾種典型脈搏波的特征點，這里以典型平脈為例。按照本文提出的方法對其作小波變換5層分解，提取第5層小波系數，如圖 5上圖所示；查找并提取所有過零點的位置，并對應到原信號上，如圖 5中圖所示；過零點位置與對應的脈搏波的特征點位置存在偏移，采用本文提出的時域差分法求出偏移量，校正偏移，可以有效地識別b點和c點的位置，如圖 5下圖所示，經過校正后，可以有效地識別脈搏波的特征點。
圖5
小波變換后的偏移校正
Figure5.
Offset correction after wavelet transform
如圖 6所示，采用本文的方法分別對實際采集到的脈搏波信號、仿真的數脈信號以及滑脈信號做特征點識別。其中計算的偏移量是每5個脈搏周期的平均值，由于硬件系統的顯示模塊，每次顯示信號的長度大約為5個脈搏周期，因此為了簡化算法并減小計算量，本文設定每5個脈搏周期進行一次偏移校正。從圖中可以看出對特征點不明顯的數脈可以有效地識別6個特征點，對于特征點被淹沒的滑脈等可以有效識別信號顯現出的特征點。由于硬件系統設計還有缺陷，實際采集到的脈搏信號f點和g點被噪聲所覆蓋，因此沒有識別出這兩個特征點，這也是本文方法在今后需要改進的地方。從醫學上來說，識別起始點b、峰值點c、降中峽d和重搏波e的位置可以得到人體心血管等系統的很多生理信息，為病例分析提供了很好的依據。
圖6
偏移校正后的結果圖
Figure6.
Results of offset correction
