如何检测（心电图）波中的模式？

Posted 2023-02-14

技术标签:

【中文标题】如何检测（心电图）波中的模式？【英文标题】：How to detect patterns in (electrocardiography) waves? 【发布时间】：2011-01-12 20:37:20 【问题描述】：

我正在尝试从心电图中读取图像并检测其中的每个主要波（P 波、QRS 复合波和 T 波）。我可以读取图像并获得矢量（如(4.2; 4.4; 4.9; 4.7; ...)）。我需要一种算法，可以遍历这个向量并检测这些波中的每一个何时开始和结束。一个例子：

如果它们始终具有相同的大小，或者如果我提前知道心电图有多少波，那就容易了。鉴于波浪：

我提取向量：

[0; 0; 20; 20; 20; 19; 18; 17; 17; 17; 17; 17; 16; 16; 16; 16; 16; 16; 16; 17; 17; 18; 19; 20; 21; 22; 23; 23; 23; 25; 25; 23; 22; 20; 19; 17; 16; 16; 14; 13; 14; 13; 13; 12; 12; 12; 12; 12; 11; 11; 10; 12; 16; 22; 31; 38; 45; 51; 47; 41; 33; 26; 21; 17; 17; 16; 16; 15; 16; 17; 17; 18; 18; 17; 18; 18; 18; 18; 18; 18; 18; 17; 17; 18; 19; 18; 18; 19; 19; 19; 19; 20; 20; 19; 20; 22; 24; 24; 25; 26; 27; 28; 29; 30; 31; 31; 31; 32; 32; 32; 31; 29; 28; 26; 24; 22; 20; 20; 19; 18; 18; 17; 17; 16; 16; 15; 15; 16; 15; 15; 15; 15; 15; 15; 15; 15; 15; 14; 15; 16; 16; 16; 16; 16; 16; 16; 16; 16; 15; 16; 15; 15; 15; 16; 16; 16; 16; 16; 16; 16; 16; 15; 16; 16; 16; 16; 16; 15; 15; 15; 15; 15; 16; 16; 17; 18; 18; 19; 19; 19; 20; 21; 22; 22; 22; 22; 21; 20; 18; 17; 17; 15; 15; 14; 14; 13; 13; 14; 13; 13; 13; 12; 12; 12; 12; 13; 18; 23; 30; 38; 47; 51; 44; 39; 31; 24; 18; 16; 15; 15; 15; 15; 15; 15; 16; 16; 16; 17; 16; 16; 17; 17; 16; 17; 17; 17; 17; 18; 18; 18; 18; 19; 19; 20; 20; 20; 20; 21; 22; 22; 24; 25; 26; 27; 28; 29; 30; 31; 32; 33; 32; 33; 33; 33; 32; 30; 28; 26; 24; 23; 23; 22; 20; 19; 19; 18; 17; 17; 18; 17; 18; 18; 17; 18; 17; 18; 18; 17; 17; 17; 17; 16; 17; 17; 17; 18; 18; 17; 17; 18; 18; 18; 19; 18; 18; 17; 18; 18; 17; 17; 17; 17; 17; 18; 17; 17; 18; 17; 17; 17; 17; 17; 17; 17; 18; 17; 17; 18; 18; 18; 20; 20; 21; 21; 22; 23; 24; 23; 23; 21; 21; 20; 18; 18; 17; 16; 14; 13; 13; 13; 13; 13; 13; 13; 13; 13; 12; 12; 12; 16; 19; 28; 36; 47; 51; 46; 40; 32; 24; 20; 18; 16; 16; 16; 16; 15; 16; 16; 16; 17; 17; 17; 18; 17; 17; 18; 18; 18; 18; 19; 18; 18; 19; 20; 20; 20; 20; 20; 21; 21; 22; 22; 23; 25; 26; 27; 29; 29; 30; 31; 32; 33; 33; 33; 34; 35; 35; 35; 0; 0; 0; 0;]

我想检测，例如：

[19 - 37] 中的 P 波。 [51 - 64] 中的 QRS 波群。等

【问题讨论】：

我知道a guy 在这个领域工作过。在这里您可以找到list of his publications。如果我没记错的话，他使用隐马尔可夫模型来可靠地检测针对已知形状训练集的波，但您可以在论文中找到更多详细信息。你已经有了很多好的答案。我很惊讶没有人建议 PhysioToolkit 中的“WFDB 软件包”，特别是 ecgpuwave 我对关于检测时间序列数据中的模式的类似问题的回答在这里 - ***.com/a/11903770/1149913 - 并包括 python 代码。我的方法是“切换自回归隐马尔可夫模型”（谷歌搜索一些相关出版物的短语）。 【参考方案1】：

我要做的第一件事就是看看那里已经有什么。事实上，这个具体问题已经得到了大量研究。以下是一些非常简单的方法的简要概述：link。

我也必须回复另一个答案。我从事信号处理和音乐信息检索方面的研究。从表面上看，这个问题确实与发病检测类似，但问题的上下文并不相同。这种类型的生物信号处理，即检测 P、QRS 和 T 相位，可以利用这些波形中每一个的特定时域特征的知识。 MIR 中的发病检测并没有，真的。（至少不可靠。）

一种适用于 QRS 检测（但不一定适用于音符起始检测）的方法是动态时间扭曲。当时域特性保持不变时，DTW 可以非常好地工作。这是一篇使用 DTW 解决此问题的简短 IEEE 论文：link。

这是一篇不错的 IEEE 杂志文章，它比较了许多方法：link。您会看到许多常见的信号处理模型已经被尝试过。略读论文，然后尝试一个您在基本水平上理解的内容。

编辑：浏览这些文章后，基于小波的方法对我来说似乎最直观。 DTW 也可以很好地工作，并且存在 DTW 模块，但小波方法对我来说似乎是最好的。其他人通过利用信号的导数来回答。我的第一个链接检查了 1990 年之前的方法，但我怀疑它们不如更现代的方法强大。

编辑：当我有机会时，我会尝试给出一个简单的解决方案，但是我认为小波在这里适合的原因是因为它们在参数化各种形状方面非常有用时间或幅度缩放。换句话说，如果您有一个具有相同重复时间形状但在不同时间尺度和幅度上的信号，小波分析仍然可以将这些形状识别为相似（粗略地说）。另请注意，我将过滤器组归为此类。类似的事情。

【讨论】：

我知道它要求很多，但我仍然希望看到一个简单的解决方案。第一个链接又失效了，但我找到了here。它是 IEEE，所以我想它应该保持不变。当然，如果有成本障碍，有类似主题的many papers。【参考方案2】：

这个难题的一部分是“onset detection”，并且已经编写了许多复杂的算法来解决这个问题。这里有更多关于onsets的信息。

下一个是Hamming Distance。该算法允许您进行模糊比较，输入是 2 个数组，输出是整数“距离”或 2 个数据集之间的差异。数字越小，2 越相似。这非常接近您的需要，但并不准确。我继续对汉明距离算法进行了一些修改以计算新的距离，它可能有一个名字，但我不知道它是什么。基本上，它将数组中每个元素之间的绝对距离相加并返回总数。这是python中的代码。

import math

def absolute_distance(a1, a2, length):
       total_distance=0
       for x in range(0,length):
               total_distance+=math.fabs(a1[x]-a2[x])
       return total_distance

print(absolute_distance([1,3,9,10],[1,3,8,11],4))

此脚本输出 2，即这 2 个数组之间的距离。

现在将这些部分放在一起。您可以使用起始检测来查找数据集中所有波的开始。然后，您可以循环通过这些位置，将每个波与样本 P-Wave 进行比较。如果您击中 QRS 波群，则距离将是最大的。如果你击中另一个 P-Wave，这个数字不会为零，但会小得多。任何 P-Wave 和任何 T-Wave 之间的距离都会非常小，但是如果您做出以下假设，这不是问题：

The distance between any p-wave and any other p-wave will be smaller than the distance between any p-wave and any t-wave.

该系列看起来像这样：pQtpQtpQt... p 波和 t 波彼此相邻，但由于该序列是可预测的，因此更易于阅读。

另一方面，这个问题可能有一个基于微积分的解决方案。然而，在我看来，曲线拟合和积分使这个问题更加混乱。我写的距离函数将找到非常相似的面积差减去两条曲线的积分。

也许可以牺牲起始计算来支持一次迭代 1 个点，从而执行 O(n) 距离计算，其中 n 是图中的点数。如果您有所有这些距离计算的列表并且知道 50 个 pQt 序列在哪里，那么您将知道 50 个最短距离不重叠所有 p 波的位置。 Bingo！为了简单起见，怎么样？但是，由于距离计算次数的增加，代价是效率损失。

【讨论】：

检测时域幅度的增加等更简单的技术通常会导致大量误报或误报。我害怕。我在我的（不是最佳的）解决方案中提出了同样的问题。是的，您的算法很有趣，但您的成功可能有限。这是一个非常复杂的问题，没有完美的解决方案。【参考方案3】：

您可以使用cross-correlation。获取每个模式的模型样本并将它们与信号相关联。您将获得相关性高的峰值。我希望使用这种提取 qrs 和 t 波的技术取得良好的效果。之后，您可以通过查找 qrs 之前的相关信号上的峰值来提取 p 波。

互相关是一种非常容易实现的算法。基本上：

x is array with your signal of length Lx
y is an array containing a sample of the signal you want to recognize of length Ly
r is the resulting correlation

for (i=0; i<Lx - Ly; i++)
  r[i] = 0;
  for (j=0; j<Ly ; j++)
    r[i] += x[i+j]*y[j];

并寻找 r 中的峰值（例如，超过阈值的值）

【讨论】：

这是一个很好的第一个尝试方法，因为波浪总是遵循特定的模式。但是对于这个问题，时间缩放和幅度缩放都可能会发生变化，因此最终，这种方法不会在受试者之间保持稳健。是的，这只是第一种方法。不健壮，但足够简单，可以尝试编码。模式匹配通常是最简单的技术，并且仍然可以提供一些结果。当然，小波肯定要好得多。【参考方案4】：

我要做的第一件事是简化数据。

不是分析绝对数据，而是分析从一个数据点到下一个数据点的变化量。

这是一个快速的单行，它将; 分隔的数据作为输入，并输出该数据的增量。

perl -0x3b -ple'( $last, $_ ) = ( $_, $_-$last )' < test.in > test.out

在您提供的数据上运行它，这是输出：

0;0;20;0;0;-1;-1;-1;0;0;0;0;-1;0;0;0;0;0;0;1;0;1 ;1;1;1;1;1;0;0;2;0;-2;-1;-2;-1;-2;-1;0;-2;-1;1;-1; 0;-1;0;0;0; 0;-1;0;-1;2;4;6;9;7;7;6;-4;-6;-8;-7;-5;-4;0;-1;0;- 1;1;1;0;1;0;-1;1;0;0;0;0;0;0;-1;0;1;1;-1;0;1;0;0;0 ;1;0;-1;1; 2;2;0;1;1;1;1;1;1;1;0;0;1;0;0;-1;-2;-1;-2;-2;-2;-2 ;0;-1;-1;0;-1;0;-1;0;-1;0;1;-1;0;0;0;0;0;0;0;0;-1; 1;1;0;0;0; 0;0;0;0;0;-1;1;-1;0;0;1;0;0;0;0;0;0;0;-1;1;0;0;0;0 ;-1;0;0;0;0;1;0;1;1;0;1;0;0;1;1;1;0;0;0;-1;-1;-2;- 1;0;-2;0; -1;0;-1;0;1;-1;0;0;-1;0;0;0;1;5;5;7;8;9;4;-7;-5;-8 ;-7;-6;-2;-1;0;0;0;0;0;1;0;0;1;-1;0;1;0;-1;1;0;0;0 ;1;0;0;0; 1;0;1;0;0;0;1;1;0;2;1;1;1;1;1;1;1;1;1;-1;1;0;0;-1; -2;-2;-2;-2;-1;0;-1;-2;-1;0;-1;-1;0;1;-1;1;0;-1;1; -1;1;0;-1; 0;0;0;-1;1;0;0;1;0;-1;0;1;0;0;1;-1;0;-1;1;0;-1;0;0 ;0;0;1;-1;0;1;-1;0;0;0;0;0;0;1;-1;0;1;0;0;2;0;1;0; 1;1;1;-1; 0;-2;0;-1;-2;0;-1;-1;-2;-1;0;0;0;0;0;0;0;0;-1;0;0; 4;3;9;8;11;4;-5;-6;-8;-8;-4;-2;-2;0;0;0;-1;1;0;0;1; 0;0;1;-1; 0;1;0;0;0;1;-1;0;1;1;0;0;0;0;1;0;1;0;1;2;1;1;2;0;1 ;1;1;1;0;0;1;1;0;0;-35;0;0;0;

_{在上述文本中插入了换行符，而这些换行符原本不存在于输出中。}

完成之后，找到 qrs 复合体就很简单了。

perl -F';' -ane'@F = map  abs($_) > 2 and $_  @F; print join ";", @F'< test.out

;;20;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;4;6;9;7;7;6;-4;-6;-8;-7;-5;-4; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;5;5;7;8;9;4;-7;-5;-8;-7;-6 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;4;3;9;8;11;4;-5;-6;-8;-8;-4; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;-35 ;;;

20 和 -35 数据点来自以 0 开头和结尾的原始数据。

要找到其他数据点，您将不得不依赖模式匹配。

如果你看第一个 p 波，你可以清楚地看到一个模式。

0;0;0;0;0;0;1;0;1;1;1;1;1;1;0;0;2;0;-2;-1;-2;-1;-2;-1;0;-2;-1;1;-1;0;-1;0;0;0;0;
#           \________ up _______/   \________ down _________/

不过，在第二个 p 浪中看到模式并不容易。这是因为第二个分散得更远了

0;0;0;1;0;1;1;0;1;0;0;1;1;1;0;0;0;-1;-1;-2;-1;0;-2;0;-1;0;-1;0;1;-1;0;0;-1;0;0;0;
#     \________ up _______/       \________________ down ________________/

第三个 p 波比其他两个更不稳定。

0;0;0;0;0;1;-1;0;1;0;0;2;0;1;0;1;1;1;-1;0;-2;0;-1;-2;0;-1;-1;-2;-1;0;0;0;0;0;
#                \_______ up ______/  \__________ down __________/

您会以与 p 波类似的方式找到 t 波。主要区别在于它们发生的时间。

这应该足以帮助您入门。

_{这两个单线只是举例，不建议日常使用。}

【讨论】：

-1：我只看到了解决问题的手动方法，用户可以通过简单地绘制数据并自己选择起点来完成。【参考方案5】：

另外两个尖峰和尖谷也是qrs复合体吗？

在我的脑海中，我认为您需要做的是计算该图在每个点的斜率。然后您还需要查看斜率变化的速度（二阶导数？？？）。如果你有一个突然的变化，那么你知道你已经达到了某种尖锐的高峰。当然，你想限制对变化的检测，所以你可能想做一些类似“如果斜率在时间间隔 T 内变化 X”之类的事情，这样你就不会拾取图表中的微小颠簸。

我已经有一段时间没有做任何数学运算了……这似乎是一道数学题 ;) 哦，我也没有做过任何类型的信号分析 :)。

只是补充一点。我认为您也可以尝试信号平均。例如，平均最后 3 或 4 个数据点。我认为您也可以通过这种方式检测突然的变化。

【讨论】：

+1 表示有趣的算法。但我觉得这个问题有点复杂。是的，另外两个峰和谷是 qrs 复合体。该图像实际上有 3 个 p 波、3 个 qrs 复合波和 3 个 t 波。这是一种有趣的方法，但如果我没有函数，我看不出如何计算二阶导数。我认为您是说要对值的变化进行评分，并选择具有高分的变化，例如某事的开始和结束，是吗？我会尝试一下，当我得到一些结果时我会发布更新。谢谢你的回答。是的，差不多。您有点得分，但您是通过计算斜率或查看幅度如何随时间变化来做到这一点的。【参考方案6】：

我不是这个特定问题的专家，但我只是从更一般的知识中浮出水面：假设您知道 QRS 复合波（或其他特征之一，但我将使用 QRS 复合波这个例子）大约发生在一个长度为 L 的固定时间段内。我想知道你是否可以将其视为一个分类问题，如下所示：

将您的信号拆分为长度为 L 的重叠窗口。每个窗口中要么包含完整的 QRS 复合波，要么不包含完整的 QRS 波群。对每个窗口进行傅里叶变换。您的特征是每个频率的信号强度。在一些人工标注的数据上训练决策树、支持向量机等。

【讨论】：

【参考方案7】：

一种很可能产生良好结果的方法是曲线拟合：

将连续波分成多个区间（最好将区间边界设置在 qrs 复合波尖峰之间的一半左右）。一次只考虑一个间隔。

定义一个模型函数，该函数可用于近似心电图曲线的所有可能变化。这并不像最初看起来那么困难。模型函数可以构造为三个函数的总和，参数分别为每个波的原点（t_）、幅度（a_）和宽度（w_）。

   f_model(t) = a_p   *  f_p  ((t-t_p  )/w_p) + 
                a_qrs *  f_qrs((t-t_qrs)/w_qrs) +
                a_t   *  f_t  ((t-t_t  )/w_t)

函数f_p(t)、f_qrs(t)、f_t(t) 是一些简单的函数，可用于对三个波中的每一个进行建模。

使用拟合算法（例如 Levenberg-Marquardt-Algorithm http://en.wikipedia.org/wiki/Levenberg%E2%80%93Marquardt_algorithm）确定每个区间数据集的拟合参数 a_p、t_p、w_p、a_qrs、t_qrs、w_qrs、a_t、t_t、w_t。

参数t_p、t_qrs和t_p是你感兴趣的。

【讨论】：

【参考方案8】：

这是一个很好的问题！我有几个想法：

Dynamic Time Warping 在这里可能是一个有趣的工具。您将为您的三个类建立“模板”，然后使用 DTW 可以看到您的模板与信号“块”之间的相关性（将信号分解为 0.5 秒位，即 0-.5 。 1-.6 .2-.7...）。我用加速度计数据进行了类似的步态分析，效果相当好。

另一种选择是组合信号处理/机器学习算法。再次将您的信号分解成“块”。再次制作“模板”（每个类需要十几个）获取每个块/模板的FFT，然后使用Naïve Bayes Classifier（或另一个ML分类器，但NB应该削减它）进行分类你的三个班级中的每一个。我也在步态数据上尝试过这个，并且能够获得超过 98% 的精度和相对复杂的信号进行召回。让我知道这是如何工作的，这是一个非常令人兴奋的问题。

【讨论】：

【参考方案9】：

“Wavelet transform”可能是相关关键字。我曾经参加过一个使用这种技术检测嘈杂心电图中不同心跳阶段的人的演讲。

就我有限的理解而言，它有点像傅立叶变换，但使用的是（缩放）副本，在你的情况下是心跳形的脉冲。

【讨论】：

【参考方案10】：

首先，标准心电图波的各种成分可能会从任何给定的图中丢失。这样的情节通常是不正常的，通常预示着某种问题，但你不能保证他们在那里。

其次，识别它们既是一门艺术，也是一门科学，尤其是在出现问题的情况下。

我的方法可能是尝试训练神经网络来识别组件。你会给它前 30 秒的数据，标准化后最低点为 0，最高点为 1.0，它将有 11 个输出。不是异常评级的输出将是最后 10 秒的权重。 0.0 表示从现在开始的 -10 秒，而 1.0 表示现在。输出将是：

最近的 P 波开始的位置最近的 P 波在哪里结束最近 P 波的异常评级，其中一个极端为“不存在”。最近 QRS 波群的起始位置最近的 QRS 波群的 Q 部分变成了 R 部分。最近的 QRS 波群的 R 部分变成了 S 部分。最近的 QRS 波群在哪里结束。最近 QRS 波群的异常评级，其中一个极端为“不存在”。最近的 T 波从哪里开始。最近的 T 波在哪里结束。最近 T 波的异常评级，其中一个极端为“不存在”。

我可能会用人们建议的一些其他类型的分析来仔细检查这一点，或者使用这些其他类型的分析以及神经网络的输出来给你答案。

当然，这种对神经网络的详细描述不应被视为规定性的。例如，我确定我不一定会选择最佳输出，我只是对它们可能是什么提出了一些想法。

【讨论】：

【参考方案11】：

小波已被证明是定位此类数据中峰值“不同大小”的最佳工具 - 小波的缩放特性使其成为此类多尺度峰值检测的理想工具。这看起来像一个非平稳信号，因此使用 DFT 不会像某些人所建议的那样是正确的工具，但如果这是一个探索性项目，您可以考虑使用信号的频谱（基本上使用自相关的 FFT 估计）信号。）

Here 是一篇很棒的论文，回顾了几种峰值检测方法 - 这是一个很好的起点。

-保罗

【讨论】：

【参考方案12】：

我还没有彻底阅读彼此的答案，但我已经扫描了它们，我注意到没有人建议查看傅里叶变换来分割这些波。

在我看来，Harmonic analysis 在数学中的应用很明确。我可能遗漏了几个微妙的点。

Discrete Fourier Transform 系数为您提供构成离散时间信号的不同正弦分量的幅度和相位，这基本上就是您想要找到的问题状态。