心电图信号的频域分析

摘要

提出了一种简便有效的、由心电图（ECG）曲线提取QRS-T复合波的新方法，从而获得P波（当房颤未发生时）或f波（当房颤时）；利用功率谱、自相关、互相关等信号处理技术，对它进行了处理和分析，特别是，比较了未房颤时和发生房颤时的相关信息的显著差别。

正文

《1 前言》

1 前言

心电图 (ECG) 是常用临床检查方法之一, 对一些疾病, 特别是对心血管病的诊断具有重要价值。 ECG信号自动分析的关键点在于QRS复合波的检测, 随之便可进行心率变异 (HVR) 分析, 以及对ECG信息进行全面研究。

在ECG信号处理中, 用于QRS复合波检测的早期方法有带通滤波、自适应滤波、非线性自适应滤波等^[1,2,3,4]。由于这类滤波方法的算法复杂, 且运算量很大, 其实际应用 (特别是临床应用) 受到限制。近年来, 随着小波技术的发展, 涌现出了基于小波变换的QRS复合波检测法^[1,5,6,7,8]。然而, 应用这类方法时需要根据所测信号的特性来选择一个合适的母波, 以期在时域以及频域内都能显示出该信号的全部特性, 而且还要选择合适的尺度, 以期小波变换后的信号突变点能够反映原信号的突变情况^[9,10], 否则便可能出现漏检 (某些QRS复合波未被检测到) 或误检 (检测到的不是QRS复合波) 。

为克服以往方法的上述弱点, 笔者提出如下的一种简便有效的、基于频数直方图的QRS-T复合波检测法。

《2 ECG曲线的P波与f波》

2 ECG曲线的P波与f波

ECG 的时域波形一般由P波、QRS 复合波、T波等波段组成, 而房颤发生时, 其中的P波则演变为f波 (房颤波) 。将 ECG 图像中的 QRS-T 复合波剔除, 留下的便是代表未房颤时的P波或房颤时的f波。

现提出一种简便有效的、基于频数直方图的 QRS-T 复合波剔除方法。

《2.1快速搜索ECG图像中的QRS-T复合波所在的区间》

2.1快速搜索ECG图像中的QRS-T复合波所在的区间

频数是指单位宽度的条块中所含的图形像素数, 频数直方图反映出图像的统计特性, 描述了图形像素在图像中的分布规律, 因此, 可以利用频数直方图所含有的信息作为图像的特征, 从而快速确定 QRS-T 复合波所在区间。ECG 的时域曲线如图1a 所示。

《图1》

图1 ECG的时域波形及频数直方图

Fig.1 ECG time domain waveform and its histogram

令:X为ECG图像矩阵 (已进行二值化处理:‘0’为白色, ‘1’为黑色) , d (j) 为每列的像素数, D为d (j) 的集合, i, j分别为像素点的行坐标、列坐标。于是,

$d (j) = \sum_{i} x (i, j), D = {d (j)} 。 (1)$ $d (j) = \sum_{i} x (i, j), D = {d (j)} 。 (1)$

i, j 的取值范围取决于所研究的ECG图像的尺寸。对于一帧图像, 一般可以取 i = 1, 2, …, 600;j = 1, 2, …, 1 200。

图1b是集合D 的图形表示, 即频数直方图, 而图1c 是图1b 的阈值化。在这里, 选择适当的阈值d₀, 派生出阈值化后的新集合

$D_{t h} = {d_{t h} (j)} 。 (2)$ $D_{t h} = {d_{t h} (j)} 。 (2)$

式 (2) 中:当 d (j) ≥ d₀ 时, d_th (j) = d (j) ;当 d (j) <d₀ 时, d_th (j) = 0。

从图1 的各图对应关系中, 可以看到, 图1c 的各个相对极大值对应着相关的R波波峰, 再令n (k) 为第k个R波波峰所在位置的列坐标, 则这类列坐标组成的集合为

$Ν = {n (k)} = {n (1), n (2), \dots n (Ν)} 。 (3)$ $Ν = {n (k)} = {n (1), n (2), \dots n (Ν)} 。 (3)$

《2.2QRS-T 复合波的检测》

2.2QRS-T 复合波的检测

设所选图像的尺寸为 M×N (行×列) , QRS-T 图形模板是M₁×N₁ 图形矩阵Q_qrst, 模板特征向量是D_qrst, 则只需依次在图1a 中、于式 (3) 所示的各列坐标附近用该模板稍事搜索, 即能找出相应的 QRS-T 复合波的位置。

搜索时所选用的判断依据是搜索区域内的差异点数e_k及图像特征向量D_k, 以及数学形态相似度。

差异点是指正被搜索区间的图像矩阵 Q_k 和模板图形矩阵Q_qrst之间的、存在着不匹配的那些点, 而差异点数 e_k 为

$e_{k} = \sum (Q_{k} x o r Q_{q r s t}) 。 (4)$ $e_{k} = \sum (Q_{k} x o r Q_{q r s t}) 。 (4)$

式 (4) 中的 xor 为逻辑运算中的‘异或’操作, 该操作是在图像矩阵Q_k 和Q_qrst中的相应矩阵元之间进行的, 而求和操作∑遍历每次搜索的小区间是对 xor 操作后形成的新矩阵的元求和。上述特征向量D_k由所搜索小区域内的各列像素数组成。

匹配规则是:在所搜索的小区域内, 符合下列式的图像即为要提取的QRS-T复合波。

$\begin{array}{l} (1 - c_{p}) D_{q r s t} \leq D_{k} \leq (1 + c_{p}) D_{q r s t}, c_{p} \in (0, 1); (5) \\ e = \min {e_{k}}_{k = 1, 2, \dots} 。 (6) \end{array}$ $\begin{array}{l} (1 - c_{p}) D_{q r s t} \leq D_{k} \leq (1 + c_{p}) D_{q r s t}, c_{p} \in (0, 1); (5) \\ e = \min {e_{k}}_{k = 1, 2, \dots} 。 (6) \end{array}$

作为例子, 讨论如何利用上述方法从ECG曲线获取其中的P波和f波。

《2.3P 波和 f 波的获取》

2.3P 波和 f 波的获取

按上述步骤, 将ECG图像中的QRS-T复合波剔除, 即可得到未房颤时的P波 (见图2a) , 或房颤时的f波 (见图2b) 。

由图2可见, P波为能量相对集中的一个小波包, 而f波则是能量分散的一组不规则振荡。

《图2》

图2 P 波和 f 波的图像

Fig.2 The curves of P wave and f wave

《3 ECG信号的处理和分析》

3 ECG信号的处理和分析

所采用的数据主要来自美国MIT-BIH数据库, 有些则是北京工业大学校医院的临床病例。其中, 发生房颤的和未房颤的各50例。

对这些数据, 从功率谱、自相关和互相关等方面进行了处理和分析。为了便于观察, 对这些函数都进行了归一化处理。

《3.1功率谱函数》

3.1功率谱函数

对所采集的ECG数据进行FFT运算, 得到其幅频特性, 据此, 即可绘出ECG的功率谱图, 见图3。

从图3可以明显看出:未房颤时功率谱的谱线分布均匀, 波谷都能回到零线;而房颤时, 谱线疏密不等, 有些谱线的波谷不能回零线。

《3.2动态频谱》

3.2动态频谱

取ECG信号每10 s为一个时段, 每个样本依次取8个时段, 剔除 QRS-T 复合波后, 留下的便是P波 (未房颤时) 或f波 (发生房颤时) , 再分别对各样本进行FFT运算, 作出其动态频谱, 如图4 所示。

对比图4a 和图4b可以看出, 房颤时和未房颤时的动态频谱有着明显的不同。房颤时f波的能量分布很集中, 浓聚在 0～5 Hz 间的一段很窄频率范围内;而未房颤时P波的能量则分布于较宽的频率范围。而且, f 波各时段的动态频谱曲线之间的差异性较大, 说明了房颤波的不规则性。

《图3》

图3 功率谱图

Fig.3 The power spectrum

《图4》

图4 动态频谱图

Fig.4 The dynamic frequency spectrum

《图5》

图5 自相关函数的图示

Fig.5 The auto-correlation function

《3.3自相关函数》

3.3自相关函数

自相关函数用于描述一个信号在不同时刻取值的相关程度, 是对信号内在联系的一种度量, 反映了该信号自身的内在联系。

令ECG信号的采样序列为 x (n) , 则其自相关函数为:

$R_{x x} (m) = \lim_{Ν \to \infty} \frac{1}{2 Ν} \sum_{n = - Ν}^{Ν} x (n) x (n - m) 。 (7)$ $R_{x x} (m) = \lim_{Ν \to \infty} \frac{1}{2 Ν} \sum_{n = - Ν}^{Ν} x (n) x (n - m) 。 (7)$

据此, 可以绘出自相关函数的示图。

图5a 和图5b分别是未房颤时和发生房颤时的自相关曲线。

比较图5a 和图5b可知, 房颤时的自相关曲线幅值急剧下降, 第一高峰R_max1 与第二高峰R_max2 的比值R_max1R_max2 ≥4.5, 而未房颤时自相关曲线幅值下降相对较慢, R_max1R_max2 ≤3。

《3.4互相关函数》

3.4互相关函数

互相关函数反映出2个信号在不同时刻取值的关联程度, 即两个信号之间的相互依赖关系。

以QRS-T的波形模板信号的采样序列作为 y (n) , ECG信号的采样序列作为 x (n) , 则互相关函数 R_xy (m) 为

$R_{x y} (m) = \lim_{Ν \to \infty} \frac{1}{2 Ν} \sum_{n = - Ν}^{Ν} x (n) y (n - m) 。 (8)$ $R_{x y} (m) = \lim_{Ν \to \infty} \frac{1}{2 Ν} \sum_{n = - Ν}^{Ν} x (n) y (n - m) 。 (8)$

据此, 可以绘出互相关函数的示图 (见图6) 。

《图6》

图6 互相关曲线

Fig.6 The cross-correlation function

从图6可以看到, ECG 信号与 QRS-T 模板之间的相关程度在各 QRS-T 复合波出现处为最大, 互相关曲线上出现了与 QRS-T 曲线近似的波形。未房颤时, QRS-T 复合波的间隔较为均匀, 各QRS-T 复合波间的波形亦较为规则。房颤发生时, QRS-T 复合波的间隔大小不等, 且于各 QRS-T 复合波间出现不规则振荡。

《4 结语》

4 结语

笔者提出了一种简便有效的、从 ECG 图像剔除 QRS-T 复合波的新方法, 从而方便地获取到未房颤时的P波或房颤时的f波。根据功率谱函数、自相关函数、互相关函数等方面的处理, 对一般难以直接识别的心电图上微小差别, 都可在频域上清楚地观察到, 从而很容易地判断是否发生了房颤。

展示更多