ECG信号读出，检测QRS，P，T 波（小波去噪，并根据检测），基于BP辨识的神经网络

转自：https://blog.csdn.net/weixin_34178244/article/details/85746737
1 ECG引进和阅读ECG信号
1）ECG介绍
详细ECG背景应用就不介绍了，大家能够參考百度 谷歌。仅仅是简单说下ECG的结构：
一个完整周期的ECG信号有 QRS P T 波组成，不同的人相应不用的波形，同一个人在不同的阶段波形也不同。我们须要依据各个波形的特点，提取出相应的特征，对不同的人进行身份识别。
2）ECG信号读取
首先须要到MIT-BIH数据库中下载ECG信号，具体的下载地址与程序读取内容介绍能够參考一下地址（讲述的非常具体）：http://blog.csdn.net/chenyusiyuan/article/details/2027887。
读代替码（基于MATLAB）例如以下：

clc; clear all;
%------ SPECIFY DATA ------------------------------------------------------
%%选择文件名称
stringname=‘111’;
%选择你要处理的信号点数
points=10000;
PATH= ‘F:\ECG\MIT-BIH database directory’; % path, where data are saved
HEADERFILE= strcat(stringname,’.hea’); % header-file in text format
ATRFILE= strcat(stringname,’.atr’); % attributes-file in binary format
DATAFILE=strcat(stringname,’.dat’); % data-file
SAMPLES2READ=points; % number of samples to be read
% in case of more than one signal:
% 2*SAMPLES2READ samples are read

%------ LOAD HEADER DATA --------------------------------------------------
fprintf(1,’\n$> WORKING ON %s …\n’, HEADERFILE);
signalh= fullfile(PATH, HEADERFILE);
fid1=fopen(signalh,‘r’);
z= fgetl(fid1);
A= sscanf(z, ‘%*s %d %d %d’,[1,3]);
nosig= A(1); % number of signals
sfreq=A(2); % sample rate of data
clear A;
for k=1:nosig
z= fgetl(fid1);
A= sscanf(z, ‘%*s %d %d %d %d %d’,[1,5]);
dformat(k)= A(1); % format; here only 212 is allowed
gain(k)= A(2); % number of integers per mV
bitres(k)= A(3); % bitresolution
zerovalue(k)= A(4); % integer value of ECG zero point
firstvalue(k)= A(5); % first integer value of signal (to test for errors)
end;
fclose(fid1);
clear A;

%------ LOAD BINARY DATA --------------------------------------------------
if dformat~= [212,212], error(‘this script does not apply binary formats different to 212.’); end;
signald= fullfile(PATH, DATAFILE); % data in format 212
fid2=fopen(signald,‘r’);
A= fread(fid2, [3, SAMPLES2READ], ‘uint8’)’; % matrix with 3 rows, each 8 bits long, = 212bit
fclose(fid2);
M2H= bitshift(A(:,2), -4);
M1H= bitand(A(:,2), 15);
PRL=bitshift(bitand(A(:,2),8),9); % sign-bit
PRR=bitshift(bitand(A(:,2),128),5); % sign-bit
M( : , 1)= bitshift(M1H,8)+ A(:,1)-PRL;
M( : , 2)= bitshift(M2H,8)+ A(:,3)-PRR;
if M(1,:) ~= firstvalue, error(‘inconsistency in the first bit values’); end;
switch nosig
case 2
M( : , 1)= (M( : , 1)- zerovalue(1))/gain(1);
M( : , 2)= (M( : , 2)- zerovalue(2))/gain(2);
TIME=(0:(SAMPLES2READ-1))/sfreq;
case 1
M( : , 1)= (M( : , 1)- zerovalue(1));
M( : , 2)= (M( : , 2)- zerovalue(1));
M=M’;
M(1)=[];
sM=size(M);
sM=sM(2)+1;
M(sM)=0;
M=M’;
M=M/gain(1);
TIME=(0:2(SAMPLES2READ)-1)/sfreq;
otherwise % this case did not appear up to now!
% here M has to be sorted!!!
disp(‘Sorting algorithm for more than 2 signals not programmed yet!’);
end;
clear A M1H M2H PRR PRL;
fprintf(1,’\n$> LOADING DATA FINISHED \n’);

%------ LOAD ATTRIBUTES DATA ----------------------------------------------
atrd= fullfile(PATH, ATRFILE); % attribute file with annotation data
fid3=fopen(atrd,‘r’);
A= fread(fid3, [2, inf], ‘uint8’)’;
fclose(fid3);
ATRTIME=[];
ANNOT=[];
sa=size(A);
saa=sa(1);
i=1;
while i<=saa
annoth=bitshift(A(i,2),-2);
if annoth59
ANNOT=[ANNOT;bitshift(A(i+3,2),-2)];
ATRTIME=[ATRTIME;A(i+2,1)+bitshift(A(i+2,2),8)+…
bitshift(A(i+1,1),16)+bitshift(A(i+1,2),24)];
i=i+3;
elseif annoth60
% nothing to do!
elseif annoth61
% nothing to do!
elseif annoth62
% nothing to do!
elseif annoth==63
hilfe=bitshift(bitand(A(i,2),3),8)+A(i,1);
hilfe=hilfe+mod(hilfe,2);
i=i+hilfe/2;
else
ATRTIME=[ATRTIME;bitshift(bitand(A(i,2),3),8)+A(i,1)];
ANNOT=[ANNOT;bitshift(A(i,2),-2)];
end;
i=i+1;
end;
ANNOT(length(ANNOT))=[]; % last line = EOF (=0)
ATRTIME(length(ATRTIME))=[]; % last line = EOF
clear A;
ATRTIME= (cumsum(ATRTIME))/sfreq;
ind= find(ATRTIME <= TIME(end));
ATRTIMED= ATRTIME(ind);
ANNOT=round(ANNOT);
ANNOTD= ANNOT(ind);

%------ DISPLAY DATA ------------------------------------------------------
figure(1); clf, box on, hold on ;grid on ;
plot(TIME, M(:,1),‘r’);
if nosig==2
plot(TIME, M(:,2),‘b’);
end;
for k=1:length(ATRTIMED)
text(ATRTIMED(k),0,num2str(ANNOTD(k)));
end;
xlim([TIME(1), TIME(end)]);
xlabel(‘Time / s’); ylabel(‘Voltage / mV’);
string=[‘ECG signal ‘,DATAFILE];
title(string);
fprintf(1,’\nKaTeX parse error: Expected 'EOF', got '\n' at position 28: … DATA FINISHED \̲n̲'); % ---------…> ALL FINISHED \n’);

以MIT-BIH数据库中111.dat 为例。

2 去除高频噪声与基线漂移
ECG读取完后，原始ECG信号含有高频噪声和基线漂移，利用小波方法能够去除对应噪声。
详细原理例如以下：将一维的ECG信号进行8层的小波分解后（MATLAB下wavedec函数，小波类型是bior2.6）得到对应的细节系数与近似系数。依据小波原理我们能够知道。1,2层的细节系数包括了大部分高频噪声，8层的近似系数包括了基线漂移。

基于此。我们将1,2层的细节系数（即高频系数置0），8成的近似系数（低频系数）置0。在对应进行小波重构，重构后我们能够明显得到去噪信号。信号无基线漂移。

以下通过图片与代码进一步解说：

小波去噪代码：（matlab）

%%%%%%%%%%%%%%%%%%%去除噪声和基线漂移%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
level=8; wavename=‘bior2.6’;
ecgdata=ECGsignalM1;
figure(2);
plot(ecgdata(1:points));grid on ;axis tight;axis([1,points,-2,5]);
title(‘原始ECG信号’);
%%%%%%%%%%进行小波变换8层
[C,L]=wavedec(ecgdata,level,wavename);
%%%%%%%提取尺度系数，
A1=appcoef(C,L,wavename,1);
A2=appcoef(C,L,wavename,2);
A3=appcoef(C,L,wavename,3);
A4=appcoef(C,L,wavename,4);
A5=appcoef(C,L,wavename,5);
A6=appcoef(C,L,wavename,6);
A7=appcoef(C,L,wavename,7);
A8=appcoef(C,L,wavename,8);
%%%%%%%提取细节系数
D1=detcoef(C,L,1);
D2=detcoef(C,L,2);
D3=detcoef(C,L,3);
D4=detcoef(C,L,4);
D5=detcoef(C,L,5);
D6=detcoef(C,L,6);
D7=detcoef(C,L,7);
D8=detcoef(C,L,8);
%%%%%%%%%%%%重构
A8=zeros(length(A8),1); %去除基线漂移,8层低频信息
RA7=idwt(A8,D8,wavename);
RA6=idwt(RA7(1:length(D7)),D7,wavename);
RA5=idwt(RA6(1:length(D6)),D6,wavename);
RA4=idwt(RA5(1:length(D5)),D5,wavename);
RA3=idwt(RA4(1:length(D4)),D4,wavename);
RA2=idwt(RA3(1:length(D3)),D3,wavename);
D2=zeros(length(D2),1); %去除高频噪声，2层高频噪声
RA1=idwt(RA2(1:length(D2)),D2,wavename);
D1=zeros(length(D1),1);%去除高频噪声，1层高频噪声
DenoisingSignal=idwt(RA1,D1,wavename);
figure(3);
plot(DenoisingSignal);
title(‘去除噪声的ECG信号’); grid on; axis tight;axis([1,points,-2,5]);
clear ecgdata;

去噪前后对照图像例如以下：
去噪前：

去噪后：

3 QRS 检測
QRS检測是处理ECG信号的基础，不管最后实现什么样的功能，QRS波的检測都是前提。所以准确的检測QRS波是特征提取的前提。我採用基于二进样条4层小波变换。在3层的细节系数中利用极大极小值方法能够非常好的检測出R波。3层细节系数的选择是基于R波在3层系数下表现的与其它噪声区别最大；详细实现例如以下：
二进样条小波滤波器： 低通滤波器：[1/4 3/4 3/4 1/4]
高通滤波器：[-1/4 -3/4 3/4 1/4]
在第3层细节系数中首先找到极大极小值对：
1）找极大值方法：找出斜率大于0的值，并赋值为1，其余为0，极大值就在序列类似1， 0这种点，即前面一个值比后面的大的值相应的位置点。
2）找极小值方法：类似极大值,找出斜率<0的值相应的位置，并赋值为1。其余的为0，极小值就在类似1,0的序列中相应的位置。即前面一个值比后面的大的值相应的位置点。
检測出的极大极小值例如以下：
3）设置阈值。提取出R波。我们能够看出。R波的值要明显大于其它位置的值，其在3层细节系数的特点也类似于此。
这样我们就能够设置一个可靠的阈值（将全部点分为4部分。求出每部分最大值的平均值T。阈值为T/3）来提取一组相邻的最大最小值对。这样最大最小值间的过0点就是相应于原始信号的R波点。

R波相应的极大极小值对例如以下：

4）补偿R波点。因为在二进样条小波变换的过程中，3层细节系数与原始信号的相应的位置有10个点的漂移。在程序中须要补偿。
（这个在程序中会给出）。

5）找Q S 波。基于R波的位置，在R波位置（在1层细节系数下）的前3个极点为Q波。在R波位置（1细节系数下）的后3个极点为S波。这样我们就将QRS波定位出来。
6）因为不同的情况，可能造成R波的漏检和错检（把T波检測为R波），我们依据相邻R波的距离进行检測漏检与错检。
当相邻R波的距离<0.4 mean(RR)平均距离时，这是错检。这样去除值小的R波。当相邻R波的距离>1.6mean(RR)时。在两个RR波间找到一个最大的极值对，定位R波。这是防止漏检。

经过上述方法，一个鲁棒性非常好的QRS检測方法就出来了。经过測试，QRS检測能达到98%。检測结果R波用红线标注，Q S 波用黑线标注。
4 T P 波检測
P T 波的检測与R波检測有非常大的相同性。仅仅只是 P T 波在4层细节系数中能够表述出更好的特性。相同依据依据极大极小值原理。能够分别检測出T P波，以及他们的起始点与终止点。即TB,TE,PB PE。详细程序我会在稍后的程序中给出。

各波段检測结果例如以下：

详细QRS T P波检查代码例如以下：

level=4;    sr=360; 
%读入ECG信号
%load ecgdata.mat;
%load ECGsignalM1.mat;
%load Rsignal.mat
mydata = DenoisingSignal;
ecgdata=mydata';
swa=zeros(4,points);%存储概貌信息
swd=zeros(4,points);%存储细节信息
signal=ecgdata(0*points+1:1*points); %取点信号

%算小波系数和尺度系数
%低通滤波器 1/4 3/4 3/4 1/4
%高通滤波器 -1/4 -3/4 3/4 1/4
%二进样条小波

for i=1:points-3
   swa(1,i+3)=1/4*signal(i+3-2^0*0)+3/4*signal(i+3-2^0*1)+3/4*signal(i+3-2^0*2)+1/4*signal(i+3-2^0*3);
   swd(1,i+3)=-1/4*signal(i+3-2^0*0)-3/4*signal(i+3-2^0*1)+3/4*signal(i+3-2^0*2)+1/4*signal(i+3-2^0*3);
end
j=2;
while j<=level
   for i=1:points-24
     swa(j,i+24)=1/4*swa(j-1,i+24-2^(j-1)*0)+3/4*swa(j-1,i+24-2^(j-1)*1)+3/4*swa(j-1,i+24-2^(j-1)*2)+1/4*swa(j-1,i+24-2^(j-1)*3);
     swd(j,i+24)=-1/4*swa(j-1,i+24-2^(j-1)*0)-3/4*swa(j-1,i+24-2^(j-1)*1)+3/4*swa(j-1,i+24-2^(j-1)*2)+1/4*swa(j-1,i+24-2^(j-1)*3);
   end
   j=j+1;
end
%画出原信号和尺度系数。小波系数
%figure(10);
%subplot(level+1,1,1);plot(ecgdata(1:points));grid on ;axis tight;
%title('ECG信号在j=1,2,3,4尺度下的尺度系数对照');
%for i=1:level
%    subplot(level+1,1,i+1);
%    plot(swa(i,:));axis tight;grid on; xlabel('time');ylabel(strcat('a  ',num2str(i)));
%end
%figure(11);
%subplot(level,1,1); plot(ecgdata(1:points)); grid on;axis tight;
%title('ECG信号及其在j=1,2,3,4尺度下的尺度系数及小波系数');
%for i=1:level
%    subplot(level+1,2,2*(i)+1);
%    plot(swa(i,:)); axis tight;grid on;xlabel('time');
%    ylabel(strcat('a   ',num2str(i)));
%    subplot(level+1,2,2*(i)+2);
%    plot(swd(i,:)); axis tight;grid on;
%    ylabel(strcat('d   ',num2str(i)));
%end

%画出原图及小波系数
%figure(12);
%subplot(level,1,1); plot(real(ecgdata(1:points)),'b'); grid on;axis tight;
%title('ECG信号及其在j=1,2,3,4尺度下的小波系数');
%for i=1:level
%    subplot(level+1,1,i+1);
%    plot(swd(i,:),'b'); axis tight;grid on;
%    ylabel(strcat('d   ',num2str(i)));
%end

%**************************************求正负极大值对**********************%
ddw=zeros(size(swd));
pddw=ddw;
nddw=ddw;
%小波系数的大于0的点
posw=swd.*(swd>0);
%斜率大于0
pdw=((posw(:,1:points-1)-posw(:,2:points))<0);
%正极大值点
pddw(:,2:points-1)=((pdw(:,1:points-2)-pdw(:,2:points-1))>0);
%小波系数小于0的点
negw=swd.*(swd<0);
ndw=((negw(:,1:points-1)-negw(:,2:points))>0);
%负极大值点
nddw(:,2:points-1)=((ndw(:,1:points-2)-ndw(:,2:points-1))>0);
%或运算
ddw=pddw|nddw;
ddw(:,1)=1;
ddw(:,points)=1;
%求出极值点的值,其它点置0
wpeak=ddw.*swd;
wpeak(:,1)=wpeak(:,1)+1e-10;
wpeak(:,points)=wpeak(:,points)+1e-10;

%画出各尺度下极值点
%figure(13);
%for i=1:level
%    subplot(level,1,i);
%    plot(wpeak(i,:)); axis tight;grid on;
%ylabel(strcat('j=   ',num2str(i)));
%end
%subplot(4,1,1);
%title('ECG信号在j=1,2,3,4尺度下的小波系数的模极大值点');

interva2=zeros(1,points);
intervaqs=zeros(1,points);
Mj1=wpeak(1,:);
Mj3=wpeak(3,:);
Mj4=wpeak(4,:);
%画出尺度3极值点
figure(14);
plot (Mj3);
%title('尺度3下小波系数的模极大值点');

posi=Mj3.*(Mj3>0);
%求正极大值的平均
thposi=(max(posi(1:round(points/4)))+max(posi(round(points/4):2*round(points/4)))+max(posi(2*round(points/4):3*round(points/4)))+max(posi(3*round(points/4):4*round(points/4))))/4;
posi=(posi>thposi/3);
nega=Mj3.*(Mj3<0);
%求负极大值的平均
thnega=(min(nega(1:round(points/4)))+min(nega(round(points/4):2*round(points/4)))+min(nega(2*round(points/4):3*round(points/4)))+min(nega(3*round(points/4):4*round(points/4))))/4;
nega=-1*(nega