Skip to content

Commit

Permalink
Add files via upload
Browse files Browse the repository at this point in the history 
A DS-based QRS detection algorithm
  • Loading branch information
@Aiwiscal
Aiwiscal committed May 16, 2018
1 parent c599851 commit 3166003
Showing 1 changed file with 190 additions and 0 deletions.
190 changes: 190 additions & 0 deletions DS_detect.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,190 @@
function [QRS_amp,QRS_ind] = DS_detect( ecg_i,gr )
%% function [QRS_amp,QRS_ind]=DS_detect(ecg_i,gr)
%% 输入
% ecg_i : 原信号，一维向量
% gr : 绘图与否，0：不绘图，1：绘图
%% 输出
% QRS_amp:QRS波振幅.
% QRS_ind:QRS波索引.
% 绘制图像.
%% 作者：刘文涵（WHliu@whu.edu.cn）
%% 版本：1.1
%% 04-24-2018
%% 代码:
if nargin < 2
gr = 1;
if nargin<1
error('The algorithm need a input:ecg_i.');
end
end
if ~isvector(ecg_i)
error('ecg_i must be a row or column vector.');
end
fs=360;
if size(ecg_i,2)<round(1.5*fs)+1
error('The algorithm need a longer input.');
end
tic,
s=ecg_i;
N=size(s,2);
ECG=s;
FIR_c1=[0.0041,0.0053,0.0068,0.0080,0.0081,0.0058,-0.0000,-0.0097,-0.0226,...
-0.0370,-0.0498,-0.0577,-0.0576,-0.0477,-0.0278,0,0.0318,0.0625,0.0867,...
0.1000,0.1000,0.0867,0.0625,0.0318,0,-0.0278,-0.0477,-0.0576,-0.0577,...
-0.0498,-0.0370,-0.0226,-0.0097,-0.0000,0.0058,0.0081,0.0080,0.0068,...
0.0053,0.0041]; % 使用fdatool设计并导出的滤波器系数,带通FIR,15~25Hz,详情使用fdatool打开DS1.fda查看
FIR_c2=[0.0070,0.0094,0.0162,0.0269,0.0405,0.0555,0.0703,0.0833,0.0928,...
0.0979,0.0979,0.0928,0.0833,0.0703,0.0555,0.0405,0.0269,0.0162,0.0094,...
0.0070]; % 使用fdatool设计并导出的滤波器系数,低通FIR,截止频率5Hz,详情使用fdatool打开DS2.fda查看

l1=size(FIR_c1,2);
ECG_l=[ones(1,l1)*ECG(1) ECG ones(1,l1)*ECG(N)]; % 数据点延拓，防止滤波边缘效应；
ECG=filter(FIR_c1,1,ECG_l); % 使用filter滤波；
ECG=ECG((l1+1):(N+l1)); % 前面延拓了数据点，这里截取有用的部分；

%% 双斜率处理
a=round(0.015*fs); % 两侧目标区间0.015~0.060s;
b=round(0.060*fs);
Ns=N-2*b; % 确保在不超过信号长度；
S_l=zeros(1,b-a+1);
S_r=zeros(1,b-a+1);
S_dmax=zeros(1,Ns);
for i=1:Ns % 对每个点双斜率处理；
for k=a:b
S_l(k-a+1)=(ECG(i+b)-ECG(i+b-k))./k;
S_r(k-a+1)=(ECG(i+b)-ECG(i+b+k))./k;
end
S_lmax=max(S_l);
S_lmin=min(S_l);
S_rmax=max(S_r);
S_rmin=min(S_r);
C1=S_rmax-S_lmin;
C2=S_lmax-S_rmin;
S_dmax(i)=max([C1 C2]);
end

%% 再次进行低通滤波，思路与上述带通滤波一致
l2=size(FIR_c2,2);
S_dmaxl=[ones(1,l2)*S_dmax(1) S_dmax ones(1,l2)*S_dmax(Ns)];
S_dmaxt=filter(FIR_c2,1,S_dmaxl);
S_dmaxt=S_dmaxt((l2+1):(Ns+l2));

%% 滑动窗口积分
w=8;wd=7;
d_l=[zeros(1,w) S_dmaxt zeros(1,w)]; % 零延拓，确保所有的点都可以进行窗口积分；?
m=zeros(1,Ns);
for n=(w+1):(Ns+w) % 滑动窗口；
m(n-w)=sum(d_l(n-w:n+w)); % 积分；
end
m_l=[ones(1,wd)*m(1) m ones(1,wd)*m(Ns)];

%% 双阈值检测与动态变化
QRS_buf1=[]; % 存储检测到的QRS波索引
AMP_buf1=[]; % 存储最近检测到的8个QRS波对应特征信号的波峰值
thr_init0=0.4;thr_lim0=0.23;
thr_init1=0.6;thr_lim1=0.3; %% 阈值变化的初始值和下限设置 ?
en=-1; % 标记波峰检出情况，高于高阈值--1，高低阈值之间--0，未检出-- -1?
thr0=thr_init0;
thr1=thr_init1;
thr1_buf=[]; %瀛樺偍闃堝?鍙樺寲鎯呭喌锛?
thr0_buf=[];
for j=8:Ns
t=1;
cri=1;
while t<=wd&&cri>0 % 检测候选波峰；
cri=((m_l(j)-m_l(j-t))>0)&&(m_l(j)-m_l(j+t)>0);
t=t+1;
end
if t==wd+1
N1=size(QRS_buf1,2); %N1:已经检测到的QRS波个数
if m_l(j)>thr1 % 高于高阈值时的处理
if N1<2 % N1小于2时直接存储；
QRS_buf1=[QRS_buf1 (j-wd)]; % j-wd 减去了滑动窗口积分带来的延迟；
AMP_buf1=[AMP_buf1 m_l(j)];
en=1;
else
dist=j-wd-QRS_buf1(N1);
if dist>0.24*fs % 检测波峰距离；
QRS_buf1=[QRS_buf1 (j-wd)];
AMP_buf1=[AMP_buf1 m_l(j)];
en=1;
else
if m_l(j)>AMP_buf1(end) % 不应期处理
QRS_buf1(end)=j-wd;
AMP_buf1(end)=m_l(j);
en=1;
end
end
end

else % 特征峰值低于高阈值

if N1<2&&m_l(j)>thr0 % 特征峰值在两阈值之间
QRS_buf1=[QRS_buf1 (j-wd)];
AMP_buf1=[AMP_buf1 m_l(j)];
en=0;
else
if m_l(j)>thr0 % 特征峰值在两阈值之间
dist_m=mean(diff(QRS_buf1));
dist=j-wd-QRS_buf1(N1);
if dist>0.24*fs && dist>0.5*dist_m % 不应期检测，并且，波峰要距离足够远（> 平均距离的一半）
QRS_buf1=[QRS_buf1 (j-wd)];
AMP_buf1=[AMP_buf1 m_l(j)];
en=0;
else
if m_l(j)>AMP_buf1(end)
QRS_buf1(end)=j-wd;
AMP_buf1(end)=m_l(j);
en=0;
end
end
else
en=-1;
end
end
end
N2=size(AMP_buf1,2);
if N2>8
AMP_buf1=AMP_buf1(2:9); % 确保只存储最近的8个特征波峰；
end
% 下面的if与博文中的公式对应
if en==1
thr1=0.7*mean(AMP_buf1);
thr0=0.25*mean(AMP_buf1);
else
if en==0
thr1=thr1-(abs(m_l(j)-mean(AMP_buf1)))/2;
thr0=0.4*m_l(j);
end
end
end
if thr1<=thr_lim1 % 确保阈值高于下限
thr1=thr_lim1;
end

if thr0<=thr_lim0
thr0=thr_lim0;
end

thr1_buf=[thr1_buf thr1];
thr0_buf=[thr0_buf thr0];
end
delay=round(l1/2)-2*w+2;
QRS_ind=QRS_buf1-delay; % 减去延迟，得到最终结果；
QRS_amp=s(QRS_ind);
toc
if gr==1 %缁樺浘
subplot(2,1,1);plot(m);axis([1 size(m,2) -0.3 1.6*max(m)]);
hold on;title('Feature signal and thresholds');grid on;
plot(QRS_buf1,m(QRS_buf1),'ro');
plot(thr1_buf,'r');
plot(thr0_buf,'k');
legend('Feature Signal','QRS Locations','Threshold1','Threshold0');
subplot(2,1,2);plot(s);%axis([1 size(s,2) min(s) 1.5*max(s)]);
xlabel('n');ylabel('Voltage / mV');
hold on;title('The result on the raw ECG');grid on;
plot(QRS_ind,QRS_amp,'ro');
legend('Raw ECG','QRS Locations');
end
end

0 comments on commit 3166003

Please sign in to comment.