CN101828916B - 心电信号处理系统 - Google Patents

Abstract

本发明涉及心电信号处理系统，所述系统基于FPGA平台，包括：预处理单元，用于对经采样的心电信号进行预处理；小波变换单元，用于对经预处理后的心电信号进行自适应提升方法的小波分解；特征提取单元，用于基于自适应提升方法定位心电信号的特征位置；信号处理控制单元，用于控制预处理单元、小波变换单元、特征提取单元的正常工作。本发明的特征提取方法准确率高、具有更强的抗干扰能力、运算量低。基于FPGA平台对心电信号进行处理，可有效提升实时信号处理能力，克服传统心电信号信号处理系统中因处理器性能限制而在心电信号信号分析处理方面的不足；还可满足便携式医疗器械在提升集成程度、降低功耗、系统小型化等的需求。

Description

心电信号处理系统

技术领域

本发明涉及心电信号处理技术，更具体地说，涉及一种心电信号处理系统。

背景技术

心电信号(electrocardiogram)是由体表电极随时间所记录的，反映心脏电活动的记录，是最早研究并应用于临床的生物电信号。ECG因为在心血管疾病的诊断、药物疗效评估等有重要临床应用，因而时至今日依然是医学界和生物医疗工程学界的研究热点。心电信号由4个基本波形所组成，分别是QRS波、P波、T波和正常情况下出现几率75％的U波，该特征划分由WillemEinthoven提出并沿用至今。

随着对ECG信号自动分析诊断的性能与功能覆盖等要求的不断提升，使得处理系统必须在对ECG信号的特征提取方面做到在保障检测准确率的同时提升信号处理的速度，并向小型化、集成化、便携式和低功耗等方向发展。

当前心电信号信号处理设备通常采用ARM、8051或MSP430等单片机做系统核心，主要由导联系统、ADC、系统核心处理器、电源模块、板载外设、LCD和存储器所构成，功能主要集中在对心电信号信号做长时间记录(通常在24小时以上)。

由于传统的便携式心电信号处理系统受到系统核心处理器性能的限制，在信号处理以及自动分析方面难以进行更深入的特征提取，如ST段、QT段检测等，并给出实时的警示信息。此外，尤其当需要进行十二导联心电信号信号的采集以及实时处理的场合，用传统的心电信号处理系统的架构就难以满足心电信号信号处理对系统性能的需求。当在系统中加入更多的外设或者处理器以解决系统处理能力受限制所带来的一系列问题时又会不可避免的受到功耗增加，系统体积增大以及重量增加等问题，给用户带来一系列的不便。

发明内容

本发明要解决的技术问题在于，针对现有技术的上述特征提取准确率低、运算量大、鲁棒性差、存储空间要求高，以及信号处理的速度低、功耗大的缺陷，提供一种心电信号处理系统。

本发明解决其技术问题所采用的技术方案是：构造一种心电信号处理系统，包括：

预处理单元，用于对经采样的心电信号进行预处理；小波变换单元，用于对经预处理后的心电信号进行自适应提升方法的小波分解；特征提取单元，用于基于自适应提升方法定位心电信号的特征位置；信号处理控制单元，用于控制预处理单元、小波变换单元、特征提取单元的正常工作。

在本发明所述的心电信号处理系统中，所述预处理单元包括工频噪声滤波器、肌电噪声滤波器和基线漂移滤波器。

在本发明所述的心电信号处理系统中，所述工频噪声滤波器包括梳状滤波器，所述梳状滤波器的转移函数为：

$H\left(z\right)=b\frac{1-z^{-N}}{1-\mathrm{\ρ}{z}^{-N}},$ $b=\frac{1+\mathrm{\ρ}}{2}$

其中，b是用于补偿滤波后信号的增益的补偿因子，ρ是影响零点间通带变化的快慢的系数，N是心电信号的采样频率与工频噪声频率的比值。

在本发明所述的心电信号处理系统中，所述肌电噪声滤波器包括多个级联的递归动求和(RRS)滤波器。

在本发明所述的心电信号处理系统中，所述肌电噪声滤波器包括4个级联的递归动求和滤波器，所述4个级联的递归动求和滤波器的转移函数为：

$H\left(z\right)={(\frac{1}{2}\×\frac{1-z^{-2}}{1-z^{-1}})}^4.$

在本发明所述的心电信号处理系统中，所述基线漂移滤波器包括巴特沃思低通滤波器，所述巴特沃思低通滤波器的转移函数为：

$H\left(z\right)=\frac{B\left(z\right)}{A\left(z\right)}=\frac{\underset{i=0}{\overset{M}{\mathrm{\Σ}}}{a}_i{z}^{-i}}{\underset{j=0}{\overset{N}{\mathrm{\Σ}}}{b}_j{z}^{-j}}$

在本发明所述的心电信号处理系统中，所述巴特沃思低通滤波器采样直接II型二阶区段(SOS)结构；所述二阶区段矩阵为：

$\mathrm{SOS}=\left[\begin{array}{cccccc}1& 2& 1& 1& -1.998& 0.998\\ 1& 2& 1& 1& -1.994& 0.994\\ 1& 2& 1& 1& -1.990& 0.990\\ 1& 2& 1& 1& -1.987& 0.987\\ 1& 2& 1& 1& -1.985& 0.985\\ 1& 2& 1& 1& -1.983& 0.983\\ 1& 2& 1& 1& -1.982& 0.982\\ 1& 1& 0& 1& -0.990& 0\end{array}\right]$

所述巴特沃思低通滤波器的每一阶的增益因子向量G为：

G＝[0.00002 0.00002 0.00002 0.00002 0.00002 0.00002 0.00454 1]

在本发明所述的心电信号处理系统中，小波变换单元包括：

分裂单元，用于对经采样的心电信号x(n)分裂为奇数序列x_o(n)和偶数序列x_e(n)，可由下式表示：x_o(n)＝x(2n+1)  x_e(n)＝x(2n)；

预测单元，用于基于所述偶数序列x_e(n)预测所述奇数序列x_o(n)，得到小波系数d(n)；可由下式表示：d(n)＝x_o(n)-P(x_e(n))，其中P为预测算子；

更新单元，用于基于所述小波系数d(n)更新所述偶数序列x_e(n)，得到尺度系数c(n)；可用下式表示：c(n)＝x_e(n)+U(x_o(n))，其中U为更新算子。

在本发明所述的心电信号处理系统中，所述预测算子P由如下公式表示：

$P\left(n\right)=\frac{1}{2}[c(n-1)+c\left(n\right)]$

所述更新算子U由如下公式表示：

$U\left(n\right)=\left\{\begin{array}{cc}d(n-1)& {\mathrm{\&Delta;}}_L<{\mathrm{\&Delta;}}_R\\ \frac{1}{2}[d(n-1)+d\left(n\right)]& {\mathrm{\&Delta;}}_L={\mathrm{\&Delta;}}_R\\ d\left(n\right)& {\mathrm{\&Delta;}}_L>-{\mathrm{\&Delta;}}_R\end{array}\right.$

Δ_L＝|c(n)-d(n-1)|  Δ_R＝|c(n)-d(n)|

在本发明所述的心电信号处理系统中，所述特征提取单元用于定位R波的位置，所述特征提取单元包括：阈值计算单元，用于依据小波变换单元进行小波分解直至达到预设尺度值后输出的小波系数d(n)计算阈值；检测单元，用于依据所述阈值检测所述心电信号中R波的位置；判断单元，用于依据尺度系数c(n)判断所述R波的位置是否正确，若否，则以预设步长调整所述阈值并由检测单元重新检测R波的位置。

在本发明所述的心电信号处理系统中，所述阈值计算单元进一步用于：将小波系数d(n)划分为多个等长区间；计算每一等长区间的极大值；计算所述多个等长区间的极大值的均值；选取所述均值的一半作为阈值。

在本发明所述的心电信号处理系统中，所述预设尺度值为4，所述等长区间为fs/4，其中fs为心电信号x(n)的采样频率；

所述检测单元进一步用于：依据所述阈值检测所述心电信号，并依据检测结果定位R波的位置并判断R波间距是否大于1/150s，若是，则由判断单元进行进一步的处理；

所述判断单元进一步用于：依据尺度系数c(n)判断R波的位置是否满足c₁(n)＜c₂(n)，c₁(n)＜c₃(n)，c₁(n)＜c₄(n)，若否，则以预设步长调整所述阈值并由所述检测单元重新检测R波的位置；其中c₁(n)为第1次小波变换的尺度系数，c₂(n)为第2次小波变换的尺度系数，c₃(n)为第3次小波变换的尺度系数，c₄(n)为第4次小波变换的尺度系数。

在本发明所述的心电信号处理系统中，所述特征提取单元还用于定位QRS波的起始点位置、终止点位置。

在本发明所述的心电信号处理系统中，所述系统还包括总线控制单元，用于实现所述系统中各个单元之间的互联。

在本发明所述的心电信号处理系统中，所述系统还包括外设单元，所述外设单元包括I²C控制器、SPI控制器、GPIO控制器和SD卡主控制器。

本发明的有益效果是，基于FPGA平台对心电信号进行处理，可实现实现心电信号信号处理系统的小型化、集成化。本发明的心电信号处理系统一方面可以有效提升便携式系统的实时信号处理能力，实现对12导联心电信号数据的并行处理，克服传统心电信号信号处理系统中因为处理器性能限制而在心电信号信号自动分析处理方面的不足；另一方面，可以满足便携式医疗器械在提升集成程度、降低功耗、系统小型化等方面的需求，实现了为使用者提供更多心脏健康情况的相关信息。

另外在本发明的心电信号处理系统，心电信号特征提取是基于自适应提升的方法来实现的，准确率jj较高、运算量较低、抗干扰能力强。

这些优点使得能够在对实时性要求高的应用场合以及嵌入式系统平台上进行准确率较高的特征提取，并有效的节约了系统资源，同时对心电信号处理系统进行功能扩展也带来了便利。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1是依据本发明一实施例的心电信号处理系统进行心电信号特征提取的方法流程示意图；

图2是基于自适应提升方法的第二代小波变换架构示意图；

图3是使用本发明一实施例的心电信号处理系统进行基于自适应提升方法的R波检测流程示意图；

图4是基于自适应提升方法和基于Marr小波的R波检测算法运行时间和检测准确率对比示意图；

图5是使用本发明一实施例的心电信号处理系统基于自适应提升方法定位的R波位置对QRS波的起始点进行定位的示意图；

图6是依据本发明一实施例的心电信号处理系统结构示意图；

图7是依据本发明另一实施例的心电信号处理系统架构示意图；

图8是图6所示的预处理单元结构示意图；

图9为改进前后梳状滤波器处理效果的对比图；

图10是小波阈值方法与RRS滤波器去噪效果对比图；

图11是中值滤波方法、EMD方法和Butterworth滤波方法分别对基线漂移进行处理的结果对比图；

图12是总线控制单元705的接口的结构示意图；

图13是I²C控制器708的结构示意图；

图14是SPI控制器706的结构示意图；

图15是GPIO控制器709的结构示意图；

图16是SD卡主控制器707的结构示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

图1是依据本发明一实施例的心电信号处理系统进行心电信号特征提取的方法流程示意图。图1所示的方法100用于心电信号R波的特征提取。

在步骤101中，基于自适应提升方法对经采样的心电信号x(n)进行多次小波变换直至达到预设尺度值；

在步骤102中，依据小波系数d(n)计算阈值；

在步骤103中，依据所述阈值检测所述心电信号中R波的位置；

在步骤104中，依据尺度系数c(n)判断所述R波的位置是否正确，若否，则进入步骤105以预设步长调整阈值并返回步骤104重新检测R波的位置。若是，则说明已正确定位R波，因此结束R波的检测。

通过方法100可以准确地定位R波的位置。

图2是基于自适应提升方法(Adaptive Lifting Scheme，ALS)的第二代小波变换架构示意图。基于提升方法的第二代小波变换包括3个步骤：分裂、预测和更新。具体方法如下：

1)分裂：将信号x(n)分裂为不相交的奇数索引序列x_o(n)和偶数索引序列x_e(n)，如式(1)所示。x_o(n)＝x(2n+1)  x_e(n)＝x(2n)

2)预测：通过预测算子P基于序列x_e(n)预测序列x_o(n)，预测值误差即小波系数d(n)。d(n)＝x_o(n)-P(x_e(n))

3)更新：通过更新算子U基于预测误差d(n)更新序列x_e(n)，所得结果即尺度系数c(n)。c(n)＝x_e(n)+U(x_o(n))

其中，在本发明一实施例中，使用的更新算子U如下式表示：

$U\left(n\right)=\left\{\begin{array}{cc}d(n-1)& {\mathrm{\&Delta;}}_L<{\mathrm{\&Delta;}}_R\\ \frac{1}{2}[d(n-1)+d\left(n\right)]& {\mathrm{\&Delta;}}_L={\mathrm{\&Delta;}}_R\\ d\left(n\right)& {\mathrm{\&Delta;}}_L>-{\mathrm{\&Delta;}}_R\end{array}\right.$

Δ_L＝|c(n)-d(n-1)|  Δ_R＝|c(n)-d(n)|

使用的预测算子P如下式表示：

$P\left(n\right)=\frac{1}{2}[c(n-1)+c\left(n\right)]$

上述使用的更新算子U和预测算子P是经实验证实，获取的检测性能较好，当然也不排除其他可能的实现方法，本发明并不限于此，只要是基于本发明的思想来实现特征的提取都属于本发明的保护范围。

图3是使用本发明一实施例的心电信号处理系统进行基于自适应提升方法的R波检测流程示意图。

首先，在步骤301中，对输入的经采样的心电信号(采样频率表示为fs)进行分裂，分裂为不相交的奇数索引序列x_o(n)和偶数索引序列x_e(n)，然后经过步骤302-304对分裂后的信号进行预测和更新(详见图2关于小波变化的描述)，直至达到预设的分解尺度值。在该实施例中，设置该分解尺度值为4，当然也可以设置为其他数值，可根据要求的精确程度以及运算速度来进行权衡，本发明并不限制于此。如果在步骤304中，判断分解已经完成，则进入下一步骤305。否则，计算进行小波变换，直至达到分解尺度4。

在步骤305中，计算阈值。将小波系数d(n)划分为多个等长区间，并选取每一等长区间的极大值，求出所选出的一系列极大值的均值，选取该均值的一半(50％)作为阈值。在本实施例中，等长区间的长度为fs/4，本发明并不限于此，只要是能够保证本方法正常运作的数值都是可以的。

在步骤306中，依据阈值进行R波检测。通常情况标准心搏数每秒低于150下，在基于标准心搏数的前提下检测R波间距是否大于1/150s，若是，则进入步骤307，若否则说明心跳很不正常。

在步骤307中，在尺度系数c(n)上基于李氏指数(李雅普洛夫指数，Lyapunov指数)判断R波位置是否满足c₁(n)＜c₂(n)，c₁(n)＜c₃(n)，c₁(n)＜c₄(n)，其中c₁(n)为第1次小波变换的尺度系数，c₂(n)为第2次小波变换的尺度系数，c₃(n)为第3次小波变换的尺度系数，c₄(n)为第4次小波变换的尺度系数。若满足，则说明已正确检测到R波的位置。若不满足，说明存在漏检、误检(错检)，这时进入步骤308，调整阈值，这里以初始阈值的10％为步长调整阈值大小并进入步骤306在指定区间重新检测R波的位置。这里的步长可依据实际情况进行选定，此处仅作为示例。所述的指定区间可以是检测出的R波间距内的区间或者包含R波间距的区间，减少了R波的检测时间。当然，也可重新依据新的阈值对整个小波系数进行检测。

在对本发明的自适应提升方法的准确率进行验证时采用了麻省理工(Massachusetts Institute of Technology)和Beth Israel医院联合建立的MIT-BIH心律失常数据库。MIT-BIH数据库中数据样本丰富、类型特征突出并有详细注释，在ECG信号自动分析处理的研究领域均使用该数据库中的样本进行实验验证。MIT-BIH心律失常数据库中共48组样本，每组样本长度30分钟，采样频率360Hz，总标准心搏数103891，并详细标明了ECG信号各个特征起始点终止点等特征信息。采用MIT-BIH数据并基于ALS方法对R波进行检测的结果统计如表1所示。其中，错检数＝漏检数+误检数。

表1

如表1中所示，由于部分样本，如105、108和228号样本受干扰较为严重，共出现错漏检99个，占到总错检数将近33％，总体上漏检数123、误检数164、总错检数287，检验的准确率达到99.724％，可以满足正常使用的需要。为了对比与其他方法性能上的优劣，本发明中实现了基于第一代小波变换的Marr小波的R波检测方法。Marr小波又被称为墨西哥草帽小波，是高斯函数的二阶导数，主要优点在于其小波变换尺度系数上的模极大值点与原始序列上模极大值点有一一对应的关系。Marr小波的母函数表达式如下式所示：

在Marr小波变换尺度4的尺度系数上，分等长区间分别求局部极大值，再对这组极大值求均值，将该均值二分之一作为阈值，求出过阈值的连续区间中极大值为R波的相应位置，再修正时移，此时与尺度4相应时延为20点，即与原始信号中R波位置有20点的延时，所采取的检错策略与基于ALS的方法相同。

基于第一代小波和第二代小波变换检测方法结果对比如表2和图4所示。图4示出了基于ALS和基于Marr小波的R波检测算法运行时间和检测准确率对比。虽然基于ALS的R波检测方法对比基于Marr小波的方法漏检多出52个，误检多出55个，总计107个，错检率高出0.106％，但检测准确率依然达到了99.724％，可以满足在嵌入式系统上对ECG信号进行实时R波检测的需要。对相同长度向量进行检测所花费的时间仅为基于Marr小波方法的11.49％，且不必进行定点化处理，大大降低了实现的难度。

表2

在定位R波后，为了便于进一步定位其他特征如P波和T波等，必须对QRS波的起始点以及终止点进行定位，具体方法可结合差分阈值的方法来实现，图5示出了使用本发明一实施例的心电信号处理系统基于自适应提升方法定位的R波位置对QRS波的起始点进行定位的示意图500，步骤描述如下：

在步骤501中，计算预设尺度值为4的尺度系数c(n)的差分信号diff(n)，可由下式表示：diff(n)＝(-2)×c(n-2)-c(n-1)+c(n+1)+2×c(n+2)；

在步骤502中，计算差分阈值diff_thresh，可由下式表示：

$\mathrm{diff}\_\mathrm{thresh}=\frac{\mathrm{thresh}\_\mathrm{param}}{16}\&times;\max ;$

其中，thresh_param为可选参数，取值为2、4、8或16，max的初始值为划分为多个等长区间的小波系数的模极大值；

在步骤503中，依据所述差分阈值diff_thresh从R波的位置向前检测QRS波的起始点位置。

在本发明另一实施例中，还可对QRS波的终止点进行定位。若是对QRS波的终止点进行定位，就要依据所述差分阈值diff_thresh从R波的位置向后检测QRS波的终止点位置。

随后在步骤504中判断是否连续检测到差分信号diff(n)大于所述差分阈值diff_thresh，若是，则可定位心电信号中QRS波的起始点；否则，进入步骤505。

若是对QRS波的终止点进行定位。步骤504是判断是否连续检测到差分信号diff(n)大于所述差分阈值diff_thresh，若是，则可定位心电信号中QRS波的终止点；否则，进入步骤505；

在步骤505中，由下式更新所述差分阈值diff_thresh，并进入步骤503重新定位QRS波的起始点；

$\max =\frac{\mathrm{first}\_\max -\max }{\mathrm{filter}\_\mathrm{param}}+\max$

其中，first_max为对应R波的幅度与对应的QRS波的起始点的幅度的差值，filter_param为可选参数，取值为2、4、8或16。

图6是依据本发明一实施例的心电信号处理系统600结构示意图。系统600是基于FPGA平台来实现的。参考图6，可知，系统600包括预处理单元601、小波变换单元602、特征提取单元603、信号处理控制单元604。

预处理单元601，用于对经采样的心电信号进行预处理；

小波变换单元602，用于对经预处理后的心电信号进行自适应提升方法的小波分解；

特征提取单元603，用于基于自适应提升方法定位心电信号的特征位置；

信号处理控制单元604，用于控制预处理单元601、小波变换单元602、特征提取单元603的正常工作。

小波变换单元602包括：分裂单元，用于对经采样的心电信号x(n)分裂为奇数序列x_o(n)和偶数序列x_e(n)；预测单元，用于基于所述偶数序列x_e(n)预测所述奇数序列x_o(n)，得到小波系数d(n)，可由下式表示：d(n)＝x_o(n)-P(x_e(n))，其中P为预测算子；更新单元，用于基于所述小波系数d(n)更新所述偶数序列x_e(n)，得到尺度系数c(n)，可用下式表示：c(n)＝x_e(n)+U(x_o(n))，其中U为更新算子。

所述预测算子P由如下公式表示：

$P\left(n\right)=\frac{1}{2}[c(n-1)+c\left(n\right)]$

所述更新算子U由如下公式表示：

$U\left(n\right)=\left\{\begin{array}{cc}d(n-1)& {\mathrm{\&Delta;}}_L<{\mathrm{\&Delta;}}_R\\ \frac{1}{2}[d(n-1)+d\left(n\right)]& {\mathrm{\&Delta;}}_L={\mathrm{\&Delta;}}_R\\ d\left(n\right)& {\mathrm{\&Delta;}}_L>-{\mathrm{\&Delta;}}_R\end{array}\right..$

Δ_L＝|c(n)-d(n-1)|  Δ_R＝|c(n)-d(n)|

特征提取单元603用于定位R波的位置，特征提取单元603包括：阈值计算单元，用于依据小波变换单元进行小波分解直至达到预设尺度值后输出的小波系数d(n)计算阈值；检测单元，用于依据所述阈值检测所述心电信号中R波的位置；判断单元，用于依据尺度系数c(n)判断所述R波的位置是否正确，若否，则以预设步长调整所述阈值并由检测单元重新检测R波的位置。

本发明的心电信号处理系统中的小波变换单元602、特征提取单元603的具体工作步骤和方法参见图1-5及其相关描述。当然本发明中的心电信号处理系统采用的小波变换和特征提取还可以采用其他的方法，也会得到较好的效果，达到本发明的目的，因此图1-图5的相关描述并不作为对本发明的限制。

图7是依据本发明另一实施例的心电信号处理系统架构700示意图。参考图7，可知，系统架构700包括信号处理控制单元701、预处理单元702、小波变换单元703、特征提取单元704。除了这些核心部件外，还包括总线控制单元705、SPI控制器706、SD卡主控制器707、I²C控制器708、GPIO控制器709，以及存储器710。

系统架构700的核心部件中，信号处理控制单元701是用于调用其他功能模块的顶层模块；预处理单元702用于对ECG信号进行预处理；小波变换单元703用于进行ALS分解；特征提取单元704用于定位ECG信号特征如QRS波、P波和T波等。关于ALS分解和特征提取的具体实现方法见上文所述(图1-图5以及相关描述)。

存储器710为内存模块，例如但不限于Block RAM。

总线控制单元705，用于各个模块间的互联和通信，遵循OpenCore组织的Wishbone总线协议。

外设模块，在图7中示出了4个外设模块：SPI控制器706、SD卡主控制器707、I²C控制器708、GPIO控制器709。SPI控制器706用于板上各个IC间的通信，与A/D转换器连接；SD卡主控制器707用于将数据保存至SD卡；I²C控制器708和GPIO控制器709用于通信和扩展系统功能。

电源模块，用于为各个外设以及系统核心供电。

系统工作的基本流程为采集到ECG信号后将数据存入指定的内存缓冲区内，当数据缓冲完成后对ECG信号首先进行预处理和ALS分解，然后进行特征提取，最后将处理完毕的数据存入SD卡中。

系统700的核心采用了Xilinx的XC2VP30，为Xilinx公司强调高性能的Virtex系列FPGA之一，该FPGA采用了130nm工艺，最高工作频率可达420MHz并集成了专用的数字信号处理处理模块DSP Slice，该FPGA的逻辑资源数量如表4.1所示^[33]。系统验证基于Xilinx所提供的XUP Virtex II Pro开发板。

系统700程序的编写采用了硬件描述语言Verilog HDL。Verilog DHL目前是IEEE标准之一，通称Verilog 2001。软件的编译工具采用了Xilinx的ISE10.1，在线调试使用了Chip Scope Pro 10.1，其中功能模块SPI控制器、SD卡主控制器参考了开源项目OpenCores的设计。

信号处理模块主要由预处理单元702、小波变换单元703、特征提取单元704三部分所组成。小波变换单元703、特征提取单元704采用上文所述的本发明的心电信号特征提取方法来实现。由于FPGA中的并发执行特性，采用了有限状态机的形式来进行编码实现，并做了定点化的处理。乘法部分则借助了XC2VP30中集成专用数字信号处理模块DSP Slice来实现，V2系列为Xilinx的FPGA中首先集成DSP Slice的产品。DSP Slice在XC2VP30为18*18bit输入、32bit输出的嵌入式乘法器，后续Xilinx的FPGA中加入了前置或后置乘加的结构，并扩展了乘法器位数进一步强化了其功能。DSP Slice在FPGA中的位置Block RAM列相邻，在使用流水线技术时结合DSP Slice的应用可以使用户应用将FPGA运行在最大时钟频率上称为可能。

下面对预处理单元702、小波变换单元703、特征提取单元704进行阐述。

对ECG信号进行预处理的目的是为了便于进行特征提取以及进一步的分析。对ECG信号进行预处理的效果对后续的特征提取效果具有至关重要的作用。ECG信号在采集过程中受到的干扰噪声主要包括工频噪声，由人体分布电容引起，频率由当地电网决定；肌电噪声，由人体肌肉活动所导致的体表电势差所引起，幅度由体表电势以及肌肉活动的状况决定；基线漂移噪声，由人体的呼吸活动引起，频率在1Hz以下，对ECG信号造成的从基线偏离的幅度通常在ECG信号的R波峰值的15％到20％之间。

图8是图6所示的预处理单元结构示意图。预处理单元601包括工频噪声滤波器801、肌电噪声滤波器802和基线漂移滤波器803。

关于工频噪声滤波器801。描述如下：

为了保证滤波效果并兼顾处理速度和便于在FPGA平台上的实现，本发明采用了梳状滤波器来处理工频噪声。梳状滤波器的原理是利用信号和自身的延时叠加，从而产生相位抵消。在ECG信号的实时处理中，为满足实时处理在速度方面的需要以及在特征提取过程中对滤波器性能要求并非很高的情况下，梳状滤波器的特点很适合用于处理工频噪声。

梳状滤波器的转移函数如下式所示：

$H\left(z\right)=\frac{1}{N}\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{z}^{-n}=\frac{1}{N}\frac{1-z^{-N}}{1-z^{-1}}$

由转移函数可知，滤波器的零点为n＝0，1，...N-1，N零点均匀分布在z平面的单位圆|z|＝1上，即单位圆被这N个零点所等分。z＝1处的极点与零点相抵消使滤波器成为FIR滤波器，具有线性相位性质。因此为了在处理工频噪声中使用梳状滤波器，就必须使采样频率f_s满足为工频噪声的整数倍N的条件。假设ECG信号的采样频率f_s＝360Hz，工频噪声频率60Hz则必须取N＝6从而使陷波的频率分别为60Hz的整数倍，因而也对工频噪声的各次谐波有抑制的作用。

上述的滤波器虽然可以有效滤除工频噪声以及其各次谐波带来的干扰，但其通带并不平缓，会导致信号失真，导致QRS波群的明显畸变，对后续的特征提取处理会带来不利的影响。为了改善滤波器的性能，本发明引入了零极点补偿以便使过零点之间的通带更为平坦。引入零极点补偿后的梳状滤波器的转移函数如下式所示。

$H\left(z\right)=b\frac{1-z^{-N}}{1-\mathrm{\&rho;}{z}^{-N}},$ $b=\frac{1+\mathrm{\&rho;}}{2}$

其中参数ρ会影响零点间通带变化的快慢，b用于补偿滤波后信号的增益，在此实施例中N＝6，ρ＝0.9，本发明并不限于此，具体取值取决于采样频率和工频噪声频率。

图9为改进前后梳状滤波器处理效果的对比图。由图9可见，由于ECG信号90％以上的能量都集中在40Hz以下的低频频段，因此确保0到60Hz之间通频带的平稳可以有效的改善滤波后导致的失真并有效的滤除工频噪声的干扰。虽然在进入零极点补偿后必须牺牲原有梳状滤波器线性相位的性质和整系数的特点，但系数ρ＝0.9，补偿因子b＝0.95，在进行定点化后运算过程中数据的动态范围有限，在后续的DSP系统以及FPGA系统中都可以较好地平衡滤波器在运算速度和性能上的需求，并未导致系统资源的浪费。

关于肌电噪声滤波器802，描述如下：

肌电噪声频率与ECG有用信号频带重叠，因此如何在确保ECG信号在滤波后失真尽可能小一直是肌电噪声滤波方法的难点。为了便于实现，提升滤波器性能并降低系统运行时的开销，采用了无乘法器结构的递归动求和(Recursive Running Sum，RRS)滤波器来处理肌电噪声。

在采样率f_s＝360Hz情况下低通滤波器截止频率选择在90Hz，为了降低实现的复杂度并节约在FPGA上实现时对延时单元的需求，取K＝2。当K＝2时滤波器的阻带衰减仅12dB，因此级联4个相同滤波器以得到足够大的阻带衰减以抑制肌电噪声，级联后的阻带衰减为级联的滤波器阻带衰减的代数和，即48dB。级联得到的滤波器转移函数如下式所示。

$H\left(z\right)={(\frac{1}{2}\&times;\frac{1-z^{-2}}{1-z^{-1}})}^4$

由于RRS滤波器具有线性相位的特点，所以在滤波后仍然可以保持ECG信号各个特征不受干扰，并且可以有效地抑制肌电噪声带来的干扰。

为了对比所使用的RRS滤波器对肌电噪声干扰的处理效果，本发明将目前研究比较深入的基于小波阈值的去噪方法与RRS滤波器的去噪方法进行对比。图10是小波阈值方法与RRS滤波器去噪效果对比图。

如图10可见，由于肌电噪声频率范围与ECG信号有重叠，基于小波的方法对肌电噪声的处理效果更好，滤波后的信号因为Coif4小波具有正则性而显得更为平滑，但对与频率较高的QRS波群则导致了R波幅度的降低；RRS滤波后因为本身具有线性相位的特点可以比较好的保持ECG信号各项特征，QRS波群失真相对较小。此外，Coif4小波分解的高低通系数向量长度均为24，进行滤波需分解到尺度5，对长度2000的ECG信号向量进行滤波所耗费的CPU时间为RRS滤波器的35倍，而系统时间为87倍。除了在处理时间上相差较大以外，在定点处理器DSP上必须对小波分解系数进行定点化而产生数据动态范围在16位以上而在计算出32位的临时结果后必须右移位，从而导致精度降低，RRS滤波器则因为整系数而不存在精度问题。

关于基线漂移滤波器803，描述如下：

ECG信号中基线漂移主要由人体呼吸活动所导致，频率在0.5到1Hz之间，对后续的特征提取会产生不利的影响。对于基线漂移的处理目前有中值滤波、小波和FIR滤波器进行高通滤波等。在嵌入式系统上实现中上述方法都不约而同地存在运算量大的问题，为了降低运算量以便于在嵌入式系统上对ECG信号进行实时的分析处理，本发明使用巴特沃思滤波器(Butterworth)低通滤波器来处理基线漂移，并与基于经验模态分解(Empirical ModeDecomposition，EMD)的方法和中值滤波的方法进行对比。

Butterworth滤波器属于IIR滤波器，最大特点在于其幅频响应中通带部分对比其他滤波器最大限度平稳，但比其他IIR滤波器在实现相同阻带性能时所需的滤波器阶数更高。

滤波的具体步骤如下：1)用低通滤波器滤除ECG有用信号；2)修正低通滤波器输出延时，其中延时取滤波器相位延时在通频带内均值1150；3)用修正延时后的滤波器输出减去原始信号得到滤波基线漂移后的结果。

所设计的Butterworth低通滤波器的截止频率1Hz，通频带0.5Hz，通带最大增益1dB，阻带最小衰减80dB，阶数15，结构上采用了直接II型(direct formII)的2阶区段(SOS)(Second Order Section)级联形式实现。SOS(Second OrderSection)矩阵和增益因子向量G分别如下面式子所示。

$\mathrm{SOS}=\left[\begin{array}{cccccc}1& 2& 1& 1& -1.998& 0.998\\ 1& 2& 1& 1& -1.994& 0.994\\ 1& 2& 1& 1& -1.990& 0.990\\ 1& 2& 1& 1& -1.987& 0.987\\ 1& 2& 1& 1& -1.985& 0.985\\ 1& 2& 1& 1& -1.983& 0.983\\ 1& 2& 1& 1& -1.982& 0.982\\ 1& 1& 0& 1& -0.990& 0\end{array}\right]$

G＝[0.00002 0.00002 0.00002 0.00002 0.00002 0.00002 0.00454 1]

其中，SOS矩阵中元素模取值范围在0.98到2之间，对实型数值取小数点后3位，而每级滤波器输出增益因子仅在0.00454时需要使用乘法器，其余均可以通过移位实现，对于16位ECG数据定点化处理过程中数据动态范围有限，算法实现简单。

下面针对本发明的基线漂移滤波方法与基于EMD和中值滤波的方法进行对比。三种方法分别对基线漂移进行处理的结果如图11所示，在相同机器上的运行时间统计如表3所示，L表示用于测试的ECG信号向量长。

图11从上到下顺序依次原始信号、中值滤波方法结果、EMD方法结果和Butterworth滤波结果，可见三种方法都可以在滤波后维持ECG信号中P波和T波等低频特征，但对于嵌入式系统平台而言基于EMD和中值滤波的方法均难以满足实时处理的需要，尤其是EMD分解不论是其高算法复杂度或者较大的运算量都使得该方法难以在目前的中低端嵌入式处理器上实现。

表3

上文关于预处理单元的工频噪声滤波器801、肌电噪声滤波器802和基线漂移滤波器803的相关描述仅作为优选实施例，还可以选择其他可以实现预处理功能的器件，本发明并不限于此。

关于图7中所示的总线控制单元705，描述如下：

总线控制单元705主要用于各个模块间的互联和通信，遵循OpenCore组织的Wishbone总线协议。Wishbone总线协议是目前各个免费硬件IP通用的总线接口规范。总线控制单元705的接口的具体实现见图12所示。所对应的信号如下所示：

●RST_I：复位信号，使所有的有限状态机返回初始状态，由顶层模块供给；

●CLK_I：时钟信号，上升沿有效，由顶层模块供给；

●ADDR_I/O：地址信号，本发明中宽度为16bit；

●DAT_I/O：数据信号，本发明中宽度16bit；

●WE_I/O：写使能信号，用于标识当前总线周期为写操作或者读操作，信号置位时表示执行写操作，相反时执行读操作；

●SEL_I/O：片选信号，用于选择进行数据操作的从设备，本文中宽度3bit；

●STB_I/O：表示设备是否处于数据操作状态；

●ACK_I/O：确认信号，当置位时表示数据操作完成；

●CYC_I/O：表示设备是否处于处于总线周期，当进行数据通信的主从设备处于工作状态时被置位；

●TAGN_I/O：用户自定义信号，如数据信号的纠错信息等。

关于图7中示出的外设模块：SPI控制器706、SD卡主控制器707、I²C控制器708、GPIO控制器709，描述如下：

I²C(Inter IC Connection Bus)是双向的半双工串行总线，主要用于板上设备间的短距离且非频繁的数据通信，具有支持多主设备、具备冲突检测和仲裁机制等特点。I2C设备的接口双线分别是串行数据线(Serial Data Line，SDA)和串行时钟线(Serial Clock Line，SCL)，最高传输速率不超过4Mbps，本发明中的I2C端口工作频率为400Kbps。

I²C控制器708的结构如图13所示。包括Wishbone总线控制器、预分频寄存器、指令寄存器、状态寄存器、传输寄存器、接收寄存器、字节指令控制器、数据I/O移位寄存器、时钟生成器以及位指令控制器。

SPI(Serial Peripheral Bus)是一种四线、同步、全双工的串行数据总线，传输速率从1MHz到70MHz。SPI设备的接口四线分别为：

●SCLK，串行时钟，由主设备供给；

●MOSI/SIMO，主设备输出/从设备输入；

●MISO/SIMO，主设备输入/从设备输出；

●SS，从设备选择信号，低电平有效，主设备供给。

SPI控制器706的结构图如图14所示。包括Wishbone总线控制器、预分频寄存器、指令寄存器、状态寄存器、传输寄存器、接收寄存器、字节指令控制器、数据I/O移位寄存器、时钟生成器以及指令控制器。

GPIO(General Purpose Input Output)是嵌入式处理器用于与外部设备连接的通用接口，输入输出方向可由用户自定义，通常用于处理器功能扩展，如用于读取输入或者输出，控制其他外部设备等，是嵌入式处理器上最为常见的外设接口。本发明中GPIO接口为8线，输入电平最高不超过5V，此处仅为示例，不作为对本发明的限制。

GPIO控制器709的架构如图15所示。包括Wishbone总线控制器、使能寄存器、控制寄存器、方向选择寄存器、传输寄存器、接收寄存器、I/O接口。

SD(Secure Digital)卡在系统中被用于保存采集到的ECG数据。SD卡的结构能保证数据的安全，且易于格式化。本发明的SD卡控制器有如下特征：

●具有缓冲描述器(Buffer Descriptor)；兼容SD主控制器规范(2.0版本)；

●支持4bit模式；可变更的读/写FIFO大小；

●具有内置CRC校验，数据线为CRC16，指令线为CRC7；

SD卡控制器由8个模块所构成，架构如图16所示。包括Wishbone总线控制器、SD控制器顶层模块、缓冲描述器、SD指令主控制器、SD数据主控制器、FIFO、SD指令串行接口、SD数据串行接口。

本发明还可基于FPGA的系统设计最终实现专用的ECG信号处理ASIC。为实现该功能，就必须对使FPGA的综合输出结构更紧凑、占用逻辑资源更少且满足信号时序要求。在对本发明的系统进行了一系列的仿真以及在线调试后，功能以及时序方面均达到了设计要求，而所占用的FPGA系统资源如表4所示。

表4

FPGA是基于可编程可配置逻辑模块(Controllable Logical Block，CLB)阵列。FPGA所具备的高性能与高灵活性等特点可以很好的满足ECG信号处理关于性能方面的需求，此外FPGA可以对数据进行并行处理的特点可以在十二导联ECG数据处理方法得到完全的发挥，这一优势是其他嵌入式处理器所难以比拟的。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种心电信号处理系统，其特征在于，所述系统基于FPGA平台，包括：

预处理单元，用于对经采样的心电信号进行预处理；

小波变换单元，用于对经预处理后的心电信号进行自适应提升方法的小波分解；

特征提取单元，用于基于自适应提升方法定位心电信号的特征位置；

信号处理控制单元，用于控制预处理单元、小波变换单元、特征提取单元的正常工作；

进一步的，所述预处理单元包括工频噪声滤波器、肌电噪声滤波器和基线漂移滤波器；

进一步的，所述工频噪声滤波器包括梳状滤波器，所述梳状滤波器的转移函数为：

$H\left(z\right)=b\frac{1-z^{-N}}{1-\mathrm{\&rho;}{z}^{-N}},b=\frac{1+\mathrm{\&rho;}}{2}$

其中，b是用于补偿滤波后信号的增益的补偿因子，ρ是影响零点间通带变化的快慢的系数，N是心电信号的采样频率与工频噪声频率的比值；

所述肌电噪声滤波器包括多个级联的递归动求和(RRS)滤波器；

所述基线漂移滤波器包括巴特沃思低通滤波器，所述巴特沃思低通滤波器的转移函数为：

$H\left(z\right)=\frac{B\left(z\right)}{A\left(z\right)}=\frac{\underset{i=0}{\overset{M}{\mathrm{\&Sigma;}}}{a}_i{z}^{-i}}{\underset{j=0}{\overset{N}{\mathrm{\&Sigma;}}}{b}_j{z}^{-j}}$

其中M表示巴特沃思低通滤波器的阶数；

所述巴特沃斯低通滤波器滤波的具体步骤如下：1)用所述巴特沃斯低通滤波器滤除ECG有用信号；2)修正所述巴特沃斯低通滤波器输出延时，其中延时取巴特沃斯低通滤波器相位延时在通频带内均值1150；3)用修正延时后的滤波器输出减去原始信号得到滤波基线漂移后的结果。

2.根据权利要求1所述的系统，其特征在于，所述肌电噪声滤波器包括4个级联的递归动求和滤波器，所述4个级联的递归动求和滤波器的转移函数为：

$H\left(z\right)={(\frac{1}{2}\&times;\frac{1-z^{-2}}{1-z^{-1}})}^4.$

3.根据权利要求1所述的系统，其特征在于，所述巴特沃思低通滤波器采样直接II型二阶区段(SOS)结构；

所述二阶区段矩阵为：

$\mathrm{SOS}=\left[\begin{array}{cccccc}1& 2& 1& 1& -1.998& 0.998\\ 1& 2& 1& 1& -1.994& 0.994\\ 1& 2& 1& 1& -1.990& 0.990\\ 1& 2& 1& 1& -1.987& 0.987\\ 1& 2& 1& 1& -1.985& 0.985\\ 1& 2& 1& 1& -1.983& 0.983\\ 1& 2& 1& 1& -1.982& 0.982\\ 1& 1& 0& 1& -0.990& 0\end{array}\right]$

所述巴特沃思低通滤波器的每一阶的增益因子向量G为：

G＝[0.00002 0.00002 0.00002 0.00002 0.00002 0.00002 0.00454 1]。

4.根据权利要求1所述的系统，其特征在于，小波变换单元包括：

分裂单元，用于对经采样的心电信号x(n)分裂为奇数序列x_o(n)和偶数序列x_e(n)，可由下式表示：

x_o(n)＝x(2n+1) x_e(n)＝x(2n)；

预测单元，用于基于所述偶数序列x_e(n)预测所述奇数序列x_o(n)，得到小波系数d(n)；可由下式表示：

d(n)＝x_o(n)-P(x_e(n))，其中P为预测算子；

更新单元，用于基于所述小波系数d(n)更新所述偶数序列x_e(n)，得到尺度系数c(n)；可用下式表示：

c(n)＝x_e(n)+U(x_o(n))，其中U为更新算子；

所述预测算子P由如下公式表示：

$P\left(n\right)=\frac{1}{2}[c(n-1)+c\left(n\right)]$

所述更新算子U由如下公式表示：

$U\left(n\right)=\left\{\begin{array}{cc}d(n-1)& {\mathrm{\&Delta;}}_L<{\mathrm{\&Delta;}}_R\\ \frac{1}{2}[d(n-1)+d\left(n\right)]& {\mathrm{\&Delta;}}_L={\mathrm{\&Delta;}}_{\text{R}}\\ d\left(n\right)& {\mathrm{\&Delta;}}_L>-{\mathrm{\&Delta;}}_R\end{array}\right.$

Δ_L＝|c(n)-d(n-1)|；Δ_R＝|c(n)-d(n)|。

5.根据权利要求4所述的系统，其特征在于，所述特征提取单元用于定位R波的位置，所述特征提取单元包括：

阈值计算单元，用于依据小波变换单元进行小波分解直至达到预设尺度值后输出的小波系数d(n)计算阈值；

检测单元，用于依据所述阈值检测所述心电信号中R波的位置；

判断单元，用于依据尺度系数c(n)判断所述R波的位置是否正确，若否，则以预设步长调整所述阈值并由检测单元重新检测R波的位置。

6.根据权利要求5所述的系统，其特征在于，所述预设尺度值为4，等长区间为fs/4，其中fs为心电信号x(n)的采样频率；

所述阈值计算单元进一步用于：将小波系数d(n)划分为多个等长区间；计算每一等长区间的极大值；计算所述多个等长区间的极大值的均值；选取所述均值的一半作为阈值；

所述检测单元进一步用于：依据所述阈值检测所述心电信号，并依据检测结果定位R波的位置并判断R波间距是否大于1/150s，若是，则由判断单元进行进一步的处理；

所述判断单元进一步用于：依据尺度系数c(n)判断R波的位置是否满足c₁(n)<c₂(n)，c₁(n)<c₃(n)，c₁(n)<c₄(n)，若否，则以预设步长调整所述阈值并由所述检测单元重新检测R波的位置；其中c₁(n)为第1次小波变换的尺度系数，c₂(n)为第2次小波变换的尺度系数，c₃(n)为第3次小波变换的尺度系数，c₄(n)为第4次小波变换的尺度系数。

7.根据权利要求1所述的系统，其特征在于，所述特征提取单元还用于定位QRS波的起始点位置、终止点位置。

8.根据权利要求1所述的系统，其特征在于，所述系统还包括：

总线控制单元，用于实现所述系统中各个单元之间的互联；

外设单元，用于对所述系统进行功能扩展；

所述外设单元包括I²C控制器、SPI控制器、GPIO控制器和SD卡主控制器。