CN104873222B - 一种基于智能移动终端的胎心率检测方法 - Google Patents

Abstract

本发明公开了一种基于智能移动终端的胎心率检测方法，其利用智能移动终端对通过音频接口输入的胎心信号进行处理，包括A/D采样、降采样、整流、滑动平均滤波、有偏估计自相关求取、无偏估计自相关求取、胎心率值判断选取等操作之后获得胎心率。本发明利用智能移动终端强大的数据处理能力，运行改进的胎心率检测算法，提高测量精度，减少反应时间，增加检测结果稳定性；同时精简了超声多普勒检测终端的硬件结构，降低了功耗和成本。同时本发明对传统自相关算法进行改进，结合无偏估计的自相关算法来加强主峰的正确判断，克服了易出现倍频的缺点，提高了算法的准确度和稳定性。

Description

一种基于智能移动终端的胎心率检测方法

技术领域

本发明属于医疗设备技术领域，具体涉及一种基于智能移动终端的胎心率检测方法。

背景技术

胎心检测是通过监测胎心率的变化来评价胎儿宫内安危的一种常用的胎儿监护方法，有助于提高胎儿分娩质量，降低胎儿畸形率和死亡率。国际上最常用的一种无创监测方法即超声多普勒法，它具有探测灵敏度高、使用方便、对胎儿影响低等优点。

自相关法是超声多普勒测量胎儿心率中最常用的方法，其利用自相关函数的周期与信号周期一致的特点，找出自相关函数的周期,即可得到胎儿心率值。信号的自相关函数是信号的延时信号与该信号本身之间的相似性度量，自相关函数具有以下的特点：

1)周期信号的自相关函数也是周期的,且两者周期相等；

2)蕴含在有用信号中噪声的自相关函数，随着延时值的增加,会出现比有用信号自相关函数更快的衰减。

当延时量T等于0或等于整数倍的信号周期，信号的自相关函数在T点处取得峰值。这样，通过寻找信号自相关函数中极大值点的位置，并计算出两个相邻极大值点的序列长度，就可以获得被噪声污染的信号中蕴含的周期性。

由于胎儿本身的心跳比较弱，再加上母体环境的干扰，超声回波信号成分比较复杂，干扰大。因此，采集的信号有时干扰较多，此时，单纯采用自相关算法得到的自相关函数，在周期检测时，容易造成周期的加倍或减半。

此外，现有的胎心仪受到外型大小、成本、功耗等因素的限制，其处理模块往往由简单的MCU构成，无法做到大数据处理，因而数据处理速度受限，精度达不到很高的要求。

发明内容

针对现有技术所存在的上述技术问题，本发明提供了一种基于智能移动终端的胎心率检测方法，依托智能移动终端强大的数据处理能力，对胎心信号进行处理得到准确、可靠的胎儿心率，特别是在信号差的情况下。

一种基于智能移动终端的胎心率检测方法，包括如下步骤：

(1)利用智能移动终端通过音频接口接收胎心仪提供的胎心信号；

(2)对所述的胎心信号进行预处理得到胎心多普勒包络信号序列；

(3)对所述的胎心多普勒包络信号序列进行自相关分析，得到有偏自相关序列和无偏自相关序列；

(4)对所述的有偏自相关序列和无偏自相关序列进行对比判断处理，以求取得到胎心率。

所述的步骤(2)至(4)均由智能移动终端内部处理器实现处理；所述的智能移动终端可以为智能手机或平板电脑。

所述的步骤(2)对胎心信号进行预处理的具体步骤如下：

2.1对所述的胎心信号进行A/D采样，得到胎心信号的采样值序列；

2.2对所述的采样值序列进行降采样处理，得到降采样值序列；

2.3对所述的降采样值序列进行绝对值处理，得到整流值序列；

2.4对所述的整流值序列进行滑动平均滤波处理，得到胎心多普勒包络信号序列。

所述的步骤2.4中通过以下公式对整流值序列进行滑动平均滤波处理：

$D\left(n\right)=\frac{1}{m}\underset{i=0}{\overset{m-1}{\mathrm{\Σ}}}C(n+i)$

其中：C(n+i)为整流值序列中的第n+i个数值，D(n)为胎心多普勒包络信号序列中的第n个信号值，n为大于0的自然数，m为滑动信号窗口的预设长度。

所述的步骤(3)对胎心多普勒包络信号序列进行自相关分析的具体步骤如下：

3.1对所述的胎心多普勒包络信号序列进行快速傅里叶变换运算，得到复数序列；

3.2对所述的复数序列进行模平方运算处理，得到对应的模平方序列；

3.3对所述的模平方序列进行快速傅里叶反变换运算，得到有偏自相关序列；

3.4根据所述的有偏自相关序列计算出对应的无偏自相关序列。

所述的步骤3.2中通过以下公式对复数序列进行模平方运算处理：

R(k)＝D*(k)D(k)＝|D(k)|²

其中：R(k)为模平方序列中的第k个数值，D(k)为复数序列中的第k个复数，D^*(k)为D(k)的共轭复数，|D(k)|为D(k)的模，k为自然数且1≤k≤N，N为胎心多普勒包络信号序列的信号值个数。

所述的步骤3.4中通过以下公式计算无偏自相关序列：

$R_{\mathrm{XX}2}\left(j\right)=\frac{N-j}{N}*R_{\mathrm{XX}1}\left(j\right)$

其中：R_xx1(j)和R_xx2(j)分别为有偏自相关序列和无偏自相关序列中的第j个数值，N为胎心多普勒包络信号序列的信号值个数，j为自然数且1≤j≤N。

所述的步骤(4)对有偏自相关序列和无偏自相关序列进行对比判断处理的具体过程如下：

4.1从有偏自相关序列中确定出两个相邻峰值之间的序列间隔长度n₁；

4.2从无偏自相关序列中确定出两个相邻峰值之间的序列间隔长度n₂；

4.3根据所述的序列间隔长度n₁和n₂，计算出对应的时间间隔T₁和T₂；

4.4根据所述的时间间隔T₁和T₂，判断求取得到所述的胎心率。

所述的步骤4.3中使所述的序列间隔长度n₁和n₂分别除以降采样的频率，得到对应的时间间隔T₁和T₂。

所述的步骤4.4中根据以下关系判断求取胎心率：

若T₁＝T₂，则胎心率F＝60/T₁；

若T₁＝2T₂，则胎心率F＝60/T₂；

若T₁＝T₂/2，则胎心率F＝60/T₁；

其他情况下，胎心率F＝60/T₁。

针对传统多普勒胎心率检测仪数据处理能力弱，测量精度低，反应慢，稳定性差同时功耗大等问题；本发明利用智能移动终端音频接口，将由超声多普勒检测终端得到的胎心信号传给智能移动终端，并利用智能移动终端强大的数据处理能力，运行改进的胎心率检测算法，提高测量精度，减少反应时间，增加检测结果稳定性；同时精简了超声多普勒检测终端的硬件结构，降低了功耗和成本。

此外，针对胎心信号的成分复杂，内包含有较多的干扰成信号，采用传统的自相关算法计算胎心率时误差大，且容易出现倍频现象的问题；本发明对传统自相关算法进行改进，结合无偏估计的自相关算法来加强主峰的正确判断，克服了易出现倍频的缺点，提高了算法的准确度和稳定性。

附图说明

图1本发明胎心率检测方法的步骤流程示意图。

图2胎心信号A/D采样序列示意图。

图3胎心信号降采样序列示意图

图4胎心信号整流滑动平均滤波序列示意图

图5胎心信号有偏估计自相关序列示意图。

图6胎心信号无偏估计自相关序列示意图。

图7为本发明与现有技术胎心率模拟测试离均差平方和的对比示意图。

具体实施方式

为了更为具体地描述本发明，下面结合附图及具体实施方式对本发明的技术方案进行详细说明。

如图1所示，本发明基于智能移动终端的胎心率检测方法，包括以下步骤：

步骤1.对胎心信号进行预处理得到胎心多普勒包络信号。

1.1对通过音频接口输入的胎心信号，进行采样率f_s1为44.1KHz的A/D采样，得到采样值序列A(n)，该信号序列如图2所示；

1.2将得到的采样值序列A(n)再次进行采样率f_s2为4.41KHz的降采样处理，得到降采样值序列B(n)，该信号序列如图3所示；

1.3对降采样值序列B(n)进行判断，如果序列值为正，则不做处理，如果序列值为负，则对该序列值做绝对值处理。通过以上整流处理之后，得到整流值序列C(n)。

1.4对整流值序列C(n)进行滑动平均滤波处理，得到胎心多普勒包络信号。即：

$D\left(n\right)=\frac{1}{m}{\mathrm{\Σ}}_{i=0}^{m-1}C(n+i),$ 其中m＝10；

具体做法为：从C(n)序列第一个值开始，依次选取C(n)序列的前10个序列值，取平均之后作为序列D(n)的第一个序列数值；之后从C(n)序列第二个值开始，再依次选取C(n)序列的10个序列值，取平均之后作为序列D(n)的第二个序列数值，依次类推，直至取完序列C(n)的所有数值，得到胎心多普勒包络信号序列D(n)，该信号序列如图4所示。

步骤2.对胎心多普勒包络信号进行自相关分析，得到处理信号。

对胎心多普勒包络信号进行有偏估计自相关处理，得到有偏估计自相关序列，继而求取无偏估计自相关序列。自相关序列的求取步骤如下：

2.1对D(n)进行FFT运算，得到D(k)；

2.2计算R(k)＝D*(k)D(k)＝|D(k)|²；

2.3最后对R(k)作IFFT运算，得到D(n)的有偏自相关序列，如图5所示；

具体实现原理为：设序列x(n)的自相关函数的表达式为：

$r_x\left(i\right)={\mathrm{\Σ}}_{n=0}^{N-1}x\left(n\right)x(n+i)$

其中：x(n)是长度为N的序列。设X(k)代表x(n)的离散傅里叶变换。即：

x(n)＝IDFT{X(k)}

$r_x\left(i\right)={\mathrm{\Σ}}_{n=0}^{N-1}x\left(n\right)x(n+i)={\mathrm{\Σ}}_{n=0}^{N-1}\left[\frac{1}{N}{\mathrm{\Σ}}_{k=0}^{N-1}X\left(k\right)e^{j\frac{2\mathrm{\π}}{N}\mathrm{kn}}\right]\left[\frac{1}{N}{\mathrm{\Σ}}_{l=0}^{N-1}X\left(l\right)e^{j\frac{2\mathrm{\π}}{N}l(n+i)}\right]---\left(1\right)$

因为序列x(n)在实际中是实序列，即x(n)＝x^*(n)。

由于：

$r_x\left(i\right)={\mathrm{\Σ}}_{n=0}^{N-1}\left[\frac{1}{N}{\mathrm{\Σ}}_{k=0}^{N-1}{X}^{*}\left(k\right)e^{-j\frac{2\mathrm{\π}}{N}\mathrm{kn}}\right]\left[\frac{1}{N}{\mathrm{\Σ}}_{l=0}^{N-1}X\left(l\right)e^{j\frac{2\mathrm{\π}}{N}l(n+i)}\right]---\left(2\right)$

$=\frac{1}{N}{\mathrm{\Σ}}_{k=0}^{N-1}{\mathrm{\Σ}}_{l=0}^{N-1}{X}^{*}\left(k\right)X\left(l\right)e^{j\frac{2\mathrm{\π}}{N}\mathrm{li}}\left[\frac{1}{N}{\mathrm{\Σ}}_{k=0}^{N-1}{e}^{j\frac{2\mathrm{\π}}{N}n(l-k)}\right]---\left(3\right)$

$\frac{1}{N}{\mathrm{\Σ}}_{k=0}^{N-1}{e}^{j\frac{2\mathrm{\π}}{N}n(l-k)}=\frac{1}{N}\frac{1-e^{j\frac{2\mathrm{\π}}{N}N(l-k)}}{1-e^{j\frac{2\mathrm{\π}}{N}(l-k)}}=\left\{\begin{array}{cc}1& l=k\\ 0& l\≠k\end{array}\right.---\left(4\right)$

可得：

$r_x\left(i\right)={\mathrm{\Σ}}_{n=0}^{N-1}x\left(n\right)x(n+i)=\frac{1}{N}{\mathrm{\Σ}}_{n=0}^{N-1}{X}^{*}\left(k\right)X\left(k\right)e^{j\frac{2\mathrm{\π}}{N}\mathrm{ki}}=\frac{1}{N}{\mathrm{\Σ}}_{n=0}^{N-1}{\left|X\left(k\right)\right|}^2{e}^{j\frac{2\mathrm{\π}}{N}\mathrm{ki}}---\left(5\right)$

由式(5)可知，自相关函数的计算可以通过计算FFT和IFFT得到。

FFT运算，即快速傅里叶变换运算，为离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，是对离散傅里叶变换的算法进行改进获得的。IFFT运算，即快速傅里叶反变换运算。

$R_{\mathrm{XX}1}\left(m\right)=\frac{1}{N-m}*{\mathrm{\Σ}}_{n=0}^{N-m-1}{X}_n*X_{n+m}---\left(6\right)$

$R_{\mathrm{XX}2}\left(m\right)=\frac{1}{N}*{\mathrm{\Σ}}_{n=0}^{N-m-1}{X}_n*X_{n+m}---\left(7\right)$

$R_{\mathrm{XX}2}\left(m\right)=\frac{N-m}{N}*\frac{1}{N-m}*{\mathrm{\Σ}}_{n=0}^{N-m-1}{X}_n*X_{n+m}=\frac{N-m}{N}*R_{\mathrm{XX}1}\left(m\right)---\left(8\right)$

其中，m＝1,2…,N；式(6)称为有偏估计自相关，式(7)称为无偏估计自相关，由式(8)进行无偏估计自相关的求取，无偏自相关序列如图6所示。

步骤3.对处理信号进行对比判断处理，求取胎儿心率。

3.1对得到的处理信号(有偏自相关序列和无偏自相关序列)进行比较，从序列中找出最大值，通过寻找最大值点的位置，并计算出两个相邻最大值点的序列间隔长度n₁和n₂；

3.2通过上述测定峰值之间的序列间隔长度，分别得到两个自相关结果序列中相邻峰值之间的时间间隔T₁和T₂；

T₁＝n₁/f_s2，T₂＝n₂/f_s2，式中f_s2为胎心信号的降采样频率。

如果上述的T₁＝T₂，则胎儿心率F＝60/T₁；

如果上述的T₁＝2T₂，则胎儿心率F＝60/T₂；

如果上述的T₁＝T₂/2，则胎儿心率F＝60/T₁。

为了验证本发明胎心率测量的准确性和稳定性，参照中华人民共和国医药行业标准超声多普勒胎儿监护仪(YY0449-2009)中胎心率测量和显示范围试验，设计实验方案。模拟胎心信号由中国计量科学研究院的超声多普勒胎儿测量仪检测装置(型号FS-3)产生。

实验通过带动钢球在充满液体的容器中上下反复运动，来模拟胎心的跳动，其心跳次数可设置。为进一步评价本发明胎心率测量精确度和稳定性，特选取市场上主流厂商京柏公司的一款家用超声多普勒胎心仪作为对比实验。实验中，设置模拟胎心率信号的频率从240次/min到60次/min，间隔15次/min，频率的误差范围±0.1次/min，在各进行30次测量。最后计算各频率测量结果的标准偏差，其结果如表1所示：

表1

由表1可知，本发明测量的标准偏差比对照组的标准偏差小一个数量级，并且在实验过程中对照组在模拟信号为120次/min的时候出现严重倍频现象，测量结果为223±2次/min；而本发明在测量过程中，未出现倍频现象。

为了评价本发明的稳定性，计算上述测量各个频率的离均差平方和DEVSQ，公式如下：

$\mathrm{DEVSQ}=\mathrm{\Σ}{(x-\stackrel{\‾}{x})}^2$

将计算结果绘制成柱状图，如图7所示；从图7中，可以看出本发明测量各个模拟胎心频率时，其离均差平方和DEVSQ在90次/min和75次/min这两个频率时与对照组很接近，其他都比对照组小很多。这说明本发明的测量胎心率结果的离散程度比对照组小，稳定性更高。

上述的对实施例的描述是为便于本技术领域的普通技术人员能理解和应用本发明。熟悉本领域技术的人员显然可以容易地对上述实施例做出各种修改，并把在此说明的一般原理应用到其他实施例中而不必经过创造性的劳动。因此，本发明不限于上述实施例，本领域技术人员根据本发明的揭示，对于本发明做出的改进和修改都应该在本发明的保护范围之内。

Claims (2)

1.一种基于智能移动终端的胎心率检测方法，包括如下步骤：

(1)利用智能移动终端通过音频接口接收胎心仪提供的胎心信号；

(2)对所述的胎心信号进行预处理得到胎心多普勒包络信号序列，具体步骤如下：

2.1对所述的胎心信号进行A/D采样，得到胎心信号的采样值序列；

2.2对所述的采样值序列进行降采样处理，得到降采样值序列；

2.3对所述的降采样值序列进行绝对值处理，得到整流值序列；

2.4通过以下公式对所述的整流值序列进行滑动平均滤波处理，得到胎心多普勒包络信号序列；

$D\left(n\right)=\frac{1}{m}\underset{i=0}{\overset{m-1}{\mathrm{\Σ}}}C(n+i)$

其中：C(n+i)为整流值序列中的第n+i个数值，D(n)为胎心多普勒包络信号序列中的第n个信号值，n为大于0的自然数，m为滑动信号窗口的预设长度；

(3)对所述的胎心多普勒包络信号序列进行自相关分析，得到有偏自相关序列和无偏自相关序列，具体步骤如下：

3.1对所述的胎心多普勒包络信号序列进行快速傅里叶变换运算，得到复数序列；

3.2通过以下公式对所述的复数序列进行模平方运算处理，得到对应的模平方序列；

R(k)＝Ｄ*(k)D(k)＝|D(k)|²

其中：R(k)为模平方序列中的第k个数值，D(k)为复数序列中的第k个复数，D^*(k)为D(k)的共轭复数，|D(k)|为D(k)的模，k为自然数且1≤k≤N，N为胎心多普勒包络信号序列的信号值个数；

3.3对所述的模平方序列进行快速傅里叶反变换运算，得到有偏自相关序列；

3.4根据所述的有偏自相关序列通过以下公式计算出对应的无偏自相关序列；

$R_{XX2}\left(j\right)=\frac{N-j}{N}*R_{XX1}\left(j\right)$

其中：R_xx1(j)和R_xx2(j)分别为有偏自相关序列和无偏自相关序列中的第j个数值，N为胎心多普勒包络信号序列的信号值个数，j为自然数且1≤j≤N；

(4)对所述的有偏自相关序列和无偏自相关序列进行对比判断处理，以求取得到胎心率，具体步骤如下：

4.1从有偏自相关序列中确定出两个相邻峰值之间的序列间隔长度n₁；

4.2从无偏自相关序列中确定出两个相邻峰值之间的序列间隔长度n₂；

4.3根据所述的序列间隔长度n₁和n₂，计算出对应的时间间隔T₁和T₂，即使所述的序列间隔长度n₁和n₂分别除以降采样的频率，得到对应的时间间隔T₁和T₂；

4.4根据所述的时间间隔T₁和T₂，通过以下关系判断求取得到所述的胎心率；

若T₁＝T₂，则胎心率F＝60/T₁；

若T₁＝2T₂，则胎心率F＝60/T₂；

若T₁＝T₂/2，则胎心率F＝60/T₁；

其他情况下，胎心率F＝60/T₁。

2.根据权利要求1所述的胎心率检测方法，其特征在于：所述的步骤(2)至(4)均由智能移动终端内部处理器实现处理；所述的智能移动终端为智能手机或平板电脑。