%本程序實現基于小波變換的QRS波的檢測   %
%鄭秀玉  深圳大學德州儀器DSPs實驗室   %
%2007年7月2日                        %

clear all;
close all;
points=1000;       level=4;    sr=360; 
%讀入ECG信號
%load ecgdata.mat;
load ecgdata.txt;
%ecgdata=mydata;
plot(ecgdata(1:points));
grid on;axis tight;axis([0,points,-2,5]);

title('ECG信號');

swa=zeros(4,points);
swd=zeros(4,points);
signal=ecgdata(0*1000+1:1*1000);

%算小波系數和尺度系數
for i=1:points-3
  swa(1,i+3)=1/4*signal(i+3-2^0*0)+3/4*signal(i+3-2^0*1)+3/4*signal(i+3-2^0*2)+1/4*signal(i+3-2^0*3);
   swd(1,i+3)=-1/4*signal(i+3-2^0*0)-3/4*signal(i+3-2^0*1)+3/4*signal(i+3-2^0*2)+1/4*signal(i+3-2^0*3);
end
j=2;
while j<=level
   for i=1:points-24
     swa(j,i+24)=1/4*swa(j-1,i+24-2^(j-1)*0)+3/4*swa(j-1,i+24-2^(j-1)*1)+3/4*swa(j-1,i+24-2^(j-1)*2)+1/4*swa(j-1,i+24-2^(j-1)*3);
     swd(j,i+24)=-1/4*swa(j-1,i+24-2^(j-1)*0)-3/4*swa(j-1,i+24-2^(j-1)*1)+3/4*swa(j-1,i+24-2^(j-1)*2)+1/4*swa(j-1,i+24-2^(j-1)*3);
   end
   j=j+1;
end
%畫出原信號和尺度系數，小波系數
figure;
subplot(level,1,1); plot(ecgdata(1:points)); grid on;axis tight;
title('ECG信號及其在j=1,2,3,4尺度下的尺度系數及小波系數');
for i=1:level
    subplot(level+1,2,2*(i)+1);
    plot(swa(i,:)); axis tight;grid on;xlabel('time');
    ylabel(strcat('a   ',num2str(i)));
    subplot(level+1,2,2*(i)+2);
    plot(swd(i,:)); axis tight;grid on;
    ylabel(strcat('d   ',num2str(i)));
end

%畫出原圖及小波系數
figure;
subplot(level,1,1); plot(real(ecgdata(1:points)),'b'); grid on;axis tight;
title('ECG信號及其在j=1,2,3,4尺度下的小波系數');
for i=1:level
    subplot(level+1,1,i+1);
    plot(swd(i,:),'b'); axis tight;grid on;
    ylabel(strcat('d   ',num2str(i)));
end

%**************************************求正負極大值對*****************************************%
ddw=zeros(size(swd));
pddw=ddw;
nddw=ddw;
%小波系數的大于0的點
posw=swd.*(swd>0);
%斜率大于0
pdw=((posw(:,1:points-1)-posw(:,2:points))<0);
%正極大值點
pddw(:,2:points-1)=((pdw(:,1:points-2)-pdw(:,2:points-1))>0);
%小波系數小于0的點
negw=swd.*(swd<0);
ndw=((negw(:,1:points-1)-negw(:,2:points))>0);
%負極大值點
nddw(:,2:points-1)=((ndw(:,1:points-2)-ndw(:,2:points-1))>0);
%或運算
ddw=pddw|nddw;
ddw(:,1)=1;
ddw(:,points)=1;
%求出極值點的值,其他點置0
wpeak=ddw.*swd;
wpeak(:,1)=wpeak(:,1)+1e-10;
wpeak(:,points)=wpeak(:,points)+1e-10;

%畫出各尺度下極值點
figure;
for i=1:level
    subplot(level,1,i);
    plot(wpeak(i,:)); axis tight;grid on;
ylabel(strcat('j=   ',num2str(i)));
end
subplot(4,1,1);
title('ECG信號在j=1,2,3,4尺度下的小波系數的模極大值點');

interva2=zeros(1,points);
intervaqs=zeros(1,points);
Mj1=wpeak(1,:);
Mj4=wpeak(3,:);

%畫出尺度3極值點
figure;
plot (Mj4);
title('尺度3下小波系數的模極大值點');

posi=Mj4.*(Mj4>0);
%求正極大值的平均
thposi=(max(posi(1:round(points/4)))+max(posi(round(points/4):2*round(points/4)))+max(posi(2*round(points/4):3*round(points/4)))+max(posi(3*round(points/4):4*round(points/4))))/4;
posi=(posi>thposi/3);
nega=Mj4.*(Mj4<0);
%求負極大值的平均
thnega=(min(nega(1:round(points/4)))+min(nega(round(points/4):2*round(points/4)))+min(nega(2*round(points/4):3*round(points/4)))+min(nega(3*round(points/4):4*round(points/4))))/4;
nega=-1*(nega<thnega/4);
%找出非0點
interva=posi+nega;
loca=find(interva);
for i=1:length(loca)-1
    if abs(loca(i)-loca(i+1))<80
       diff(i)=interva(loca(i))-interva(loca(i+1));
    else
       diff(i)=0;
    end
end
%找出極值對
loca2=find(diff==-2);
%負極大值點
interva2(loca(loca2(1:length(loca2))))=interva(loca(loca2(1:length(loca2))));
%正極大值點
interva2(loca(loca2(1:length(loca2))+1))=interva(loca(loca2(1:length(loca2))+1));
intervaqs(1:points-10)=interva2(11:points);
count=zeros(1,1);
count2=zeros(1,1);
count3=zeros(1,1);
mark1=0;
mark2=0;
mark3=0;
i=1;
j=1;
Rnum=0;
%*************************求正負極值對過零，即R波峰值，并檢測出QRS波起點及終點*******************%
while i<points
    if interva2(i)==-1
       mark1=i;
       i=i+1;
       while(i<points&interva2(i)==0)
          i=i+1;
       end
       mark2=i;
%求極大值對的過零點
       mark3= round((abs(Mj4(mark2))*mark1+mark2*abs(Mj4(mark1)))/(abs(Mj4(mark2))+abs(Mj4(mark1))));
%R波極大值點
       R_result(j)=mark3-10;
       count(mark3-10)=1;
%求出QRS波起點
       kqs=mark3-10;
       markq=0;
     while (kqs>1)&&( markq< 3)
         if Mj1(kqs)~=0
            markq=markq+1;
         end
         kqs= kqs -1;
     end
  count2(kqs)=-1;
  
%求出QRS波終點  
  kqs=mark3-10;
  marks=0;
  while (kqs<points)&&( marks<2)
      if Mj1(kqs)~=0
         marks=marks+1;
      end
      kqs= kqs+1;
  end
  count3(kqs)=-1;
  i=i+60;
  j=j+1;
  Rnum=Rnum+1;
 end
i=i+1;
end
%************************刪除多檢點，補償漏檢點**************************%
num2=1;
while(num2~=0)
   num2=0;
%j=3,過零點
   R=find(count);
%過零點間隔
   R_R=R(2:length(R))-R(1:length(R)-1);
   RRmean=mean(R_R);
%當兩個R波間隔小于0.4RRmean時,去掉值小的R波
for i=2:length(R)
    if (R(i)-R(i-1))<=0.4*RRmean
        num2=num2+1;
        if signal(R(i))>signal(R(i-1))
            count(R(i-1))=0;
        else
            count(R(i))=0;
        end
    end
end
end

num1=2;
while(num1>0)
   num1=num1-1;
   R=find(count);
   R_R=R(2:length(R))-R(1:length(R)-1);
   RRmean=mean(R_R);
%當發現R波間隔大于1.6RRmean時,減小閾值,在這一段檢測R波
for i=2:length(R)
    if (R(i)-R(i-1))>1.6*RRmean
        Mjadjust=wpeak(4,R(i-1)+80:R(i)-80);
        points2=(R(i)-80)-(R(i-1)+80)+1;
%求正極大值點
        adjustposi=Mjadjust.*(Mjadjust>0);
        adjustposi=(adjustposi>thposi/4);
%求負極大值點
        adjustnega=Mjadjust.*(Mjadjust<0);
        adjustnega=-1*(adjustnega<thnega/5);
%或運算
        interva4=adjustposi+adjustnega;
%找出非0點
        loca3=find(interva4);
        diff2=interva4(loca3(1:length(loca3)-1))-interva4(loca3(2:length(loca3)));
%如果有極大值對,找出極大值對
        loca4=find(diff2==-2);
        interva3=zeros(points2,1)';
        for j=1:length(loca4)
           interva3(loca3(loca4(j)))=interva4(loca3(loca4(j)));
           interva3(loca3(loca4(j)+1))=interva4(loca3(loca4(j)+1));
        end
        mark4=0;
        mark5=0;
        mark6=0;
    while j<points2
         if interva3(j)==-1;
            mark4=j;
            j=j+1;
            while(j<points2&interva3(j)==0)
                 j=j+1;
            end
            mark5=j;
%求過零點
            mark6= round((abs(Mjadjust(mark5))*mark4+mark5*abs(Mjadjust(mark4)))/(abs(Mjadjust(mark5))+abs(Mjadjust(mark4))));
            count(R(i-1)+80+mark6-10)=1;
            j=j+60;
         end
         j=j+1;
     end
    end
 end
end
%畫出原圖及標出檢測結果
figure;
plot(ecgdata(0*1000+1:1*1000)),grid on,axis tight,axis([1,points,-2,5]);
title('ECG信號的R波峰值及QRS波波段');
hold on
plot(count,'r');
plot(count2,'k');
plot(count3,'k');
for i=1:Rnum
    if R_result(i)==0;
        break
    end
    plot(R_result(i),ecgdata(R_result(i)),'bo','MarkerSize',10,'MarkerEdgeColor','g');
end
hold off