function rozmereni(poradove_cislo_signalu,VS,kolik)
% funkce provede podroblé rozměření jednoho cyklu signálu zadaného počátky
% a konci jednotlivých vln z databáze TESPRO a cyklu vybraného jako referenční
%
% poradove_cislo_signalu ... pořadové číslo signálu v databázi CSE
% VS...vstupní signál načtený z databáze CSE
% kolik...kolik má daný signál svodů
%
% příklad spuštění: rozmereni(20,x,12)
close all
clc

%% AUTOMATICKÉ PODROBNÉ ROZMĚŘENÍ JEDNOHO CYKLU ZA POMOCI BODŮ V DATABÁZI
%% TESPRO

TESPRO = xlsread('testro_excel.xls');  % načtení databáze TESPRO z MS Excel

pozorovani=poradove_cislo_signalu;  % nastavení parametru

% automatické načítámí bodů z databáze TESPRO
PonsetT = TESPRO(pozorovani,3);
PoffsetT = TESPRO(pozorovani,4);
QRSonsetT = TESPRO(pozorovani,5);
QRSoffsetT = TESPRO(pozorovani,6);
ToffsetT = TESPRO(pozorovani,7);

% parciální rozměření vzdálenosti významných bodů EKG záznamu získaných z databáze TESPRO
PQintervalT = PoffsetT - PonsetT;
PQsegmentT = QRSonsetT - PoffsetT;
delkaQRST = QRSoffsetT - QRSonsetT;
QTintervalT = ToffsetT - QRSonsetT;
STintervalT = ToffsetT - QRSoffsetT;

% hlavicka výpisu
vysledky(1,:)={'svod' 'P1' 'P2' 'PN' 'delkaP' 'delkaPN'...
    'Q' 'q' 'delkaQ' 'R' 'r' 'delkaR1' 'R2'...
    'r2' 'delkaR2' 'S' 's' 'delkaS' 'T1' 'T2'...
    'TN' 'delkaT' 'delkaTN'};


for svod=1:kolik
vstupni_signalT=VS(svod,TESPRO(pozorovani,2):(TESPRO(pozorovani,2)+TESPRO(pozorovani,7)));

%% určení hodnoty izolinie a její odečtení od vstupního signálu
hodnota1 = floor(PQsegmentT/2);  % určení délky daného cyklu
hodnota2 = PoffsetT+hodnota1;  % určení pozice středového bodu
hodnota_izolinie1 = floor(sum(vstupni_signalT(hodnota2-10:hodnota2+10))/21);  % váhovaná hodnota izilinie pro danný cyklus

ma=max(vstupni_signalT(hodnota2-10:hodnota2+10));
mi=min(vstupni_signalT(hodnota2-10:hodnota2+10));
hodnota_izolinieT=hodnota_izolinie1+mi-floor(ma/3); % hodnota izolinie váhovaná možným zkreslením 

signalT=vstupni_signalT-hodnota_izolinie1;
vyslekdy_pomocna=[];
%% detekce bodů
% T...označení prvku detekovaného z databáze TESPRO

[P1,P2,PN,delkaP,delkaPN,Q,q,delkaQ,R,r,delkaR1,R2,r2,delkaR2,S,s,delkaS,T1,T2,TN,delkaT,delkaTN]= detekce_bodu(signalT,hodnota_izolinieT,PonsetT,PoffsetT,QRSonsetT,QRSoffsetT,ToffsetT);

vysledky(svod+1,:)={svod P1 P2 PN delkaP delkaPN Q q delkaQ R r delkaR1 R2 r2 delkaR2 S s delkaS T1 T2 TN delkaT delkaTN};

end

% export výsledků
xlswrite('vysledky.xls',vysledky,'A:W',1);


%% NALEZENÍ A ROZMĚŘENÍ REPREZENTATICVNHO CYKLU V ZÁZNAMU EKG

%% funkce nalezení významných bodů EKG - Ing. Martin Vítek, PhD. 
[QRS,QRSonset,QRSoffset,T,Toffset,P,Ponset,Poffset] = EKGrozmereni_global(VS,'c',500,'n');

%% hlednání homogenního cyklu v seřazených po sobě jdoucích pozicích jednotlivých významných bodů nalezených v EKG signálu
QRSonset=floor(QRSonset);
QRSoffset=floor(QRSoffset);
Ponset=floor(Ponset);
Poffset=floor(Poffset);
Toffset=floor(Toffset);
pomocna=[QRSonset, QRSoffset, Toffset, Ponset, Poffset];
serazeny=sort(pomocna);  % pozice jednotlivých významých bodů seřazené vzestupně
b=zeros(1,length(serazeny));   % pomocný vektor detekovanýc bodů

PQinterval=[];
PQsegment=[];
delkaQRS=[];
QTinterval=[];
STinterval=[];

% cyklus prochází jednotlivé vektory a hledá shodu se seřazeným vektorem
% poloh významných bodů
c=0;
for aa=1:length(serazeny)
    if NajdiVPoli(serazeny(aa), Ponset)
        b(aa)=1;
    end
    if b(aa)==0
        if NajdiVPoli(serazeny(aa), Poffset)
            b(aa)=2;
        end
    end    
    if b(aa)==0
        if NajdiVPoli(serazeny(aa), QRSonset)
            b(aa)=3;
        end
    end    
    if b(aa)==0
        if NajdiVPoli(serazeny(aa), QRSoffset)
            b(aa)=4;
        end
    end
    
    if b(aa)==0
        if NajdiVPoli(serazeny(aa), Toffset)
            b(aa)=5;    % pokud jsou před Tofset všechny důležité dody, je provedeno parciální rozměření vzdáleností mezi body
            if (aa>4) && b(aa-4)==1 && b(aa-3)==2 && b(aa-2)==3 && b(aa-1)==4 
                c=c+1; 
                PQinterval(c) = serazeny(aa-3) - serazeny(aa-4);  % parciální rozměření pro analýzu mediánu vzdáleností
                PQsegment(c) = serazeny(aa-2) - serazeny(aa-3);
                delkaQRS(c) = serazeny(aa-1) - serazeny(aa-2);
                QTinterval(c)= serazeny(aa) - serazeny(aa-2);
                STinterval(c) = serazeny(aa) - serazeny(aa-1);  
                Ponset_p(c)=serazeny(aa-4);
                Poffset_p(c)=serazeny(aa-3);
                QRSonset_p(c)=serazeny(aa-2);
                QRSoffset_p(c)=serazeny(aa-1);
                Toffset_p(c)=serazeny(aa);
            end  
        end
    end
 end

%% určení mediánové vzdálenosti pro odstanění vlivu odchlujících hodnot vzdálednotí
PQint=median(PQinterval);
PQseg=median(PQsegment);
dQRS=median(delkaQRS);
QTint=median(QTinterval);
STint=median(STinterval);

%% nalezení reprezentativního cyklu 
pomocna=1000;

% hledání nejpodobněšího cyklu, za pomoci odchylek od ideální hodnoty
for ee=1:c 
   mi1=PQinterval(ee)-PQint;
   mi2=PQsegment(ee)-PQseg; 
   mi3=delkaQRS(ee)- dQRS; 
   mi4=QTinterval(ee)-QTint;
   mi5=STinterval(ee)-STint;
   soucet(ee)=sqrt(abs(mi1)^2+abs(mi2)^2+abs(mi3)^2+abs(mi4)^2+abs(mi5)^2);  % součet odchylek od medianove hodnoty 
   if soucet(ee)<=pomocna  % čím menší je odchylka, tím ideálnější by měl být cylkus
       pomocna=soucet(ee);
       Ponset=Ponset_p(ee);
       Poffset=Poffset_p(ee);
       QRSonset=QRSonset_p(ee);
       QRSoffset=QRSoffset_p(ee);
       Toffset=Toffset_p(ee);
   end
end

% podrobné rozměření vybraného cyklu signálu
for svod=1:kolik
    
signal=VS(svod,:);

% odečtení izolinie os signálu
hodnota1 = floor((Poffset-Ponset)/2);  % určení délky daného cyklu
hodnota2 = Poffset+hodnota1;  % určení pozice středového bodu
hodnota_izolinie = floor(sum(signal(hodnota2-5:hodnota2+5))/11);
signal=signal-hodnota_izolinie;

% provedení rozměření
[P1,P2,PN,delkaP,delkaPN,Q,q,delkaQ,R,r,delkaR1,R2,r2,delkaR2,S,s,delkaS,T1,T2,TN,delkaT,delkaTN]= detekce_bodu(signal,hodnota_izolinie,Ponset,Poffset,QRSonset,QRSoffset,Toffset);

vysledky_reprezentacni_cyklus(svod+1,:)={svod P1 P2 PN delkaP delkaPN Q q delkaQ R r delkaR1 R2 r2 delkaR2 S s delkaS T1 T2 TN delkaT delkaTN};
end

% halavička výstupu rozměření
vysledky(1,:)={'svod' 'P1' 'P2' 'PN' 'delkaP' 'delkaPN'...
    'Q' 'q' 'delkaQ' 'R' 'r' 'delkaR1' 'R2'...
    'r2' 'delkaR2' 'S' 's' 'delkaS' 'T1' 'T2'...
    'TN' 'delkaT' 'delkaTN'};

%% zhodnocení rozdílu mezi počátky a konci vln nalezenými a z databáze
% TESPRO
nalezene_cykly={' ','TESPRO','výběr'; 'Ponset',PonsetT,Ponset;'Pofset',PoffsetT,Poffset;...
    'QRSonset',QRSonsetT,QRSonset;'QRSoffset',QRSoffsetT,QRSoffset;'Toffset',ToffsetT,Toffset};

%% export výsledků do MS Excel
xlswrite('vysledky_reprezentativni_cyklus.xls',vysledky_reprezentacni_cyklus,'A:W',1);
xlswrite('nalezene_parametry_cyklu.xls',nalezene_cykly);

disp('hotovo'); % konec analýzy
end

function [P1,P2,PN,delkaP,delkaPN,Q,q,delkaQ,R,r,delkaR1,R2,r2,delkaR2,S,s,delkaS,T1,T2,TN,delkaT,delkaTN]= detekce_bodu(signal,hodnota_izolinie,Ponset,Poffset,QRSonset,QRSoffset,Toffset);
% funkce provede rozčlenení cyklu na jednotlivé funkční úseky a jejich rozměžení s následným vykreslením detekovaných bodů 

%% detekce bodů v komplexu QRS
QRS=signal(QRSonset:QRSoffset);
[Maxi1poziceQRS,Mini1poziceQRS,Maxi2poziceQRS,Mini2poziceQRS,delkaMaxima1QRS,delkaMaxima2QRS,delkaMinima1QRS,delkaMinima2QRS]=HledaniVrcholu(QRS,hodnota_izolinie); 

%% detekce vrcholu P a T vlny
signal_pomocny=[signal(1:QRSonset-1), linspace(signal(QRSonset),signal(QRSoffset),length(QRS)), signal(QRSoffset+1:end)];

% filrace 
h = fir1(40,30/(500/2),'low') ; % koeficient jmenovatele prenosove fuknce              
vysledek = filtfilt(h,1,signal_pomocny);

% P vlna
P=vysledek(Ponset:Poffset);
[Maxi1poziceP,Mini1poziceP,Maxi2poziceP,Mini2poziceP,delkaMaxima1P,delkaMaxima2P,delkaMinima1P,delkaMinima2P]=HledaniVrcholu(P,hodnota_izolinie);

% vlna T
T=vysledek(QRSoffset:Toffset);
[Maxi1poziceT,Mini1poziceT,Maxi2poziceT,Mini2poziceT,delkaMaxima1T,delkaMaxima2T,delkaMinima1T,delkaMinima2T]=HledaniVrcholu(T,hodnota_izolinie);

%% vykreslení signálu a získaných parametrů
signaly=signal(Ponset:Toffset);

figure
plot(signaly+hodnota_izolinie,'k');
hold on
plot(vysledek(Ponset:Toffset)+hodnota_izolinie,'g');
plot(linspace(hodnota_izolinie,hodnota_izolinie,length(signaly)),'b')
plot(linspace(hodnota_izolinie+25,hodnota_izolinie+25,length(signaly)),'r')
plot(linspace(hodnota_izolinie-25,hodnota_izolinie-25,length(signaly)),'r')

title('originální a filtrovaný signál')
legend('originální','filtrovaný','izolinie')
xlabel([' úsek signálu v rozmezí cyklu ',num2str(Ponset),' : ',num2str(Toffset), ' na globální ose']);
ylabel('U [\muV]');

% vykreslení detekovaných vrcholů
P=P+hodnota_izolinie;
QRS=QRS+hodnota_izolinie;
T=T+hodnota_izolinie;
vykresleni_maxim_y=[P(Maxi1poziceP+1),P(Maxi2poziceP+1),QRS(Maxi1poziceQRS+1),QRS(Maxi2poziceQRS+1),...
   T(Maxi1poziceT+1)];
vykresleni_maxim_x=[(Maxi1poziceP),Maxi2poziceP,Maxi1poziceQRS+QRSonset-Ponset,Maxi2poziceQRS+QRSonset-Ponset,...
   Maxi1poziceT+QRSoffset-Ponset,];
plot(vykresleni_maxim_x,vykresleni_maxim_y,'oc');

vykresleni_minim_y=[P(Mini1poziceP+1),QRS(Mini1poziceQRS+1),QRS(Mini2poziceQRS+1),...
   T(Mini1poziceT+1)];
vykresleni_minim_x=[(Mini1poziceP),(Mini1poziceQRS+QRSonset-Ponset),(Mini2poziceQRS+QRSonset-Ponset),...
   (Mini1poziceT+QRSoffset-Ponset)];
plot(vykresleni_minim_x,vykresleni_minim_y,'or');
hold off

%% hodnocení výsledků pro rozdělení do jednotlivých tříd  a převedení délek do časového měřítka
R=0;    r=0;        P1=0;         TN=0;     delkaP=0;   delkaQ=0;
R2=0;   r2=0;       P2=0;                   delkaPN=0;  delkaS=0;
Q=0;    q=0;        T1=0;                   delkaR1=0;  delkaT=0;
S=0;    s=0;                                delkaR2=0;  delkaTN=0;

% pro rozděnelí kmitů R a r
if Maxi1poziceQRS~=0 
    if QRS(Maxi1poziceQRS)>5
        R=abs(QRS(Maxi1poziceQRS));
    else
        r=abs(QRS(Maxi1poziceQRS));        
    end
    delkaR1=PrevodJednotek(delkaMaxima1QRS);
else R=0; r=0; delkaR1=0;
end

if Maxi2poziceQRS~=0 
    if QRS(Maxi2poziceQRS)>5
        R2=abs(QRS(Maxi2poziceQRS));
    else
        r2=abs(QRS(Maxi2poziceQRS));
    end
    delkaR2=PrevodJednotek(delkaMaxima2QRS);
else R2=0; r2=0; delkaR2=0;
end

% pro vrhol Q a q
if Mini1poziceQRS~=0 && ( Maxi1poziceQRS>Mini1poziceQRS || Maxi2poziceQRS>Mini1poziceQRS )
    if abs(QRS(Mini1poziceQRS))>5
        Q=abs(QRS(Mini1poziceQRS));
    else
        q=abs(QRS(Mini1poziceQRS));
    end
    delkaQ=PrevodJednotek(delkaMinima1QRS);
elseif Mini2poziceQRS~=0 && ( Maxi1poziceQRS>Mini2poziceQRS || Maxi2poziceQRS>Mini2poziceQRS )
    if abs(QRS(Mini2poziceQRS))>5
        Q=abs(QRS(Mini2poziceQRS));
    else
        q=abs(QRS(Mini2poziceQRS));
    end
    delkaQ=PrevodJednotek(delkaMinima2QRS);
else Q=0; q=0; delkaQ=0;
end

% pro vrhol S a s
if Mini1poziceQRS~=0 && (Maxi1poziceQRS<Mini1poziceQRS || Maxi2poziceQRS<Mini1poziceQRS)
    if abs(QRS(Mini1poziceQRS))>5
        S=abs(QRS(Mini1poziceQRS));
    else
        s=abs(QRS(Mini1poziceQRS));
    end
    delkaS=PrevodJednotek(delkaMinima1QRS);
elseif  Mini2poziceQRS~=0 && (Maxi1poziceQRS<Mini2poziceQRS || Maxi2poziceQRS<Mini2poziceQRS)
    if abs(QRS(Mini2poziceQRS))>5
        S=abs(QRS(Mini2poziceQRS));
    else
        s=abs(QRS(Mini2poziceQRS));
    end
    delkaS=PrevodJednotek(delkaMinima2QRS);
else S=0; s=0; delkaS=0;
end

%pro rozdělení P
if Maxi1poziceP~=0 
        P1=P(Maxi1poziceP);
        delkaP=PrevodJednotek(delkaMaxima1P);
else P1=0; delkaP=0;
end

if Maxi2poziceP~=0 
        P2=abs(P(Maxi2poziceP));
else P2=0; 
end

% pro rozdělení P negativni
if Mini1poziceP~=0 
        PN=abs(P(Mini1poziceP));
        delkaPN=PrevodJednotek(delkaMinima1P);
else PN=0; 
end

% pro rozdělení T vrcholů
if Maxi1poziceT~=0 
        T1=abs(T(Maxi1poziceT));
        delkaT=PrevodJednotek(delkaMaxima1T);
else T1=0;
end

% pro rozdělení T vrcholů
if Maxi2poziceT~=0 
        T2=T(Maxi2poziceT);
else T2=0;
end

% pro rozdělení T negativni
if Mini1poziceT~=0 
        TN=abs(T(Mini1poziceT));
        delkaTN=PrevodJednotek(delkaMinima1T);
else TN=0;
end

end

function [Maxi1pozice,Mini1pozice,Maxi2pozice,Mini2pozice,delkaMaxima1,delkaMaxima2,delkaMinima1,delkaMinima2]=HledaniVrcholu(signal,hodnota_izolinie)
% funkce provede podrovné rozměření vstupního úseku signálu

% průchody nulovou izolinií
p=0;
pozice = [];
for k=1:length(signal)-1
   if sign(signal(k+1)) ~= sign(signal(k))
      p=p+1;
      pozice(p)=k;  
   end
end
pozice=[0 pozice k k+1];  % počátek a konec zaručují paspoň nějaký průchod osou

% první diference
for i = 1:length(signal)-3
    if signal(i)==signal(i+1)
        if signal(i)==signal(i+2)
            diference(i) = signal(i+3) - signal(i);  % diference pokud by trend signálu pokračoval
        else
           diference(i) = signal(i+2) - signal(i);  % diference při lineárním trendu signálu
        end
    else
        diference(i) = signal(i+1) - signal(i);  % diference  
    end
end

MaxiIntervalu=[];
MinIntervalu=[]; 
p=0;
q=0;

% určení vrcholů podle podmínek analýzy průběhu funce
for i=2:length(diference)-1   
    if diference(i-1)>0 && diference(i+1)<0 
        if signal(i) > 25   % pokud by sestoupení na izilinii nebylo provedeno již dříve, 
            p=p+1;          % je možnost využít: hodnota_izoline+25
            MaxiIntervalu(p)=i;
        end
    elseif diference(i-1)<0 && diference(i+1)>0 
            if signal(i) < -25
                q=q+1;
                MinIntervalu(q)=i;
            end  
    end
end

%% hledání hlavních vrcholů 
po=0;
qo=0;
Mini1=0;   Maxi1pozice=0;
Maxi1=0;   Mini1pozice=0;
Mini2=0;   Maxi2pozice=0;
Maxi2=0;   Mini2pozice=0;      
PrahMinima=-25;
PrahMaxima=25;


for i=1:2
while po+1 <= length(MinIntervalu) && length(MinIntervalu)~=0
        po=po+1;
        if signal(MinIntervalu(po)) < Mini1    
            Mini1=signal(MinIntervalu(po));
            Mini1pozice=MinIntervalu(po)-1;  
        else
            if signal(MinIntervalu(po)) < Mini2 && signal(MinIntervalu(po))~=Mini1 && abs(Mini1pozice-MinIntervalu(po))>10            
                Mini2=signal(MinIntervalu(po));
                Mini2pozice = MinIntervalu(po)-1; % zohledňujeme vliv diference na detekci
            end
        end       
end

while qo+1 <= length(MaxiIntervalu) && length(MaxiIntervalu)~=0 
        qo=qo+1;
        if signal(MaxiIntervalu(qo)) > Maxi1
            Maxi1=signal(MaxiIntervalu(qo));
            Maxi1pozice=MaxiIntervalu(qo);
        else
            if signal(MaxiIntervalu(qo)) > Maxi2 && signal(MaxiIntervalu(qo))~=Maxi1 && abs(Maxi1pozice-MaxiIntervalu(qo))>10  % se zadaným pásem možného výskytu extrému
                  Maxi2 = signal(MaxiIntervalu(qo));  
                  Maxi2pozice=MaxiIntervalu(qo); % zohledňujeme vliv diference na detekci
            end
        end
end
end

% hledání vzdálenodtí kmitů mezi průchody nulovou linií jednotlivých hledanýchvrcholů
delkaMinima1=0;    delkaMaxima1=0;
delkaMinima2=0;    delkaMaxima2=0;

for n=1:length(pozice)-1
    if pozice(n)<Mini1pozice && Mini1pozice<pozice(n+1) 
        delkaMinima1=pozice(n+1)-pozice(n);
        delkaMinima1 = PrevodJednotek(delkaMinima1);
    end    
    if pozice(n)<Mini2pozice && Mini2pozice<pozice(n+1) 
        delkaMinima2=pozice(n+1)-pozice(n);
        delkaMinima2 = PrevodJednotek(delkaMinima2);
    end
    if pozice(n)<Maxi1pozice && Maxi1pozice<pozice(n+1)
        delkaMaxima1=pozice(n+1)-pozice(n);
        delkaMaxima1 = PrevodJednotek(delkaMaxima1);
    end
    if pozice(n)<Maxi2pozice && Maxi2pozice<pozice(n+1)
        delkaMaxima2=pozice(n+1)-pozice(n);
        delkaMaxima2 = PrevodJednotek(delkaMaxima2);
    end
end
end

function cas = PrevodJednotek(delka)
% funkce převede vzdálenost vzorků na časovou vzdálenost

d=(1/500)*delka;
cas=d*1000;

% možný převdod pro přímí výstup
%if cas>1000
%    cas=d;
%   disp([num2str(cas) 's']);
%else
%   disp([num2str(cas) 'ms']);
%end
end

function nalezeno = NajdiVPoli(prvek, pole)
% funkce detekuje shodu zkoumaného prvku s jedním prvkem z vektoru prvků
    nalezeno=false;
    ab=1;
    while (ab<=length(pole)) && (nalezeno == false)
        if prvek == pole(ab)
            nalezeno=true;
        end
        ab=ab+1;
    end
end

%--------------------------------------------------------------------------
% FUNKCE PRO NALEZENÍ SPECIFICKÝCH BODŮ V EKG SIGNÁLU
%--------------------------------------------------------------------------
% Ing. Martin Vítek - 30.3.2011
%--------------------------------------------------------------------------
function [QRS_glob,QRSon_glob,QRSoff_glob,T_glob,Toff_glob,P_glob,Pon_glob,Poff_glob] = EKGrozmereni_global(vstup_sig,typ,fvz,zobraz)
%%------ NAPOVEDA ---------------------------------------------------------

% vstupni parametry
% vstup_sig : vstupni signal ve formatu *.mat,
%             jde obecne o matici, kde radky jsou svody
% typ : 'c' nebo 'k' (c = clovek, k = kralik), typ EKG
% fvz : pouzita vzorkovaci frekvence v Hz
% zobraz : 'a' nebo 'n' (a = ano, n = ne), zobrazi pozice detekovanych bodu

% vystupni parametry
% pozice vyznamnych bodu spolecne pro vsechny svody
% postupne: QRS_glob,QRSon_glob,QRSoff_glob,T_glob,Toff_glob,P_glob,Pon_glob,Poff_glob

% priklad spusteni
% [QRS,QRSonset,QRSoffset,T,Toffset,P,Ponset,Poffset] =
% EKGrozmereni_global(x,'c',500,'a')

%% program
svodu = size(vstup_sig,1);

for i = 1:svodu
    
    [QRS,QRSon,QRSoff,T,Toff,P,Pon,Poff] = EKGrozmereni_local(vstup_sig(i,:),typ,fvz,'n');
    
    Pon_loc{i} = Pon;
    P_loc{i} = P;
    Poff_loc{i} = Poff;
    QRSon_loc{i} = QRSon;
    QRS_loc{i} = QRS;
    QRSoff_loc{i} = QRSoff;
    T_loc{i} = T;
    Toff_loc{i} = Toff;
end

tol = 50;

Pon_glob = vyber(Pon_loc,tol,svodu,1);
P_glob = vyber(P_loc,tol,svodu,0);
Poff_glob = vyber(Poff_loc,tol,svodu,2);
QRSon_glob = vyber(QRSon_loc,tol,svodu,3);
QRS_glob = vyber(QRS_loc,tol,svodu,0);
QRSoff_glob = vyber(QRSoff_loc,tol,svodu,4);
T_glob = vyber(T_loc,tol,svodu,0);
Toff_glob = vyber(Toff_loc,tol,svodu,5);

end

%% pomocné funkce

function [QRS,QRSonset,QRSoffset,T,Toffset,P,Ponset,Poffset] = EKGrozmereni_local(vstup_sig,typ,fvz,zobraz)
%% parametry programu / menit jen v pripade spatne detekce QRS
if      typ == 'c'
        pasma = [15 41];    % zakladni pasmo pro cloveka
        konst(1) = 1.55;    % prahova konstanta pro R
        konst(2) = 1.55;    % prahova konstanta pro TP
        konst(3) = 0.1;     % minimalni vzdalenost mezi dvema QRS komplexy v sekundach  
        
elseif  typ == 'k'
        pasma = [5 14];     % zakladni pasmo pro kralika
        konst(1) = 1.85;    % prahova konstanta pro R
        konst(2) = 1.85;    % prahova konstanta pro TP
        konst(3) = 0.03;    % minimalni vzdalenost mezi dvema QRS komplexy v sekundach  
else
        error('Zvolte typ bud c (clovek), nebo k (kralik)')
end

pasma = pasma*fvz/500;
vlnka = 'bior1.5';          % typ vlnky

parametry = {typ,fvz,konst};

%% transformace pro detekci QRS

band(1,:) = transformace(vstup_sig,pasma(1),vlnka);

%% pruchody

[nuly_R, extremy_R] = pruchody(band(1,:));
prah_R = konst(1)*std(band(1,:),1);

%% detekce QRS

[QRS,QRSpol,RRmed] = QRSdetect(nuly_R,extremy_R,prah_R,parametry);

[QRSonset] = QRSonsetdetect(nuly_R,extremy_R,prah_R,RRmed,QRSpol,band(1,:),parametry);

[QRSoffset] = QRSoffsetdetect(nuly_R,extremy_R,prah_R,RRmed,QRSpol,band(1,:));

%% odstraneni QRS komplexu

[inter_sig] = interpolace(vstup_sig,QRSonset,QRSoffset);

%% transformace pro detekci P a T

band(2,:) = transformace(inter_sig,pasma(2),vlnka);

%% pruchody

[nuly_TP, extremy_TP] = pruchody(band(2,:));
prah_TP = konst(2)*std(band(2,:),1);

%% detekce T

[T,Tpol] = Tdetect(nuly_TP,extremy_TP,prah_TP,QRSonset,QRSoffset,RRmed);

[Toffset] = Tenddetect(nuly_TP,extremy_TP,RRmed,QRSonset,T,Tpol,band(2,:));

%% detekce P

[P,Ppol] = Pdetect(nuly_TP,extremy_TP,prah_TP,QRSonset,QRSoffset,RRmed,1);

[Ponset] = Ponsetdetect(nuly_TP,extremy_TP,RRmed,QRSoffset,P,Ppol,band(2,:));

[P,Ppol] = Pdetect(nuly_TP,extremy_TP,prah_TP,QRSonset,QRSoffset,RRmed,2);

[Poffset] = Poffsetdetect(nuly_TP,extremy_TP,RRmed,QRSonset,P,Ppol,band(2,:));

%% zobrazeni

if zobraz == 'a'
scrsz = get(0,'ScreenSize');
width = 3*scrsz(3)/4;
height = scrsz(4)/2;
figure('Position',[(scrsz(3)-width)/2 (scrsz(4)-height)/2 width height])

plot(vstup_sig,'Color','b')
hold on;

Y = axis;
od = Y(4)-0.98*(abs(Y(4)-Y(3)));
do = Y(4)-0.02*(abs(Y(4)-Y(3)));

if ~isempty(QRSonset)
    line([QRSonset;QRSonset],[od;do],'Color','k','LineStyle','--','LineWidth',2)
end
if ~isempty(QRSoffset)
    line([QRSoffset;QRSoffset],[od;do],'Color','k','LineStyle','--','LineWidth',2)
end
if ~isempty(Ponset)
    line([Ponset;Ponset],[od;do],'Color','r','LineStyle','--','LineWidth',2)
end
if ~isempty(Poffset)
    line([Poffset;Poffset],[od;do],'Color','r','LineStyle','--','LineWidth',2)
end
if ~isempty(Toffset)
    line([Toffset;Toffset],[od;do],'Color','g','LineStyle','--','LineWidth',2)
end

end

end
function skupinky = seskup2(QRSm, vzdalenost)
% function SESKUP2 je zalozena na shlukove analyze
% 
%     skupinky = seskup2(QRSm, vzdalenost)
%        QRSm - bunkove pole, co bunka to svod
%             - kazda bunka obsahuje pozice nadetekovanych bodu v danym svodu
%        vzdalenost - minimalni vzdalenost (ve vzorcich), 
%                     pro rozdeleni na samostatne shluky
%        skupinky - bunkove pole, co bunka to bod
%                 - kazda bunka obsahuje blizke detekce ze vsech (nekterych) svodu (i vicenasobne)
%
%

% 2010-04-27, Jan Hrubes

m = sort(cell2mat(QRSm));
% cas = 0.1; % s (odhad casu mezi vlnama) 
% fvz = 500; % Hz
% vzdalenost = cas * fvz; % pocet vzorku
if length(m)>=2
    Y = pdist(m'); % vypocet vzdalenosti
    Z = linkage(Y); % shlukovani
    T = cluster(Z,'cutoff', vzdalenost,'criterion','distance')'; % rozdeleni do shluku
    [uT,oT] = unique(T,'first'); % zjisteni poctu shluku
    [~,pos] = sort(oT); % zjisteni poradi shluku
    skupinky = cell(1,length(uT)); % a nasazeni do bunkoveho pole
    for i = uT
        skupinky{i} = sort(m(T==pos(i)));
    end
else
    skupinky = {m};
end
end
function pozice = vyber2(rada,typ,dat)
% if length(rada{1})>=dat/2
% subskupinky = seskup2(rada,15);
% x = [];
% for i = 1:length(subskupinky)
%     x(i) = length(subskupinky{1,i});
% end
% y = find(max(x)==x);
% z = subskupinky{y};
% pozice = z(end);
% else
%     pozice =[];
% end
if length(rada)>=dat/2
    switch typ
        case 0
            pozice = median(rada);
        case 1
            %pozice = rada(round(1*length(rada)/4));
            pozice = median(rada);
        case 2
            % pozice = rada(round(3*length(rada)/4));
            pozice = median(rada);
        case 3
            % pozice = rada(round(1*length(rada)/4));
            pozice = median(rada);
        case 4
            % pozice = rada(round(3*length(rada)/4));
            pozice = median(rada);
        case 5
            % pozice = rada(round(5*length(rada)/9));
            pozice = median(rada);
    end
else
    pozice =[];
end
end
function [glo] = vyber(loc,tol,dat,typ)

    pom = seskup2(loc,tol);
    glo = [];
    for j = 1:length(pom)
        k = vyber2(pom{1,j},typ,dat);
        glo = [glo k];
    end
    %glo{i,1} =  pom2;


end
function [band] = transformace(vstup_sig,pasmo,vlnka)
%------ NAPOVEDA ----------------------------------------------------------

% funkce provadi rychlou CWT s vyuzitim konvoluce

%------ VYTVORENI FILTRU --------------------------------------------------

if  pasmo == 15 && strcmp(vlnka,'bior1.5')
    imp_char = [-7.60565599708887e-16,-1.85659367047585e-11,-4.45140058022962e-09,8.97251273307321e-08,-6.53831413699163e-07,-3.58815479937255e-06,1.33324726230374e-05,1.72677727375859e-05,-8.13034241591086e-05,-0.000146486039457173,-0.000158209850004651,0.000145613639643926,0.00107927082419663,0.00121693755225996,0.000132862465243511,-0.00200069372513942,-0.00517203607827470,-0.00658411822754586,-0.00590073070315290,-0.00683068849394393,-0.00668885498518011,0.000415044740183573,0.0138419336827928,0.0364237444725993,0.0550911264396067,0.0562097611661096,0.0464790430912854,0.0309942198664056,-0.00410399118501448,-0.0576403173228762,-0.135672219588872,-0.238644260830236,-0.303562256437038,-0.331380919672457,-0.340130389822616,-0.308264798217232,-0.247467189706591,-0.152279038815097,-0.00840605215956051,0.139800722880035,0.239634586842784,0.303380284390641,0.338943983730319,0.332878660953073,0.307990683019033,0.244958260784066,0.148316651348171,0.0643972223407127,0.00884067750687203,-0.0285538987166462,-0.0452044560796249,-0.0554507879029005,-0.0560798652309547,-0.0389972661784107,-0.0156749561683949,-0.00145430879176865,0.00624697761669970,0.00702698923766024,0.00588666001736941,0.00655645951063974,0.00538761083014166,0.00235368499262149,3.87916588733341e-05,-0.00115460182689205,-0.00116549674890849,-0.000213031569019047,0.000152489701845874,0.000145536489723662,9.46800031922444e-05,-1.18267930921447e-05,-1.69189414038318e-05,3.42459175725230e-06,1.09439007460723e-06,-1.76864434296067e-07,1.25969867651953e-08,2.08973938848534e-10];
elseif pasmo == 41 && strcmp(vlnka,'bior1.5')    
    imp_char = [-1.25742756689685e-15,-1.26881392334892e-13,-8.73922309057134e-12,-2.88645843074840e-10,2.97511947557746e-10,5.26727722869422e-09,-3.85925438746758e-08,-4.98158678247319e-08,1.38768100942093e-07,4.17137830605845e-07,-1.05643472903079e-07,-1.38865815307360e-06,-2.13690043263204e-06,-2.09377032496627e-06,-1.97023584679570e-06,2.39748851998860e-06,1.23486742999098e-05,1.94488878778454e-05,1.51784277959329e-05,2.52125057878260e-06,-1.79340650284360e-05,-4.84485880654922e-05,-7.97602833560836e-05,-9.47492497252575e-05,-8.66161406134489e-05,-8.45116725592607e-05,-9.83820978497574e-05,-0.000105001790035676,-4.77225496081356e-05,4.82895244543950e-05,0.000204124685821793,0.000433427530170947,0.000699583234529610,0.000834758415229956,0.000826072378541680,0.000712739837390467,0.000523355623062375,0.000238704928453196,-0.000142000883153561,-0.000582574208477941,-0.00111678631390131,-0.00170671966780629,-0.00255294069225966,-0.00321453155356668,-0.00371366380068779,-0.00402910409422038,-0.00406701277460427,-0.00377511700943481,-0.00356401342128039,-0.00355623518565073,-0.00367967556443057,-0.00398034003646380,-0.00454347381780396,-0.00474231981668435,-0.00411038585644414,-0.00270336142960599,-0.000953120375467694,0.00115742563169297,0.00386087844623606,0.00730828699222139,0.0112316524688249,0.0160393029098479,0.0218181037465998,0.0275732902310478,0.0317314832059487,0.0344470636794224,0.0356636029416735,0.0351534851762619,0.0332361373995689,0.0312156560190402,0.0288557167904198,0.0260202584551251,0.0227187280935774,0.0187504961487219,0.0126907456233335,0.00427873444741241,-0.00573890534619091,-0.0167496344934771,-0.0290146064179105,-0.0438370354098923,-0.0606339181897053,-0.0796755004643571,-0.0969583674084762,-0.124726758030298,-0.146417967009399,-0.164232250289230,-0.178378314311040,-0.189240541420099,-0.195938838583113,-0.200515494199318,-0.203806466540297,-0.205403214959108,-0.204504484593999,-0.201416165839947,-0.195558086652470,-0.186386706649258,-0.173906028449381,-0.159464825486535,-0.143307411178697,-0.123984058705305,-0.100204426861308,-0.0730368340364784,-0.0415330650695232,-0.00537037651685421,0.0319251040679111,0.0646461208167127,0.0929789916496002,0.117701454228338,0.138363090846845,0.155101757772216,0.170043987858716,0.183166383615305,0.193342817648000,0.200095236739690,0.203853723034840,0.205403880218856,0.204460255706795,0.201587133006526,0.197317641128100,0.191524183895262,0.181784529058429,0.168540447735490,0.151769347058811,0.126448539430783,0.108151635804853,0.0854955007483108,0.0656785442975078,0.0483923171104386,0.0328834482711410,0.0200122447567570,0.00871412569850927,-0.00156949938283585,-0.0105801216638661,-0.0172817797255345,-0.0217481857604045,-0.0251200124310823,-0.0281250518503231,-0.0305923260386450,-0.0326849045749701,-0.0346764028796881,-0.0357244909098922,-0.0349314750193850,-0.0326458100138954,-0.0288894088423525,-0.0235464705580775,-0.0175487556090831,-0.0124859514843156,-0.00834430845645248,-0.00476978783852583,-0.00183243434276977,0.000400430378799387,0.00224235561300576,0.00377814994283247,0.00464850489835909,0.00468088476109749,0.00411175320161333,0.00374607536966372,0.00357149738534964,0.00355166657327058,0.00368394086355591,0.00401466535109187,0.00406948640488193,0.00382800248803475,0.00336561452306181,0.00266532787104792,0.00200986318495701,0.00129236369301020,0.000717493087728296,0.000261858890739376,-0.000141753035798926,-0.000452970416298559,-0.000669020848803373,-0.000801000985317976,-0.000847059653789262,-0.000752557340832607,-0.000512580097557159,-0.000259523835261330,-8.50711939328661e-05,2.44545414162085e-05,9.48245104263474e-05,0.000104460291201567,8.64898764976676e-05,8.42755622739342e-05,9.43739310153585e-05,8.59290885945942e-05,5.83333604731448e-05,2.51867213902288e-05,2.53599482507695e-06,-1.25248535268426e-05,-1.91860853303255e-05,-1.54475077203722e-05,-4.54267790716980e-06,1.20319401187173e-06,2.27685253125787e-06,2.07732346221462e-06,1.76458893852726e-06,3.92036692800090e-07,-3.51672021222431e-07,-2.60408465211960e-07,3.64562732130277e-08,4.61779802462585e-08,2.74283112596608e-09,-3.35123314764449e-09,6.55357819646668e-10,-4.61460502880443e-11,-8.09747680208346e-13];
else
    impulz = [zeros(1,pasmo*5) 1 zeros(1,pasmo*5)];
    odezva = cwt(impulz,pasmo,vlnka);
    imp_char = odezva(find(odezva));
end
%------ KONVOLUCE ---------------------------------------------------------

zp = (length(imp_char)-1)/2;
zac(1:ceil(zp)) = vstup_sig(1);
kon(1:ceil(zp)) = vstup_sig(end);

vstup_sig = [zac vstup_sig kon];

band = conv(vstup_sig,imp_char);

%------ ELIMINACE ZPOZDENI ------------------------------------------------

band = band(ceil(zp)+1:end-fix(zp));
band = band(ceil(zp)+1:end-ceil(zp));

end
function [signal] = interpolace(vstup_sig,QRSonset,QRSoffset)

if min(QRSonset)<=max(QRSoffset)
    
QRSoffset(QRSoffset<QRSonset(1))=[];
QRSonset(QRSonset>QRSoffset(end))=[];

delka = min(length(QRSonset),length(QRSoffset));
QRSonset = QRSonset(1:delka);
QRSoffset = QRSoffset(1:delka);

if length(QRSoffset)~=length(QRSonset)
    error('nestejny pocet QRS on a off')
end



%------ NALEZENI A VYJMUTI QRS --------------------------------------------

buffer = [];
pozice = 1:length(vstup_sig);
for i = 1:length(QRSonset)
    buffer = [buffer pozice(QRSonset(i):QRSoffset(i))];
end

pozice(buffer) = [];

%------ INTERPOLACE QRS----------------------------------------------------

signal = interp1(pozice,vstup_sig(pozice),1:length(vstup_sig),'linear');

%--------------------------------------------------------------------------
else
    signal = vstup_sig;
end
end
function [u,v] = pruchody(x)

%------ NALEZENÍ PRŮCHODŮ NULOU -------------------------------------------

counter = 1;

for i = 1:length(x)-1
    if x(i)<=0 && x(i+1)> 0
       u(counter) = i;
       counter = counter + 1;
    elseif x(i)>=0 && x(i+1)< 0
       u(counter) = i;
       counter = counter + 1;
    end
end

% u = [1 u length(x)];

%------ NALEZENÍ MAXIM MEZI NULAMI ----------------------------------------

v = zeros(1,length(u)-1);

for i = 1:length(u)-1;
    if  x(u(i))>0
        v(i) = min(x(u(i):u(i+1)-1));
    elseif  x(u(i))==0
        if  x(u(i+1))< 0
            v(i) = min(x(u(i):u(i+1)-1));
        else
            v(i) = max(x(u(i):u(i+1)-1));
        end
    else
        v(i) = max(x(u(i):u(i+1)-1));
    end
end
end
function [QRS,QRSpol,RRmed] = QRSdetect(nuly,hladiny,prah,parametry)

%------ ZÁKLADNÍ ALGORITMUS -----------------------------------------------

counter = 1;
for i = 2:length(nuly)-1
    if  hladiny(i)>prah && hladiny(i-1)<(prah*(-1)) || ...  % 1 podminka
        hladiny(i-1)>prah && hladiny(i)<(prah*(-1));        % 2 podminka
        QRSpol(counter) = i;
        counter = counter+1;
    end
end
QRS = nuly(QRSpol);

%------ ELIMINACE MULTIVLN ------------------------------------------------

elim = diff(QRS);
kde = find(elim < parametry{3}(3)*parametry{2});
QRS(kde+1) = [];
QRSpol(kde+1) = [];

%------ RR INTERVAL -------------------------------------------------------

RR = diff(QRS);
RRmed = median(RR);

%------ DODATEČNÉ HLEDÁNÍ -------------------------------------------------

missing = find(RR>1.6*RRmed);
for j = 1:length(missing)
    missing_area = find(nuly>QRS(missing(j))+0.3*RR(missing(j)) & nuly<QRS(missing(j)+1)-0.3*RR(missing(j)));
    for i = missing_area(1):missing_area(end)
        if  abs(hladiny(i-1)-hladiny(i)) > 2*prah;
            QRS(end+1) = nuly(i);
            QRSpol(end+1) = i;
        end
    end
end

%------ SEŘAZENÍ KOMPLEXŮ A JEJICH POLOH ----------------------------------

QRS = sort(QRS);
QRSpol = sort(QRSpol);
end
function [QRSonset] = QRSonsetdetect(nuly,hladiny,prah,RRmed,QRSpol,bands,parametry)

%------ PREALOKACE --------------------------------------------------------

QRSonsetpol = zeros(1,length(QRSpol));
pozice = zeros(1,length(QRSpol));

%------ PARAMETRY ---------------------------------------------------------

const1 = 0.07;                  % oboje naraz
const2 = 11;                    % 
const3 = 59*parametry{2}/500;   % 

% const1 = 0.14;                  % 12
% const2 = 11;                    % 
% const3 = 59*parametry{2}/500;   % 

% const1 = 0.18;                  % xyz
% const2 = 11;
% const3 = 57*parametry{2}/500;

%------ ALGORITMUS --------------------------------------------------------

for i = 2:length(QRSpol)
    podminka = 0;
    while podminka ~= -1
        podminka = podminka + 1;
        xxx = nuly(QRSpol(i)) - nuly(QRSpol(i)- (podminka-1));
        xxxx = nuly(QRSpol(i)- (podminka-1)) - nuly(QRSpol(i)- podminka);
        xxxx2 = abs(hladiny(QRSpol(i)-podminka));
        if  xxxx2 < prah*const1 || ...              % 1
            xxx > RRmed/const2 || ...               % 2
            xxxx2/xxxx < prah/const3                % 3
            QRSonsetpol(i) = QRSpol(i) - (podminka - 1);
            podminka = -1;
        end
        if  podminka == -1
            x = find(abs(bands(1,nuly(QRSonsetpol(i)):nuly(QRSonsetpol(i)+1)))>const1*prah); % 0.15*
            d = diff(x);
            c = find(d==1);
            x(c+1) = [];
            pozice(i) = x(end);
        end
    end
end
QRSonsetpol(1) = [];
pozice(1) = [];
QRSonset = nuly(QRSonsetpol)+pozice;
end
function [QRSoffset] = QRSoffsetdetect(nuly,hladiny,prah,RRmed,QRSpol,bands)

%------ PREALOKACE --------------------------------------------------------

QRSoffsetpol = zeros(1,length(QRSpol)-1);
pozice = zeros(1,length(QRSpol)-1);

%------ PARAMETRY ---------------------------------------------------------

const1 = 0.18;        % oboje
const2 = 6;
const3 = 1;
const4 = 0.75;

% const1 = 0.23;      % 12
% const2 = 6;         
% const3 = 1;
% const4 = 0.6;
% 
% const1 = 0.21;      % xyz
% const2 = 5;
% const3 = 3;
% const4 = 0.75;
%------ ALGORITMUS --------------------------------------------------------

for i = 1:length(QRSpol)-1
    podminka = 0;
    while podminka ~= -1
        podminka = podminka + 1;
        xxxx = abs(hladiny(QRSpol(i)+podminka-1));
        xxxx2 = abs(hladiny(QRSpol(i)+podminka));
        xxx = nuly(QRSpol(i) + (podminka - 1)) - nuly(QRSpol(i));
        xxx2 = nuly(QRSpol(i) + podminka) - nuly(QRSpol(i));
        if  (xxxx < prah*const1 || xxx > RRmed/const2) && ...   
            (xxxx2 < const3*xxxx || xxxx2 < prah*const4 || xxx2 > RRmed/const2)
            QRSoffsetpol(i) = QRSpol(i) + (podminka - 1);
            podminka = -1;
        end
        if  podminka == -1
            x = find(abs(bands(1,nuly(QRSoffsetpol(i)-1):nuly(QRSoffsetpol(i))))>const1*prah); % 0.28*
        if  length(x)==1
            pozice(i)= x;
        else
            if x(1) ==1
               x(1) = [];
            end
            d = diff(x);
            c = find(d == 1);
            xx = x;
            xx(c) = [];
            if  length(xx) ==2
                pom = max(abs(bands(1,nuly(QRSoffsetpol(i)-1):nuly(QRSoffsetpol(i)-1)+xx(1))));
                pom2 = max(abs(bands(1,nuly(QRSoffsetpol(i)-1)+xx(1)+1:nuly(QRSoffsetpol(i)-1)+xx(2))));
                if  pom2 >= 3*pom
                    pozice(i) = xx(2);
                else
                    pozice(i) = xx(1);
                end
            else
                pozice(i) = xx(1);
            end
        end
        end
    end
end
QRSoffset = nuly(QRSoffsetpol-1)+pozice;
end
function [T,Tpol] = Tdetect(nuly,hladiny,prah,QRSonset,QRSoffset,RRmed)

%------ PREALOKACE --------------------------------------------------------

T = zeros(1,length(QRSonset));
Tpol = zeros(1,length(QRSonset));

%------ PARAMETRY ---------------------------------------------------------
const1 = 0.41;       % vzdalenost v procentech RRmed    % oboje
const2 = 0.14;      % vzdalenost v procentech RRmed
const3 = 0.1;       % krok prahu
prah = 0.07*prah;

% const1 = 0.40;       % vzdalenost v procentech RRmed    % 12
% const2 = 0.15;      % vzdalenost v procentech RRmed
% const3 = 0.1;       % krok prahu                        
% prah = 0.07*prah;   % prepocet prahu
% 
% const1 = 0.45;       % vzdalenost v procentech RRmed    % xyz
% const2 = 0.14;      % vzdalenost v procentech RRmed
% const3 = 0.1;       % krok prahu
% prah = 0.05*prah;   % prepocet prahu

%------ ALGORITMUS --------------------------------------------------------

for i = 1:length(QRSonset)
    waves = [];
    Vpol = [];
    delka = QRSonset(i) - QRSoffset(i);
    if  delka > const1*RRmed
        waves = find(nuly > QRSoffset(i)+const2*RRmed & nuly < QRSoffset(i)+const1*RRmed);
    end
  if ~isempty(waves)
      k = 0;
    while isempty(Vpol)
        counter = 1;
        k = k+const3;
    for j = 1:length(waves)
        if  abs(hladiny(waves(j)-1)) > prah*(1-k) && abs(hladiny(waves(j))) > prah*(1-k)
            Vpol(counter) = j;
            counter = counter + 1;
        end
    end
    
    end
    
    if  ~isempty(Vpol)
        T(i) = nuly(waves(Vpol(1)));
        Tpol(i) = waves(Vpol(1));
    end
  end
end

%------ NALEZENÍ T VLN NA NENULOVÝCH POZICÍCH -----------------------------

pom = find(T);
T = T(pom);
Tpol = Tpol(pom);
end
function [Tend] = Tenddetect(nuly,hladiny,RRmed,QRSonset,T,Tpol,band)

%------ PREALOKACE --------------------------------------------------------

Tendpol = zeros(1,length(Tpol));
pozice = zeros(1,length(Tpol));

%------ PARAMETRY ---------------------------------------------------------
const1 = 0.7;       % oboje
const2 = 0.25;
const3 = 0.15;

% const1 = 0.7;       % 12
% const2 = 0.25;
% const3 = 0.24;
% 
% const1 = 0.7;       % xyz
% const2 = 0.18;
% const3 = 0.35;
%------ ALGORITMUS --------------------------------------------------------

for i = 1:length(Tpol)
podminka = 1;
velikost = min(abs(hladiny(Tpol(i)-1:Tpol(i))));
    while podminka ~= -1
        podminka = podminka + 1;
        xxx = abs(hladiny(Tpol(i)+podminka-1));
        xx = find(QRSonset > T(i),1);
    if  xxx < velikost*const1 || nuly(Tpol(i)+podminka)> QRSonset(xx) || (nuly(Tpol(i)+podminka)-nuly(Tpol(i))>const2*RRmed)...
        || podminka == 2
        Tendpol(i) = Tpol(i) + (podminka-1);
        podminka = -1;
    end
    if  podminka == -1
        x = find(abs(band(nuly(Tendpol(i)-1):nuly(Tendpol(i))))>const3*abs(hladiny(Tendpol(i)-1)));
        if  length(x)==1
            pozice(i)= x;
        else
            if x(1) ==1
               x(1) = [];
            end
            d = diff(x);
            c = find(d == 1);
            xx = x;
            xx(c) = [];
            pozice(i) = xx(1);
        end
    end
    end
end

Tend = nuly(Tendpol-1)+pozice;
end
function [P,Ppol] = Pdetect(nuly,hladiny,prah,QRSonset,QRSoffset,RRmed,OnOff)

%------ PREALOKACE --------------------------------------------------------

P = zeros(1,length(QRSonset));
Ppol = zeros(1,length(QRSonset));

%------ PARAMETRY ---------------------------------------------------------
% oboje
const2 = 0.59;       % vzdalenost v procentech RRmed    % 59
const3 = 0.035;       % vzdalenost v procentech RRmed
const4 = 0.05;       % krok prahu
prah = 0.25*prah;    % prepocet prahu

% % 12
% const2 = 0.57;       % vzdalenost v procentech RRmed    % 0.54
% const3 = 0.04;       % vzdalenost v procentech RRmed    % 0.05
% const4 = 0.05;       % krok prahu
% prah = 0.25*prah;    % prepocet prahu                   % 0.26
% 
% % xyz
% const2 = 0.54;       % vzdalenost v procentech RRmed
% const3 = 0.05;       % vzdalenost v procentech RRmed
% const4 = 0.05;       % krok prahu
% prah = 0.26*prah;    % prepocet prahu
%------ ALGORITMUS --------------------------------------------------------

for i = 1:length(QRSonset)
    waves = [];
    Vpol = [];
    delka = QRSonset(i) - QRSoffset(i);
    if  delka > (const2+const3)*RRmed
        waves = find(nuly > QRSoffset(i)+const2*RRmed & nuly < QRSonset(i)-const3*RRmed);
    end
  if ~isempty(waves)
      k = 0;
    while isempty(Vpol)
        counter = 1;
        k = k+const4;
    for j = 1:length(waves)
        if  abs(hladiny(waves(j)-1)) > prah*(1-k) && abs(hladiny(waves(j))) > prah*(1-k)
            Vpol(counter) = j;
            counter = counter + 1;
        end
    end
    
    end
    
    if  ~isempty(Vpol)
        switch(OnOff)
            case 2
                P(i) = nuly(waves(Vpol(ceil(length(Vpol)/2))));
                Ppol(i) = waves(Vpol(ceil(length(Vpol)/2)));
            case 1
                P(i) = nuly(waves(Vpol(end)));
                Ppol(i) = waves(Vpol(end));
        end
    end
  end
end

%------ NALEZENÍ P VLN NA NENULOVÝCH POZICÍCH -----------------------------

pom = find(P);
P = P(pom);
Ppol = Ppol(pom);
end
function [Ponset] = Ponsetdetect(nuly,hladiny,RRmed,QRSoffset,P,Ppol,band)

%------ PREALOKACE --------------------------------------------------------

Ponsetpol = zeros(1,length(Ppol));
pozice = zeros(1,length(Ppol));

%------ PARAMETRY ---------------------------------------------------------

const1 = 0.70;      % oboje
const2 = 0.27;      %
const3 = 0.4;       %

% const1 = 0.70;      % 12
% const2 = 0.30;      %
% const3 = 0.63;      %
% 
% const1 = 0.70;      % xyz
% const2 = 0.30;      % 
% const3 = 0.61;      % 
%------ ALGORITMUS --------------------------------------------------------

for i = 1:length(Ppol)
podminka = 0;
velikost = min(abs(hladiny(Ppol(i)-1:Ppol(i))));
    while podminka ~= -1
        podminka = podminka + 1;
        xxx = abs(hladiny(Ppol(i)-podminka));
        xx = find(QRSoffset < P(i),1,'last');
    if  xxx < velikost*const1 || nuly(Ppol(i)-podminka) < QRSoffset(xx) || nuly(Ppol(i))-nuly(Ppol(i)-podminka)>const2*RRmed ...
        || podminka == 3
        Ponsetpol(i) = Ppol(i) - (podminka-1);
        podminka = -1;
    end
    if  podminka == -1
        x = find(abs(band(nuly(Ponsetpol(i)):nuly(Ponsetpol(i)+1)))>const3*abs(hladiny(Ponsetpol(i)))); % 0.15*
        d = diff(x);
        c = find(d==1);
        x(c+1) = [];
        pozice(i) = x(end);
    end
    end
end

Ponset = nuly(Ponsetpol)+pozice;

end
function [Poffset] = Poffsetdetect(nuly,hladiny,RRmed,QRSonset,P,Ppol,band)

%------ PREALOKACE --------------------------------------------------------

Poffsetpol = zeros(1,length(Ppol));
pozice = zeros(1,length(Ppol));

%------ PARAMETRY ---------------------------------------------------------
const1 = 0.75;        % oboje 0.75
const2 = 0.25;         %  0.25
const3 = 0.6;         % 0.6

% const1 = 0.75;      % 12
% const2 = 0.2;       % 
% const3 = 0.73;      %
% 
% const1 = 0.7;       % xyz
% const2 = 0.14;      % 
% const3 = 0.39;      % 

%------ ALGORITMUS --------------------------------------------------------

for i = 1:length(Ppol)
podminka = 1; % 0
velikost = min(abs(hladiny(Ppol(i)-1:Ppol(i))));
    while podminka ~= -1
        podminka = podminka + 1;
        xxx = abs(hladiny(Ppol(i)+podminka-1));
        xx = find(QRSonset > P(i),1);
    if  xxx < velikost*const1 || nuly(Ppol(i)+podminka)> QRSonset(xx)-2 || (nuly(Ppol(i)+podminka)-nuly(Ppol(i))>const2*RRmed)...
        || podminka == 1
        Poffsetpol(i) = Ppol(i) + (podminka-1);
        podminka = -1;
    end
    if  podminka == -1
        x = find(abs(band(nuly(Poffsetpol(i)-1):nuly(Poffsetpol(i))))>const3*abs(hladiny(Poffsetpol(i)-1))); % 0.15*
        if  length(x)==1
            pozice(i)= x;
        else
            if x(1) ==1
               x(1) = [];
            end
            d = diff(x);
            c = find(d == 1);
            xx = x;
            xx(c) = [];
            pozice(i) = xx(1);
        end
    end
    end
end

Poffset = nuly(Poffsetpol-1)+pozice;
end