使用 matlab 心电数据捕获心率

心电数据采集到了，如何计算心率成了一个问题，不过这个问题很早就有人解决了，但是网上能查找到的资料少之又少，代码也很少公开。这里博主就提供一个使用matlab程序计算心率的方法，只是提供了一个简单的思路和程序，此方法和程序里面的参数可能只适合博主自己做的心电器，未必适合其他心电数据，具体情况具体分析，欢迎大家来讨论。

1、数据导入！

第一步就是数据导入了，博主这里有一个转换好的txt数据文件，直接用load指令就能导入了。

data=load('ecg.txt');

然后可以使用plot指令查看一下波形

plot(data)

有可能这个数据是反的，这时候则需要添加一个“-”号，放大数据如下图所示：

2、数据预处理

这一步很关键，观察心电波形，R波峰很明显，通常也是使用探测R波峰的方式捕获心率的，但是上图所测量到的T波也比较大，如果不处理的话很可能会因为基线漂移造成捕获失败，程序将R波和T波混淆。

从信号成分分析，R波的信号成分在20Hz以上，T波成分一般在10Hz以下，这样我们就可以使用滤波器，将T波的信号衰减，突出R波的信号。同时使用高通滤波，去除基线偏移。这样我们就需要设计一个带通滤波器了，带宽在20Hz-30Hz，20Hz是为了去除基线和衰减T波信号，30Hz是保证心电信号保留的状态下，去除一些抖动噪声，特别是50hz的工频干扰。由于博主的数据采样率是500Hz，所以设计滤波器的时候需要按照500hz采样率来设计。这里使用的是32阶FIR滤波器，使用hamming窗。

Fs = sampleRate;           % Sampling Frequency
N    = 32;          % Order
Fc1  = 20;           % First Cutoff Frequency
Fc2  = 30;          % Second Cutoff Frequency
flag = 'scale';     % Sampling Flag
win = hamming(N+1); % Create the window vector for the design algorithm.
b  = fir1(N, [Fc1 Fc2]/(Fs/2), 'bandpass', win, flag); % Calculate the coefficients using the FIR1 function.
filteredA = filtfilt(b,1,data);

滤波效果如下，是不是基线很平稳，R波很突出，T波基本衰减很小，这样捕获R波就很方便了。

3、设置捕获条件！

接下来就是捕获R波，博主使用的是matlab自带的 findpeaks 函数，该函数是寻找极大值，同时可以限定条件，比如峰值阈值，极大值间隔等等。这样就很方便了，在一段时间内（窗口为1分钟），各个R峰值的阈值一般在信号的70%（这个值可以自己调试，看你心电数据了）以上，通过获取某段时间（窗口为1分钟）的最小值和最大值，来设置 'MinPeakHeight 捕获条件阈值。其次就是极大值间隔，一般正常人的心率不会超过120次每分钟，findpeaks 里面的 'MinPeakDistance' 可以设置最小间隔就是0.5s了。

获取某分钟的心电数据（500指的是采样率）：

time_watch =26;
RSF=filteredA(time_watch*60*500:(time_watch+1)*60*500);

幅度阈值设置（这里设置的是70%以上才可以捕获极大值）：

Threshold = (max(RSF) - min(RSF))*0.7 + min(RSF);

接着是间隔设置（这里设置的是0.5s的最小间隔）：

[R_pks,R_locs]=findpeaks(RSF,500,'MinPeakDistance',0.5,'MinPeakHeight',Threshold);

接着就是显示是否捕获到R峰了：

plot(RSF);
title(['第',num2str(time_watch),' 分钟 EEG'])
hold on
plot(R_locs*500,R_pks,'x')

运行效果如图，效果还是不错的：