使用小波变换和 ST 段检测的 ECG 特征提取（Matlab）

引言

有关ECG信号和导联的理解，请参考http://en.wikipedia.org/wiki/Electrocardiography。传统上，这类ECG信号由ECG采集设备采集，这些设备会生成导联输出的打印件。心脏病专家分析数据以检查信号的异常或正常。但近年来，自动ECG处理受到了极大的关注。https://codeproject.org.cn/KB/cpp/ecg_dsp.aspx提供了一个关于通过廉价硬件将ECG信号采集和滤波到PC的精彩概述。您可以从http://www.physionet.org/physiobank/database/mitdb/下载各种疾病的ECG信号样本。现在，主要的问题是如何开发一个系统来提取ECG信号的特征，以便这些特征可以用于自动疾病诊断。在本文中，我们将讨论一种从ECG信号中提取特征的技术，并进一步分析ST段的抬高和压低，这是心肌缺血的症状。您可以从http://circ.ahajournals.org/content/110/17/2721.full了解更多关于心血管异常及其与ECG峰值相关性的信息。

背景

我们将详细讨论处理来自MIT-BIH数据库且格式为.mat的ECG信号的算法。对于当前的分析，我们考虑正常窦性心律和ST段抬高信号。最后，使用阈值检查信号的正常性。

Using the Code

首先选择一个.mat格式的文件名并加载文件。

[fname path]=uigetfile('*.mat');
fname=strcat(path,fname);
load(fname );

在信号前后附加100个零，以消除窗口在查找峰值位置时跨越信号边界的可能性。

z=zeros(100,1);

A=[z;A;z];

执行小波分解。小波分解的过程会对信号进行降采样。这本质上意味着以远低于原始信号的频率采样。因此，细节被减少，QRS复合体被保留。

[c,l]=wavedec(s,4,'db4');

提取变换后的系数

ca1=appcoef(c,l,'db4',1);

ca2=appcoef(c,l,'db4',2);

ca3=appcoef(c,l,'db4',3);

ca4=appcoef(c,l,'db4',4);

如果绘制系数图，您会注意到频带被分离，ca1、ca2、ca3和ca4是更干净的信号。但由于降采样，它们的样本数量会少于实际信号。您可以看到第一个信号类似于实际信号，但样本数量恰好是原始信号的四分之一，因为信号被分解了4个级别。第2级样本数量是第1级的正好一半，第3级样本数量是第2级的正好一半。由于样本数量的减少，这类信号也称为降采样信号。

很明显，第2级分解的数据是无噪声的。因此，我们将此信号视为理想的ECG信号，从中必须检测到QRS。但第一个R峰位于第3级分解信号中，大约在第40个样本处，而在原始信号中，它位于第260个位置。因此，一旦在第3级重建信号中检测到R峰，就必须在实际信号中进行交叉验证。

在降采样信号中检测R峰

首先找到大于原始信号最大值60%的值。这些值几乎都是R峰。由于分解后的信号是无噪声信号，首先需要在无噪声信号中检测到第一个R峰。但请记住，最终目标是在原始信号中检测峰值。原始信号中的样本值将与分解后的信号不同。因此，我们的策略是先在降采样信号中检测R峰，然后交叉验证原始信号中的这些点。

令y1为分解后的信号。

m1=max(y1)*.60;

P=find(y1>=m1);

因此，P现在是满足上述条件的点的集合。如果您仔细观察信号，R峰不是一个单一的脉冲峰，因此同一峰值中可能有多个点满足条件。需要记住的一点是，在500Hz采样信号中，P到R的位置不会低于350个样本。在第4级分解顺序中，这个值约为20。所以我们首先会移除距离太近的R位置. 

P1=P; 

P2=[]; 

last=P1(1);

P2=[P2 last];

for(i=2:1:length(P1))

    if(P1(i)>(last+10))

% In this step we find R peaks which are atleast 10 samples apart

        last=P1(i);

        P2=[P2 last];

    end

end

现在变量P2代表降采样信号中R峰的位置. 

在原始信号中检测R峰 

搜索信号y1中大于该值m1的所有位置。它们是R位置。但在我们继续之前，您必须知道，Rt中的R位置至少是相同点的实际R位置的1/4。因此，我们将首先通过乘以4来映射检测到的位置到原始信号. 

P3=P2*4;

%Multiply the current location with 4 to get the actual scale.

您必须记住的另一个重要事情是，降采样信号中的R位置永远不会在原始信号的4倍尺度上。降采样过程总是会偏离信号的位置。因此，我们需要在从参考R点（P3获得）的+-20个样本的窗口内搜索原始信号中的最大值。

Rloc=[];

for( i=1:1:length(P3))

  range= [P3(i)-20:P3(i)+20]

%Search within a window of +-20 samples in the original signal with reference to up scaled  R locations detected in downsampled signal.

    m=max(A(range));

    l=find(A(range)==m);

    pos=range(l);

    Rloc=[Rloc pos];

end

Ramp=A(Rloc);

现在很清楚，Ramp和Rloc代表了原始尺度的R峰幅度和位置。让我们看看在波形中对它们的标记。

 

参照R峰检测其他峰

从R峰向前和向后遍历，搜索最小值和最大值，它们分别是P、Q、T、S峰。所以循环遍历Rloc并搜索其他峰值. 

首先，如果您观察波形，会非常清楚，从R位置开始，如果您选择一个Rloc-100到Rloc-50的窗口并找到最大值，那么这个最大值就是P峰. 

X=Rloc;

y1=A;

for(i=1:1:1) % If you have a 12 lead data than, for(i=1:1:12)

for(j=1:1:length(X))

a=Rloc(i,j)-100:Rloc(i,j)-10;

    m=max(y1(a));

    b=find(y1(a)==m);

    b=b(1);

    b=a(b);

    Ploc(i,j)=b;

    Pamp(i,j)=m;

%The minima in the Window of   Rloc-100 to Rloc-10  is essentially the Q peak. 

    %% Q  Detection

    a=Rloc(i,j)-50:Rloc(i,j)-10;

    m=min(y1(a));

    b=find(y1(a)==m);

    b=b(1);

    b=a(b);

    Qloc(i,j)=b;

    Qamp(i,j)=m;

%With similar logic you  can detect the S and T peaks.

    %% S  Detection

    a=Rloc(i,j)+5:Rloc(i,j)+50;

    m=min(y1(a));

    b=find(y1(a)==m);

    b=b(1);

    b=a(b);

    Sloc(i,j)=b;

    Samp(i,j)=m;

    %% T Peak

    a=Rloc(i,j)+25:Rloc(i,j)+100;

    m=max(y1(a));

    b=find(y1(a)==m);

    b=b(1);

    b=a(b);

    Tloc(i,j)=b;

    Tamp(i,j)=m;

end

end

一旦所有峰值都正确检测到，您就可以找到每个导联的零交叉点作为起始点和结束点。ST段可以从S-Offset和T-Onset计算得出. 

兴趣点 

您也可以考虑在使用Symlet或任何其他滤波技术之前清理ECG信号.