基于短時傅里葉變換和盲分離的胎兒心率檢測方法

【作 者】王旭，蔡坤

1 廣東工業大學自動化學院，廣州市，510006

2 華南農業大學電子工程學院，廣州市，510642

基于短時傅里葉變換和盲分離的胎兒心率檢測方法

【作 者】王旭1，蔡坤2

1 廣東工業大學自動化學院，廣州市，510006

2 華南農業大學電子工程學院，廣州市，510642

胎兒心率的變化是循環系統和中樞神經系統機能調節的表現，在圍產期對胎兒進行胎心率檢測具有重要的意義，基于此該文提出了一種基于短時傅里葉變換和盲分離的胎兒心率檢測方法。首先對混疊心電信號進行預處理，然后運用小波變換技術分離出含有噪聲的胎兒心電信號，再對它進行短時傅里葉變換和盲分離，之后計算它們的相關系數，最后選取與原信號相關性最強的一個獨立分量進行峰值檢測并計算出胎兒瞬時心率。實驗結果表明：該方法能夠提高胎兒峰值(R波)檢測率，而且在信噪比較低的情況下，它對胎兒峰值(R波)定位具有較高的準確性。

胎兒心電信號；瞬時心率；短時傅里葉變換；盲信號分離

0 引言

胎兒心電圖是一種檢測胎兒在子宮內健康狀況的重要方法，能夠反映胎兒心臟的心率，而胎兒心率的變化是循環系統和中樞神經系統機能調節的表現，因此在圍產期對胎兒進行胎心率檢測可以了解胎兒在子宮內的健康狀況[1]。然而由于采集的胎兒心電信號幅度小，信號的噪聲比低且大都混雜有幅度大、分布廣的噪聲干擾，例如常見的工頻干擾，呼吸、肌電引起的噪聲等，它給醫學診斷帶來極大的困難。在時域和頻域上，這些干擾信號混疊在一起對胎兒心電信號的瞬時心率計算造成很大的影響[2]。因此，研究如何準確、便捷，有效地測量圍產期胎心電的瞬時心率具有及其重要的實用價值和臨床意義。

目前已經有許多提取胎兒心電信號(FECG)方法的報道。如匹配濾波法[3]，它易于實現，結構簡單，但是計算的標準胎兒心電周期不準，對胎兒心電的正確識別率較低。如神經網絡算法[4]，它提取胎兒性能較好，但是此算法需要較長的訓練時間，收斂速度慢，不適合臨床應用。如自適應濾波法[5-6]，它可以在沒有或者只有很少信號先驗統計知識的情況下提取胎兒心電，該方法用宮底電極輸入作為參考信號，基本上沒有胎兒心電，這樣既可以消除母親心電又能盡可能抵消肌電噪聲，但是電極之間的延時會影響到算法的準確率。此外，自適應濾波算法很大程度上依賴參考信號與理想信號的獨立性及兩信號的相關性，這會影響到自適應噪聲抵消的效果。如小波分析法[7]，它提取到的胎兒心電信號準確性較高，效果較理想，但是對于不同的數據參數值會有很大的改變，計算量大。

前面有學者運用小波變換技術[8]，提出基于小波變換的胎兒心電信號提取方法，本文的算法是在此方法的基礎上進行研究和改進，改進之處在于：(1)運用小波變換技術從混疊信號中分離出含有噪聲的胎兒心電信號，去除了混疊信號中較大的母親心電信號干擾；(2)短時傅里葉變換中窗函數具有特征識別的作用，合適的窗函數可以較好地識別胎兒心電的QRS波，準確地計算胎兒心電速率；(3)運用盲分離理論中的特征提取的特點可以更好地提取QRS波，提高準確率。

1 基于短時傅里葉變換和盲分離的胎兒心電速率檢測方法

基于短時傅里葉變換和盲分離的胎兒心電速率檢測方法流程圖如圖1所示。

圖1 方法流程圖Fig.1 Method flow chart

1.1 信號預處理

梳狀濾波器可以在保持信號帶寬不變的情況下，使其在0 Hz、50 Hz及其高次諧波處有很窄的阻帶，這樣就可以運用梳狀濾波器消除50 Hz及其高次諧波的工頻干擾和頻率較低的基線漂移。

1.2 小波變換分離胎兒心電信號

小波變換可以等效為一組濾波器，信號通過一個分解高通濾波器和分解低通濾波器，對應輸出高頻分量和低頻分量，即成為細節分量和近似分量。小波重構就是將分解之后的近似信號與細節信號疊加得到原始信號[9]。由于db小波函數雖然不具有對稱性，但具有緊支集正交性，能夠比較準確地進行信號重構。前面有學者運用db2小波基進行四尺度小波變換得到較好的處理效果[10]，此文采用db2小波基進行小波變換分離胎兒心電信號。源心電信號進行小波變換得到近似信號和細節信號，在每個尺度信號中進行峰值檢測和閾值處理，得到只含有母親心電信息的各尺度信號，然后利用小波重構得到母親心電信號，再用源心電信號減去母親心電信號即可得到含有噪聲的胎兒心電信號。

1.3 信號的短時傅里葉變換

短時傅里葉變換（STFT）其主要思想是將信號加窗，將加窗后的信號進行傅里葉變換，加窗后使得變換為時間t附近很小時間上的局部譜，窗函數可以根據t的位置變化在整個時間軸上平移，利用窗函數可以得到任意位置附近的時間段頻譜，在時間域上實現信息的局域化，突出局部化的信息[11]。

設信號為x(t)，t∈(-∞，+∞)，分析窗函數為s(t)，則非平穩信號x(t)的連續短時傅里葉變換定義[12]為：

t 時刻處短時傅里葉變換的計算過程如下：

(1) 將分析窗w(τ)從時間零處平移到時間t處，得到w(τ-t)；

(2) 利用平移后的分析窗w(τ-t)對原信號做加窗截斷處理，得到短時信號xi(τ)=x(τ)w(τ-t)；

(3) 對短時信號xi(τ)進行傅里葉變換分析得到傅里葉頻譜[13]。

我們的目的在于定位胎兒心電的R波，那么就需要選取一種合適且穩定的窗函數。根據胎兒心電R波和各種窗函數的性質，經過對幾種窗函數的實驗和比較，選定窗函數為漢明窗，設窗口長度為P，P代表一個胎兒心電周期的數據點數，設計過程為：

首先在[-π, π]區間上均勻產生P個數據點，保存為向量T=[t1, t2, ..., tp]，漢明窗函數為：

式中a+b=1，n=M+1，a，b為調節函數的參數，M為窗口中QRS的長度，再利用漢明窗函數得到向量X=[x1, x2, ..., xn]，對X在向量T上進行DTFT計算得到Z=[z1, z2, ..., zn]，取Z的幅值作為所需窗的值Y=[y1, y2, ..., yn]；最后調節參數a, b, M，使得窗函數主瓣的長度約等于QRS波的長度，滿足此條件的窗函數即為本文中所需要的窗函數，如圖2所示。

1.3 盲分離算法

盲分離是一種多維信號處理方法，它指在未知原信號以及混合模型也未知的情況下，僅從觀測信號中恢復出源信號各個獨立分量的過程。設源信號為n個相互獨立的源信號s=(s1, s2, ..., sn)T，通過一個線性瞬時矩陣A得到的是m個混合信號x=(x1, x2, ..., xn)T，矩陣表示為x=As，其中s∈Rn×1為源信號，A∈Rm×n為未知混疊矩陣，x為觀測到的混疊信號。盲分離的目的就是在源信號s和混合矩陣A未知的情況下，從混合信號x中分離出y，即找到一個分離矩陣B=(b1, b2, ..., bm)，使得y =Bx，y=(y1, y2, ..., yn)T為源信號s的估計。變換模型如圖3所示。

圖2 窗函數Fig.2 Window function

圖3 盲分離模型Fig.3 Model of blind separation

其中，s為短時傅里葉變換之后的信號，s和線性混合矩陣A都是未知的，B為待求的分離矩陣。源信號的估計矢量y即為獨立分量分析后得到的多個獨立分量，已知的觀測矢量x為不同頻率的胎心電和噪聲的混合信號。y中的某一行數據為分離得到的胎兒心電信號的估計。本文盲分離算法采用FastICA算法[14]。

1.4 胎兒心電信號分量的提取

胎兒心電信號的提取就是從盲分離得到的多路獨立分量中提取出與胎兒心電信號最相近的一路信號。將盲分離得到的多路信號降維至m路信號，將m路信號與y=(y1, y2, ..., yn)T進行相關系數的計算，選取相關性最大的獨立分量作為以QRS波為主要信息的胎兒心電信號yx。

1.5 胎兒瞬時心率的計算

首先對胎兒心電信號采用峰值檢測[15]，得到該曲線所有的峰值點，然后對所有峰值點依次進行前向差分，再根據瞬時心率公式求得胎兒的瞬時心率。瞬時心率的定義如下：

其中，△t為每兩相鄰峰值的時間間隔，所求的υ即為瞬時心率。

2 仿真實驗與結果分析

基于短時傅里葉變換和盲分離的胎兒心率檢測方法在MATLAB 2012b上編程實現，在Windows 7環境中運行。

2.1 算法對比

我們現在來討論一種基于小波變換的心電R波檢測算法，選取標準心電數據庫中的某一例信號進行實驗，實現過程如下所示：

(1) 對原始數據利用梳狀濾波進行預處理，數據長度為3 584個；

(2) 進行小波變換，利用db2小波進行多尺度分解得到第四層細節信號S4；

(3) 利用峰值檢測函數進行峰值檢測，并設定母親心電信號閾值，檢測出母親心電信號；

(4) 細節信號S4減去母親心電信號得到含有噪聲胎兒心電信號；

(5) 再進行峰值檢測，設定胎兒心電信號閾值，定位胎兒心電QRS波；

(6) 最后，利用瞬時心率公式計算胎兒心電速率。

對MIT-BIH數據庫的10例信號用基于小波變換的心電R波檢測算法和本文所用的方法進行實驗，每例數據長度為40 960。實驗結果對比如圖4所示。

圖4 兩種算法檢測結果對比圖（實線：小波變換算法，虛線：本文使用的算法）Fig.4 The comparison chart of two algorithm detection result（the solid line: the algorithm of wavelet transform, the dotted line: the algorithm used by the paper）

由圖4結果對比分析得出如下結論：(1)本文使用的算法檢測率整體高于小波變換算法的檢測率；(2)當胎兒心電與母親心電重合時，兩種算法都會產生漏值點，而本文所使用的算法因進行特征提取，降低了漏值率；(3)對比檢測率，胎兒心電信號幅度較小即信噪比較低時，小波變換算法無法檢測和判斷，而本文所使用的算法，極大地提高了檢測率。

2.2 臨床實驗

本文的心電數據來自于MIT-BIH數據庫，實際處理過程如下所示：

第一步，梳狀濾波消除50 Hz及其高次諧波的工頻干擾以及基線漂移得到信號Sl，然后利用小波變換技術分離出母親心電信號Sm，再將信號Sl減去母親心電信號Sm得到含有噪聲的胎兒心電信號Sf，見圖5。

圖5 混合心電信號，母親心電信號，包含噪聲的胎兒心電信號Fig.5 Mixed ECG, mother ECG, fetal ECG with noise

第二步，首先根據窗函數設定方法得到合適的窗函數，然后對信號Sf進行加窗處理，再進行傅里葉變換，得到信號其中三路信號見圖6。

圖6 短時傅里葉變換后的第二、三、五路信號Fig.6 The second, third, fifth signal after STFT

第三步，利用FastICA算法對信號sstft進行盲分離，得到信號Sn=(s1, s2, ..., sn)T，降維，再計算信號Sn與信號Sf的相關系數，選取相關性最大的獨立分量作為以QRS波為主要信息的胎兒心電信號sx，見圖7。

圖7 盲分離之后的三路信號Fig.7 Three signals after blind separate

第四步，對信號Sx進行峰值檢測，然后依次進行前向差分，再根據瞬時心率公式求得胎兒的瞬時心率，見圖8。

圖8 以QRS波為主要信息的胎兒心電信號，胎兒心電瞬時心率Fig.8 Fetal ECG signal that QRS wave as the main information, instantaneous heart rate of fetal ECG signal

3 結論

本文提出了基于短時傅里葉變換和盲分離的胎兒心率檢測方法，此方法的目的在于定位胎兒心電QRS波，從而計算胎兒瞬時心率。實驗結果表明，運用小波變換技術去除較大的干擾信號母親心電信號，短時傅里葉變換設計合適的窗函數識別QRS波，FastICA盲分離算法提取QRS波的特征，不僅提高峰值(R波)的檢測率，而且在低信噪比的條件下能夠準確地識別QRS波。

[1] 賈文娟. 胎兒心電信號檢測方法及監護系統的研究[D].北京:北京工業大學, 2011.

[2] 陳壽齊, 沈越泓, 許魁. 噪聲背景下的胎兒心電提取[J]. 數據采集與處理, 2010, 25(2): 292.

[3] 申麗巖. 非入侵式胎兒心電信號提取方法研究及其DSP系統實現[D]. 北京: 北京工業大學, 2006: 2-3.

[4] Agarwal N,Prasad DV, Swarnalatha R. Extraction of fetal electrocardiographic signals using neural network[C]//WCB 2010, ICBME and APBiomec, 2010: 1350-1353.

[5] 劉森, 周禮杲, 楊福生. 應用自適應噪聲抵消系統作胎兒心電信號處理以實現胎兒監護[J]. 中國生物醫學工程學報, 1985, 4(4): 220-229.

[6] 李章勇, 謝正祥, 牛永紅. 自適應干擾對消技術提取胎兒心電的可視化仿真實現[J]. 上海生物醫學工程, 2001, 22(4): 13-17.

[7] 賈文娟, 吳水才, 白燕萍, 等. 小波變換模極大值算法用于胎兒心電信號提取的研究[J]. 醫療衛生裝備, 2010, 32(12): 14-17.

[8] 李卓妮. 成人及胎兒心電信號R波檢測算法的研究[D]. 廣州: 華南理工大學, 2013: 46-53.

[9] 張德豐.MATLAB小波分析[M]. 北京: 清華大學出版社, 2009.

[10] 嚴文鴻, 蔣寧. 基于小波變換和匹配濾波的胎兒心電信號R波檢測[J]. 中國醫療器械雜志, 2015, 39(5): 318-320.

[11] 祁才君. 數字信號處理技術的算法分析與應用[M]. 北京: 機械工業出版社, 2005.

[12] 皇甫堪, 陳建文. 現代數字信號處理[M].北京: 電子工業出版社, 2003.

[13] 馮愛玲. 基于短時傅里葉變換的胎心率檢測算法與實現[D]. 廣州: 廣東工業大學, 2014.

[14] 楊福生, 洪波.獨立分量分析的原理與應用[M]. 北京:清華大學出版社, 2006.

[15] 于宣福, 嚴文鴻, 劉輝. 一種基于香農包絡的胎兒瞬時心率檢測方法[J]. 電腦編程技巧與維護, 2015, (11): 75-76.

Detection of Heart Rate of Fetal ECG Based on STFT and BSS

【 Writers 】WANG Xu1, CAI Kun2
1 School of Automation, Guangdong University of Technology, Guangzhou, 510006
2 School of Electronic Engineering, South China Agricultural University, Guangzhou, 510642

fetal ECG, instantaneous heart rate, short-time Fourier transform, blind source separation

TH772.2

A

10.3969/j.issn.1671-7104.2016.01.006

1671-7104(2016)01-0022-05

2015-10-21

王旭，E-mail: 954582137@qq.com

【 Abstract 】Changes in heart rate of fetal is function regulating performance of the circulatory system and the central nervous system , it is significant to detect heart rate of fetus in perinatal fetal. This paper puts forward the fetal heart rate detection method based on short time Fourier transform and blind source separation. First of all, the mixed ECG signal was preprocessed, and then the wavelet transform technique was used to separate the fetal ECG signal with noise from mixed ECG signal, after that, the short-time Fourier transform and the blind separation were carried on it, and then calculated the correlation coefficient of it, Finally, An independent component that it has strongest correlation with the original signal was selected to make FECG peak detection and calculated the fetal instantaneous heart rate. The experimental results show that the method can improve the detection rate of the FECG peak (R), and it has high accuracy in fixing peak(R) location in the case of low signal-noise ratio.