gusucode.com > 能检测心电特征点PQRST，用matlab编写了界面，并能计算识别率 > ecg/PQRST_detect.m

    function [R Q S P T X1 Y4 Y1 se Qstart Send]=PQRST_detect(X,sfreq)   %差分阈值法特征点检测
%确定QRS波群
for n=2:3599
    Y1(n)=X(n+1)-X(n-1);          %一阶导数
end
for n=3:3598
    Y2(n)=X(n+2)-2*X(n)+X(n-2);   %二阶导数
end
max1=max(Y1);
max2=max(Y2);
for n=3:3598
    Y3(n)=max1*Y1(n)+max2*Y2(n);
end
a=0.25;
max3=max(Y3);
for i=3:3598
    if (abs(Y3(i))>max3*a)
       Y4(i)=1;
    else
       Y4(i)=0;
    end
end


%对Y4的处理1，去毛刺
for i=4:3597
    if (Y4(i)==1&&Y4(i-1)==0&&Y4(i+1)==0)
        Y4(i)=0;
    end
end


%对Y4的处理2，使QRS波群范围内一阶导数全为1，无“零点”，最终确定QRS波群
L=280;                          %窗口大小，根据两个R峰之间的点数来确定
for i=3:3542
    if (Y4(i)==0&&Y4(i-1)==1)
        for j=1:L/5
            if Y4(i+j)==1
               Y4(i)=1;
               break;
            end
        end
    end
end
for i=3598:-1:3542
    if Y4(i)==1
        for j=3543:i
            Y4(j)=1;
        end
        break;
    end
end


%确定R峰
B=[];
for i=4:3597                  %找出QRS波群的起始点
    if (Y4(i)==1&&Y4(i-1)==0)
        B=[B,i];
    elseif (Y4(i)==1&&Y4(i+1)==0)
        B=[B,i];
    end
end
s=size(B);
R=[];
for j=1:2:s(2)
     a=B(j);
     b=B(j+1);
    [maxR,r0]=max(X(1,a:b));  %找出QRS波群内幅值最大的一点，即为R峰
     r=a+r0-1;                %得出的最大值的位置是以a开头向后数（R0-1）位
     R=[R,r]; 
end


%求峰域值
dsum=[];
for n=2:3599
    d=Y1(n).^2;
    dsum=[dsum,d];
end
dert=sum(dsum)/3598;  %方差
se=26*dert;


%R峰起点
Rnum=size(R);
Rstart=[];
for i=1:Rnum(2)
    s1=0;     
    for j=(R(i)-1):-1:1
        if abs(Y1(j))>se
            s1=j;        %R峰左边斜率首先超过se线的点为起点s1
            break;
        end
    end
    for m=s1:-1:1
        if (abs(Y1(m))<se)     %第一个斜率穿过se线的即为R峰起点
            Rstart=[Rstart,m];
            break;
        end
    end 
end


%Q峰位置
Q=[];
for i=1:Rnum(2)
    for j=Rstart(i):-1:1
        if (Rstart(i)-j)<0.04*sfreq                             %由R峰起点向前0.04s内（因为Q波的宽度小于0.04s）
            if (abs(X(j))>=abs(X(j-1))&&abs(X(j))>=abs(X(j+1))) %第一个绝对值极大值即为Q峰 
                Q=[Q,j];
                break;
            end
        end
    end
end


%R峰终点
Rend=[];
for i=1:Rnum(2)
    s2=0;
    c1=0;
    c2=0;
    for j=(R(i)+1):3598
        if abs(Y1(j))>se  %R峰右边斜率首先超过se线的点为起点s2
            s2=j;
            break;
        end
    end
    for k=(s2+1):3598
        if (abs(Y1(k))>se&&abs(Y1(k+1))<se) %第一个自下而上穿越-se线的即为R峰终点
            c2=k;
            break;
        end
    end
    Rend=[Rend,c2];
end


%S峰位置
S=[];
len=round(0.15*sfreq);  %由R峰起点向后0.15s内
for i=1:Rnum(2)
    a=Rend(i);
    [MIN s0]=min(X(1,(a:a+len)));
    s=a+s0-1;
    S=[S,s];
end


%确定S峰终点
Send=[];
for i=1:Rnum(2)
    for j=S(i):3600
        if (abs(Y1(j))>=abs(Y1(j+1))&&abs(Y1(j))>=abs(Y1(j-1)))  %S峰之后第一个斜率的绝对值极大值点即为S峰终点
            Send=[Send,j];
            break;
        end
    end
end 


%确定Q峰起点
Qstart=[];
for i=1:Rnum(2)
    for j=Q(i):-1:1
        if (abs(Y1(j))>=abs(Y1(j+1))&&abs(Y1(j))>=abs(Y1(j-1))) %Q峰之前第一个斜率的绝对值极大值点即为Q峰起点
            Qstart=[Qstart,j];
            break;
        end
    end
end 


%心电去噪
for i=5:3596
    X1(i)=(X(i-4)+2*X(i-3)+3*X(i-2)+4*X(i-1)+5*X(i)+4*X(i+1)+3*X(i+2)+2*X(i+3)+X(i+4))/25;
end
X1(1)=X(1);
X1(2)=X(2);
X1(3)=X(3);
X1(4)=X(4);
X1(3597)=X(3597);
X1(3598)=X(3598);
X1(3599)=X(3599);
X1(3600)=X(3600);


%确定窗的宽度
wsum=0;
for i=2:Rnum(2)
    wsum=wsum+R(i)-R(i-1);
end
RRavi=wsum/(Rnum(2)-1);
RRav=RRavi/sfreq;
if RRav>0.7
    bw=0.03*sfreq;
    ew=0.4*sfreq;
else 
    bw=0.02*sfreq;
    ew=0.4*RRav*sfreq;
end


%确定T峰
wstart=round(Send+bw); %S峰之后的窗
wend=round(Send+ew);
T=[];
for i=1:Rnum(2)
    [tmax,t0]=max(X1(1,wstart(i):wend(i)));  %窗内心电信号的最大值即为T峰
    t=wstart(i)+t0-1;
    T=[T,t];
end


%确定P峰
w1start=round(Qstart-ew); %Q峰之前的窗
w1end=round(Qstart-bw);
if w1start(1)<0
    w1start(1)=1;
end
P=[];
for i=1:Rnum(2)
    [pmax,p0]=max(X1(1,w1start(i):w1end(i)));  %窗内心电信号的最大值即为P峰
    p=w1start(i)+p0-1;
    P=[P,p];
end
end