detect.asv

function [ fmg,fmgvt,fmgvp,Rpeak,Qlocatins,Slocatins,Qbegin,Send,Ppeak,Pbegin,Pend,Tbegin,Tend,Tpeak,NOofR,NOofT ] = detect(ecgsignal )
%参数说明：
%fmg:QRS的腐蚀结果
%fmgvt:T腐蚀结果
%fmgvp:P腐蚀结果
%Rpeak：R点位置，Qlocations:Q点位置，Slocations:S点位置，Qbegin：Q波起点，Send:S波终点
%Ppeak：p波，Pbegin：P波起点，Pend：P波终点
%Tpeak:T波，Tbegin:T波起点，Tend：T波终点
%ecgsignal:输入心电信号

% a=load('F:\examplejson.txt');
% b = resample(a,500,446);
%[tm,sig]=rdsamp('E:\PHY\CCChallenge\set-a\1045434',[]);
coefs = cwt(ecgsignal,38,'db3');
hi = ecgsignal;
%mask
%直方图定阈值
jidazhi = coefs(find(diff(sign(diff(coefs.*(coefs>0))))==-2)+1);
[counts,centers]=hist(jidazhi);
sum1 = 0;
for i =1:10
    sum1 =sum1 + counts(i)*centers(i);
end
mcent = sum1/sum(counts);
%直方图定阈值

%初步候选区
for i =1:length(coefs)
    if abs(coefs(i))>mcent;
        mp(i)=1;
    else
        mp(i)=0;
    end
end
%初步候选区

%腐蚀出新的候选区
s=ones(1,50);
fmg = enrosin( s,mp )+ones(1,length(mp));
vqrs=hi.*fmg';
clear s
%QRS dection
diffmg = diff(fmg);
st = find(diffmg==1);
en = find(diffmg(st(1)+2:end)==-1)+st(1)+1;
%腐蚀出新的候选区 st起点 en终点

%qrs方向判定
qrsflag = findjizhidianfangxiang(coefs(st(1):en(1)));
if qrsflag ==1
    disp('qrs波正向')
else
    vqrs = -vqrs;
    coefs = -coefs;
    disp('qrs波反向')
end
%qrs方向判定


lengthtouse = min(length(st),length(en));
fmginterval = en(1:lengthtouse)-st(1:lengthtouse);
Rpeak = [];
Qlocatins = [];
Slocatins = [];
Send = [];
Rs = [];
Re = [];
Qbegin = [];
%遍历每一个候选区域
for i = 1:lengthtouse
    if en(i)-st(i) < 0.5*mean(fmginterval) || en(i)-st(i) < 80
        for mj = st(i):en(i)+1
            fmg(mj)=0;
        end
        st(i)=0;
        en(i)=0;
        continue;
    end
    
    %定R
    tmparray = vqrs(st(i):en(i));
    tmparraypos = tmparray.*(tmparray>0);
    tmparrayneg = tmparray.*(tmparray<0);
    [maxpos,posloc] = max(tmparraypos);
    [mineg,negpoc] = min(tmparrayneg);
    Rpeak =[Rpeak, posloc(length(posloc))+st(i)-1];
    %定R
    
    
    %记录检出R波区域的起点和终点
    Rs = [Rs, st(i)];
    Re = [Re, en(i)];
   % tempq = vqrs(st(i):Rpeak(i));
   
   %定Q
     distancesq = zeros(1,1);
    for j = st(i):posloc(length(posloc))+st(i)-1
        distancesq(j-st(i)+1) = dist([st(i),vqrs(st(i))],[posloc(length(posloc))+st(i)-1,vqrs(posloc(length(posloc))+st(i)-1)],[j,vqrs(j)]);
    end
    [maxValue,location] = max(distancesq);
    Qlocatins = [Qlocatins,location(end)+st(i)-1];
    %定Q
    
    %定Q begin
    Qbeginarray = coefs(st(i):location(end)+st(i)-1);
    [minValues,minlocs] = min(Qbeginarray);
    Qbegin = [Qbegin, st(i)+minlocs(end)-1];
    %定Q begin
    
    %定S
    distancesS = zeros(1,1);
    for s = posloc(end)+st(i)-1:en(i)
        distancesS(s-(posloc(end)+st(i)-1)+1) = dist([posloc(end)+st(i)-1,vqrs(posloc(length(posloc))+st(i)-1)],[en(i),vqrs(en(i))],[s,vqrs(s)]);
    end
    [maxValues,locationS] = max(distancesS);
    tmps = vqrs( posloc(end)+st(i)-1:en(i));
    jixiaozhilocs = find(diff(sign(diff(tmps)))==2)+1;
    flags = xifiny(locationS(end),jixiaozhilocs);
    if flags == 1
        Slocatins =[Slocatins, locationS(end)+posloc(length(posloc))+st(i)-1-1];
        Sendarray = coefs(locationS(end)+posloc(end)+st(i)-1-1:en(i));
        %[minvaluesofS,minlocofs] = min(Sendarray);
        %Send =[Send,minlocofs(end)+locationS(length(locationS))+posloc(length(posloc))+st(i)-1-1-1];
    else
        jixiaoarray = vqrs( locationS(end)+posloc(end)+st(i)-1-1:en(i));
        jixiaozhilocfin = find(diff(newsign(diff(jixiaoarray)))==2)+1;
        if length(jixiaozhilocfin)==0
            Slocatins = [Slocatins,en(i)];
            %Send =[Send,en(i)];
        else
             Slocatins =[Slocatins, locationS(end)+posloc(end)+st(i)-1-1+jixiaozhilocfin(1)-1];
%              Sendarray = coefs( locationS(end)+posloc(length(posloc))+st(i)-1-1+jixiaozhilocfin(1)-1:en(i));
%              [minvaluesofS,minlocofs] = min(Sendarray);
%              Send =[Send,minlocofs(end)+locationS(end)+posloc(length(posloc))+st(i)-1-1+jixiaozhilocfin(1)-1-1];
        end
    end
    %定S： Slocatins
    
    %定Send
    tmpp = vqrs(Slocatins(end):en(i))';
    fdf = findlocations(vqrs(Qbegin(end)),tmpp);
    Send = [Send,Slocatins(length(Slocatins))+fdf(1)];
    %定Send
    
    clear tmparray
    clear distancesq
    clear distancesS
    clear Sendarray
    clear Qbeginarray
end

%因为不满足条件的st，en候选框被置零，所以这里取非零
 st = st((st~=0));
 en = en((en~=0));
 
%校验
%删除多检点
flag = 1;

while(flag==1)
    flag = 0;
    jianyanR = ones(1,length(Rpeak));
    meanRR = mean(diff(Rpeak));
    for i = 2:length(Rpeak)
        if Rpeak(i)-Rpeak(i-1)<0.6*meanRR;
            flag = flag + 1;
            if coefs(Rpeak(i))>coefs(Rpeak(i-1))
                jianyanR(i-1)=0;
            else
                jianyanR(i)=0;
            end
        end
    end
    if flag == 0
        break;
    end
   NewRR=jianyanR.*Rpeak;
   mm = Rs;
   nn = Re;
   Rs = Rs.*jianyanR;
   Re = Re.*jianyanR;
   newQlocatins =  Qlocatins.*jianyanR;
   newSlocatins = Slocatins.*jianyanR;
   newSend = Send.*jianyanR;
   newQbegins = Qbegin.*jianyanR;
   Rpeak = NewRR((NewRR~=0));
   Qlocatins = newQlocatins((newQlocatins~=0));
   Slocatins = newSlocatins((newSlocatins~=0));
   Send = newSend((newSend~=0));
   Qbegin = newQbegins((newQbegins~=0));
   slocs = find(Rs==0);
   for i = 1:length(slocs)
       for j = mm(slocs(i)):nn(slocs(i))
           fmg(j)=0;
       end
   end
   Rs =  Rs((Rs~=0));
   Re =  Re((Re~=0));
   flag = 1;
   clear jianyanR
   clear NewRR
   clear mm
   clear nn
end

%校验
NOofR = [];
for i = 1:length(Rpeak)
    NOofR(i) = i;
end


%T wave detection
%清除QRS对应的候选框所对应的信号
vt = hi-hi.*fmg';
%变换系数
coefvt = cwt(vt,98,'db3');
%变换系数

%直方图阈值
jidazhiofT = coefvt(find(diff(sign(diff(coefvt.*(coefvt>0))))==-2)+1);
[countvt,centervt]=hist(jidazhiofT);
sum2 = 0;
for i =1:length(centervt)
sum2 =sum2 + countvt(i)*centervt(i);
end
mcentvt = sum2/sum(countvt);
%直方图阈值

%初步候选区
for i =1:length(coefvt)
    if abs(coefvt(i))>mcentvt;
        mpvt(i)=1;
    else
        mpvt(i)=0;
    end
end
%初步候选区

%腐蚀出新的候选区
s=ones(1,50);
fmgvt = enrosin( s,mpvt )+ones(1,length(mpvt));
nifmg = 1 - fmg;
fmgvt = nifmg.*fmgvt;
clear s
ecgvt = vt.*fmgvt';
diffmgvt = diff(fmgvt);
st = find(diffmgvt==1);
en = find(diffmgvt(st(1)+2:end)==-1)+st(1)+1;
%腐蚀出新的候选区 增强版fmgvt = nifmg.*fmgvt;


%T波方向检测
Tflag = findjizhidianfangxiang(coefvt(st(1):en(1)));
if Tflag ==1
    disp('t波正向')
else
    ecgvt = -ecgvt;
    coefvt = -coefvt;
    disp('t波反向')
end
%T波方向检测 Tflag，若反了，就翻转信号


Tpeak = [];
Tbegin = [];
Tend = [];
NOofT = [];
i = 1;
%遍历T的每一个候选框
while(i<min(length(st),length(en)))
    
    %找候选框中最大值，早期没做翻转动作，做麻烦了
    tmparray = ecgvt(st(i):en(i));
    tmparraypos = tmparray.*(tmparray>0);
    tmparrayneg = tmparray.*(tmparray<0);
    [maxpos,posloc] = max(tmparraypos);
    [mineg,negpoc] = min(tmparrayneg);
    Tpeak = [Tpeak,posloc(length(posloc))+st(i)-1];
    %找候选框中最大值
    
    %定Tbegin
    Tbeginarray = ecgvt(st(i):posloc(length(posloc))+st(i)-1);
    [tmpbegin,beginloc] = min(Tbeginarray);
    Tbegin = [Tbegin,beginloc(length(beginloc))+st(i)-1];
    %定Tbegin
    
    %定T end
    Tendarray = ecgvt(posloc(length(posloc))+st(i)-1:en(i));
    [tmpend,endloc] = min(Tendarray);
    Tend = [Tend,endloc(length(endloc))+posloc(length(posloc))+st(i)-2];
    %定T end
    
    %以R为基准，去掉错误的候选框
    
    index1 =  findlocations( st(i),Rpeak );
    NOofT = [NOofT index1-1 ];
    index2 =  findlocations( st(i+1),Rpeak );
    if (index1 == index2)
        for j = st(i+1):en(i+1)+3
            fmgvt(j)=0;
        end
        i = i + 2;
    else
        i = i + 1;
    end
    %以R为基准，去掉错误的候选框
    
    clear tmparray
    clear tmparraypos
    clear tmparrayneg
    clear posloc
    clear negpoc Tbeginarray Tendarray endloc beginloc tmparray tmparraypos tmparrayneg
end
clear st
clear en
clear s


%P wvae detection
%置零T的候选框
ecgvts = hi.*fmgvt';
ecgvp = hi - ecgvts;
%置零T的候选框

%根据Tpeak和Qbegin继续置零信号
pidenticator = ones(1,length(hi));
for i = 1:min(length(Tpeak),length(Qlocatins))
    index = findlocations(Qlocatins(i),Tpeak);
    if index > length(Tpeak)
        for ds = Tpeak(end):length(ecgvp)
            ecgvp(ds) = 0;
        end
    else
        for j = Qlocatins(i)-5:Tpeak(index)
            ecgvp(j)=0;
            pidenticator(j)=0;
        end
    end
end
%根据Tpeak和Qbegin继续置零信号

%变换系数
coefvp = cwt(ecgvp ,78,'db3');
%变换系数

%直方图确定阈值，直方图统计变换系数的极大值
jidazhivp = coefvp(find(diff(sign(diff(coefvp.*(coefvp>0))))==-2)+1);
[countsvp,centersvp]=hist(jidazhivp);
sum3 = 0;
for i =1:10
    sum3 =sum3 + countsvp(i)*centersvp(i);
end
mcent = sum3/sum(countsvp);
%直方图确定阈值

%初步候选区域
for i =1:length(coefvp)
    if abs(coefvp(i))>mcent;
        mpvp(i)=1;
    else
        mpvp(i)=0;
    end
end
%初步候选区域

%腐蚀出新的候选区域
s=ones(1,50);
fmgvp = enrosin( s,mpvp )+ones(1,length(mpvp));
Ppeak = [];
Pbegin = [];
Pend = [];
NOofP = [];
diffmgvp = diff(fmgvp);
st = find(diffmgvp==1);
if(isempty(st))
    return;
end
en = find(diffmgvp(st(1)+2:end)==-1)+st(1)+1;
%腐蚀出新的候选区域

%判断p波正反，若反则倒过来
%P波方向判定
Pflag = findjizhidianfangxiang(coefvp(st(1):en(1)));
if Pflag ==1
    disp('p波正向')
else
    ecgvp = -ecgvp;
    coefvp = -coefvp;
    disp('p波反向')
end
%判断p波正反，若反则倒过来

i = 1;
%遍历每一个p波候选框
while(i<min(length(st),length(en)))
    %去掉异常的p波候选框
    if abs(st(i)-en(i))<10
        i = i + 1;
        continue;
    end
    %找到st(i)在Rpeak序列中的位置，若在R1和R2之间，则返回2
    index1 =  findlocations( st(i),Rpeak );
    index2 =  findlocations( st(i+1),Rpeak ); 
    if (index1 == index2)
        for j = st(i):en(i)+3
            fmgvp(j)=0;
        end
        i = i + 1;
        continue;
    end
    NOofP = [NOofP index1];
    %去掉异常的p波候选框
    
    %找到st(i)之后最近的Qbeigin点
    tselection = Qbegin.*(Qbegin>st(i));
    indexd = find(tselection);
    %找到st(i)之后最近的Qbeigin点
    
    %框定st(i)到Qbegin之前5个点之间的框为精确候选框
    tmparray = ecgvp(st(i):Qbegin(indexd(1))-5);
    %框定st(i)到Qbegin之前5个点之间的框为精确候选框
    
    %取候选框内非零点对应的横坐标
    %tmparraypos = tmparray.*(tmparray>0);
    %tmparrayneg = tmparray.*(tmparray<0);
    tpatt =find(tmparray);
    if length(tpatt)<10
        i = i+1;
        continue;
    end
     %取候选框内非零点对应的横坐标tpatt
   
    %取p波对应的CWT系数
    xxx = coefvp(st(i)+tpatt(2):st(i)+tpatt(end));
    %取p波对应的CWT系数
    
    %取极大值
    jidazhiloc = find(diff(sign(diff(xxx.*(xxx>0))))==-2)+1;
    if(length(jidazhiloc)==0)
        [ss,jidazhiloc]=max(xxx.*(xxx>0));
    end
    %取极大值位置jidazhiloc
    
    %再次精确候选区域
    houxuanp = ecgvp(jidazhiloc(end)+st(i)+tpatt(2)-4:jidazhiloc(end)+st(i)+tpatt(2)+4);
    [ps,locp] = max(houxuanp);
    %再次精确候选区域信号:houxuanp，locp为最大值的位置
    
    %P点波峰
    Ppeak = [Ppeak,jidazhiloc(end)+st(i)+tpatt(2)-4+locp(end)-1];
    
    %P点波峰
    
    %时域曲线和变换系数曲线左交点为Pbegin
    sss = coefvp(jidazhiloc(end)+st(i)+tpatt(2)-4+locp(end)-1:-1:st(i)+tpatt(2)-1);
    uuu = ecgvp(jidazhiloc(end)+st(i)+tpatt(2)-4+locp(end)-1:-1:st(i)+tpatt(2)-1);
    locbe = findcrosspoint(uuu,sss);
    Pbegin = [Pbegin,jidazhiloc(end)+st(i)+tpatt(2)-4+locp(end)-1-locbe+1];
    %时域曲线和变换系数曲线左交点为Pbegin
    
    %时域曲线和变换系数曲线右交点为Pend
    eee = coefvp(jidazhiloc(end)+st(i)+tpatt(2)-4+locp(end)-1:st(i)+tpatt(end)-1);
    vvv = ecgvp(jidazhiloc(end)+st(i)+tpatt(2)-4+locp(end)-1:st(i)+tpatt(end)-1);
    locen = findcrosspoint(vvv,eee);
    
    Pend = [Pend, jidazhiloc(end)+st(i)+tpatt(2)-4+locp(end)-1+locen-1];
    %时域曲线和变换系数曲线右交点为Pend
    
    
    i = i + 1;
    
    clear tmparray tmparraypos tmparrayneg posloc negpoc Pbeginarray Pendarray sss eee uuu vvv
end
figure
plot(hi)
hold on
plot(Ppeak,hi(Ppeak),'*')
hold on
plot(Rpeak,hi(Rpeak),'*')
hold on
plot(Tpeak,hi(Tpeak),'*')
hold on
plot(Qlocatins,hi(Qlocatins),'o')
hold on
plot(Slocatins,hi(Slocatins),'o')
hold on
plot(Qbegin,hi(Qbegin),'o')
hold on
plot(Send,hi(Send),'o')
hold on
plot(Tend,hi(Tend),'o')
hold on
plot(Tbegin,hi(Tbegin),'o')
hold on
plot(Pbegin,hi(Pbegin),'*')
hold on
plot(Pend,hi(Pend),'*')
end