【作 者】王旭，蔡坤

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.