【生物医学信号处理及其MATLAB应用】Chapter 5 初步应用：ECG信号的时间特征检测和提取

5.1 ECG信号产生原理

​ 心电信号产生原理：心脏在收缩之前先有电激动，约在0.02~0.07s后才有机械收缩。心脏电激动产生动作电流。人体是良好的容积导体，心脏动作电流可以传导身体各部，在两个体表部位放置点击，连接心电图机，获得的心电活动曲线，即心电图（ECG）。下图为典型的心电信号，包含R/P/T等波，这些波段的潜伏期和幅度能一定程度反映心脏生理特征。因此首先需要能检测和提取出这些事件。

5.2 心电R波尖峰的检测

5.2.1 阈值法检测心电R波尖峰

1. 对信号进行扫描，找到其中的峰值。

2. 取一个阈值，阈值的设定可以很灵活。一般与最大值，平均值有关，具体应视情况而定，通过试错获得。

3. 对于所有大于阈值的峰值点作为检测到的R波尖峰。

4. 由生理基础可以知道，R波间隔是相对稳定的。可以通过检测峰值点的间隔，去除那些较高或较小的间隔伪迹。

clear;clc;close all;
% 读入数据
fnam = input('Enter the ECG file name:', 's');
fid = fopen(fnam);
ecg = fscanf(fid, '%f');
fclose('all');

% 对数据进行加 0 处理，防止 R 波出现在数据最前或最后
g = [zeros(10,1);ecg;zeros(10,1)];
% 用最大值归一化数据
g = g/max(g);
fs = 200;  % 采样率 200Hz
N = length(g);
t = 0:1/fs:(N-1)/fs;

% 阈值法
rate = 0.7;  % 阈值大小
th = max(g)*rate;
a = find(g > th);  % 获得大于阈值的数据
len_a = length(a);
rwave = [];
m = 5;              %窗口宽度
% 比较附近的连续数据，获得局部最大值的横坐标(对应时间)和纵坐标(对应幅度)
for i = 1:len_a
    if g(a(i)) == max(g((a(i)-m):(a(i)+m)))
        rwave = [rwave;a(i) g(a(i))];
    end
end
%间隔判断：比较两相邻峰的距离，剔除伪迹和其他波段的干扰
j = 1;
while j < length(rwave)
    if (rwave(j+1, 1) - rwave(j, 1))/fs < 0.4
        if rwave(j, 2) < rwave(j+1, 2)
           rwave(j, :) = [];
        else
           rwave(j+1, :) = [];
        end
    else
        j = j + 1;
    end
end
% 画信号的 R 波检测图
figure;
plot(t, g);
hold on;
plot((rwave(:, 1) - 1)/fs, rwave(:, 2), '*');

% 计算R-R间隔和心率的平均值
R_R = mean(diff(rwave(:, 1)))/fs
heart_rate = 60/R_R

5.2.2 Pan-Tompkins法检测心电R波尖峰

1. 将信号分别通过给定的低通滤波器、高通滤波器

2. 对滤波后的信号求一阶导数

3. 对求导之后的信号进行平方运算

4. 将信号通过滑动窗口进行积分，选取窗口长度为N

5. 应用阈值法检测经过前四步处理之后的心电信号R波尖峰

clc;clear;close all
% 读入数据
ECG3 = load('...\ECG3.dat');
ECG4 = load('...\ECG4.dat');
ECG5 = load('...\ECG5.dat');
ECG6 = load('...\ECG6.dat');

fs = 200;
N = 4000;
n = 1:4000;
% 原始信号时域图
figure;subplot(4, 1, 1);
plot(n/fs, ECG3);
title('四个心电信号');
ylabel('uV');xlabel('时间/s');
subplot(4, 1, 2);
plot(n/fs, ECG4);
ylabel('uV');xlabel('时间/s');
subplot(4, 1, 3);
plot(n/fs, ECG5);
ylabel('uV');xlabel('时间/s');
subplot(4, 1, 4);
plot(n/fs, ECG6);
ylabel('uV');xlabel('时间/s');

% 加窗并滤波
N1 = 32;
Fc1 = 11;
Fc2 = 50;
flag = 'scale';
win = hamming(N1 + 1);
b = fir1(N1, [Fc1 Fc2]/(fs/2), 'bandpass', win, flag);
ECG3 = filtfilt(b, 1, ECG3);
ECG4 = filtfilt(b, 1, ECG4);
ECG5 = filtfilt(b, 1, ECG5);
ECG6 = filtfilt(b, 1, ECG6);
figure;
subplot(4, 1, 1);
plot(n/fs, ECG3);
title("四个信号加窗滤波后");
ylabel('uV');xlabel('时间/s');
subplot(4, 1, 2);
plot(n/fs, ECG4);
ylabel('uV');xlabel('时间/s');
subplot(4, 1, 3);
plot(n/fs, ECG5);
ylabel('uV');xlabel('时间/s');
subplot(4, 1, 4);
plot(n/fs, ECG6);
ylabel('uV');xlabel('时间/s');

ecg = ECG3;
% 求一阶导数（diff）、平方
ecg = diff(ecg).^2;
% 滑动窗口积分
N2 = 30;
temp2 = [];
for i = 1:length(ecg) - N2
    temp2 = [temp2;trapz(ecg(i:i + N2 - 1))];
end
ecg = temp2;
figure;
plot(ecg);
ylabel('uV');xlabel('时间/s');

% 阈值法
ecg = [zeros(10, 1);ecg;zeros(10, 1)];
rate = 0.7;  % 阈值大小
th = max(ecg)*rate;
a = find(ecg > th);  % 获得大于阈值的数据
len_a = length(a);
rwave = [];
rwave1 = [];
m = 5;  % 窗口宽度
% 比较附近的连续数据，获得局部最大值的横坐标(对应时间)和纵坐标(对应幅度)
for i = 1:len_a
    if ecg(a(i)) == max(ecg((a(i)-m):(a(i)+m)))
        rwave = [rwave;a(i) ecg(a(i))];
    end
end
% 间隔判断：比较两相邻峰的距离，剔除伪迹和其他波段的干扰
j = 1;
while j < length(rwave)
    if (rwave(j+1, 1) - rwave(j, 1))/fs < 0.4
        if rwave(j, 2) < rwave(j+1, 2)
           rwave(j, :) = [];
        else
           rwave(j+1, :) = [];
        end
    else
        j = j + 1;
    end
end

threshold = (max(ecg) - min(ecg))*0.2 + min(ecg);
[rwave1(:,1), rwave1(:,2)] = findpeaks(ecg, 200, 'MinPeakDistance', 0.25, 'MinPeakHeight', threshold);
heartrate = length(rwave1)*60*fs/N
t = 0:1/fs:(length(ecg) - 1)/fs;
figure;
plot(t, ecg);hold on
plot(rwave1(:, 2), rwave1(:, 1), '*');
xlabel('时间/s');title('积分后信号')