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

Question

我正在尝试从心电图中读取图像并检测其中的每个主要波（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;]

我想检测，例如：

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

Answer 1

The first thing that I would do is see what is already out there. Indeed, this specific problem has already been heavily researched. Here is a brief overview of some really simple methods: link.

I must respond to another answer, as well. I do research in signal processing and music information retrieval. On the surface, this problem does appear similar to onset detection, but the problem context is not the same. This type of biological signal processing, i.e., detection of the P, QRS, and T phases, can exploit knowledge of specific time-domain characteristics of each of these waveforms. Onset detection in MIR doesn't, really. (Not reliably, at least.)

One approach that would work well for QRS detection (but not necessarily for note onset detection) is dynamic time warping. When time-domain characteristics remain invariant, DTW can work remarkably well. Here is a short IEEE paper that uses DTW for this problem: link.

This is a nice IEEE magazine article that compares many methods: link. You'll see that many common signal processing models have been tried. Skim the paper, and try one that you understand at a basic level.

EDIT: After browsing these articles, a wavelet-based approach seems most intuitive to me. DTW will work well, too, and there exist DTW modules out there, but the wavelet approach seems best to me. Someone else answered by exploiting derivatives of the signal. My first link examines methods from before 1990 that do that, but I suspect that they are not as robust as more modern methods.

EDIT: I'll try to give a simple solution when I get the chance, but the reason why I think wavelets are suited here are because they are useful at parameterizing a wide variety of shapes regardless of time or amplitude scaling. In other words, if you have a signal with the same repeated temporal shape but at varying time scales and amplitudes, wavelet analysis can still recognize these shapes as being similar (roughly speaking). Also note that I am sort of lumping filter banks into this category. Similar things.

Answer 2

这个难题的一部分是“发病检测”，并且已经编写了许多复杂的算法来解决这个问题。这是有关发病的更多信息。

下一块是汉明距离。该算法允许您进行模糊比较，输入是 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，即这两个数组之间的距离。

现在把这些碎片放在一起。您可以使用起始检测来查找数据集中所有波的开始。然后，您可以循环通过这些位置，将每个波与样本 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 波的所有位置重叠。 答对了！ 为了简单起见，那是怎么回事？然而，权衡是由于距离计算数量增加而导致效率损失。

Answer 3

您可以使用互相关。获取每个模式的模型样本并将它们与信号相关联。您将获得相关性高的峰值。我希望使用这种提取 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 中的峰值（例如，超过阈值的值）

Answer 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 波。主要区别在于它们何时发生。

---

这应该足以让您入门。

_{这两个单眼线仅作为示例，不建议日常使用。}

Answer 5

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

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

Answer 6

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

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

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

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

Answer 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 是您感兴趣的参数。

Answer 8

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

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

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

Answer 9

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

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

Answer 10

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

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

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

1. 最近的 P 波开始的地方
2. 最近的 P 波在哪里结束
3. 最近 P 波的异常评级，其中一个极端是“不存在”。
4. 最近的 QRS 波群开始的位置
5. 最近的 QRS 波群的 Q 部分变成了 R 部分。
6. 最近的 QRS 波群的 R 部分变成了 S 部分。
7. 最近的 QRS 波群在哪里结束。
8. 最近 QRS 波群的异常评级，其中一个极端是“不存在”。
9. 最近的 T 波开始的地方。
10. 最近的 T 波在哪里结束。
11. 最近 T 波的异常评级，其中一个极端为“不存在”。

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

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

Answer 11

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

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

-保罗

Answer 12

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

对我来说，这似乎是谐波分析在数学中的明确应用。可能有几个我可能遗漏的细微之处。

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

我可能在这里遗漏了一些东西......