CN105870549A - 一种分数阶90°移相器及其qrs检波方法 - Google Patents

Abstract

本发明公开了一种分数阶90°移相器及其QRS检波方法，属于数字信号处理技术领域。基于Grünwald‑Letnikov分数阶微积分，定义分数阶中心差分算子，其本质为一维分数阶90°移相器，其相频特性与传统一阶导数相同，恒为90°相移，由于其幅频特性为频率的分数阶次幂函数，可以通过调节分数阶次在检测精度与抗噪性之间取得一定平衡，将心电信号x(n)通过所设计的分数阶90°移相器，进行R波检测，在已确定R波的基础上，位于R波之前的极小值点为Q波，位于R波之后的极小值点为S波。本发明计算简便，对受噪严重、波形异常的心电信号具有较好的检波效果。

Description

一种分数阶90°移相器及其QRS检波方法

技术领域

本发明公开了一种分数阶90°移相器及其QRS检波方法，属于数字信号处理技术领域。

背景技术

一个正常的心电信号是由一系列的波形组成的，主要包含P波、QRS复合波、T波、ST段以及U波。QRS复合波的检测是心电信号波形检测中的首要问题。精确可靠的QRS波检测是诊断心电疾患的首要依据，是其他波形诊断的重要参考。QRS复合波的检测主要是基于其形态、频率特性与其他波形的差异来实现的。近年来，QRS复合波检测手段主要有：阈值法，小波分析法，模板匹配法，神经网络法等等。

QRS复合波检测的首要问题是波峰值点的检测，即R波的检测。由于R波具有幅度大和斜率高的特点，因此最为常见的R波检测方法有幅度法、面积法和斜率法。幅度法和面积法在遇到非常规的心电信号，如高T波、高P波时容易出错；斜率法则优于这两种方法。斜率法，又称差分阈值法，是将滤波后的信号进行差分运算，再采用预先给定的阈值进行判别。另一种阈值法则考虑心电信号QRS复合波的宽度检测，当检测到QRS时输出一个脉冲信号，并与某个事先设定的阈值进行比较，从而检测出QRS复合波。总的来说，阈值法设计思路简单，计算量小，易于工程实现，但同时阈值的设置存在随机性，精度不能得到一致的保证。

小波变换可以将信号分解成不同尺度下的子信号，属于局部化分析方法。它研究信号频谱和小波变换尺度之间的关系，通过将心电信号中不同波形和干扰源在频谱上进行对比，从频谱的分布上检测出QRS波。小波变换通过反复实验，学习并优化阈值的选择，运用多标准判定，保证了较高的检测精度。小波变换的另一个重要特点是抗噪性强，减少漏检和误检的概率，其缺点是计算量大，不适宜实时工程应用。

模板匹配法将QRS波群、T波等近似认为是单独的模板，通过幅值归一化处理，信号中不同的幅频分量对应各个模板的能量分布，通过不同幅频分量的模板之间的匹配来检测出QRS波群，将对QRS波群的检测转化为模板间的匹配鉴别。该方法对幅值判断和频率判断进行了统一，优点是抗干扰能力强，精度高，缺点是耗时长，计算量大。

神经网络法首先对信号进去降噪处理，得到噪声较小段信号的QRS复合波模板作为参考模板。然后按次序在心电信号中找到匹配度较高的心电信号QRS复合波。神经网络的特点是判别速度快，但需要较长的训练时间。

发明内容

发明目的：为了克服现有技术中存在的不足，本发明提供一种分数阶90°移相器及其QRS检波方法，该移相器相频特性与传统一阶导数相同，恒为90°相移，由于其幅频特性为频率的分数阶次幂函数，可以通过调节分数阶次在检测精度与抗噪性之间取得一定平衡。

为实现上述目的，本发明采用的技术方案为：一种分数阶90°移相器，包括中心差分Grünwald-Letnikov微积分算子所述中心差分Grünwald-Letnikov微积分算子包括前向Grünwald-Letnikov微积分算子后向Grünwald-Letnikov微积分算子以及比例运算器；所述前向Grünwald-Letnikov微积分算子滤波器后向Grünwald-Letnikov微积分算子滤波器相互并联在一起之后再与比例运算器串联；所述输入信号分别通过前向Grünwald-Letnikov微积分算子滤波器后向Grünwald-Letnikov微积分算子滤波器进行滤波，滤波后的信号通过比例运算器调节后得到输出信号。

优选的：所述中心差分Grünwald-Letnikov微积分算子滤波器与前向Grünwald-Letnikov微积分算子滤波器和后向Grünwald-Letnikov微积分算子滤波器的关系为：

$D_d^v=\frac{1}{2sin\left(\frac{v\mathrm{\π}}{2}\right)}(D_l^v-D_r^v)$

其中，v≠2n，n∈Z。

优选的：所述输出信号为：

Y(e^iω)＝X(e^iω)·isgn(ω)|ω|^v

其中，v为分数阶微积分滤波器的阶次，取任意非偶实数；X(e^iω)为输入信号；i为虚数单位，ω为信号频率。

优选的：所述中心Grünwald-Letnikov微积分算子滤波器前向Grünwald-Letnikov微积分算子滤波器以及后向Grünwald-Letnikov微积分算子滤波器均通过差分方程的形式来构建。

一种基于分数阶90°移相器的QRS检波方法，将心电信号x(n)通过所设计的分数阶90°移相器，进行R波检测，在已确定R波的基础上，位于R波之前的极小值点为Q波，位于R波之后的极小值点为S波。

优选的：将心电信号x(n)通过所设计的分数阶90°移相器，利用阀值法滤除心电信号中的高频小幅毛刺，然后进行R波检测。

优选的：所述R波检测，其步骤为：1、区间法选出备选的正波和负波；2、排除法剔除不符合要求的极值点，选择相邻两个极值点中正确的极值点；3、配对法进一步排除T波和P波；4、坐标修正法确定原始信号中R波的正确坐标。

优选的：所述QRS波形检测，Q波和S波检测：位于R波之前的极小值点为Q波，之后的极小值点为S波。

本发明提供的分数阶90°移相器及其QRS波形检测，相比现有技术，具有以下有益效果：

(1)分数阶90°移相器，其相频特性与传统一阶导数相同，恒为90°相移，然而由于其幅频特性为频率的分数阶次幂函数，可以通过调节分数阶次在检测精度与抗噪性之间取得一定平衡。

(2)设计方法简便，通过对信号进行分数阶微积分运算，分别实现对信号的分数阶前向滤波和后向滤波，通过并联设计最终得到分数阶90°移相器，设计方法简便且算法效率高。

(3)设计灵活性大，通过调节分数阶微积分的阶次，可以依据实际需要的设计指标选择合适的阶次。

综上所述：将本发明基于分数阶90°移相器的检波方法应用于心电信号QRS波形检测，可取得较好的检测效果。设计方法简便，设计灵活性大，通过调节分数阶微积分的阶次，可以依据实际需要的设计指标选择合适的阶次。本发明的分数阶90°移相器的检波方法可应用于非平稳信号处理，尤其适用于批信号处理领域。

附图说明

图1为本发明分数阶90°移相器的流程图；

图2为卷积模板曲线图，其中v分别取-0.5和0.5两种情况，m＝100；

图3为本发明分数阶90°移相器的QRS波形检测流程图；

图4为分数阶90°移相器作用于输入信号前后的变化图；

图5为阈值处理前后信号对比图；

图6为检测R波第一步，采用区间法选出备选的正波和负波后的曲线图；

图7为检测R波第二步，采用排除法剔除不符合要求的极值点后的曲线图；

图8为检测R波第三步，采用配对法进一步排除T波和P波后的曲线图；

图9为检测R波第四步，采用坐标修正法确定原始信号中R波的正确坐标的曲线图；

图10为美国麻省理工学院MIT-BIH心电数据库第106号数据检测出的QRS波群。

具体实施方式

本发明所解决的技术问题在于获得一种比现有的基于一阶差分检波方法更理想的分数阶信号处理算法，所设计的滤波器其本质为一维分数阶90°移相器，其相频特性与传统一阶导数相同，恒为90°相移，由于其幅频特性为频率的分数阶次幂函数，可以通过调节分数阶次在检测精度与抗噪性之间取得一定平衡。在具体说明本发明内容之前，有必要对本说明书所用符号涵义进行简要说明。

根据一阶左导数的定义，可得到：

$D_{l}f\left(t\right)=\underset{h\→0}{\lim }\frac{f\left(t\right)-f(t-h)}{h}$

其中，f(t)为输入信号。

将定义为一阶n重左导数，可得到：

$D_l^{n}f\left(t\right)=\underset{h\→0}{\lim }\frac{1}{h^n}\underset{k=0}{\overset{n}{\mathrm{\Σ}}}{(-1)}^k\left(\begin{array}{c}n\\ k\end{array}\right)f(t-\mathrm{kh})$

其中，为二项式系数。

当n取分数时，该二项式系数可替换为其中v为任意非整数，由此我们定义v阶前向Grünwald-Letnikov微积分算子用符号表示，对函数f(t)进行v阶前向Grünwald-Letnikov微积分，得到下式：

$D_l^{v}f\left(t\right)=\underset{h\→0}{\lim }\frac{1}{h^n}\underset{k=0}{\overset{\[\frac{t-a}{h}\]}{\Σ}}{(-1)}^k\left(\begin{array}{c}v\\ k\end{array}\right)f(t-kh)$

其中，t和a分别是前向Grünwald-Letnikov微积分算子的上限和下限，该微积分算子将n阶微分与n阶积分运算统一用一个微积分算子来表示，且将整数域的运算推广到非整数域。当v＞0时，上式表示微分运算，当v＜0时，上式表示积分运算，微积分算子由此得名，信号进行v阶前向Grünwald-Letnikov微积分的傅里叶变换表示为：

$F\[D_l^{v}f\left(t\right)\]={\left(i\mathrm{\ω}\right)}^{v}f\left(\mathrm{\ω}\right)$

因此，v阶前向Grünwald-Letnikov微积分算子的频域函数为：

$H_l^v\left(e^{iw}\right)={\left(i\mathrm{\ω}\right)}^v={\left|\mathrm{\ω}\right|}^v{e}^{i\frac{v\mathrm{\π}}{2}\mathrm{sgn}\left(\mathrm{\ω}\right)}$

对应的幅频和相频响应分别为：

根据一阶右导数的定义，可得到：

$D_{r}f\left(t\right)=\underset{h\→0}{\lim }\frac{f\left(t\right)-f(t+h)}{h}$

其中f(t)为输入信号。

将定义为一阶n重右导数，可得到：

$D_r^{v}f\left(t\right)=\underset{h\→0}{\lim }\frac{1}{h^v}\underset{k=0}{\overset{\[\frac{b-t}{h}\]}{\Σ}}{(-1)}^k\left(\begin{array}{c}v\\ k\end{array}\right)f(t-kh)$

其中，为二项式系数。

当n取分数时，该二项式系数可替换为其中v为任意非整数。由此我们定义v阶后向Grünwald-Letnikov微积分算子用符号表示。对函数f(t)进行v阶后向Grünwald-Letnikov微积分，得到下式：

$D_r^{v}f\left(t\right)=\underset{h\→0}{\lim }\frac{1}{h^v}\underset{k=0}{\overset{\[\frac{b-t}{h}\]}{\Σ}}{(-1)}^k\left(\begin{array}{c}v\\ k\end{array}\right)f(t-kh)$

其中，b和t分别是后向Grünwald-Letnikov微积分算子的上限和下限。信号进行v阶后向Grünwald-Letnikov微积分的傅里叶变换表示为：

$F\[D_r^{v}f\left(t\right)\]={(-i\mathrm{\ω})}^{v}f\left(\mathrm{\ω}\right)$

因此，v阶后向Grünwald-Letnikov微积分算子的频域函数为：

$H_r^v\left(e^{iw}\right)={(-i\mathrm{\ω})}^v={\left|\mathrm{\ω}\right|}^v{e}^{-i\frac{v\mathrm{\π}}{2}\mathrm{sgn}\left(\mathrm{\ω}\right)}$

对应的幅频和相频响应分别为：

前向Grünwald-Letnikov微积分算子和后向Grünwald-Letnikov微积分算子可被视为两种特殊的滤波器，其幅值响应随信号频率和分数阶微积分阶次变化而变化，信号通过前向Grünwald-Letnikov微积分算子相角滞后而信号通过后向Grünwald-Letnikov微积分算子相角超前

下面结合附图对本发明作更进一步的说明。

本发明设计的分数阶90°移相器，如图1所示，包括中心Grünwald-Letnikov微积分算子滤波器所述中心Grünwald-Letnikov微积分算子滤波器包括前向Grünwald-Letnikov微积分算子滤波器后向Grünwald-Letnikov微积分算子滤波器以及比例运算器，所述前向Grünwald-Letnikov微积分算子滤波器后向Grünwald-Letnikov微积分算子滤波器相互并联在一起之后再与比例运算器串联；所述输入信号分别通过前向Grünwald-Letnikov微积分算子滤波器和后向Grünwald-Letnikov微积分算子滤波器进行滤波，滤波后的信号通过比例运算器调节后得到输出信号。

用数学表达式来说明与的关系，如下式

$D_d^v=\frac{1}{2sin\left(\frac{v\mathrm{\π}}{2}\right)}(D_l^v-D_r^v)$

其中v≠2n，n∈Z。

因此，对函数f(t)进行v阶中心Grünwald-Letnikov微积分可以表示为：

$D_d^{v}f\left(t\right)=\frac{1}{2\sin \left(\frac{v\mathrm{\π}}{2}\right)}\underset{h\→0}{\lim }\frac{1}{h^v}\underset{k=0}{\overset{N}{\Σ}}{(-1)}^k\left(\begin{array}{c}v\\ k\end{array}\right)\[f\left(t-kh\right)-f(t+kh)\]$

其中

假设将中心Grünwald-Letnikov微积分算子应用于一个有限长的输入信号序列f(t)，相应的傅里叶变换为

$F\[D_d^{v}f\left(t\right)\]=\frac{1}{2sin\left(\frac{v\mathrm{\π}}{2}\right)}\[{\left(i\mathrm{\ω}\right)}^v-{(-i\mathrm{\ω})}^v\]F\left(\mathrm{\ω}\right)=i\mathrm{sgn}\left(\mathrm{\ω}\right){\left|\mathrm{\ω}\right|}^v$

因此，中心Grünwald-Letnikov微积分算子的频率函数为

$H_d^v\left(\mathrm{\ω}\right)=i\mathrm{sgn}\left(\mathrm{\ω}\right){\left|\mathrm{\ω}\right|}^v$

对应的幅频特性和相频特性分别为：

具体实施时，由于数字信号处理是基于离散信号进行的，因此我们需将上述时域函数进行离散化操作，即用差分方程的形式来构建数字信号滤波器。取h＝1，且当或足够大的时候，可省略公式中的取极限这一步骤。前向Grünwald-Letnikov微积分算子后向Grünwald-Letnikov微积分算子和中心Grünwald-Letnikov微积分算子的差分近似定义如下：

${\mathrm{\Δ}}_l^{v}f\left(t\right)=\underset{i=1}{\overset{N}{\Σ}}{(-1)}^k\left(\begin{array}{c}v\\ k\end{array}\right)f(t-k)$

${\mathrm{\Δ}}_r^{v}f\left(t\right)=\underset{k=0}{\overset{N}{\Σ}}{(-1)}^k\left(\begin{array}{c}v\\ k\end{array}\right)f(t+k)$

其中v≠2n，n∈Z

令

为了计算出前向Grünwald-Letnikov微积分算子的差分方程，取卷积模板模板长度为m+1。那么，前向Grünwald-Letnikov微积分算子的差分方程可写为

${\mathrm{\Δ}}_l^{v}f\left(t\right)=f\left(t\right)*M_l^v$

同理，为了计算出后向Grünwald-Letnikov微积分算子的差分方程，取卷积模板模板长度为m+1。那么，后向Grünwald-Letnikov微积分算子的差分方程可写为

当中心差分算子应用在有限长度的信号f(t)上时，差分方程形式为

${\mathrm{\Δ}}_d^{v}f\left(t\right)=\frac{1}{2sin\left(\frac{v\mathrm{\π}}{2}\right)}\[{\mathrm{\Δ}}_l^{v}f\left(t\right)-{\mathrm{\Δ}}_r^{v}f\left(t\right)\]=f\left(t\right)*M_d^v\left(t\right)$

其中，卷积模板长度为2m+1。

图2所示为卷积模板曲线图，其中v分别取0.5和-0.5两种情况，m＝100，即卷积模板长度为201，两条卷积模板曲线奇对称且在两端趋于0，表现了90°移相的通性，且保证了适度的计算复杂度。

图3为QRS波形检测流程图，首先，我们将原始心电信号x(n)通过所设计的分数阶90°移相器得到信号y(n)作为后续实验对象，见图4，接着，利用常见的阈值法滤除心电信号中的高频小幅毛躁，使得滤波后的输出信号较平稳的区间更加光滑，便于后续R波检测，图5是阈值作用后的心电信号，两条水平线分别为正阈值与负阈值。阈值法作用后得到的信号，实验中记为信号xdf。

接着，第三个环节通过4个步骤得到R波的正确坐标。

图6为检测R波的第一步，用区间法选出备选的正波和负波。将信号等距离划分区间，逐个区间检索出正、负极值。每个区间包含50个采样点(由于信号采样频率为360Hz，因此50个采样点对应的区间时长大约为0.14s)。如果区间内有非零的极值点(绝对值)，则把这些极值点选出来并记录坐标。R-R峰值间隔医学统计上大约为0.4-1.2秒之间，每个R-R峰间，大约有3-9个搜索区间不等。图6(b)中，圆圈表示的是区间法确定的正负极值点。可以看到，通过区间法初步筛选出来的极值点很多，且由于区间起始点的随机性和区间大小的设定具有先决性，检测出的“假”极值点很多。

图7为第二步，用排除法剔除不符合要求的极值点，选择相邻两个的极值点中正确的极值点。测量相邻两个非零正极值点间的距离。选择出所有相邻距离小于区间长度的幅值较小的点，并删除。

图8为第三步，配对法进行进一步排除T波和P波。为了防止偶尔出现的高幅值T波和P波，需要使用配对法进行排除，找到正确的正负极值对。测量两个相邻的正负极值点(正负极值点放在一起考虑)。有两种情况：1)若距离小于等于区间长度，即50个采样点或0.14秒，且正极值点出现在前，负极值点在后，则两个极值点为一个正确的正-负配对。2)若距离大于区间长度，两个极值点不构成一对配对，需要删除。图8(b)中显示了将错误归为R波的T波排除后的效果。

图9为第四步，坐标修正法确定原始信号中R波的正确坐标。一旦我们在差分信号中找到R波，搜索将鉴于此时确定的信号坐标，回到原始信号，将此坐标作为新的搜索中心点，向此坐标前后拓宽一定的长度，得到新的搜索窗，窗口大小大概取在125ms～250ms。通过这个步骤的处理，弥补移相和滤波过程中造成的信号点R波相位的偏移。图9(b)展示的是在原始信号x(n)中使用本方法检测出来的R波，图中x指原始心电信号，R指检测出来的R波。

图10为美国麻省理工学院MIT-BIH心电数据库第106号数据检测出的QRS波群，由于Q波和S波分别是位于R波前后的负向极值点波，且距离R波较近，QRS波群的长度一般为0.08～0.12s(即29-44个采样点之间)，因此以R波为中心，前后各22个采样点作为检测窗口。位于R波之前的极小值点为Q波，之后的极小值点为S波。

Claims (7)

1.一种分数阶90°移相器，其特征在于：包括中心Grünwald-Letnikov微积分算子滤波器所述中心Grünwald-Letnikov微积分算子包括前向Grünwald-Letnikov微积分算子和后向Grünwald-Letnikov微积分算子所述输入信号分别通过前向Grünwald-Letnikov微积分算子和后向Grünwald-Letnikov微积分算子进行滤波，滤波后的信号通过比例运算器调节后得到输出信号。

2.根据权利要求1所述的分数阶90°移相器，其特征在于：所述中心Grünwald-Letnikov微积分算子滤波器与前向Grünwald-Letnikov微积分算子滤波器以及后向Grünwald-Letnikov微积分算子滤波器的关系为：

$D_d^v=\frac{1}{2sin\left(\frac{v\mathrm{\π}}{2}\right)}(D_l^v-D_r^v)$

其中，v≠2n，n∈Z。

3.根据权利要求2所述的分数阶90°移相器，其特征在于：所述输出信号为：

Y(e^iω)＝X(e^iω)·isgn(ω)|ω|^v

其中，v为分数阶微积分滤波器的阶次，取任意非偶实数；X(e^iω)为输入信号；i为虚数单位，ω为信号频率。

4.根据权利要求3所述的分数阶90°移相器，其特征在于：所述中心Grünwald-Letnikov微积分算子滤波器前向Grünwald-Letnikov微积分算子滤波器以及后向Grünwald-Letnikov微积分算子滤波器均通过差分方程的形式来构建。

5.一种基于权利要求1所述的分数阶90°移相器的QRS检波方法，其特征在于：将心电信号x(n)通过所设计的分数阶90°移相器，进行R波检测，在已确定R波的基础上，位于R波之前的极小值点为Q波，位于R波之后的极小值点为S波。

6.根据权利要求5所述的QRS检波方法，其特征在于：所述R波检测的步骤为：1、区间法选出备选的正波和负波；2、排除法剔除不符合要求的极值点，选择相邻两个极值点中正确的极值点；3、配对法进一步排除T波和P波；4、坐标修正法确定原始信号中R波的正确坐标。

7.根据权利要求5所述的QRS检波方法，其特征在于：所述QRS波形检测、Q波和S波检测：位于R波之前的极小值点为Q波，之后的极小值点为S波。