基于离散小波变换(DWT)的心电信号伪影去除及心电信号PQRST波检测（MATLAB R2018）

心电信号是心脏神经，肌肉组织电化学活动的表现形式。这些电化学活动使心脏内部产生一系列非常协调的电刺激脉冲，分别使心房，心室的肌肉细胞兴奋，从而有节律的舒张和收缩。这些生物电活动在体表的不同部位形成不同的电位差变化，就是心电信号。其他生理电活动，如脑电（Electroencephalogram，EEG），肌电（Electromyography，EMG）等都是如此。

心脏主要由心肌组成四腔室结构，上部分为两心房，下部分为两心室。心脏收缩时

的心电活动称为除极，心脏舒张时的心电活动称为复极。一个心动周期开始于窦房结，它是心脏的最高起搏点（也叫一级起搏点）。它出发的激励指令经结间束首先传给房室结（也称第二级起搏点）。房室结向下发出一通传导路，称为房室束，它位于室间隔内。房室束往下又不断发左右两个束支，越分越细，最后分别形成互相交织的网状结构，称为普肯耶纤维，最后终止于心肌内。心电传导变化复杂，呈混沌状，但其有序结果通过周身组织传遍全身，使身体各部位出现有规律的，各向异性的电位变化。一个周期的心电信号如下图所示：

将测量电极按照一定规则放置在人体表面的特定位置，会监测到心电信号的变化曲线，不同的导联方式得到的心电信号波形会有所不同，将这些变化曲线记录下来，就是临床诊断中用到的心电图。在一般情况下，心电图由一系列波群组成，各个波段反映不同阶段心电变化，它们是P波，QRS波群，T波和可能出现的一小段U波。标准心电图如下所示：

心电信号是心脏物理规律的微弱电反映，它有如下特点：

（1）信号微弱，容易被噪声淹没。心电信号幅度为毫伏级，约为0.5~5mV，典型值为1mV。

（2）具有低频和能量集中的特性。心电信号是一种低频信号，频率范围在0.05Hz~100Hz。其能量主要集中在0.25Hz~40Hz。

（3）随机性强，不同人的心电信号都不一样，甚至有通过心音进行人员身份识别的技术手段，但是标准心电可以通过函数模拟产生。

（4）准周期性。心电信号是一种典型的准周期信号，各个周期的信号虽然不完全相同但它们的特征相同。

鉴于此，提出一种基于离散小波变换(DWT)的心电信号伪影去除及心电信号PQRST波检测方法，运行环境为MATLAB R2018，心电信号伪影去除代码如下：

clear all
close all
clc

Fs = 360; % Sampling Frequency
Fnotch = 0.67; % Notch Frequency
BW = 5; % Bandwidth
Apass = 1; % Bandwidth Attenuation
[b, a] = iirnotch (Fnotch/ (Fs/2), BW/(Fs/2), Apass);
Hd = dfilt.df2 (b, a);

load ('100m.mat');
ecgsig = val/200;
t = 0:length(ecgsig)-1;
tx = t./Fs;

subplot (4, 1, 1), plot(tx,ecgsig), title ('ECG Signal with baseline wander'), grid on
y0=filter (Hd, ecgsig);
subplot (4, 1, 2), plot(tx,y0), title ('ECG signal with low-frequency noise (baseline wander) Removed'), grid on

Fnotch = 50; % Notch Frequency
BW = 50; % Bandwidth
Apass = 1; % Bandwidth Attenuation
[b, a] = iirnotch (Fnotch/ (Fs/2), BW/ (Fs/2), Apass);
Hd1 = dfilt.df2 (b, a);
y1=filter (Hd1, y0);
subplot (4, 1, 3), plot (tx,y1), title ('ECG signal with power line noise Removed'), grid on

d = fdesign.lowpass('Fp,Fst,Ap,Ast',0.4,0.5,1,80);
Hd2 = design(d,'equiripple');
y2 = filter(Hd2,y1);
subplot(4,1,4)
plot(tx,y2),title('ECG Signal with high frequency noise removed'), grid on

wt = modwt(y2,4,'sym4');
wtrec = zeros(size(wt));
wtrec(3:4,:) = wt(3:4,:);

y3 = imodwt(wtrec,'sym4');
y3 = abs(y3).^2;
avg = mean(y3);
[Rpeaks,locs] = findpeaks(y3,t,'MinPeakHeight',8*avg,'MinPeakDist',50);
nohb = length(locs);
timelimit = length(ecgsig)/Fs;
hbpermin = (nohb*60)/timelimit;
disp(strcat('Heart Rate = ',num2str(hbpermin)))

figure
plot(t,y3)
grid on
xlim([0,length(ecgsig)])
hold on
plot(locs,Rpeaks,'^r');
％完整代码：mbd.pub/o/bread/ZZWZk5py
xlabel('samples'), title(strcat('Rpeaks found and Heart Rate : ',num2str(hbpermin)))