% function [NR,RRIM,XINLV,NQRSON,NQRSEN,ESNQRSON,ESNQRSEN]=findQRS(ecg,IBW)
 function NR=findQRS(ecg,IBW)
%输入为 【一维数组】ECG信号,【逻辑数值】是否去除基线漂移/补偿基准电压
%输出为【一维数组】R峰值点位置、【数值】平均RR间隔、【一维数组】QRS波群起始点位置、【一维数组】终止点位置、【一维逻辑数组】能否找到起始点、【一维逻辑数组】能否找到终止点

%by zghzgh1779，北京邮电大学，一边开着matlab一边用纸笔计算总觉得哪里不对

%读取ecg信号，tm是采样点，ecg是采样频率360hz的离散数据信号
N=length(ecg)-15;
ecgN=ecg(1:N+15);

%————————

%多孔算法，高通和低通滤波器（补偿时移
ecgH1=filter([1/8,3/8,3/8,1/8],[1,0],ecgN);%一层低通,到下一级 长度4
ecgG1=filter([2,-2],1,ecgN);%一层高通，尺度1空间 长2
ecgH2=filter([1/8,0,3/8,0,3/8,0,1/8],[1,0,0],ecgH1);%2层低通，到下一级 长7
ecgG2=filter([2,0,-2],1,ecgH1);%2层高通，尺度2空间 长3
ecgG2(1:N)=ecgG2(3:N+2);%去群时延长度2
ecgH3=filter([1/8,0,0,0,3/8,0,0,0,3/8,0,0,0,1/8],[1,0,0,0,0],ecgH2);%3层低通，到下一级
ecgG3=filter([2,0,0,0,-2],1,ecgH2);%3层低通，尺度3 长5
ecgG3(1:N)=ecgG3(7:N+6);%去群时延长度6
ecgG4=filter([2,0,0,0,0,0,0,0,-2],1,ecgH3);%4层低通，尺度4 长9
ecgG4(1:N)=ecgG4(15:N+14);%去群时延长度14



%计算R峰阈值
TfQRS1=rms(ecgG1(1:N));%阈值是对应尺度空间的1到N样本点的均方根值的比例项
TfQRS2=rms(ecgG2(1:N));
TfQRS3=1.5*rms(ecgG3(1:N));
TfQRS4=0.5*rms(ecgG4(1:N));

%计算模极大值序列，对正负幅值分开处理
[~,NQRS4]=findpeaks(abs(ecgG4(1:N)),'MinPeakHeight',TfQRS4);%4层尺度空间，取过阈值峰值，得4层模极大值序列

%处理正极大值序列
NQRS4P=NQRS4(ecgG4(NQRS4)>0);%取4层尺度空间正模极大值位置序列
Max4P=length(NQRS4P);%记录NQRS4P长度，也就是处理正模极大值序列总长

%处理负极大值序列
NQRS4N=NQRS4(ecgG4(NQRS4)<0);%取4层尺度空间负模极大值位置序列
Max4N=length(NQRS4N);%记录NQRS4P长度，也就是处理负模极大值序列总长

%用4层求3层正序列
NQRS3P=[];
for i=1:Max4P
if NQRS4P(i)-10>0;%避免第一个样本点位置小于5导致findpeaks函数错误
    FPS=NQRS4P(i)-10;
else
    FPS=1;
end
[modQRS3P,NQRS3TP]=findpeaks(ecgG3(FPS:NQRS4P(i)+10),'MinPeakHeight',TfQRS3); %在NQRS4P序列每个样本点左右5个样本点内搜索正极大值，并把备选极大值丢进临时组
%modQRS3为正小波系数峰值，NQRS3TP为位置,注意NQRS3tP是12345，不是从NQRS4-5开始算
EM3P=isempty(NQRS3TP);%判断是否是空数组，em3P=1则空数组
Max3TP=length(NQRS3TP);%取临时数据组长度
[MAXecg3P,LOCecg3P]=max(modQRS3P);%取最大值
modQRS3P(LOCecg3P)=0;%并删掉对应位置最大值
    if  EM3P==1; %若没找到最大值
        NQRS3P(i)=0;%把NQRS3P序列对应元素设置为0
    elseif MAXecg3P>1.2*modQRS3P(1:Max3TP); %最大的模极大值点大于1.2倍其余极值点
        NQRS3P(i)=FPS-1+NQRS3TP(LOCecg3P);%则把最大值处的3层尺度空间上的位置丢进序列
    else %取距离NQRS4P（i）最近的模极大值点
        NQRS3TnoNP=FPS-1+NQRS3TP(1:Max3TP)-NQRS4P(i);%NQRS3TP序列值减去N，算序列距离
        [~,LOCNQRS3TnoNP]=min(abs(NQRS3TnoNP));%取绝对值，判断距离最小，Loc就是nqrs3tP数组里最近的序号
        NQRS3P(i)=FPS-1+NQRS3TP(LOCNQRS3TnoNP(1));%取nqrs3tP第loc个序列信息存进3尺度序列数组
    end
end

%用4层求3层负序列
NQRS3N=[];
for i=1:Max4N
if NQRS4N(i)-10>0;%避免第一个样本点位置小于5导致findpeaks函数错误
    FPS=NQRS4N(i)-10;
else
    FPS=1;
end
[modQRS3N,NQRS3TN]=findpeaks(-ecgG3(FPS:NQRS4N(i)+10),'MinPeakHeight',TfQRS3); %在NQRS4N序列每个样本点左右5个样本点内搜索负极大值，并把备选极大值丢进临时组
%modQRS3N为负小波系数峰值，NQRS3TN为位置,注意NQRS3tN是12345，不是从NQRSN4-5开始算
EM3N=isempty(NQRS3TN);%判断是否是空数组，em3N=1则空数组
Max3TN=length(NQRS3TN);%取临时数据组长度
[MAXecg3N,LOCecg3N]=max(modQRS3N);%取最大值
modQRS3N(LOCecg3N)=0;%并删掉对应位置最大值
    if  EM3N==1; %若没找到最大值
        NQRS3N(i)=0;%把NQRS3N序列对应元素设置为0
    elseif MAXecg3N>1.2*modQRS3N(1:Max3TN); %最大的模极大值点大于1.2倍其余极值点
        NQRS3N(i)=FPS-1+NQRS3TN(LOCecg3N);%则把最大值处的3层尺度空间上的位置丢进序列
    else %取距离NQRS4N（i）最近的模极大值点
        NQRS3TnoNN=NQRS3TN(1:Max3TN)-NQRS4N(i);%NQRS3TP序列值减去N，算序列距离
        [~,LOCNQRS3TnoNN]=min(abs(NQRS3TnoNN));%取绝对值，判断距离最小，Loc就是nqrs3tN数组里最近的序号
        NQRS3N(i)=FPS-1+NQRS3TN(LOCNQRS3TnoNN(1));%取nqrs3tP第loc个序列信息存进3尺度序列数组
    end
end

% 在3层过滤异常极大值
% 第一次处理正极大序列
NQRS3P2=NQRS3P;%备份3层正极大值序列
for i=1:Max4P;
    if NQRS3P(i)<=0
    else
        NQRS3PR=NQRS3P(abs(NQRS3P-NQRS3P(i))<46);%取45个样本点之内的重复样本点
        NQRS3PR(NQRS3PR==0)=[];%避免已置0点干扰
        if length(NQRS3PR)>1 && ecgG3(NQRS3P(i))<max(ecgG3(NQRS3PR));%若左右45范围内不止原样本点一个正极大值，且原样本点不是最大值
            NQRS3P(i)=0;%把此样本点当做无用正极大并置0
        end
    end
end
%处理负极大序列
for i=1:Max4P;
    NQRS3NR=NQRS3N(abs(NQRS3N-NQRS3P(i))<46);%取正极大值点45个样本点之内的负极大值点
    NQRS3NR(NQRS3NR==0)=[];%避免已置0点干扰
    NQRS3NRLEN=length(NQRS3NR);
    if NQRS3NRLEN>1;%如果只有一个负极大值点则不作处理
        AAA=ecgG3(NQRS3NR);
        BBB=NQRS3NR-NQRS3P(i);
        BBB=BBB';
        QRS3AL=abs(AAA./BBB);%取幅值/距离作为判断条件
        [MQRS3AL,NQRS3AL]=max(QRS3AL);%找出最大值
        QRS3AL(NQRS3AL)=0;
        if MQRS3AL>1.2*QRS3AL;%如果最大值大于1.2倍其余
            NQRS3NR(NQRS3AL)=0;%此点置0，其余点应滤去
        elseif NQRS3NR-NQRS3P(i)>0;%如果都右侧
                NQRS3NR(1)=0;%右侧第一个点置0
        elseif NQRS3NR-NQRS3P(i)<0;%如果都左侧
                NQRS3NR(end)=0;%左侧最后一个点置0
        else
            NQRS3NRD=NQRS3NR-NQRS3P(i);%找右侧最近的
            NQRS3NRD(NQRS3NRD<0)=45;
            [~,NNQRS3NRD]=min(NQRS3NRD);
            NQRS3NR(NNQRS3NRD)=0;%此点置0
        end
        for k=1:Max4N
            IMNQRS3NR=ismember(NQRS3N(k),NQRS3NR);
           if IMNQRS3NR==1;
               NQRS3N(k)=0;%滤掉所有NQRSNR数组中不为0的样本点
           end
        end
    end
end
%再次过滤正极大序列
for i=1:Max4N
    NQRS3PR2=NQRS3P2(abs(NQRS3P2-NQRS3N(i))<46);%取负极大值点45个样本点内正极大值点
    NQRS3PR2(NQRS3PR2==0)=[];%避免已置0点干扰
    NQRS3NR2LEN=length(NQRS3PR2);
    if NQRS3NR2LEN>1;%如果只有一个正极大值则不作处理
        AAA2=ecgG3(NQRS3PR2);
        BBB2=NQRS3PR2-NQRS3N(i);
        BBB2=BBB2';
        QRS3AL2=abs(AAA2./BBB2);%取幅值/距离作为判断条件
        [MQRS3AL2,NQRS3AL2]=max(QRS3AL2);%找出最大值
        QRS3AL2(NQRS3AL2)=0;
        if MQRS3AL2>1.2*QRS3AL2;%如果最大值大于1.2倍其余
            NQRS3PR2(NQRS3AL2)=0;%此点置0，其余点应滤去
        elseif NQRS3PR2-NQRS3N(i)>0;%如果都右侧
                NQRS3PR2(1)=0;%右侧第一个点置0
        elseif NQRS3PR2-NQRS3N(i)<0;%如果都左侧
                NQRS3PR2(end)=0;%左侧最后一个点置0
        else
            NQRS3PRD=NQRS3N(i)-NQRS3PR2;%找左侧最近的
            NQRS3PRD(NQRS3PRD<0)=45;
            [~,NNQRS3PRD]=min(NQRS3PRD);
            NQRS3PR2(NNQRS3PRD)=0;%此点置0
        end
        for k=1:Max4P
            IMNQRS3PR2=ismember(NQRS3P2(k),NQRS3PR2);
           if IMNQRS3PR2==1;
               NQRS3P2(k)=0;%滤掉所有在NQRSPR2数组中不为0的样本点位置
           end
        end
    end
end
NQRS3P=NQRS3P2;%使用第二次处理的正极大值序列


%用3层求2层正序列
NQRS2P=[];
for i=1:Max4P
if NQRS3P(i)-15>0;%避免第一个样本点位置小于5导致findpeaks函数错误
    FPS=NQRS3P(i)-15;
else
    FPS=1;
end
    if NQRS3P(i)==0;%若上一级序列为0，这一级序列置0
        NQRS2P(i)=0;
    else
    [modQRS2P,NQRS2TP]=findpeaks(ecgG2(FPS:min(NQRS3P(i)+15,length(ecgG2))),'MinPeakHeight',TfQRS2); %在NQRS3P序列每个样本点左右5个样本点内搜索极大值，并把备选极大值丢进临时组
    %modQRS2P为模绝对值峰值，NQRS2TP为位置
    EM2P=isempty(NQRS2TP);%判断是否是空数组，em2P=1则空数组
    Max2TP=length(NQRS2TP);%取临时数据组长度
    [MAXecg2P,LOCecg2P]=max(modQRS2P);%取最大值
    modQRS2P(LOCecg2P)=0;%并删掉对应位置最大值
        if  EM2P==1; %若没找到最大值
            NQRS2P(i)=0;%把NQRS2P序列对应元素设置为0
        elseif MAXecg2P>1.2*modQRS2P(1:Max2TP); %模极大值点大于1.2倍其余
            NQRS2P(i)=FPS-1+NQRS2TP(LOCecg2P);%则把最大值处的2层尺度空间上的位置丢进序列
        else %取距离NQRS3P（i）最近的模极大值点
            NQRS2TnoNP=FPS-1+NQRS2TP(1:Max2TP)-NQRS3P(i);%NQRS2TP序列值减去N，算序列距离
            [~,LOCNQRS2TnoNP]=min(abs(NQRS2TnoNP));%取绝对值，判断距离最小，Loc就是nqrs2tP数组里最近的序号
            NQRS2P(i)=FPS-1+NQRS2TP(LOCNQRS2TnoNP(1));%取nqrs2tP第loc个序列信息存进2尺度序列数组
        end
    end
end

%用2层求1层正序列
NQRS1P=[];
for i=1:Max4P
if NQRS2P(i)-15>0;%避免第一个样本点位置小于5导致findpeaks函数错误
    FPS=NQRS2P(i)-15;
else
    FPS=1;
end
    if NQRS2P(i)==0;
        NQRS1P(i)=0;
    else
    [modQRS1P,NQRS1TP]=findpeaks(ecgG1(FPS:min(NQRS2P(i)+15,length(ecgG1))),'MinPeakHeight',TfQRS1); %在NQRS2P序列每个样本点左右5个样本点内搜索极大值，并把备选极大值丢进临时组
    %modQRS2P为模绝对值峰值，NQRS2TP为位置
    EM1P=isempty(NQRS1TP);%判断是否是空数组，em1P=1则空数组
    Max1TP=length(NQRS1TP);%取临时数据组长度
    [MAXecg1P,LOCecg1P]=max(modQRS1P);%取最大值
    modQRS1P(LOCecg1P)=0;%并删掉对应位置最大值
        if  EM1P==1; %若没找到最大值
            NQRS1P(i)=0;%把NQRS2序列对应元素设置为0
        elseif MAXecg1P>1.3*modQRS1P(1:Max1TP); %模极大值点大于1.3倍其余
            NQRS1P(i)=FPS-1+NQRS1TP(LOCecg1P);%则把最大值处的2层尺度空间上的位置丢进序列
        else %取距离NQRS2P（i）最近的模极大值点
            NQRS1TnoNP=FPS-1+NQRS1TP(1:Max1TP)-NQRS2P(i);%NQRS2T序列值减去N，算序列距离
            [~,LOCNQRS1TnoNP]=min(abs(NQRS1TnoNP));%取绝对值，判断距离最小，Loc就是nqrs2t数组里最近的序号
            NQRS1P(i)=FPS-1+NQRS1TP(LOCNQRS1TnoNP(1));%取nqrs2t第loc个序列信息存进2尺度序列数组
        end
    end
end

%用3层求2层负序列
NQRS2N=[];
for i=1:Max4N
if NQRS3N(i)-15>0;%避免第一个样本点位置小于5导致findpeaks函数错误
    FPS=NQRS3N(i)-15;
else
    FPS=1;
end
    if NQRS3N(i)==0;%若上一级序列为0，这一级序列置0
        NQRS2N(i)=0;
    else
    [modQRS2N,NQRS2TN]=findpeaks(-ecgG2(FPS:min(NQRS3N(i)+15,length(ecgG2))),'MinPeakHeight',TfQRS2); %在NQRS3N序列每个样本点左右5个样本点内搜索极大值，并把备选极大值丢进临时组
    %modQRS2N为模绝对值峰值，NQRS2TN为位置
    EM2N=isempty(NQRS2TN);%判断是否是空数组，em2N=1则空数组
    Max2TN=length(NQRS2TN);%取临时数据组长度
    [MAXecg2N,LOCecg2N]=max(modQRS2N);%取最大值
    modQRS2N(LOCecg2N)=0;%并删掉对应位置最大值
        if  EM2N==1; %若没找到最大值
            NQRS2N(i)=0;%把NQRS2N序列对应元素设置为0
        elseif MAXecg2N>1.2*modQRS2N(1:Max2TN); %模极大值点大于1.2倍其余
            NQRS2N(i)=FPS-1+NQRS2TN(LOCecg2N);%则把最大值处的2层尺度空间上的位置丢进序列
        else %取距离NQRS3N（i）最近的模极大值点
            NQRS2TnoNN=FPS-1+NQRS2TN(1:Max2TN)-NQRS3N(i);%NQRS2TN序列值减去N，算序列距离
            [~,LOCNQRS2TnoNN]=min(abs(NQRS2TnoNN));%取绝对值，判断距离最小，Loc就是nqrs2tN数组里最近的序号
            NQRS2N(i)=FPS-1+NQRS2TN(LOCNQRS2TnoNN(1));%取nqrs2tN第loc个序列信息存进2尺度序列数组
        end
    end
end

%用2层求1层负序列
NQRS1N=[];
for i=1:Max4N
if NQRS2N(i)-15>0;%避免第一个样本点位置小于5导致findpeaks函数错误
    FPS=NQRS2N(i)-15;
else
    FPS=1;
end
    if NQRS2N(i)==0;
        NQRS1N(i)=0;
    else
        [modQRS1N,NQRS1TN]=findpeaks(-ecgG1(FPS:min(NQRS2N(i)+15,length(ecgG1))),'MinPeakHeight',TfQRS1); %在NQRS2N序列每个样本点左右5个样本点内搜索极大值，并把备选极大值丢进临时组
    %modQRS2N为模绝对值峰值，NQRS2TN为位置
    EM1N=isempty(NQRS1TN);%判断是否是空数组，em1N=1则空数组
    Max1TN=length(NQRS1TN);%取临时数据组长度
    [MAXecg1N,LOCecg1N]=max(modQRS1N);%取最大值
    modQRS1N(LOCecg1N)=0;%并删掉对应位置最大值
        if  EM1N==1; %若没找到最大值
            NQRS1N(i)=0;%把NQRS2序列对应元素设置为0
        elseif MAXecg1N>1.3*modQRS1N(1:Max1TN); %模极大值点大于1.3倍其余
            NQRS1N(i)=FPS-1+NQRS1TN(LOCecg1N);%则把最大值处的2层尺度空间上的位置丢进序列
        else %取距离NQRS2N（i）最近的模极大值点
            NQRS1TnoNN=FPS-1+NQRS1TN(1:Max1TN)-NQRS2N(i);%NQRS1TN序列值减去N，算序列距离
            [~,LOCNQRS1TnoNN]=min(abs(NQRS1TnoNN));%取绝对值，判断距离最小，Loc就是nqrs1tN数组里最近的序号
            NQRS1N(i)=FPS-1+NQRS1TN(LOCNQRS1TnoNN(1));%取nqrs2t第loc个序列信息存进1尺度序列数组
        end
    end
end

%去除孤立极大值

NQRS1P(NQRS1P==0)=[];
NQRS1N(NQRS1N==0)=[];
Max1P=length(NQRS1P);
Max1N=length(NQRS1N);
for i=1:Max1P
    ISON=0;%判断是否孤立极大值
    for k=1:Max1N
        if abs(NQRS1P(i)-NQRS1N(k))<45
            ISON=1;
        end
    end
    if ISON==0;
        NQRS1P(i)=0;
    end
end
NQRS1P(NQRS1P==0)=[];
Max1P=length(NQRS1P);
for i=1:Max1N
    ISON=0;%判断是否孤立极大值
    for k=1:Max1P;
        if abs(NQRS1N(i)-NQRS1P(k))<45
            ISON=1;
        end
    end
    if ISON==0;
        NQRS1N(i)=0;
    end
end
NQRS1N(NQRS1N==0)=[];

%找零交点
NQRS=[NQRS1P,NQRS1N]; %把1尺度空间最后得到的序列数保存
NQRS(NQRS==0)=[]; %去除序列数组中为0的点（孤立极大、被滤去点等）
NQRS=sort(NQRS,1,'ascend');%升序排列
NQRS=unique(NQRS);%去重复点
MaxNQRS=length(NQRS);%取长度
NR=[];%R峰序列
i=1;
for i=1:(MaxNQRS-1);
    if NQRS(i+1)-NQRS(i)>46 %45个样本点内没有下一个极值，置0
        NR(i)=0;
    elseif i>1 && NQRS(i)-NR(i-1)<46 %与上一R峰在45样本点内，置0
        NR(i)=0;
    else %否则求该点与下一个极值点之间的零交点
        if NQRS(i+1)>length(ecgG1)
            LOCNR=pickzero(ecgG1(NQRS(i):end));
        else
            LOCNR=pickzero(ecgG1(NQRS(i):NQRS(i+1)));%找ecgG1信号过零点
        end
        if(isempty(LOCNR));LOCNR=0;end
        if ecgG1(NQRS(i))>0;
            NR(i)=NQRS(i)-1+LOCNR(1);%先正后负，把第一个零点保存进R峰序列
        else
            NR(i)=NQRS(i)-1+LOCNR(end);%先负后正，把最后一个零点保存进R峰序列
        end
    end 
end
NR(NR==0)=[];%去掉NR中的零点

%在原信号空间内标记极大值点
if IBW==1
LENR=length(NR);
    for i=1:LENR;
        if ecgN(NR(i))>0;%找左右两侧同符号的极值点作为峰值
            [~,LOCMAXNR]=max(ecgN(NR(i)-5:NR(i)+5));
            NR(i)=NR(i)-6+LOCMAXNR;
        else
            [~,LOCMINNR]=min(ecgN(NR(i)-5:NR(i)+5));
            NR(i)=NR(i)-6+LOCMINNR;
        end
    end
end

%【说明】此处NR数组即标定的R峰的位置，元素数值代表样本点在时间上的位置（也即tm数组中的第几个元素），单位为1/360秒

%求RR间期
% LENR=length(NR);
% NR1=NR(1:LENR-1);
% NR2=NR(2:LENR);
% RRI=NR2-NR1;%此数组即RR间隔
% RRIM=mean(RRI);%RR间隔均值
% RRI(RRI>1.5*RRIM)=[];
% RRIM=mean(RRI);
% %【说明】RRIM数值即平均RR间隔，可以用来计算心率
% 
% %计算每分钟心率
% XINLV=60*360./RRIM;
% XINLV=round(XINLV);
% 
% %拆分极值点对,按奇偶项分别保存
% NPRE=NQRS(1:2:end);%极值点对前者
% NPOST=NQRS(2:2:end);%极值点对后者
% LENR=length(NR);%R峰点和极值点对个数
% for i=1:LENR
%     if NPRE(i)-15>0;%避免第一个样本点搜索窗溢出
%         FNPRE=NPRE(i)-15;
%     else
%         FNPRE=1;
%     end
%     TFQRSPRE=0.03*abs(ecgG2(NPRE(i)));%特征峰pre阈值
%     %处理起始点序列
%     [~,NPRETP]=findpeaks(ecgG2(FNPRE:NPRE(i)-1),'MinPeakHeight',TFQRSPRE); %第二层尺度空间npre前15个样本点内找过阈值的局部最大模
%     [~,NPRETN]=findpeaks(-ecgG2(FNPRE:NPRE(i)-1),'MinPeakHeight',TFQRSPRE); %分正负找峰值点
%     NPRET=[NPRETP:NPRETN];
%     NPRET=sort(NPRET,1,'ascend');%升序排列
%     EM1=isempty(NPRET);
%     if EM1==1;%如果找不到符合条件的特征峰值点
%         NPRETNO=pickzero(ecgG2(FNPRE:NPRE(i)-1));%如果没有，取最近的QRS2过零点作为起始点【阈值没有
%         NSPRE(i)=0;%特征峰值点置0
%         if isempty(NPRETNO)==1;%如果搜索窗内找不到过零点
%             [~,NPRETNO]=min(abs(ecgG2(FNPRE:NPRE(i)-1)));%把最小值标记为起始点
%             ESNQRSON(i)=0;%标记起始点无法找到
%         else
%             ESNQRSON(i)=1;%标记起始点可以找到
%         end
%         NQRSON(i)=FNPRE-1+NPRETNO(end);%最后一个最近
%     else %如果找到符合的特征峰值点
%         NPRET=FNPRE-1+NPRET;
%         NPRET(ecgG2(NPRET)*ecgG2(NPRE(i))>0)=[];%去掉与NPRE点同号的样本点，留异号的
%         EM2=isempty(NPRET);
%         if EM2==1;
%             NPRETNO=pickzero(ecgG2(FNPRE:NPRE(i)-1));%如果没有异号峰值点，取最近的QRS2过零点作为起始点
%             NSPRE(i)=0;%特征峰值点置0
%             if isempty(NPRETNO)==1;%如果搜索窗内找不到过零点
%                 [~,NPRETNO]=min(abs(ecgG2(FNPRE:NPRE(i)-1)));%把最小值标记为起始点
%                 ESNQRSON(i)=0;%标记起始点无法找到
%             else
%                 ESNQRSON(i)=1;%标记起始点可以找到
%             end
%             NQRSON(i)=FNPRE-1+NPRETNO(end);
%         else
%             NSPRE(i)=NPRET(end);%距离NPRE最近的局部峰值点被标记为特征峰值点,
%             NSPRET=pickzero(ecgG2(FNPRE:NSPRE(i)-1));%从特征峰值向前15个样本点寻找过零点
%             if isempty(NSPRET)==1;%如果搜索窗内找不到过零点
%                 [~,NSPRET]=min(abs(ecgG2(FNPRE:NPRE(i)-1)));%把最小值标记为起始点
%                 ESNQRSON(i)=0;%标记起始点无法找到
%             else
%                 ESNQRSON(i)=1;%标记起始点可以找到
%             end
%             NQRSON(i)=FNPRE-1+NSPRET(end);%标记最近的过零点为起始点
%         end
%     end
%     
%     if NPOST(i)+30>N;
%         LNPOST=N;
%     else
%         LNPOST=NPOST(i)+30;
%     end
% %     if NPOST(i)>length(ecgG2)
% %     else
%     TFQRSPOST=0.03*abs(ecgG2(NPOST(i)));%特征峰POST阈值
% %     end    %A00181出错
%     %处理终止点序列
%     [~,NPOSTTP]=findpeaks(ecgG2(NPOST(i)+1:LNPOST),'MinPeakHeight',TFQRSPOST); %第二层尺度空间npost后30个样本点内找过阈值的局部最大模
%     [~,NPOSTTN]=findpeaks(-ecgG2(NPOST(i)+1:LNPOST),'MinPeakHeight',TFQRSPOST); %第二层尺度空间npost后30个样本点内找过阈值的局部最大模
%     NPOSTT=[NPOSTTP:NPOSTTN];
%     NPOSTT=sort(NPOSTT,1,'ascend');%升序排列
%     EM1=isempty(NPOSTT);
%     if EM1==1;%如果找不到符合条件的特征峰值点
%         NPOSTTNO=pickzero(ecgG2(NPOST(i)+1:LNPOST));%如果没有，取最近的QRS2过零点作为起始点【阈值没有
%         NSPOST(i)=0;%特征峰值点置0
%         if isempty(NPOSTTNO)==1;%如果搜索窗内找不到过零点
%             [~,NPOSTTNO]=min(abs(ecgG2(FNPRE:NPRE(i)-1)));%把最小值标记为终止点
%             ESNQRSEN(i)=0;%标记终止点无法找到
%         else
%             ESNQRSEN(i)=1;%标记终止点可以找到
%         end
%         NQRSEN(i)=NPOST(i)+NPOSTTNO(1);%第一个最近
%     else %如果找到符合的特征峰值点
%         NPOSTT=NPOST(i)+NPOSTT;
%         NPOSTT(ecgG2(NPOSTT)*ecgG2(NPOST(i))>0)=[];%去掉与NPost点同号的样本点，留异号的
%         EM2=isempty(NPOSTT);
%         if EM2==1;
%             NPOSTTNO=pickzero(ecgG2(NPOST(i)+1:LNPOST));%如果没有，取最近的QRS2过零点作为起始点
%             NSPOST(i)=0;%特征峰值点置0
%             if isempty(NPOSTTNO)==1;%如果搜索窗内找不到过零点
%                 [~,NPOSTTNO]=min(abs(ecgG2(FNPRE:NPRE(i)-1)));%把最小值标记为终止点
%                 ESNQRSEN(i)=0;%标记终止点无法找到
%             else
%                 ESNQRSEN(i)=1;%标记终止点可以找到
%             end
%             NQRSEN(i)=NPOST(i)+NPOSTTNO(1);
%         else
%             NSPOST(i)=NPOSTT(1);%距离Npost最近的局部峰值点被标记为特征峰值点
%             NSPOSTT=pickzero(ecgG2(NSPOST(i)+1:LNPOST));%从特征峰值向后30个样本点寻找过零点
%             if isempty(NSPOSTT)==1;%如果搜索窗内找不到过零点
%                 [~,NSPOSTT]=min(abs(ecgG2(FNPRE:NPRE(i)-1)));%把最小值标记为终止点
%                 ESNQRSEN(i)=0;%标记终止点无法找到
%             else
%                 ESNQRSEN(i)=1;%标记终止点可以找到
%             end
%             NQRSEN(i)=NSPOST(i)+NSPOSTT(1);%标记最近的过零点为终止点
%         end
%     end
end
%【说明】NQRSON是起始点序列，NQRSEN是终止点序列，ESNQRSON数组标志第i个峰值前起始点能不能正常找到，1为可以找到，0为无法正常找到，选择了一个可能的接近点代替
