参考链接：ECG ×AI: 机器/深度学习的ECG应用入门（3）_Aiwiscal的博客-CSDN博客_ecg 深度学习

function [QRS_amp,QRS_ind] = DS_detect( ecg_i,gr )
%% function [QRS_amp,QRS_ind]=DS_detect(ecg_i,gr)
%% 输入
% ecg_i : 原信号，一维向量;1*N形式
% 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;
size(ecg_i,2)
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