心脏标测信号的频域相关分析方法与流程

本发明涉及一种心脏标测信号的频域相关（F相关）分析方法。

背景技术：

心脏标测技术是利用多电极点采集心脏电信号的一种技术手段，目前标测位点呈现 “大数量、高精度”的发展趋势，从而更细致地研究和分析心脏的电活动情况。通过分析各标测信号在心律失常期间的变化以及它们之间的关联与差异，研究者可以细致深入地理解心律失常的电生理机制，从中寻觅与心律失常的触发和维持机制紧密相连的关键因素。

节律分析是目前新兴的一种心脏标测信号分析方向，旨在研究标测信号所对应心肌细胞群参与心脏电活动的节律性（也即周期性）的状态。有关主导频峰（主导节律）的分析方法在节律分析中占有重要的位置。主导频峰以傅里叶变换为基础，计算标测信号的功率谱，并以功率谱中的最大峰值作为计算结果，反映心肌细胞群的主导节律。目前，主导频峰（Dominant frequency，简写DF）在寻找房颤异位兴奋源和指导射频消融上面已经发挥着一定的用途。但当信号的规整程度较差时，心脏内可能存在多个异位兴奋源，主导频峰的结果的准确度不高。而且，主导频峰仅能反映某处信号的主导节律，不能反映两处不同位置信号节律之间的关系或者同一位置不同时间段信号之间的关联。

相关分析是基于相关函数的一种信号相似度的分析方法，它描述两路信号在不同时刻的取值之间的相互关联程度，能够度量两个信号之间波形的相似程度以及与相对时间位移之间的关系。心脏不同部位处标测的信号波形各有不同，如果直接研究不同部位处信号的相关分析结果，那么易受到波形因素的影响。而且，相关分析不能反映两路信号频域上的关联程度。尽管相干分析可以用于分析两路信号在频域上的关联程度，但是其应用的局限性较大（对于不同起源位置的信号，各处频率相干系数均较低），很难用于分析多起源点的复杂心律失常信号（如房颤）。

技术实现要素：

为了克服普通相关分析只能反映时域波形的关联程度，不能反映频域节律之间的关系；同时为了克服频域相干分析用于复杂心律失常下的心电信号分析的局限性，本发明的目的在于提出一种全新的用于评价节律相关程度的频域相关分析方法——F相关。此方法以节律分析为基础，将相关分析从时域延伸至频域，研究两路或多路信号间节律的相似程度，可用于分析房颤等复杂心律失常下的心电信号，提供潜在的不易察觉的电生理信息，以供科学研究参考。

本发明提出的心脏标测信号的频域相关分析方法（F相关），通过一系列信号预处理步骤，最终得到名为“最大相关系数”和“最大相关频率”两个评价指标，用于评价两路节律信号在频域上的关联程度。具体步骤如下：

(1)利用心脏标测技术采集心电信号，所述心电信号包括体表标测信号、心内膜标测信号与心外膜标测信号；

(2)将步骤（1）采集到的心电信号通过带通滤波器进行带通滤波，控制带通滤波器的通带频率范围为20-100Hz；

(3)对步骤（2）得到的经过带通滤波的心电信号取绝对值；

(4)计算步骤（3）得到的心电信号的功率谱密度G(f) ，得到频谱波形，所述频谱波形包括直流成分、主导频峰和谐波成分；

(5)剔除步骤（4）得到的直流成分，只保留步骤（4）得到的功率谱密度G(f) 中频率范围为0.2-20Hz的部分，得到频率范围为0.2-20Hz的低频段功率谱密度G'(f) ，所述低频段功率谱密度中含有主导频峰和有限个谐波成分的频率；

(6)任意选取i路和j路两路心电信号，通过步骤（1）至（5）得到步骤（5）中定义的i路、j路的频率范围为0.2-20Hz的低频段功率谱密度G_i'(f) 和G_j'(f) ，计算G_i'(f) 和G_j'(f)的互相关，公式如下：

其中：ω 和f 都是频率自变量，FC 代表频域互相关，FC_ij(f)表示i路、j路的频率范围为0.2-20Hz的低频段功率谱密度G_i'(f) 和G_j'(f) 的频域互相关；

(7)通过步骤（6）的互相关计算得到最大相关系数F-C_max以及最大相关频率F-f_max，分析和研究步骤（6）中的互相关结果；

最大相关系数F-C_max：反映了两路心电信号预处理后0.2-20Hz的频谱波形的最大相似程度，在0.2-20Hz内的频谱波形中，主导频峰（DF）是最大的峰（若主导频峰<10Hz还会有主导频峰的谐波成分），F-C_max则正好反映的就是由这两个DF主导的波形的最大相似程度，当两路心电信号的节律越相似，F-C_max则越大（越接近1）；

最大相关频率F-f_max：主要反映了两个主导频峰之间的频率差ΔDF。

本发明提出的方法，首先是基于节律（特别是主导节律）信号，其次是在频域上再进行相关分析。上述步骤中，步骤(2)和(3)是为了更好地得到信号的主导节律以及谐波成分。

步骤(4)中所得信号的功率谱密度，其中最大值是0Hz（步骤(3)取绝对值后直流成分增大），次大的就是主导频峰（DF）。主导频峰DF反映了信号的主导节律。

步骤(5)中G'(f)的峰值就是主导频峰，它剔除了直流成分的影响，并且包含了主导频峰和有限个谐波成分。

步骤(6)中，将G_i'(f) 和G_j'(f) 当作两路新的信号来处理，计算它们的互相关，得出的互相关函数是以频率f 作为自变量，而不再是以时间t 作为自变量。由于G'(f) 的波形主要由主导频峰和其谐波决定，所以这种互相关能够反映节律之间的关系，特别是两路信号的主导频峰之间的关联。

本发明中提出的F相关方法，主要有两个分析指标：最大相关系数，记作F-C_max；以及该最大系数处的频率——最大相关频率F-f_max。

最大相关系数F-C_max：反映了两个信号预处理后0.2-20Hz的频谱波形的最大相似程度。在0.2-20Hz内的频谱中，主导频峰DF是最大的峰，若DF<10Hz还会有DF的谐波成分，所以0.2-20Hz的频谱形状跟DF有着密切的关联。而F-C_max则正好反映的就是由这两个DF主导的波形的最大相似程度，两路信号的节律越相似，F-C_max则越大（越接近1）。

最大相关频率F-f_max：主要反映了两个主导频峰之间的频率差ΔDF。由于0.2-20Hz内的频谱波形在很大程度上由主导频峰决定，所以，最大相关频率F-f_max的值能反映这两个最大峰（DF）的频率之差。

本发明具有以下有益效果：

1. 本发明的F相关分析方法，避免了不同标测部位波形不一致对结果的影响，是将相关分析从时域延伸至频域的一种方法。同时它克服了频域的相干分析难以用于复杂心律失常信号分析的问题。

2. F相关反映了两处信号的节律关联程度，能够用于时间维度——同一标测位点处前后不同时间段的两段信号之间的节律关联性。基于时间维度的应用，可以用于窦性心律监测、易颤性分析等等。

3. 本发明的F相关分析方法能够用于空间维度，即研究同一时间不同的标测位点处信号之间的节律关联性。基于空间的应用，结合伪彩图，可以研究不同部位节律的关联程度。

4. 本方法的应用范围可推广至所有的电生理节律信号的量化研究，如脑电和肌电等。

附图说明

图1是心房标测电信号的示意图（窦性心律信号（a）和房颤信号（b））；

图2是窦性心律信号的预处理和频谱示意图；其中：（a）为某一标测位点处的一段5s时长窦性心律信号，（b）为信号（a）的频谱，（c）为信号（a）通过20-100Hz的带通滤波器后得到的波形，（d）为（c）的频谱，（e）为（c）取绝对值得到的波形，（f）为（e）的频谱。

图3是窦性心律情况下，同一时刻不同标测位置信号的相关分析、相干分析以及F相关示意图；其中：（a）为某一通道的一段5s时长窦性心律信号，（b）为同一时刻另一通道的一段5s时长窦性心律信号，（c）为将信号（a）和（b）进行传统的相关运算得到的波形，（d）为将信号（a）和（b）进行传统的相干运算得到的波形，（e）为将信号（a）和（b）按照本发明方法的步骤（1）至步骤（4）的流程（依次进行20-100Hz的带通滤波器、取绝对值、计算功率谱密度）得到的频谱波形，（f）为将（e）中的两路结果按照本发明方法的步骤（5）和步骤（6）的流程得到的F相关结果。

图4是窦性心律情况下，同一标测位置不同时刻信号的相关分析、相干分析以及F相关示意图；其中：（a）为某一通道的一段5s时长窦性心律信号，（b）为同一通道另一时刻的一段5s时长窦性心律信号，（c）为将信号（a）和（b）进行传统的相关运算得到的波形，（d）为将信号（a）和（b）进行传统的相干运算得到的波形，（e）为将信号（a）和（b）按照本发明方法的步骤（1）至步骤（4）的流程（依次进行20-100Hz的带通滤波器、取绝对值、计算功率谱密度）得到的频谱波形，（f）为将（e）中的两路结果按照本发明方法的步骤（5）和步骤（6）的流程得到的F相关结果。

图5是房颤情况下，同一时刻不同标测位置信号的相关分析、相干分析以及F相关示意图；其中：（a）为某一通道的一段5s时长房颤信号，（b）为同一时刻另一通道的一段5s时长房颤信号，（c）为将信号（a）和（b）进行传统的相关运算得到的波形，（d）为将信号（a）和（b）进行传统的相干运算得到的波形，（e）为将信号（a）和（b）按照本发明方法的步骤（1）至步骤（4）的流程（依次进行20-100Hz的带通滤波器、取绝对值、计算功率谱密度）得到的频谱波形，（f）为将（e）中的两路结果按照本发明方法的步骤（5）和步骤（6）的流程得到的F相关结果。

图6是房颤情况下，同一标测位置不同时刻信号的相关分析、相干分析以及F相关示意图；其中：（a）为某一通道的一段5s时长房颤信号，（b）为同一通道另一时刻的一段5s时长房颤信号，（c）为将信号（a）和（b）进行传统的相关运算得到的波形，（d）为将信号（a）和（b）进行传统的相干运算得到的波形，（e）为将信号（a）和（b）按照本发明方法的步骤（1）至步骤（4）的流程（依次进行20-100Hz的带通滤波器、取绝对值、计算功率谱密度）得到的频谱波形，（f）为将（e）中的两路结果按照本发明方法的步骤（5）和步骤（6）的流程得到的F相关结果。

图7是F相关的两个指标F-C_max和F-f_max示意图。

图8是F相关在时间维度上的应用举例示意图。

图9是F相关在空间维度上的应用举例示意图；（a）为最大相关系数F-C_max矩阵的图像，（b）为最大相关频率F-f_max矩阵的图像。

具体实施方式

下面结合附图和实施例对本发明方法及其应用做进一步说明。

实施例1：将本发明的F相关方法分别用于窦性心律标测信号（不同部位节律一致）和房颤信号（存在多个节律，不同部位节律可能不一致），并与传统的相关分析和相干分析进行比较。本实施例采用心外膜标测信号，工作流程如下：

（1）.利用心外膜标测技术采集犬的心房外膜标测信号（数字信号）。如图1所示，图1(a)是某一标测位点处的窦性心律信号，图1(b)是另一标测位点的房颤信号。图1(a)与图1(b)不是同一时间采集的。

（2）. 将步骤（1）采集到的窦性心律信号进行20-100Hz的带通滤波。如图2所示，图2(a)是某一标测位点处的一段5s时长窦性心律信号，图2(b)是其频谱；将图2(a)中信号通过20-100Hz的带通滤波器后，得到图2(c)所示的信号，其频谱如图2(d)所示。

（3）. 将带通滤波后的信号（图2(c)）取绝对值，得到如图2(e)所示的信号。

（4）. 计算步骤（3）得到的信号的功率谱密度，得到如图2(f)所示的频谱波形。本实例是采用MATLAB软件中的快速傅立叶变换（FFT）来估计信号的功率谱。为了方便观察，图2(f)中只显示了0-40Hz的频谱波形，其中直流成分最大，次大的是主导频峰（椭圆圈出），谐波成分也很明显。在进行后面的F相关分析时只取类似于图2(f)中的0.2-20Hz这一段频谱波形。

（5）.取两路信号按照步骤（1）至步骤（4）的流程得到各自在0.2-20Hz范围内的频谱波形。如图3所示，图3(a)和图3(b)分别是同一时刻不同通道（不同标测位置）的窦性心律信号。对图3(a)和(b)所示的信号按照步骤（1）至步骤（4）的流程依次进行20-100Hz的带通滤波器、取绝对值、计算功率谱密度，得到图3(e)所示的频谱波形。图3 (e)中的每路波形的最大峰就是反映该信号主导节律的主导频峰（已剔除直流成分）。

（6）. 将步骤（5）得到的两路信号0.2-20Hz范围内的功率谱密度当作两路新的信号，即：将图3 (e)中实线信号的0.2-20Hz范围内低频段功率谱密度记作G_i'(f)；将图3 (e)中虚线信号的0.2-20Hz范围内低频段功率谱密度记作G_j'(f) ，计算G_i'(f) 和G_j'(f) 的互相关（F相关），结果如图3(f)所示。计算公式如下：

其中ω 和f 都是频率自变量，FC_ij(f) （这里，FC 即代表频域互相关）表示G_i'(f) 和G_j'(f) 的互相关。

（7）. 分析和研究步骤（6）得到的互相关结果，两个分析指标为最大相关系数F-C_max，以及取得该最大系数处的频率——最大相关频率F-f_max。图3(f)中F相关的最大相关系数F-C_max是接近1，表明两路信号具有很大的关联性，并且是在F-f_max=0Hz处取得最大值，这表明二者的节律（主要指主导频峰DF）是一致的（即节律差异为0Hz）。

下面将本发明的方法与传统的相关分析和相干分析进行比较。

对图3(a)和(b)的两路信号（同一时刻不同标测位置）分别做传统的相关分析和相干分析，结果如图3中的图3(c)和图3(d)所示。可以发现，相关分析的系数较低（参见图3(c)）。理论上，二者起源于同一起搏点窦房结，相关系数应该较高，但是由于采集部位不同导致的信号波形形态差异，使得传统的相关分析不能很好地反映这一电生理现象。由于这两路信号起源于同一起搏点窦房结，因此传统的相干系数在0-200Hz内相对较高（参见图3(d)）。

图4则是同一通道不同时刻的两路窦性心律信号的分析结果，计算方法同图3的介绍。这里，两路信号的节律理论上是一致的，波形也类似。但是由于每个激动波的细节部分会有差异导致相关系数不是很高，相干系数也不高，容易得到两路信号节律不一致的判断结果。但是，此时F相关的系数依旧是接近1，并且在0Hz处取得最大值。结果表明，在这种情况下采用F相关分析方法可得到二者节律一致的结果，这一结果与临床电生理的共识是一致的。

图5所示是同一时刻不同通道的房颤信号，由于是不同位点标测所得信号，两路信号的波形和主导节律都不一样，信号x 的节律是5s激动50次，即10Hz，信号y 的节律是5s激动42次，即8.4Hz。此时传统相关分析的相关系数已经很低（参见图5(c)），最大值为0.13，相干系数也低且无规律（参见图5(d)）。但对于本发明的F相关分析而言，在10-8.4=1.6Hz处取得了F相关的最大相关系数（参见图5(f)），正确地反映了两路信号的节律差异。

图6所示，则是同一通道不同时刻的房颤信号，二者节律一致，波形也类似（相似程度比窦性心律时差），传统相关系数依然很低（参见图6(c)），而F相关则在0Hz附近处取得F相关的最大相关系数0.83（图6(f)），正确地反映了二者的节律关系。

通过对图3至图6的两类窦性心律信号和房颤信号的分析，可以知道，传统的相关分析和相干分析，所得的关联系数较低（尤其是房颤时这两个系数下降非常多），难以判断信号间的节律关联程度。而对于本发明的F相关分析方法而言，在窦性心律和房颤时均能够很好地反映两路房颤信号的节律差异情况（最大相关频率F-f_max）及其关联程度（最大相关系数F-C_max）。所以，本发明在分析心脏标测信号的节律上面有着比较明显的作用，这个作用是传统的相关分析和相干分析所不具有的。

图7是用于解释本发明的两个指标F-C_max和F-f_max的示意图，它是图6(f)局部放大的结果。

实施例2：图8是将本发明的F相关（计算步骤同实施例1）用于128路标测信号时间维度上的实施例。在该实例中，同时采集心房外表面128点的电信号，并以一定的采样频率连续采集。将t₁时刻开始采集的信号与t₁+Δt时刻开始采集的信号的每一个对应通道实施本发明的F相关（如取Δt=10s，信号时长为5s，就是将1号通道t₁时刻开始采集的一段5s时长的信号与1号通道t₁+10s时刻开始采集的另一段5s时长的信号实施F相关；依次实施到所有128路信号），得到F相关的两个指标F-C_max和F-f_max。图8是取Δt=10s时的结果示意图，采集信号为房颤发生前的窦性心律信号。图8中的最大F相关系数F-C_max以实线表示，最大F相关系数处的频率F-f_max用虚线表示。F-f_max在某些通道处的值为0，而在另一些通道不为0，表明仅仅部分通道信号的节律发生了改变，而不是整体一致变动（不同步），同时F-C_max的整体值不高也提示前后不同时间段信号的节奏可能发生了紊乱。这与正常窦性心律时F-C_max普遍处于0.9-1.0之间，且F-f_max从1号到128号通道都保持同一个稳定值有着较大区别。这里，图8提示着心房信号的节律正发生着不正常的变化。对比分析连续多幅的时间维度结果，可以有助于追踪和捕捉异常心律失常的发生和变化。

实施例3：图9是将本发明的F相关（计算步骤同实施例1）用于128路标测信号空间维度上的实施例。该实例取自房颤发生时的128路同步采集的心房心外膜信号，时长5s。将每两路信号均实施本发明的F相关，并分别得出相应的两个指标F-C_max和F-f_max。这样，F-C_max和F-f_max的结果在两两F相关之后形成128*128的矩阵，将矩阵用图像表示，如图9所示。对比分析多幅不同时刻的图像或由多幅连续图像构成的视频，可以有助于分析和研究心律失常的发生和维持机制。