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

【作 者】王旭，蔡坤

1 廣東工業(yè)大學(xué)自動(dòng)化學(xué)院，廣州市，510006

2 華南農(nóng)業(yè)大學(xué)電子工程學(xué)院，廣州市，510642

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

【作 者】王旭1，蔡坤2

1 廣東工業(yè)大學(xué)自動(dòng)化學(xué)院，廣州市，510006

2 華南農(nóng)業(yè)大學(xué)電子工程學(xué)院，廣州市，510642

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

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

0 引言

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

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

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

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

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

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

1.1 信號預(yù)處理

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

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

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

1.3 信號的短時(shí)傅里葉變換

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

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

t 時(shí)刻處短時(shí)傅里葉變換的計(jì)算過程如下：

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

(2) 利用平移后的分析窗w(τ-t)對原信號做加窗截?cái)嗵幚?，得到短時(shí)信號xi(τ)=x(τ)w(τ-t)；

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

我們的目的在于定位胎兒心電的R波，那么就需要選取一種合適且穩(wěn)定的窗函數(shù)。根據(jù)胎兒心電R波和各種窗函數(shù)的性質(zhì)，經(jīng)過對幾種窗函數(shù)的實(shí)驗(yàn)和比較，選定窗函數(shù)為漢明窗，設(shè)窗口長度為P，P代表一個(gè)胎兒心電周期的數(shù)據(jù)點(diǎn)數(shù)，設(shè)計(jì)過程為：

首先在[-π, π]區(qū)間上均勻產(chǎn)生P個(gè)數(shù)據(jù)點(diǎn)，保存為向量T=[t1, t2, ..., tp]，漢明窗函數(shù)為：

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

1.3 盲分離算法

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

圖2 窗函數(shù)Fig.2 Window function

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

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

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

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

1.5 胎兒瞬時(shí)心率的計(jì)算

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

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

2 仿真實(shí)驗(yàn)與結(jié)果分析

基于短時(shí)傅里葉變換和盲分離的胎兒心率檢測方法在MATLAB 2012b上編程實(shí)現(xiàn)，在Windows 7環(huán)境中運(yùn)行。

2.1 算法對比

我們現(xiàn)在來討論一種基于小波變換的心電R波檢測算法，選取標(biāo)準(zhǔn)心電數(shù)據(jù)庫中的某一例信號進(jìn)行實(shí)驗(yàn)，實(shí)現(xiàn)過程如下所示：

(1) 對原始數(shù)據(jù)利用梳狀濾波進(jìn)行預(yù)處理，數(shù)據(jù)長度為3 584個(gè)；

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

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

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

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

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

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

圖4 兩種算法檢測結(jié)果對比圖（實(shí)線：小波變換算法，虛線：本文使用的算法）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結(jié)果對比分析得出如下結(jié)論：(1)本文使用的算法檢測率整體高于小波變換算法的檢測率；(2)當(dāng)胎兒心電與母親心電重合時(shí)，兩種算法都會產(chǎn)生漏值點(diǎn)，而本文所使用的算法因進(jìn)行特征提取，降低了漏值率；(3)對比檢測率，胎兒心電信號幅度較小即信噪比較低時(shí)，小波變換算法無法檢測和判斷，而本文所使用的算法，極大地提高了檢測率。

2.2 臨床實(shí)驗(yàn)

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

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

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

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

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

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

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

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

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

3 結(jié)論

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

[1] 賈文娟. 胎兒心電信號檢測方法及監(jiān)護(hù)系統(tǒng)的研究[D].北京:北京工業(yè)大學(xué), 2011.

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

[3] 申麗巖. 非入侵式胎兒心電信號提取方法研究及其DSP系統(tǒng)實(shí)現(xiàn)[D]. 北京: 北京工業(yè)大學(xué), 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] 劉森, 周禮杲, 楊福生. 應(yīng)用自適應(yīng)噪聲抵消系統(tǒng)作胎兒心電信號處理以實(shí)現(xiàn)胎兒監(jiān)護(hù)[J]. 中國生物醫(yī)學(xué)工程學(xué)報(bào), 1985, 4(4): 220-229.

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

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

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

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

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

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

[12] 皇甫堪, 陳建文. 現(xiàn)代數(shù)字信號處理[M].北京: 電子工業(yè)出版社, 2003.

[13] 馮愛玲. 基于短時(shí)傅里葉變換的胎心率檢測算法與實(shí)現(xiàn)[D]. 廣州: 廣東工業(yè)大學(xué), 2014.

[14] 楊福生, 洪波.獨(dú)立分量分析的原理與應(yīng)用[M]. 北京:清華大學(xué)出版社, 2006.

[15] 于宣福, 嚴(yán)文鴻, 劉輝. 一種基于香農(nóng)包絡(luò)的胎兒瞬時(shí)心率檢測方法[J]. 電腦編程技巧與維護(hù), 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.