心沖擊(BCG)信號是反映心臟機械運動的生理信號,測量中無需在受試者身體表面貼附電極,能實現無感覺生理監護。但BCG信號微弱,易受到干擾,測量時經常被淹沒在噪聲中。為了有效識別BCG信號,提出一種基于時頻聯合分布和經驗模態分解(EMD)的BCG信號降噪方法。該方法先建立BCG信號的自適應最優核,然后在時頻平面內提取BCG信號分量,最后根據EMD原理對BCG信號分量進行濾波,從而實現BCG信號降噪。仿真研究表明,該方法克服了EMD處理在不同時間含有相同或相似頻率成分信號時的不足,所提出方法實現了BCG信號降噪,可以有效還原BCG信號特征。
引用本文： 楊丹, 徐彬, 葉琳琳, 金晶晶. 心臟心沖擊信號降噪方法研究. 生物醫學工程學雜志, 2014, 31(6): 1368-1372. doi: 10.7507/1001-5515.20140259 復制
引言
心臟泵血活動會使人體產生微小震動，將這種微小震動記錄下來即稱為心沖擊(ballistocardiogram,BCG)信號。正常的BCG與心跳同步，具有重復性，主要包含H、I、J、K、M和N波。J波是BCG周期中幅值最大點，估計相鄰周期的J波間隔，即為心動周期；確定同周期IJ幅值，即為心臟泵血的收縮力大小，如圖 1所示。
圖1
一段正常的BCG信號及其特征
Figure1.
A segment of normal BCG signal and its features
相對于心電圖(electrocardiogram,ECG)信號，BCG信號可在無約束情況下從人體表面測得，其測量系統也可以設計在人們的日常生活家具中如床[1]、椅子[2]、人體秤[3]，因此，BCG信號適合日常生活的無感覺人體健康監測。
基于BCG信號的人體健康監測方法分兩個步驟實現。一是檢測BCG信號，二是提取BCG信號中所包含的生理參數或信息，進行分析評估。然而BCG信號自身十分微弱，且易受到呼吸、身體運動、工頻噪聲的干擾[4-5]，導致直接測量得到的BCG信號常常淹沒在噪聲中，無法獲得其包含的心率、呼吸等生理特征信息。因此，BCG信號降噪處理對于識別其特征參數是非常必要的。目前BCG信號降噪主要采用小波變換方法，其效果取決于小波閾值參數選擇[6]；而從信號分解角度，小波變換中基函數是固定的，不能匹配信號的全部。近年來，出現了一種數據驅動的非線性、非平穩信號處理方法--經驗模態分解(empirical mode decomposition,EMD)。它的基函數是一系列幅度和頻率均可變的正弦余弦函數，可通過信號自適應獲得，其特性依賴于信號本身，能取得很好的分解效果，被廣泛應用到語音[7]、圖像[8]以及ECG信號處理中[9]。
BCG信號是多周期非平穩隨機信號，隨心臟活動出現周期性變化，在不同的時間段有相似的頻率成分出現，直接利用EMD進行信號降噪，不能實現信號成分中預期的本征模態函數(intrinsic mode function,IMF)。為了有效識別用于人體健康監測的多周期BCG信號，本文提出一種基于時頻聯合分布和EMD的BCG信號降噪方法。先對BCG信號進行自適應最優核的時頻分布計算，然后在時頻平面內自適應提取以心動周期為基準的BCG信號分量，再對其進行EMD，最后按照時間順序連接重構低頻成分，實現多周期BCG信號降噪。
1 心臟BCG信號的采集
BCG信號采集裝置由傳感器電路和信號調理電路組成。傳感器電路嵌入改裝的家庭人體秤中，受試者靜止地站于人體秤上，心臟泵血時引起的身體微小震動即為BCG信號。BCG信號通過傳感器電路耦合后，再經差分前置放大電路、多級放大電路、隔直電路、陷波電路、帶通濾波電路，最終得到測量的BCG信號。采集電路的總體構圖如圖 2所示。其中差分放大電路抑制共模信號，隔直電路抑制直流噪聲，陷波電路抑制工頻干擾，帶通電路濾除測量中引發的干擾信號，多級放大實現BCG信號幅值放大。
圖2
BCG信號的采集電路
Figure2.
BCG signals recording circuit
2 提出的BCG信號降噪方法
2.1 基于EMD的降噪
對于給定信號x(t)，經EMD分解可表示為
| $x\left( t \right)=\sum\limits_{i=1}^{N}{IM{{F}_{i}}\left( t \right)+{{r}_{N}}\left( t \right)},$ |
其中rN(t)為趨勢項，IMF1(t)，…，IMFN(t)為信號的IMF。它們分別包含了信號從高到低不同頻率段的成分。基于EMD信號降噪的原理是對信號EMD分解后的多階IMF函數進行篩選，消除含有噪聲成分的IMF，保留含有信號成分的IMF并重構，完成降噪處理，其原理如圖 3所示。EMD降噪的關鍵IMF階數D的選取，通常從第D個IMF到最后一個IMF再加上趨勢項rN(t)即為所估計的低頻信號(t), 而
| $x\left( t \right)=\sum\limits_{i=1}^{N}{IM{{F}_{k}}\left( t \right)+{{r}_{N}}\left( t \right)}$ |
圖3
信號的EMD降噪原理
Figure3.
De-noising theory using EMD
對于非平穩信號更加完備的EMD分解可表達為
| $x\left( t \right)=\sum\limits_{k=1}^{M}{{{a}_{k}}\left( t \right)rect}\left( \frac{t-{{\tau }_{k}}}{{{b}_{k}}} \right)cos({{\theta }_{k}}\left( t \right))+r\left( t \right),$ |
其中ak(t)是IMFk(t)的包絡線，θk(t)是IMFk(t)的相位，rect((t-τk)/bk)是寬度為bk、中心為τk的矩形函數。信號x(t)可通過cos(θk(t))和rect((t-τk)/bk)在時頻平面上分別對時間和頻率進行定位。如果信號由一系列分量(分量定義的是信號的時頻平面上能量集中的區域)疊加而成，當每個分量相對于時間軸連續且頻率軸可分，EMD才能將其分解為多個IMF。BCG信號具備周期性，且在不同時間段具有相同成分，因此，要利用EMD分解多周期BCG信號，需提取BCG信號分量，再對其進行EMD分解，從而實現降噪處理。
2.2 基于時頻分布的BCG信號分量提取
2.2.1 核函數的定義及最優化條件
自適應最優核函數g(θ,τ)，其極坐標形式為
| $g\left( r,\Phi \right)=exp\left( -\frac{{{r}^{2}}}{2{{\sigma }^{2}}\left( \Phi \right)} \right),$ |
其中0≤g(r,Φ)≤1，r=為極徑。
對于給定信號及其模糊函數，最優核函數g(r,Φ)應滿足優化問題：
| $maxg\int_{0}^{2\pi }{\int_{0}^{\infty }{{{\left| AF\left( r,\Phi \right)g\left( r,\Phi \right) \right|}^{2}}rdrd\Phi }}\text{ }$ |
約束條件為
| $\frac{1}{4}{{\pi }^{2}}\int_{0}^{2\pi }{\int_{0}^{\infty }{{{\left| g\left( r,\Phi \right) \right|}^{2}}rdrd\Phi \le \alpha ~,}}$ |
其中AF(r,Φ)為待分析信號模糊函數的極坐標表示，α≥0為核函數的能量體積。
2.2.2 BCG信號最優核函數的時頻聯合分布實現
設BCG信號模糊函數的極坐標AF(r,Φ)，其中0≤Φ≤π，0≤r≤，L為BCG信號采樣點數。令r=pΔr，p=0，1，…，P-1,Φ=qΔΦ,q=0,1,…,Q-1,P、Q分別表示半徑和極角方向的采樣點數，Δr、ΔΦ表示采樣間隔。
| $\left\{ \begin{align} & \Delta r=2\sqrt{\pi /L},P=L/\sqrt{2} \\ & \Delta \Phi =\pi /L,Q=L \\ \end{align} \right.$ |
此時BCG信號的自適應核函數為
| ${{g}_{d}}\left( p,q \right)=exp\left( -\frac{{{(p\Delta r)}^{2}}}{2{{\sigma }^{2}}(q\Delta \Phi )~} \right)$ |
最優核函數(5)及約束條件(6)可寫為
| $max\sigma \sum\limits_{p=0}^{P-1}{\sum\limits_{q=0}^{Q-1}{p}}\left| A{{F}_{d}}\left( p,q \right)exp\left( -\frac{{{(p\Delta r)}^{2}}}{2{{\sigma }^{2}}(q\Delta \Phi )} \right) \right|$ |
| $\sum\limits_{q=0}^{Q-1}{{{\sigma }^{2}}(q\Delta \Phi )\le 4{{\pi }^{2}}\alpha }$ |
令
| $f=\sum\limits_{p=0}^{P-1}{\sum\limits_{q=0}^{Q-1}{p}}{{\left| A{{F}_{d}}\left( p,q \right)exp\left( -\frac{{{(p\Delta r)}^{2}}}{2{{\sigma }_{q}}^{2}} \right) \right|}^{2}}$ |
則
| $\sigma \left( k+1 \right)=\sigma \left( k \right)+\mu \left( k \right)\nabla f\left( k \right),$ |
其中μ(k)為迭代步長，?f(k)是函數f在σ(k)的梯度函數，其表達式如(10)和(11)。
| $\nabla f\left( k \right)={{\left[ \frac{\partial f}{\partial {{\sigma }_{0}}\left( k \right)},\ldots ,\frac{\partial f}{\partial {{\sigma }_{Q-1}}\left( k \right)} \right]}^{T}}$ |
| $\frac{\partial f}{\partial {{\sigma }_{q}}\left( k \right)}=\frac{2{{(\Delta r)}^{2}}}{{{\sigma }_{q}}^{3}\left( k \right)}\sum\limits_{p=0}^{P-1}{\sum\limits_{q=0}^{Q-1}{{{p}^{3}}}}{{\left| A{{F}_{d}}\left( p,q \right) \right|}^{2}}exp\left( -\frac{{{(p\Delta r)}^{2}}}{{{\sigma }_{q}}^{2}\left( k \right)} \right)$ |
對式(9)進行迭代，滿足約束條件(8)時，迭代終止，得到最優的σq,計算最優核函數。將BCG信號模糊函數AF(r,Φ)與最優核函數相乘，根據Cohen類分布的計算方法即可得到基于BCG信號的自適應最優核的時頻聯合分布。
2.2.3 基于時頻分布的BCG信號分量提取
信號的時頻聯合分布反映了信號能量的分布，由在不同時間具有相似分量的信號的時域及其時頻分布可知，時頻分布中存在能將時頻平面進行分段的極小值點，對其進行標記，可在時頻平面提取信號的分量，具體實現步驟：
(1)對BCG信號的自適應最優核函數時頻聯合分布，按照時間軸進行歸一化處理。
(2)提取時頻分布的頻率軸上的極大值，將所有極大值進行曲線擬合，得到時頻分布的包絡線。
(3)計算包絡線的局部極大值，擬合曲線的極大值smax(i),i=1,2,…,nmax。
(4)計算包絡線的局部極小值，擬合曲線的極小值smin(j),j=1,2,…,nmin。設定最小值數量的閾值N，重復步驟(4)，當nmin＞N時停止最小值的篩選。
(5)計算相鄰兩個smin的極大值點s′max(k),k=1,2,…,(nmin-1)，，并計算其中的最小值，記為sth，如果smax(i)≥sth，則記為s′min(l)，l=1,2,…,n′min。
(6)將s′min(l),l=1,2,…,n′min，對應到時域信號中，提取BCG信號分量。
2.3 EMD的IMF階數確定
根據EMD原理每個IMF的均值都為零，如果信號存在非零趨勢項，會對序數較大的模式分量產生影響，定義：
| ${{x}_{aver}}\left( t \right)=\frac{\sum\limits_{j=1}^{k}{\sum\limits_{i=1}^{M}{IM{{F}_{j}}({{t}_{i}})}}}{M},$ |
其中k為IMF序數，取值k∈[1,N]，i是采樣時刻，M是采樣點數。本文根據xaver(t)隨k變化的情況確定IMF重建階數，取使xaver(t)＞0時快速偏離零點的k值作為D，再用式(2)重構，獲得信號分量的低頻成分，最后按照時間順序連接重構后的低頻成分，實現BCG信號的降噪。
3 實驗研究
3.1 BCG信號采集
本實驗中利用型號為EB9301的香山人體電子秤，將其稱重傳感器改裝成傳感器電路，受試者站在改裝后的人體秤上，將傳感器的輸出端接在信號調理電路上，用示波器測量受試者在安靜狀態下的BCG信號。為了驗證測量BCG信號的準確性，實驗中同步采集了ECG信號，如圖 4所示。可以看出BCG信號中包含大量噪聲信號，但波形與ECG信號同步，呈現一定的規律性變化。
圖4
同步采集BCG和ECG信號的波形
Figure4.
Measuring BCG & ECG signal synchronously
3.2 基于時頻聯合分布的BCG信號分量提取
計算α=1時的BCG信號自適應最優核函數時頻聯合分布的包絡線，設閾值N=10，求取包絡線的極小值s′min和s′max極大值，所得的極小值將BCG信號分為與心動周期同步的6段，如圖 5所示。
圖5
BCG信號的自適應最優核函數時頻分布包絡線
Figure5.
Adaptive optimal kernel time-frequency distribution envelope of BCG signal
3.3 BCG信號分量的EMD分解中IMF階數確定
計算BCG信號分量IMF的xaver(t)與分解IMF個數k的關系，如圖 6所示。選擇使xaver(t)快速偏離零點的值作為D：對于分量1和分量5，在k=4時xaver(t)迅速偏離零點，其余分量則在k=3時xaver(t)迅速偏離原點。根據k值對各分量進行重構，消除噪聲成分。
圖6
BCG信號各分量的xaver(t)與k的關系曲線
Figure6.
Relation curves between xaver(t) and k of each BCG signal component
3.4 BCG信號降噪
圖 7顯示BCG信號降噪前后的波形，從中可以看到進行降噪處理后，BCG信號特征得到明顯的保留。計算降噪前后BCG信號的功率譜密度，如圖 8所示。由圖 8可知，濾波前后的BCG信號功率譜，在高頻部分平滑很多，且在低頻部分沒有明顯衰減，說明在去除噪聲的同時很好地保留了BCG信號的低頻特征。
圖7
降噪前后的BCG信號
Figure7.
BCG signals before and after de-noising
圖8
降噪前后BCG信號的PSD
Figure8.
PSD of BCG signals before and after adaptive filtering
3.5 降噪BCG的特征分析
利用實驗室的BCG信號檢測平臺對10名健康受試者(5名男性和5名女性)進行同步BCG信號及ECG信號測量，利用上述算法對BCG信號進行降噪處理，并提取J波及ECG信號的R波，計算其JJ間期、RR間期、心動周期數據如表 1所示。可以看出BCG信號檢測的心動周期與ECG信號檢測的心動周期基本一致，進一步說明了BCG信號對心動周期檢測具有合理性。
表1
健康受試者的BCG信號和ECG信號特征分析
Table1.
Normal person’s BCG and ECG features analysis
4 結論
為了有效地識別BCG信號特征，本文提出一種基于時頻聯合分布和EMD的BCG信號降噪方法。該方法先利用自適應最優核函數得到BCG信號的時頻分布，根據時頻分布與時域信號的對應關系，得到以心動周期為基準的信號分量，再應用EMD的濾波器原理去除了原始BCG信號中的噪聲成分。該方法克服了EMD處理在不同時間含有相同或相似頻率成分信號時的不足，實現了多周期BCG信號的降噪，取得了良好的仿真效果。
引言
心臟泵血活動會使人體產生微小震動，將這種微小震動記錄下來即稱為心沖擊(ballistocardiogram,BCG)信號。正常的BCG與心跳同步，具有重復性，主要包含H、I、J、K、M和N波。J波是BCG周期中幅值最大點，估計相鄰周期的J波間隔，即為心動周期；確定同周期IJ幅值，即為心臟泵血的收縮力大小，如圖 1所示。
圖1
一段正常的BCG信號及其特征
Figure1.
A segment of normal BCG signal and its features
相對于心電圖(electrocardiogram,ECG)信號，BCG信號可在無約束情況下從人體表面測得，其測量系統也可以設計在人們的日常生活家具中如床[1]、椅子[2]、人體秤[3]，因此，BCG信號適合日常生活的無感覺人體健康監測。
基于BCG信號的人體健康監測方法分兩個步驟實現。一是檢測BCG信號，二是提取BCG信號中所包含的生理參數或信息，進行分析評估。然而BCG信號自身十分微弱，且易受到呼吸、身體運動、工頻噪聲的干擾[4-5]，導致直接測量得到的BCG信號常常淹沒在噪聲中，無法獲得其包含的心率、呼吸等生理特征信息。因此，BCG信號降噪處理對于識別其特征參數是非常必要的。目前BCG信號降噪主要采用小波變換方法，其效果取決于小波閾值參數選擇[6]；而從信號分解角度，小波變換中基函數是固定的，不能匹配信號的全部。近年來，出現了一種數據驅動的非線性、非平穩信號處理方法--經驗模態分解(empirical mode decomposition,EMD)。它的基函數是一系列幅度和頻率均可變的正弦余弦函數，可通過信號自適應獲得，其特性依賴于信號本身，能取得很好的分解效果，被廣泛應用到語音[7]、圖像[8]以及ECG信號處理中[9]。
BCG信號是多周期非平穩隨機信號，隨心臟活動出現周期性變化，在不同的時間段有相似的頻率成分出現，直接利用EMD進行信號降噪，不能實現信號成分中預期的本征模態函數(intrinsic mode function,IMF)。為了有效識別用于人體健康監測的多周期BCG信號，本文提出一種基于時頻聯合分布和EMD的BCG信號降噪方法。先對BCG信號進行自適應最優核的時頻分布計算，然后在時頻平面內自適應提取以心動周期為基準的BCG信號分量，再對其進行EMD，最后按照時間順序連接重構低頻成分，實現多周期BCG信號降噪。
1 心臟BCG信號的采集
BCG信號采集裝置由傳感器電路和信號調理電路組成。傳感器電路嵌入改裝的家庭人體秤中，受試者靜止地站于人體秤上，心臟泵血時引起的身體微小震動即為BCG信號。BCG信號通過傳感器電路耦合后，再經差分前置放大電路、多級放大電路、隔直電路、陷波電路、帶通濾波電路，最終得到測量的BCG信號。采集電路的總體構圖如圖 2所示。其中差分放大電路抑制共模信號，隔直電路抑制直流噪聲，陷波電路抑制工頻干擾，帶通電路濾除測量中引發的干擾信號，多級放大實現BCG信號幅值放大。
圖2
BCG信號的采集電路
Figure2.
BCG signals recording circuit
2 提出的BCG信號降噪方法
2.1 基于EMD的降噪
對于給定信號x(t)，經EMD分解可表示為
| $x\left( t \right)=\sum\limits_{i=1}^{N}{IM{{F}_{i}}\left( t \right)+{{r}_{N}}\left( t \right)},$ |
其中rN(t)為趨勢項，IMF1(t)，…，IMFN(t)為信號的IMF。它們分別包含了信號從高到低不同頻率段的成分。基于EMD信號降噪的原理是對信號EMD分解后的多階IMF函數進行篩選，消除含有噪聲成分的IMF，保留含有信號成分的IMF并重構，完成降噪處理，其原理如圖 3所示。EMD降噪的關鍵IMF階數D的選取，通常從第D個IMF到最后一個IMF再加上趨勢項rN(t)即為所估計的低頻信號(t), 而
| $x\left( t \right)=\sum\limits_{i=1}^{N}{IM{{F}_{k}}\left( t \right)+{{r}_{N}}\left( t \right)}$ |
圖3
信號的EMD降噪原理
Figure3.
De-noising theory using EMD
對于非平穩信號更加完備的EMD分解可表達為
| $x\left( t \right)=\sum\limits_{k=1}^{M}{{{a}_{k}}\left( t \right)rect}\left( \frac{t-{{\tau }_{k}}}{{{b}_{k}}} \right)cos({{\theta }_{k}}\left( t \right))+r\left( t \right),$ |
其中ak(t)是IMFk(t)的包絡線，θk(t)是IMFk(t)的相位，rect((t-τk)/bk)是寬度為bk、中心為τk的矩形函數。信號x(t)可通過cos(θk(t))和rect((t-τk)/bk)在時頻平面上分別對時間和頻率進行定位。如果信號由一系列分量(分量定義的是信號的時頻平面上能量集中的區域)疊加而成，當每個分量相對于時間軸連續且頻率軸可分，EMD才能將其分解為多個IMF。BCG信號具備周期性，且在不同時間段具有相同成分，因此，要利用EMD分解多周期BCG信號，需提取BCG信號分量，再對其進行EMD分解，從而實現降噪處理。
2.2 基于時頻分布的BCG信號分量提取
2.2.1 核函數的定義及最優化條件
自適應最優核函數g(θ,τ)，其極坐標形式為
| $g\left( r,\Phi \right)=exp\left( -\frac{{{r}^{2}}}{2{{\sigma }^{2}}\left( \Phi \right)} \right),$ |
其中0≤g(r,Φ)≤1，r=為極徑。
對于給定信號及其模糊函數，最優核函數g(r,Φ)應滿足優化問題：
| $maxg\int_{0}^{2\pi }{\int_{0}^{\infty }{{{\left| AF\left( r,\Phi \right)g\left( r,\Phi \right) \right|}^{2}}rdrd\Phi }}\text{ }$ |
約束條件為
| $\frac{1}{4}{{\pi }^{2}}\int_{0}^{2\pi }{\int_{0}^{\infty }{{{\left| g\left( r,\Phi \right) \right|}^{2}}rdrd\Phi \le \alpha ~,}}$ |
其中AF(r,Φ)為待分析信號模糊函數的極坐標表示，α≥0為核函數的能量體積。
2.2.2 BCG信號最優核函數的時頻聯合分布實現
設BCG信號模糊函數的極坐標AF(r,Φ)，其中0≤Φ≤π，0≤r≤，L為BCG信號采樣點數。令r=pΔr，p=0，1，…，P-1,Φ=qΔΦ,q=0,1,…,Q-1,P、Q分別表示半徑和極角方向的采樣點數，Δr、ΔΦ表示采樣間隔。
| $\left\{ \begin{align} & \Delta r=2\sqrt{\pi /L},P=L/\sqrt{2} \\ & \Delta \Phi =\pi /L,Q=L \\ \end{align} \right.$ |
此時BCG信號的自適應核函數為
| ${{g}_{d}}\left( p,q \right)=exp\left( -\frac{{{(p\Delta r)}^{2}}}{2{{\sigma }^{2}}(q\Delta \Phi )~} \right)$ |
最優核函數(5)及約束條件(6)可寫為
| $max\sigma \sum\limits_{p=0}^{P-1}{\sum\limits_{q=0}^{Q-1}{p}}\left| A{{F}_{d}}\left( p,q \right)exp\left( -\frac{{{(p\Delta r)}^{2}}}{2{{\sigma }^{2}}(q\Delta \Phi )} \right) \right|$ |
| $\sum\limits_{q=0}^{Q-1}{{{\sigma }^{2}}(q\Delta \Phi )\le 4{{\pi }^{2}}\alpha }$ |
令
| $f=\sum\limits_{p=0}^{P-1}{\sum\limits_{q=0}^{Q-1}{p}}{{\left| A{{F}_{d}}\left( p,q \right)exp\left( -\frac{{{(p\Delta r)}^{2}}}{2{{\sigma }_{q}}^{2}} \right) \right|}^{2}}$ |
則
| $\sigma \left( k+1 \right)=\sigma \left( k \right)+\mu \left( k \right)\nabla f\left( k \right),$ |
其中μ(k)為迭代步長，?f(k)是函數f在σ(k)的梯度函數，其表達式如(10)和(11)。
| $\nabla f\left( k \right)={{\left[ \frac{\partial f}{\partial {{\sigma }_{0}}\left( k \right)},\ldots ,\frac{\partial f}{\partial {{\sigma }_{Q-1}}\left( k \right)} \right]}^{T}}$ |
| $\frac{\partial f}{\partial {{\sigma }_{q}}\left( k \right)}=\frac{2{{(\Delta r)}^{2}}}{{{\sigma }_{q}}^{3}\left( k \right)}\sum\limits_{p=0}^{P-1}{\sum\limits_{q=0}^{Q-1}{{{p}^{3}}}}{{\left| A{{F}_{d}}\left( p,q \right) \right|}^{2}}exp\left( -\frac{{{(p\Delta r)}^{2}}}{{{\sigma }_{q}}^{2}\left( k \right)} \right)$ |
對式(9)進行迭代，滿足約束條件(8)時，迭代終止，得到最優的σq,計算最優核函數。將BCG信號模糊函數AF(r,Φ)與最優核函數相乘，根據Cohen類分布的計算方法即可得到基于BCG信號的自適應最優核的時頻聯合分布。
2.2.3 基于時頻分布的BCG信號分量提取
信號的時頻聯合分布反映了信號能量的分布，由在不同時間具有相似分量的信號的時域及其時頻分布可知，時頻分布中存在能將時頻平面進行分段的極小值點，對其進行標記，可在時頻平面提取信號的分量，具體實現步驟：
(1)對BCG信號的自適應最優核函數時頻聯合分布，按照時間軸進行歸一化處理。
(2)提取時頻分布的頻率軸上的極大值，將所有極大值進行曲線擬合，得到時頻分布的包絡線。
(3)計算包絡線的局部極大值，擬合曲線的極大值smax(i),i=1,2,…,nmax。
(4)計算包絡線的局部極小值，擬合曲線的極小值smin(j),j=1,2,…,nmin。設定最小值數量的閾值N，重復步驟(4)，當nmin＞N時停止最小值的篩選。
(5)計算相鄰兩個smin的極大值點s′max(k),k=1,2,…,(nmin-1)，，并計算其中的最小值，記為sth，如果smax(i)≥sth，則記為s′min(l)，l=1,2,…,n′min。
(6)將s′min(l),l=1,2,…,n′min，對應到時域信號中，提取BCG信號分量。
2.3 EMD的IMF階數確定
根據EMD原理每個IMF的均值都為零，如果信號存在非零趨勢項，會對序數較大的模式分量產生影響，定義：
| ${{x}_{aver}}\left( t \right)=\frac{\sum\limits_{j=1}^{k}{\sum\limits_{i=1}^{M}{IM{{F}_{j}}({{t}_{i}})}}}{M},$ |
其中k為IMF序數，取值k∈[1,N]，i是采樣時刻，M是采樣點數。本文根據xaver(t)隨k變化的情況確定IMF重建階數，取使xaver(t)＞0時快速偏離零點的k值作為D，再用式(2)重構，獲得信號分量的低頻成分，最后按照時間順序連接重構后的低頻成分，實現BCG信號的降噪。
3 實驗研究
3.1 BCG信號采集
本實驗中利用型號為EB9301的香山人體電子秤，將其稱重傳感器改裝成傳感器電路，受試者站在改裝后的人體秤上，將傳感器的輸出端接在信號調理電路上，用示波器測量受試者在安靜狀態下的BCG信號。為了驗證測量BCG信號的準確性，實驗中同步采集了ECG信號，如圖 4所示。可以看出BCG信號中包含大量噪聲信號，但波形與ECG信號同步，呈現一定的規律性變化。
圖4
同步采集BCG和ECG信號的波形
Figure4.
Measuring BCG & ECG signal synchronously
3.2 基于時頻聯合分布的BCG信號分量提取
計算α=1時的BCG信號自適應最優核函數時頻聯合分布的包絡線，設閾值N=10，求取包絡線的極小值s′min和s′max極大值，所得的極小值將BCG信號分為與心動周期同步的6段，如圖 5所示。
圖5
BCG信號的自適應最優核函數時頻分布包絡線
Figure5.
Adaptive optimal kernel time-frequency distribution envelope of BCG signal
3.3 BCG信號分量的EMD分解中IMF階數確定
計算BCG信號分量IMF的xaver(t)與分解IMF個數k的關系，如圖 6所示。選擇使xaver(t)快速偏離零點的值作為D：對于分量1和分量5，在k=4時xaver(t)迅速偏離零點，其余分量則在k=3時xaver(t)迅速偏離原點。根據k值對各分量進行重構，消除噪聲成分。
圖6
BCG信號各分量的xaver(t)與k的關系曲線
Figure6.
Relation curves between xaver(t) and k of each BCG signal component
3.4 BCG信號降噪
圖 7顯示BCG信號降噪前后的波形，從中可以看到進行降噪處理后，BCG信號特征得到明顯的保留。計算降噪前后BCG信號的功率譜密度，如圖 8所示。由圖 8可知，濾波前后的BCG信號功率譜，在高頻部分平滑很多，且在低頻部分沒有明顯衰減，說明在去除噪聲的同時很好地保留了BCG信號的低頻特征。
圖7
降噪前后的BCG信號
Figure7.
BCG signals before and after de-noising
圖8
降噪前后BCG信號的PSD
Figure8.
PSD of BCG signals before and after adaptive filtering
3.5 降噪BCG的特征分析
利用實驗室的BCG信號檢測平臺對10名健康受試者(5名男性和5名女性)進行同步BCG信號及ECG信號測量，利用上述算法對BCG信號進行降噪處理，并提取J波及ECG信號的R波，計算其JJ間期、RR間期、心動周期數據如表 1所示。可以看出BCG信號檢測的心動周期與ECG信號檢測的心動周期基本一致，進一步說明了BCG信號對心動周期檢測具有合理性。
表1
健康受試者的BCG信號和ECG信號特征分析
Table1.
Normal person’s BCG and ECG features analysis
4 結論
為了有效地識別BCG信號特征，本文提出一種基于時頻聯合分布和EMD的BCG信號降噪方法。該方法先利用自適應最優核函數得到BCG信號的時頻分布，根據時頻分布與時域信號的對應關系，得到以心動周期為基準的信號分量，再應用EMD的濾波器原理去除了原始BCG信號中的噪聲成分。該方法克服了EMD處理在不同時間含有相同或相似頻率成分信號時的不足，實現了多周期BCG信號的降噪，取得了良好的仿真效果。
