在Matlab中使用DWT检测心电图QRS峰值和心率

QRS结合了典型心电图上的三个偏转（Q、R和S）。它对应于人类心脏的左右心室的去极化和大心室肌肉的收缩。

在数值和功能分析中，离散小波变换（DWT）是任何小波被离散采样的小波变换。

离散小波变换有许多工程、数学和计算机科学的应用。最值得注意的是，它被用于信号编码，以更加冗余的形式表示离散的信号，通常作为数据压缩的前提条件。

sym4 小波与QRS相似，适用于QRS检测。因此，这个过程可以帮助诊断各种心脏疾病。本教程将研究如何利用心电图数据库获得这些心电图的峰值和检测率。从这个方法，我们可以得到心率。

前提条件

要通过本教程，你需要。

安装[MATLAB]。
对[MATLAB]基础知识的正确理解。

QRS复合物

正如我们前面所说，它是典型心电图信号中三个偏转（Q、R和S）的组合。

image of qrs

其中：
P是第一个偏转。
Q是对基线的第一次负偏转。
R是对基线的最高正偏转。
S是对基线的第二次负偏转。

正常QRS的振幅为5至30毫米，持续时间为0.06至0.12秒。QRS复合体的宽度、振幅和形状有助于诊断室性心律失常、传导异常、心室肥大、心肌梗塞、电解质重排和其他疾病状态。

请注意，QRS复合体并不总是有所有三个QRS。它可以有各种形状，如下图所示。

qrs complex configuration

PhysioNet上的心电图数据库

在本教程中，我们使用来自MIT-BIH心律失常的信号，以及从PhysioNet下载的ECG-ID数据库。PhysioNet上的每个ECG信号都有以下三个文件。

*.atr:参考注释。
*.dat。数据文件（信号）。
*.hea: 头文件。

然而，Matlab不能读取这些文件，因此我们必须将我们的ECG转换为.mat 。为了做到这一点，我们使用PhysioNet ATM。ATM银行的界面如下所示。

ATM interface

你可以在输入中通过点击下拉箭头来选择你的数据库。注意，这里有所有的PhysioNet ecg数据库。

Show all the databases in the dropdown

您可以选择记录、信号、注释、输出长度、时间格式和数据格式，因为它们都有选项。当你到达工具箱部分时，你也可以选择你的选项，当你选择plot waveforms ，你会有如下所示的波形图。

waveform

因为我们需要在Matlab中读取它，所以我们导出它。要做到这一点，我们选择export signal as .mat ，然后在工具箱上下载它。

exporting

downloading

由于我们只需要信号，我们下载.mat 文件。

使用symlet4小波进行心电图信号分析

sym4 小波与QRS复合体相似。这就是为什么它是QRS检测的首选。为了明确这一点，看一下提取的QRS复合体和扩张的sym4 小波的图像并进行比较。

comparison

正如你所看到的，ecg的QRS复合物与sym4 小波的形状相当相似。这就是为什么sym4 小波总是被优先用于心电图信号分析的原因。

建议基于DWT的QRS检测

以下是基本的心电图信号，如果我们仔细观察，我们可以找到具有特定频率贡献的标记区域。

image of freq distribution

f1：代表高频噪声，有一定的频率f1。
f2：是QRS，有f2的频率贡献。
f3：心电图中缓慢变化的内容，有频率贡献f3。

这三个频率之间的关系将是f1>f2>f3 。

我们的目标是保留所有的R峰，消除所有其他的频率。为了清楚起见，我们说我们要消除f1 和f3 ，但保留f2 。

这就是所谓的带通滤波。你可以在小波变换的帮助下实现它。小波变换将相同频段的信号分组。因此，你可以通过消除一些频段来实现带通滤波。

这种带通滤波可以通过消除心电图信号的一些低尺度（高频）和高尺度（低频）的小波系数来实现。为此，使用未定额的小波变换来获得小波系数。

什么是未截断的小波变换？

那么，在一个正常的mra 小波中，转换信号在每一个分解级别后都被降频为2，通过这种方式，它的大小在每一个分解级别都会减少。

因此，在一个未减值的小波中，信号的长度保持不变。使用sym4 ，对一个心电图信号进行4级分解，如下图所示。

sym4 decomposition

第一个图是心电图信号。d's 是心电图信号每一级的详细系数。a4 是第四级的近似系数。

我们将通过去除系数a4 来获得带通滤波，因为它将不被考虑--类似地，我们去除d1 和d2 。

我们不考虑它的原因是，它是一个近似的系数。它承载了所有的低频细节。不考虑d1 和d2 ，因为它们包含了信号的高频细节。考虑d2 和d4 ，以重建或实现带通滤波的信号。

我们只考虑d3 和d4 ，并采取反小波变换的方式得到以下信号。

signal

在标准峰值检测算法的帮助下，我们可以定位这些R峰。同时，你可以找到给定时间间隔内的总R峰数，从而找到心率。

例如，假设我们有一个10秒的心电图信号，并且R-peaks的总数有一些数值，那么我们可以找到一分钟内的R-peaks的数量，代表每分钟的心跳，这就是心率。

用Matlab代码从心电图信号中获得QRS峰和心率

第一步是输入我们的信号。用户应该输入信号，所以Matlab应该要求输入信号。为了让Matlab允许用户从文件夹中选择信号，我们使用uigetfile 函数。这个函数考虑到了路径和文件名。

%program to get QRS peaks and heart rate from ecg signal

[filename, pathname]=uigetfile('*.*', 'select your ecg signal');
filewithpath = strcat(pathname, filename);

接下来，我们需要信号的采样频率。这些采样频率是在数据库中定义的。我们使用input 函数，因为用户定义了采样频率。这个函数会读取用户的输入。在这之后，使用load 函数加载数据。

Fs = input('Enter sampling rate:');
ecg = load(filename);   %reading ecg signal

之后，我们对振幅进行标准化处理。它是通过将ecg值除以增益来完成的。这个增益值也是在数据库中给出的。我们还使用函数length 来获得信号的长度，这个函数将信号作为输入。

这个长度有助于确定信号所需的时间。

ecgsig = (ecg.val)./200; %Normalize gain
t = 1:length(ecgsig);  %No. of samples
tx = t./Fs;  %Getting time vector

接下来，我们需要使用sym4 ，计算4级小波变换的未截断的小波。这可以确保信号的长度保持不变。为了计算这个，我们使用modwt 函数。这个函数取ecg signal 和sym4 4级。

wt = modwt(ecgsig,4, 'sym4');  %4-level undecimated DWT using sym4
wtrec = zeros(size(wt));

如前所述，我们的小波变换有5行，即 $a\_n、d\_4、d\_3、d\_2$ 和 $d\_1$ 。我们不需要近似的高频系数 $d\_1$ 和 $d\_2$ 。所以我们提取 $d\_3$ 和 $d\_4$ ，也就是第三和第四行。

我们使用要提取的行作为未减值的DWT的参数。

wtrec(3:4, :)= wt(3:4,:);   %extracting only d3 and d4 coefficients

通过进行反离散小波变换（IDWT），我们将得到一个只保留了r峰的信号。逆向小波变换在执行完小波变换后将信号恢复到原始形式。

在Matlab中，我们使用imodwt 函数来做IDWT，参数为信号与提取的部分wtrec 。

y=imodwt(wtrec, 'sym4');  %IDWT with only d3 and d4.
y=abs(y).^2;  %magnitude square

然后我们找到信号的平均值。在寻找信号的峰值时，它将被用作阈值。寻找平均值是通过使用mean 函数完成的。

avg = mean(y);   %getting average of y^2 as threshold

%Finding peaks
[Rpeaks, locs] = findpeaks(y,t, 'MinPeakHeight', 8*avg, 'MinPeakDistance', 50);

find peaks 是信号处理工具箱中的一个变量，用于寻找峰值。我们将最小的峰值距离作为，以避免在峰值相近的情况下出现错误检测。这可能是由于不适当的过滤造成的。，给出了R峰的位置。50 locs

让我们根据信号的长度来确定R峰的位置。它代表了拍子的数量。然后将这个节拍数转换为每分钟的节拍数。

nohb = length(locs);   %No. of beats
timelimit = length(ecgsig)/Fs;  %getting the time function of the signal
hbpermin = (nohb*60)/timelimit;   %Getting Beat per minute.
disp(strcat('HeartRate= ', num2str(hbpermin)))  %displaying the heartrate

绘制正常心电图信号与时间的对比图，这样我们就能看到其中的差别。

%displaying ecg signal and detected R-peaks
subplot(211)
plot(tx, ecgsig);

xlim([0, timelimit])
grid on
xlabel('seconds')
ylabel('ECG signal')

此外，将过滤后的信号与检测到的峰值一起绘制出来。

subplot(212)
plot(t,y)
grid on
xlim([0, length(ecgsig)]);
hold on
plot(locs, Rpeaks, 'ro')
xlabel('samples')
title(strcat('R peaks found and heartrate: ', num2str(hbpermin)))  % displayes the heartrate.

当我们执行我们的程序时，我们会有以下的输出。

output

结论

使用离散小波变换，心电图、QRS和心率检测变得更加容易。Matlab是小波分析的最佳软件。正如我们所看到的，这些变换已经完成并以内置的形式存在。因此，它使执行操作变得容易。

另外，Matlab除了有内置的变换形式外，还有其他的内置功能，帮助分析信号。这些函数包括length ，用于获取长度。此外，ecg信号的数据库与Matlab兼容，因为它提供了下载Matlab文件的选项。

如何在Matlab中使用DWT检测心电图QRS峰值和心率