clc; clear all; close all;
%%AUTOMATICKE ROZMERENI EKG SIGNALU --- EVA OLŠANSKÁ
%%(2013)
%% K fungování je třeba mít ve složce s tímto programem soubory všech signálů a referenční pozic
%% Ke spuštení celeho rozmereni stací stlacit run
%% Otevreni signalu z databaze CSE
for signal=1:125
    if  signal == 70|| signal == 67 % odstraneni nerozmerych signalu, signaly nerozmerene u vlny P, jsou
        signal= signal +1;            %odecitany az u vypoctu odchylek
    end
if length(num2str(signal))==1
    load (['CSE_MO1_12\MO1_00',num2str(signal),'_12.mat'])  
end
if length(num2str(signal))==2
    load (['CSE_MO1_12\MO1_0',num2str(signal),'_12.mat'])
end
 if length(num2str(signal))==3
    load (['CSE_MO1_12\MO1_',num2str(signal),'_12.mat'])   
 end

QRSbg_vysledny=[];QRSend_vysledny=[];

fvz = 500; %vzorkovaci frekvence
n=97; %pocet vzorku impulsni charakteritiky filtru 
if signal == 65 % osetreni extremnich nefyziolog. kmitu z duvodu selh.det.pro tento sig.
    x(:,1:400)= zeros;
end
%rozdeleni do jednotlivych svodu
I = x(1,:);II = x(2,:);III = x(3,:);avr = x(4,:);avf = x(5,:);avl = x(6,:);
V1 = x(7,:);V2 = x(8,:);V3 = x(9,:);V4 = x(10,:);V5 = x(11,:);V6 = x(12,:);
I_orig = x(1,:);II_orig = x(2,:);III_orig = x(3,:);avr_orig = x(4,:);avf_orig  = x(5,:);avl_orig  = x(6,:);
V1_orig  = x(7,:);V2_orig  = x(8,:);V3_orig  = x(9,:);V4_orig  = x(10,:);V5_orig  = x(11,:);V6_orig  = x(12,:);

%nastaveni pasmove propusti --->ECGPP
b=fir1(n,[(0.8/(fvz/2)) (24/(fvz/2))]);

I = filtfilt (b,1,I);II = filtfilt (b,1,II);III = filtfilt(b,1,III);avr = filtfilt (b,1,avr);avl = filtfilt (b,1,avl);avf = filtfilt (b,1,avf);
V1= filtfilt(b,1,V1);V2= filtfilt (b,1,V2);V3= filtfilt (b,1,V3);V4= filtfilt(b,1,V4);V5= filtfilt (b,1,V5);V6= filtfilt(b,1,V6);

%derivace --->ECGDER
Id = diff(I);IId = diff(II);IIId = diff(III);avld = diff(avl);avrd = diff(avr);avfd = diff(avf);
V1d = diff(V1);V2d = diff(V2);V3d = diff(V3);V4d = diff(V4);V5d = diff(V5);V6d = diff(V6);

%matice filtrovaneho a derivovaneho signalu 
ecgder_bezumoc = [Id;IId;IIId;avld;avrd;avfd;V1d;V2d;V3d;V4d;V5d;V6d];

%umocnění na druhou
Id2 = (Id).^2;IId2 = (IId).^2;IIId2 = (IIId).^2;avrd2 = (avrd).^2;avld2 = (avld).^2;avfd2 = (avfd).^2;
V1d2 = (V1d).^2;V2d2 = (V2d).^2;V3d2 = (V3d).^2;V4d2 = (V4d).^2;V5d2 = (V5d).^2;V6d2 = (V6d).^2;

%INTEGRACE - blok slouzici k vyhlazeni signalu
delka_okna=47; %fvz = 500...1s...95ms = 47,5N
zp=(delka_okna-1)/2;   % zpozdeni
pocet_svodu = 12;
ecgder=[Id2; IId2; IIId2; avrd2 ;avld2 ;avfd2 ;V1d2 ;V2d2;V3d2;V4d2;V5d2;V6d2];
d=length(Id2);
 for j = 1:pocet_svodu
for i=1:d-2*zp     %prumerovani
    ecgder_prum(j,i+zp)=(sum(ecgder(j,i:i+delka_okna-1))/delka_okna);
   end
 end
 
%DETEKCE R VLNY
cellR=cell(12,1);
velikostiR=cell(12,1);
prah=[];
for k =1:pocet_svodu
    
% % k vytvoreni prvniho prahu pro detekci
 R=(max(ecgder_prum(k,(1:1000)))+(max(ecgder_prum(k,1000:2000)))+(max(ecgder_prum(k,2000:3000)))+(max(ecgder_prum(k,3000:4000)))+(max(ecgder_prum(k,4000:end))))/5;
 
 prah1 = 0.5;
 prah =R*prah1; min_prah = 0.1 *R;

vektor_polohR = []; pomocna = [];ind=0;krok = 1; % inicializace
delka_okna = 120;                   % delka okna pro hledani maxim
ref = 0.20 *fvz;                     
pomocna=ecgder_prum(k,:);
pomocna = [pomocna , zeros(1,delka_okna)]; % prodlouzeni

%hledani nadprahove hodnoty
while krok<=length(pomocna)
      if pomocna (krok)> prah 
        kus = pomocna(krok:krok+delka_okna-1);   
        [hodnota,poloha] = max( kus );      % ...hledani maxima v dece okna
        vektor_polohR = [vektor_polohR, krok+poloha];          
        ind = 0;
        krok=krok+poloha+ref;               % přidání refrakterní doby

      end
       
      if length(vektor_polohR) >= 2 && krok > fvz  && ind < 1      % nekontroluje se pokud RR = 0 
         RR = vektor_polohR(end) - vektor_polohR (end-1);             % délka RR intervalu
          if krok-vektor_polohR(end)>=1.66*RR              % pokud se behem 1.66 nasobku nevyskytne R 
            prah=0.2*prah;                       % prah se snizi na 20% 
            if prah <= min_prah
            	prah = min_prah;
            end
            ind = ind + 1;                       % indikace ke snizeni prahu
            krok=vektor_polohR(end)+ref;                   
          end
      end
    krok=krok+1;
end
             
cellR{k,1}=vektor_polohR;

end


% ZPRESNENI R, v tomto bloku dochazi k hledani pruchodu nulovu hladinou
% obklopenou max a min v urcitem okne, v pripade nalezeni vice pruchodu se
% bere ten nejblizsi
cellR_upr=cell(12,1);vektor1=[];vektor2=[];
delka_okna=30;

for k = 1:pocet_svodu 
vysledne_polohyR =[];
for q = cellR{k,1}  %hledani bodu prochazejicich nulovou hladinou mezi maximem a minimem
    if   q+delka_okna < 4999 && q - delka_okna > 0
    okno = ecgder_bezumoc(k,q-delka_okna:q+delka_okna); 
[hodnota1,poloha1] = max(okno); 
[hodnota2,poloha2] = min(okno); 
poloha1 = poloha1+q-delka_okna;
poloha2 = poloha2+q-delka_okna;
   if poloha1 > poloha2

     for i= poloha2:poloha1

     if ((ecgder_bezumoc(k,i)>0 && ecgder_bezumoc(k,i+1)<0) || (ecgder_bezumoc(k,i)<0 && ecgder_bezumoc(k,i+1)>0))

    vysledna_poloha =(i); 
     end
     end
   end
    if poloha1 < poloha2
    
    for i= poloha1:poloha2

    if ((ecgder_bezumoc(k,i)>0 && ecgder_bezumoc(k,i+1)<0) || (ecgder_bezumoc(k,i)<0 && ecgder_bezumoc(k,i+1)>0))
    
 vysledna_poloha =(i); 
    end
    end
    end
   end         
vysledne_polohyR = [vysledne_polohyR vysledna_poloha];
end       
    
  cellR_upr{k,1}=vysledne_polohyR;
end
% --------------------------------------------------------------------------------------

%DETEKCE Q (bohuzel zbytecne slozity cyklus, nicmene funkci, tak jsem jej nemela potrebu menit)
%v okne pred vlnou R se detekuje pruchod nulou
delka_okna= 65;
hodnota = -80;
cellQ=cell(12,1);
velikostQ= [];
vektor_polohQ= [];
velikostiQ=cell(12,1);
j=1;pom= 0;
for k = 1:pocet_svodu 
    pomocna=cellR_upr{k,1};

    if length(pomocna)<=1 %osetreni v pripade selhani deteke R v celem svodu
        k= k+1;
        pomocna=cellR_upr{k,1}; 
    end
     if length(pomocna)<=1 %osetreni v pripade selhani deteke R v celem svodu
        k= k+1;
        pomocna=cellR_upr{k,1}; 
     end
      
     
 if ecgder_bezumoc(k,pomocna(1,2)-5)<0 
  
for q= cellR_upr{k,1}    % hledani pruchodu nulovou hladinou odpovidajicich vlne Q pred R vlnou
     if q -delka_okna >0
         okno= q-delka_okna:q-10;
          okno = fliplr(okno);
  for i = okno
     
    
    if ecgder_bezumoc(k,i)>hodnota 
        hodnota = ecgder_bezumoc(k,i);
    end
    if ecgder_bezumoc(k,i+1)>0 &&  pom <1        
    vektor_polohQ(j)=find(ecgder_bezumoc(k,:)==hodnota);
     
    velikostQ(j)=hodnota; 
     j=j+1; hodnota = -50;pom = 2;

    end   
  end
       pom=0;hodnota = -80;
     cellQ{k,1}=vektor_polohQ; 
     velikostiQ{k,1}= velikostQ;
end

 end
 j=1;
 end

  pomocna=cellR_upr{k,1};
 if ecgder_bezumoc(k,pomocna(1,2)-5)>0 
  
for q= cellR_upr{k,1}    
   hodnota=50 ;  
   if q-delka_okna > 0          
         okno= q-delka_okna:q-10;
          okno = fliplr(okno);
  for i = okno
 
    if ecgder_bezumoc(k,i)<hodnota && ecgder_bezumoc(k,i)>0
          hodnota = ecgder_bezumoc(k,i);
    end
    if ecgder_bezumoc(k,i+1)<0 &&  pom <1    &&  hodnota ~= 50;  
    vektor_polohQ(j)=find(ecgder_bezumoc(k,:)==hodnota);
     
    velikostQ(j)=hodnota; 
     j=j+1; hodnota = 50;pom = 2;

    end   
  end

       pom=0;hodnota = -80;
     cellQ{k,1}=vektor_polohQ;
       velikostiQ{k,1}= velikostQ;
   end
end
j=1;
 end
end
% --------------------------------------------------------------------
%DETEKCE S
%v okne po vlne R se detekuje pruchod nulou
delka_okna= 50;
cellS=cell(12,1);
vysledne_polohyS =[];
k=1;pom=0;refrakt =10;vysledna_poloha = [];

for k = 1:pocet_svodu 
vysledne_polohyS =[];

 for q = cellR_upr{k,1}  %opet hledani pruchodu nulovou hladinou, odpovidajicih vlne S, lezicich za vlnou R
    
    if   q + delka_okna+refrakt < 4998 
    okno = q+refrakt:q+refrakt+ delka_okna ;
    for i = okno
     if ((ecgder_bezumoc(k,i)>0 && ecgder_bezumoc(k,i+1)<0) || (ecgder_bezumoc(k,i)<0 && ecgder_bezumoc(k,i+1)>0)) && pom==0

    vysledna_poloha =(i ); pom=1 ;
    
     end  
    end     
     pom=0;
    end 
vysledne_polohyS = [vysledne_polohyS vysledna_poloha];vysledna_poloha = [];
  end    
    
  cellS{k,1}=vysledne_polohyS;
end
 
    
%detekce zacatku komplexu QRS-> nejdrive zjisteni maximalniho sklonu

cellQbg=cell(12,1);
velikostiQbg_pk=cell(12,1);

for k =1:pocet_svodu
hodnoty=[]; vektor_polohQbg=[];delka_okna= 15;
for q= cellQ{k,1}
    if q>delka_okna
    okno = ecgder_bezumoc(k,q-delka_okna:q); 
[hodnota,poloha] = max(abs(okno ));  % nalezeni maximalni strmosti pík pb
vektor_polohQbg= [vektor_polohQbg, (q-delka_okna)+poloha]; 
hodnoty = [hodnoty hodnota];
    end 
end
  cellQbg{k,1}=vektor_polohQbg;  velikostiQbg_pk{k,1}=hodnoty;
    end
%SAMOTNA DETEKCE ZACATKU KOMPLEXU QRS
% vydeleni pbg konstantou a nalezeni odpovidajici urovne signalu
cellQbg_upr=cell(12,1);
vektor_upr_polohQbg=[];
vektor_upr_polohyQbg=[];
 konst=1.41;% nastaveni konstanty pro detekci QRSb
for k =1:pocet_svodu
    
   for q= cellQbg{k,1}
       if q-15> 0 
       okno = ecgder_bezumoc(k, q-15:q-1); % nalezeni konecnych pozic vydelenim pb urcitou konstantou
       okno =fliplr(okno); %pretoceni vektoru
    vektor_upr_polohQbg =find(abs(okno) < abs(ecgder_bezumoc(k,q))/(konst), 1);%beru prvni prvek, ktery bude mensi nez dany |prah|
   vektor_upr_polohyQbg=[ vektor_upr_polohyQbg  (q-1-vektor_upr_polohQbg)];
    kontrola = isempty (vektor_upr_polohQbg) ;
   if kontrola == 1 ; 
        vektor_upr_polohyQbg=[ vektor_upr_polohyQbg  (q)];
   end
       end
       end
   cellQbg_upr{k,1}=vektor_upr_polohyQbg;
   vektor_upr_polohQbg=[];vektor_upr_polohyQbg=[];
 
end

% ------------------------------------------------------------------------

%DETEKCE QRS END
% hledani maximalni strmosti v signalu za vlnou S, pík pend
cellSend=cell(12,1);
velikostiSend_pk=cell(12,1);

for k =1:pocet_svodu
hodnoty=[]; vektor_polohSend=[];delka_okna= 25;

for q= cellS{k,1}
    if q+delka_okna < 4999
    okno = ecgder_bezumoc(k,q:q+delka_okna); 
[hodnota,poloha] = max(abs(okno ));   %nalezeni max strmosti
vektor_polohSend= [vektor_polohSend, q+poloha]; 
hodnoty = [hodnoty hodnota];
    end 
end
  cellSend{k,1}=vektor_polohSend;  velikostiSend_pk{k,1}=hodnoty;
end

%SAMOTNA DETEKCE KONCE KOMPLEXU QRS
%vydeleni pend konstantou a nalezeni odpovidajici urovne signalu
cellSend_upr=cell(12,1);
vektor_upr_polohSend=[];
vektor_upr_polohySend=[];
 konst=1.32 ;% nastaveni konstanty pro detekci QRSend
for k =1:pocet_svodu
    
   for q= cellSend{k,1}
       if q+11< 4999
       okno = ecgder_bezumoc(k, q:q+11);
    vektor_upr_polohSend =find(abs(okno) < abs(ecgder_bezumoc(k,q))/(konst), 1);%beru prvni prvek, ktery bude mensi nez dany |prah|
   vektor_upr_polohySend=[ vektor_upr_polohySend  (q+vektor_upr_polohSend)];
       end
     kontrola = isempty (vektor_upr_polohSend) ;
   if kontrola == 1 ; 
        vektor_upr_polohySend=[ vektor_upr_polohySend  (q)];
   end
   end
   cellSend_upr{k,1}=vektor_upr_polohySend;
   vektor_upr_polohSend=[];vektor_upr_polohySend=[];
 
end
% nacteni referencnich pozic
 load poziceEKGvCSE2

% %SHLUKOVA ANALÝZA PRO ZJIŠTĚNÍ VÝSLEDNÝCH VEKTORŮ QRS BEGINING

% 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)


vzdalenost = 30;
m = sort(cell2mat(cellQbg_upr'));
% cas = 0.1; % s (odhad casu mezi vlnami) 
% 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
pom2 = [];pom=skupinky;

    for i = 1:length(pom)
        if length(pom{1,i})>length(cellQbg_upr')/2
           pom2 = [pom2 round(median(pom{1,i}))];
        end
    end
QRSbg_vysledny = pom2;

%SHLUKOVA ANALÝZA PRO QRS END
vzdalenost = 30;
m = sort(cell2mat(cellSend_upr'));
% 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
pom2 = [];pom=skupinky;
   
    for i = 1:length(pom)
        if length(pom{1,i})>length(cellSend_upr')/2
           pom2 = [pom2 round(median(pom{1,i}))];
        end
    end
QRSend_vysledny = pom2;

QRSbg_vysledny_vsechny{signal,:} =QRSbg_vysledny;
QRSend_vysledny_vsechny{signal,:} =QRSend_vysledny;

% ------------------------------------------------------------------------
% filtrace signálu pomoci DP, k detekci vln P a T

poc_vz = 71;
mezni_dolni_kmitocet=14;
bb = fir1(poc_vz,mezni_dolni_kmitocet/(fvz/2),'low');
DERFI=[];

for i = 1:12
    DERFI(i,:) = filtfilt(bb,1,ecgder_bezumoc(i,:));
    
end

% zpresneni Q, kvuli filtraci DP
% dochazi ke zpresneni pozice vlny Q, tak aby odpovidala presne prochodu
% nulovou hladinou
cellQ_upr=cell(12,1);

delka_okna=50;
delka_okna2=8;
vektor1=[];
vektor2=[];


for k = 1:pocet_svodu 
vysledne_polohyQ =[]; vysledna_poloha =[];
for q = cellQ{k,1}  
    if   q+delka_okna < 4999 && q - delka_okna > 0
    okno = DERFI(k,q-delka_okna:q+delka_okna2); 
[hodnota1,poloha1] = max(okno); 
[hodnota2,poloha2] = min(okno); 
poloha1 = poloha1+q-delka_okna;
poloha2 = poloha2+q-delka_okna;
   if poloha1 > poloha2

     for i= poloha2:poloha1

     if ((DERFI(k,i)>0 && DERFI(k,i+1)<0) || (DERFI(k,i)<0 && DERFI(k,i+1)>0))

    vysledna_poloha =(i); 
     end
     end
   end
    if poloha1 < poloha2
    
    for i= poloha1:poloha2

    if ((DERFI(k,i)>0 && DERFI(k,i+1)<0) || (DERFI(k,i)<0 && DERFI(k,i+1)>0))
    
 vysledna_poloha =(i); 
    end
    end
    end
   end         
vysledne_polohyQ = [vysledne_polohyQ vysledna_poloha];vysledna_poloha=[];
end       
    
  cellQ_upr{k,1}=vysledne_polohyQ;vysledna_poloha=[];vysledne_polohyQ=[];
end

% --------------------------------------------------------------------------
% DETEKCE VLNY P

cellP=cell(12,1);velikostP= [];velikostiP=cell(12,1);vektor_polohP=[];
vektor1=[];vektor2=[];

% nastaveni delky oken a refrakternich fazi podle tepove frekvence
% nasledne nalezeni pruchodu nulou mezi maximem a minimem v okne
 if   length(cellR{10,:}) >= 17 
refrakt = 14; delka_okna=80;
end
  if length(cellR{10,:}) == 13|| length(cellR{10,:}) == 16|| length(cellR{10,:}) == 12|| length(cellR{10,:}) == 11|| length(cellR{10,:}) == 14|| length(cellR{10,:}) == 15
refrakt = 24; delka_okna=100;
  end
  if length(cellR{10,:}) <= 10
refrakt = 25; delka_okna=130;
 end

 
for k = 1:pocet_svodu 
vysledne_polohyP =[];
for q = cellQ_upr{k,1}  
    if   q - delka_okna > 0
    okno = DERFI(k,q-delka_okna:q-refrakt); 
[hodnota1,poloha1] = max(okno); 
[hodnota2,poloha2] = min(okno); 
poloha1 = poloha1+q-delka_okna;
poloha2 = poloha2+q-delka_okna;
   if poloha1 > poloha2

     for i= poloha2:poloha1

     if ((DERFI(k,i)>0 && DERFI(k,i+1)<0) || (DERFI(k,i)<0 && DERFI(k,i+1)>0))

    vysledna_poloha =(i); 
     end
     end
   end
    if poloha1 < poloha2
    
    for i= poloha1:poloha2

    if ((DERFI(k,i)>0 && DERFI(k,i+1)<0) || (DERFI(k,i)<0 && DERFI(k,i+1)>0))
    
 vysledna_poloha =(i); 
    end
    end
    end
   end         
vysledne_polohyP = [vysledne_polohyP vysledna_poloha];
end       
    
  cellP{k,1}=vysledne_polohyP;
end
 

%samotna detekce Pbg nejdrive zjisteni maximalniho sklonu
cellPbg=cell(12,1);velikostiPbg_pk=cell(12,1);vektor_polohPbg =[];

for k =1:pocet_svodu
hodnoty=[]; vektor_polohPbg=[];delka_okna= 45;
for q= cellP{k,1}
    if q>delka_okna
    okno = DERFI(k,q-delka_okna:q); 
[hodnota,poloha] = max(abs(okno ));  
vektor_polohPbg= [vektor_polohPbg, (q-delka_okna)+poloha]; 
hodnoty = [hodnoty hodnota];
    end 
end
  cellPbg{k,1}=vektor_polohPbg;  velikostiPbg_pk{k,1}=hodnoty;
    end
%SAMOTNA DETEKCE ZACATKU VLNY P
%vydeleni maximalni strmosti konstantou
cellPbg_upr=cell(12,1);
vektor_upr_polohPbg=[];
vektor_upr_polohyPbg=[];
 konst=1.3;% nastaveni konstanty pro detekci Pb
for k =1:pocet_svodu
    
   for q= cellPbg{k,1}
       if q-15> 0
       okno = DERFI(k, q-15:q-1);
       okno =fliplr(okno);% pretoceni vektoru
    vektor_upr_polohPbg =find(abs(okno) < abs(DERFI(k,q))/(konst), 1);%beru prvni prvek, ktery bude mensi nez dany |prah|
   vektor_upr_polohyPbg=[ vektor_upr_polohyPbg  (q-1-vektor_upr_polohPbg)];
     kontrola = isempty (vektor_upr_polohPbg) ;
   if kontrola == 1 ; 
       vektor_upr_polohyPbg=[ vektor_upr_polohyPbg  (q)];
   end
       end
       end
   cellPbg_upr{k,1}=vektor_upr_polohyPbg;
   vektor_upr_polohPbg=[];vektor_upr_polohyPbg=[];
 
end
% % %SHLUKOVA ANALÝZA PRO ZJIŠTĚNÍ VÝSLEDNÝCH VEKTORŮ P ZACATEK
% 
vzdalenost = 45;
m = sort(cell2mat(cellPbg_upr'));
% cas = 0.1; % s (odhad casu mezi vlnami) 
% 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
pom2 = [];pom=skupinky;

    for i = 1:length(pom)
        if length(pom{1,i})>length(cellPbg_upr')/2
           pom2 = [pom2 round(median(pom{1,i}))];
        end
    end
    
Pbg_vysledny = pom2;    
Pbg_vysledny_vsechny{signal,:} =Pbg_vysledny;

%---------------------------------------------------------------------------

%DETEKCE P END
%nalezeni maximalni strmosti v danem okne po vlne P
cellPend=cell(12,1);
velikostiPend_pk=cell(12,1);

for k =1:pocet_svodu
hodnoty=[]; vektor_polohPend=[];delka_okna= 25;

for q= cellP{k,1}
    okno = DERFI(k,q:q+delka_okna); 
[hodnota,poloha] = max(abs(okno ));  
vektor_polohPend= [vektor_polohPend, q+poloha]; 
hodnoty = [hodnoty hodnota];
end      
  cellPend{k,1}=vektor_polohPend;  velikostiPend_pk{k,1}=hodnoty;
end

%SAMOTNA DETEKCE KONCE VLNY P
%vyydeleni maximalni strmosti konstantou a nalezeni odpovidajiciho signalu
cellPend_upr=cell(12,1);
vektor_upr_polohPend=[];
vektor_upr_polohyPend=[];
 konst=1.4;% nastaveni konstanty pro detekci Pend
for k =1:pocet_svodu
    
   for q= cellPend{k,1}
       okno = DERFI(k, q:q+35);
    vektor_upr_polohPend =find(abs(okno) < abs(DERFI(k,q))/(konst), 1);%beru prvni prvek, ktery bude mensi nez dany |prah|
   vektor_upr_polohyPend=[ vektor_upr_polohyPend  (q+vektor_upr_polohPend)];
    kontrola = isempty (vektor_upr_polohPend) ;
   if kontrola == 1 ; 
       vektor_upr_polohyPend=[ vektor_upr_polohyPend  (q)];
   end
   end
   cellPend_upr{k,1}=vektor_upr_polohyPend;
   vektor_upr_polohPend=[];vektor_upr_polohyPend=[];
 
end
%SHLUKOVA ANALÝZA PRO P END
vzdalenost = 45;
m = sort(cell2mat(cellPend_upr'));
% 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
pom2 = [];pom=skupinky;
   
    for i = 1:length(pom)
        if length(pom{1,i})>length(cellPend_upr')/2
           pom2 = [pom2 round(median(pom{1,i}))];
        end
    end
    
Pend_vysledny = pom2;
Pend_vysledny_vsechny{signal,:} =Pend_vysledny;
  

% zpresneni S, kvuli filtraci DP
cellS_upr=cell(12,1);

delka_okna=60;
delka_okna2=10;
vektor1=[];
vektor2=[];


for k = 1:pocet_svodu 
vysledne_polohyS =[]; vysledna_poloha =[];
for q = cellS{k,1}  
       if   q+delka_okna < 4998 && q - delka_okna2 > 0
    okno = DERFI(k,q-delka_okna2:q+delka_okna); 
[hodnota1,poloha1] = max(okno); 
[hodnota2,poloha2] = min(okno); 
poloha1 = poloha1+q-delka_okna2;
poloha2 = poloha2+q-delka_okna2;
   if poloha1 > poloha2

     for i= poloha2:poloha1

     if ((DERFI(k,i)>0 && DERFI(k,i+1)<0) || (DERFI(k,i)<0 && DERFI(k,i+1)>0))

    vysledna_poloha =(i); 
     end
     end
   end
    if poloha1 < poloha2
    
    for i= poloha1:poloha2

    if ((DERFI(k,i)>0 && DERFI(k,i+1)<0) || (DERFI(k,i)<0 && DERFI(k,i+1)>0))
    
 vysledna_poloha =(i); 
    end
    end
    end
   end         
vysledne_polohyS = [vysledne_polohyS vysledna_poloha];
end       
    
  cellS_upr{k,1}=vysledne_polohyS;
end

%---------------------------------------------------------------------------
% DETEKCE VLNY T

cellT=cell(12,1);velikostT= [];velikostiT=cell(12,1);vektor_polohT=[];
vektor1=[];vektor2=[];
% nastaveni oken a refrakternich fazi podle tepove frekvence
 if   length(cellR{10,:}) >= 17 
refrakt = 45; delka_okna=70;
 end
 if length(cellR{10,:}) >= 11 && length(cellR{10,:}) <=16
refrakt = 53; delka_okna=95;
 end
  if length(cellR{10,:}) <= 10
refrakt = 55; delka_okna=130;
 end

 %hledani pruchodu nulou v okne
for k = 1:pocet_svodu 
vysledne_polohyT =[];
for q = cellS_upr{k,1}  
    if   q + delka_okna+refrakt < 4998
    okno = DERFI(k,q+refrakt:q+refrakt+ delka_okna); 
[hodnota1,poloha1] = max(okno); 
[hodnota2,poloha2] = min(okno); 
poloha1 = poloha1+q+refrakt;
poloha2 = poloha2+q+refrakt;
   if poloha1 > poloha2

     for i= poloha2:poloha1

     if ((DERFI(k,i)>0 && DERFI(k,i+1)<0) || (DERFI(k,i)<0 && DERFI(k,i+1)>0))

    vysledna_poloha =(i); 
     end
     end
   end
    if poloha1 < poloha2
    
    for i= poloha1:poloha2

    if ((DERFI(k,i)>0 && DERFI(k,i+1)<0) || (DERFI(k,i)<0 && DERFI(k,i+1)>0))
    
 vysledna_poloha =(i); 
    end
    end
    end
   end         
vysledne_polohyT = [vysledne_polohyT vysledna_poloha];vysledna_poloha = [];
end       
    
  cellT{k,1}=vysledne_polohyT;
end
 

%DETEKCE T END
cellTend=cell(12,1);
velikostiTend_pk=cell(12,1);
% hledani maximalni strmosti signalu
for k =1:pocet_svodu
hodnoty=[]; vektor_polohTend=[];delka_okna= 30;

for q= cellT{k,1}
    if q+delka_okna  < 4999 % osetreni presahnuti delky okna delku signalu
         if DERFI(k,q-2)<0 
    okno = DERFI(k,q:q+delka_okna); 
[hodnota,poloha] = max(okno );  
vektor_polohTend= [vektor_polohTend, q+poloha]; 
hodnoty = [hodnoty hodnota];
         end 
          if DERFI(k,q-2)>0 
    okno = DERFI(k,q:q+delka_okna); 
[hodnota,poloha] = min(okno );  
vektor_polohTend= [vektor_polohTend, q+poloha]; 
hodnoty = [hodnoty hodnota];
         end 
    end
end
  cellTend{k,1}=vektor_polohTend;  velikostiTend_pk{k,1}=hodnoty;
end

%SAMOTNA DETEKCE KONCE VLNY T
%vydeleni urovne signalu konstantou
cellTend_upr=cell(12,1);
vektor_upr_polohTend=[];
vektor_upr_polohyTend=[];
 konst=1.75;% nastaveni konstanty pro detekci Tend
for k =1:pocet_svodu
    
   for q= cellTend{k,1}
   if q+20 < 4999
       okno = DERFI(k, q:q+20);
    vektor_upr_polohTend =find(abs(okno) < abs(DERFI(k,q))/(konst), 1);%beru prvni prvek, ktery bude mensi nez dany |prah|
   vektor_upr_polohyTend=[ vektor_upr_polohyTend  (q+vektor_upr_polohTend)];
    kontrola = isempty (vektor_upr_polohTend) ;
   if kontrola == 1 ; 
       vektor_upr_polohyTend=[ vektor_upr_polohyTend  (q)];
   end
   end
   
   end
   cellTend_upr{k,1}=vektor_upr_polohyTend;
   vektor_upr_polohTend=[];vektor_upr_polohyTend=[];
 
end
%SHLUKOVA ANALÝZA PRO T END
vzdalenost = 65;
m = sort(cell2mat(cellTend_upr'));
% 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
pom2 = [];pom=skupinky;
   
    for i = 1:length(pom)
        if length(pom{1,i})>=length(cellTend_upr')/2
           pom2 = [pom2 round(median(pom{1,i}))];
        end
    end
    
Tend_vysledny = pom2;
Tend_vysledny_vsechny{signal,:} =Tend_vysledny;
    
    end


%------------ VYSLEDNE VYPOCTY--------------------------------
% urceni odhylky pro QRS zacatek
hledana_pozice=[];
skutecna_pozice=[];
for odch=1:125
    if  odch == 70|| odch == 67
        odch= odch +1;
    end 
    pomocna=QRSbg_vysledny_vsechny{odch,1};
   
   hledana_pozice = find(poziceEKGvCSE2(odch,3)-50 < pomocna & pomocna < poziceEKGvCSE2(odch,3)+40 );
   skutecna_pozice(odch)= pomocna(hledana_pozice);
   
   odchylka_QRS_zacatek(odch) = skutecna_pozice(odch) - poziceEKGvCSE2(odch,3);
end
%kontrola spravnych detekci pro vyrazeni z pocitani
% v pripade ze odchylka od referencni pozice je prilis vysoka, tak je
% signal zarazen do chybne detekce, nezapocitava se do smerodatne odchylky,
% ale snizi hodnotu senzitivity
chybna_detekce_zacatek=0;
odchylka_QRS_zacatek_upr =[];
casova_hranice=18;
for i = 1: length(odchylka_QRS_zacatek)
if abs(odchylka_QRS_zacatek(i))<= casova_hranice
    odchylka_QRS_zacatek_up=odchylka_QRS_zacatek(i);
     odchylka_QRS_zacatek_upr=[odchylka_QRS_zacatek_upr,odchylka_QRS_zacatek_up];
else abs(odchylka_QRS_zacatek(i))> casova_hranice;
    chybna_detekce_zacatek = chybna_detekce_zacatek+1;
end
end
odchylka_QRS_zacatek_upr(67)=[] ; %odstraneni nerozmerenych signalu
odchylka_QRS_zacatek_upr(69)=[] ;%odstraneni nerozmereneho signalu, signalu 70, jehoz cislo je zde 69, kvuli odecteni predchoziho

%senzitivita
Senzitivita_zacatek_QRS=length(odchylka_QRS_zacatek_upr)/123 %pocet sig. ktere se zapocitavaji do spravnych det/celk.pocet 

prumerna_odchylka_QRS_zacatek = mean(odchylka_QRS_zacatek_upr);
casova_prumerna_odchylka_QRS_zacatek = prumerna_odchylka_QRS_zacatek *(1/fvz);
smerodatna_zacatek = std( odchylka_QRS_zacatek_upr);
casova_odch_zacatek=smerodatna_zacatek *(1/fvz);
disp(['smerodatna odchylka QRS zacatek ' num2str(casova_odch_zacatek*1000) ' ms'])
disp(['prumerna odchylka QRS zacatek ' num2str(casova_prumerna_odchylka_QRS_zacatek*1000) ' ms'])

% urceni odhylky pro QRS konec
hledana_pozice=[];
skutecna_pozice=[];
for odch=1:125
    if  odch == 70|| odch == 67
        odch= odch +1;
    end 
    pomocna=QRSend_vysledny_vsechny{odch,1};
   
   hledana_pozice = find(poziceEKGvCSE2(odch,4)-50 < pomocna & pomocna < poziceEKGvCSE2(odch,4)+40 );
   skutecna_pozice(odch)= pomocna(hledana_pozice);
   
   odchylka_QRS_konec(odch) = skutecna_pozice(odch) - poziceEKGvCSE2(odch,4);
end
%kontrola spravnych detekci pro vyrazeni z pocitani
% v pripade ze odchylka od referencni pozice je prilis vysoka, tak je
% signal zarazen do chybne detekce, nezapocitava se do smerodatne odchylky,
% ale snizi hodnotu senzitivity
chybna_detekce=0;
odchylka_QRS_konec_upr =[];
casova_hranice=16;
for i = 1: length(odchylka_QRS_konec)
if abs(odchylka_QRS_konec(i))<= casova_hranice
    odchylka_QRS_konec_up=odchylka_QRS_konec(i);
     odchylka_QRS_konec_upr=[odchylka_QRS_konec_upr,odchylka_QRS_konec_up];
else abs(odchylka_QRS_konec(i))> casova_hranice;
    chybna_detekce = chybna_detekce+1;
end
end
odchylka_QRS_konec_upr(67)=[] ; %odstraneni nerozmerenych signalu
odchylka_QRS_konec_upr(69)=[] ;%odstraneni nerozmerenych signalu 70
%senzitivita
Senzitivita_konec_QRS=length(odchylka_QRS_konec_upr)/123
prumerna_odchylka = mean(odchylka_QRS_konec_upr);
casova_prumerna_odchylka = prumerna_odchylka *(1/fvz);
smerodatna_konec = std( odchylka_QRS_konec_upr);
casova_odch_konec=smerodatna_konec *(1/fvz);
disp(['smerodatna odchylka QRS konec ' num2str(casova_odch_konec*1000) ' ms'])
disp(['prumerna odchylka QRS konec ' num2str(casova_prumerna_odchylka*1000) ' ms'])


% % urceni odhylky pro P konec

hledana_pozice=[];
skutecna_pozice=[];
odchylka_P_konec=[];
for odch=1:125
    if   odch == 120 ||odch == 111 || odch == 100 ||odch == 92 ||odch == 89 ||odch == 57 ||odch == 52 ||odch == 45 ||odch == 18 || odch == 10 ||odch  == 70|| odch == 67 % odstraneni nerozmerych signalu
        odch= odch +1;
    end
    if odch == 93;
          odch= odch +1;
    end
    pomocna=Pend_vysledny_vsechny{odch,1};
   
   hledana_pozice = find(poziceEKGvCSE2(odch,2)-40 < pomocna & pomocna < poziceEKGvCSE2(odch,2)+40 );
   skutecna_pozice(odch)= pomocna(hledana_pozice);
   
   odchylka_P_konec(odch) = skutecna_pozice(odch) - poziceEKGvCSE2(odch,2);
end

odchylka_P_konec_pro_kont =odchylka_P_konec ;
odchylka_P_konec(10)=[] ; %odstraneni nerozmerenych signalu, cislo odstranovaneho signalu se pak vzdy snizi oproti signalu orig., dojde k posunu
odchylka_P_konec(17)=[] ;%odstraneni nerozmerenych signalu 
odchylka_P_konec(43)=[] ;%odstraneni nerozmerenych signalu 
odchylka_P_konec(49)=[] ; %odstraneni nerozmerenych signalu 
odchylka_P_konec(53)=[] ;%odstraneni nerozmerenych signalu
odchylka_P_konec(62)=[] ; %odstraneni nerozmerenych signalu
odchylka_P_konec(64)=[] ;%odstraneni nerozmerenych signalu
odchylka_P_konec(82)=[] ;%odstraneni nerozmerenych signal
odchylka_P_konec(64)=[] ;%odstraneni nerozmerenych signalu
odchylka_P_konec(83)=[] ; %odstraneni nerozmerenych signalu
odchylka_P_konec(90)=[] ;%odstraneni nerozmerenych signalu
odchylka_P_konec(100)=[] ;%odstraneni nerozmerenych signal
odchylka_P_konec(108)=[] ;%odstraneni nerozmerenych signalu

% kontrola spravnych detekci pro vyrazeni z pocitani
% v pripade ze odchylka od referencni pozice je prilis vysoka, tak je
% signal zarazen do chybne detekce, nezapocitava se do smerodatne odchylky,
% ale snizi hodnotu senzitivity
chybna_detekce_konec=0;
odchylka_P_konec_upr =[];
casova_hranice=25;

for i = 1: length(odchylka_P_konec)
if abs(odchylka_P_konec(i))<= casova_hranice
    odchylka_P_konec_up=odchylka_P_konec(i);
     odchylka_P_konec_upr=[odchylka_P_konec_upr,odchylka_P_konec_up];
else abs(odchylka_P_konec(i))> casova_hranice;
    chybna_detekce_konec = chybna_detekce_konec+1;
end
end
%senzitivita
Senzitivita_konec_P=length(odchylka_P_konec_upr)/112
prumerna_odchylka = mean(odchylka_P_konec_upr);
casova_prumerna_odchylka = prumerna_odchylka *(1/fvz);
smerodatna_konec = std( odchylka_P_konec_upr);
casova_odch_konec=smerodatna_konec *(1/fvz);
disp(['smerodatna odchylka P konec ' num2str(casova_odch_konec*1000) ' ms'])
disp(['prumerna odchylka P konec ' num2str(casova_prumerna_odchylka*1000) ' ms'])





% urceni odhylky pro P zacatek

hledana_pozice=[];
skutecna_pozice=[];
for odch=1:125
    if   odch == 120 ||odch == 111 || odch == 100 ||odch == 92 ||odch == 89 ||odch == 57 ||odch == 52 ||odch == 45 ||odch == 18 || odch == 10 ||odch  == 70|| odch == 67 % odstraneni nerozmerych signalu
        odch= odch +1;
    end
    if odch == 93;
          odch= odch +1;
    end
    pomocna=Pbg_vysledny_vsechny{odch,1};
   
   hledana_pozice = find(poziceEKGvCSE2(odch,1)-40 < pomocna & pomocna < poziceEKGvCSE2(odch,1)+40 );
   skutecna_pozice(odch)= pomocna(hledana_pozice);
   
   odchylka_P_zacatek(odch) = skutecna_pozice(odch) - poziceEKGvCSE2(odch,1);
end

odchylka_P_zacatek_pro_kont =odchylka_P_zacatek ;
odchylka_P_zacatek(10)=[] ; %odstraneni nerozmerenych signalu
odchylka_P_zacatek(17)=[] ;%odstraneni nerozmerenych signalu 
odchylka_P_zacatek(43)=[] ;%odstraneni nerozmerenych signalu 
odchylka_P_zacatek(49)=[] ; %odstraneni nerozmerenych signalu 
odchylka_P_zacatek(53)=[] ;%odstraneni nerozmerenych signalu
odchylka_P_zacatek(62)=[] ; %odstraneni nerozmerenych signalu
odchylka_P_zacatek(64)=[] ;%odstraneni nerozmerenych signalu
odchylka_P_zacatek(82)=[] ;%odstraneni nerozmerenych signal
odchylka_P_zacatek(64)=[] ;%odstraneni nerozmerenych signalu
odchylka_P_zacatek(83)=[] ; %odstraneni nerozmerenych signalu
odchylka_P_zacatek(90)=[] ;%odstraneni nerozmerenych signalu
odchylka_P_zacatek(100)=[] ;%odstraneni nerozmerenych signal
odchylka_P_zacatek(108)=[] ;%odstraneni nerozmerenych signal

% kontrola spravnych detekci pro vyrazeni z pocitani
% v pripade ze odchylka od referencni pozice je prilis vysoka, tak je
% signal zarazen do chybne detekce, nezapocitava se do smerodatne odchylky,
% ale snizi hodnotu senzitivity
chybna_detekce_zacatek=0;
odchylka_P_zacatek_upr =[];
casova_hranice=25;

for i = 1: length(odchylka_P_zacatek)
if abs(odchylka_P_zacatek(i))<= casova_hranice
    odchylka_P_zacatek_up=odchylka_P_zacatek(i);
     odchylka_P_zacatek_upr=[odchylka_P_zacatek_upr,odchylka_P_zacatek_up];
else abs(odchylka_P_zacatek(i))> casova_hranice;
    chybna_detekce_zacatek = chybna_detekce_zacatek+1;
end
end
% senzitivita
Senzitivita_zacatek_P=length(odchylka_P_zacatek_upr)/112

prumerna_odchylka_P_zacatek = mean(odchylka_P_zacatek_upr);
casova_prumerna_odchylka_P_zacatek = prumerna_odchylka_P_zacatek *(1/fvz);
smerodatna_zacatekP = std( odchylka_P_zacatek_upr);
casova_odch_zacatekP=smerodatna_zacatekP *(1/fvz);
disp(['smerodatna odchylka P zacatek ' num2str(casova_odch_zacatekP*1000) ' ms'])
disp(['prumerna odchylka P zacatek ' num2str(casova_prumerna_odchylka_P_zacatek*1000) ' ms'])


% % urceni odhylky pro T konec

hledana_pozice=[];
skutecna_pozice=[];
odchylka_T_konec=[];
for odch=1:125
    if   odch  == 70|| odch == 67 % odstraneni nerozmerych signalu
        odch= odch +1;
    end
   
    pomocna=Tend_vysledny_vsechny{odch,1};
   
   hledana_pozice = find(poziceEKGvCSE2(odch,5)-40 < pomocna & pomocna < poziceEKGvCSE2(odch,5)+40 );
   skutecna_pozice(odch)= pomocna(hledana_pozice);
   
   odchylka_T_konec(odch) = skutecna_pozice(odch) - poziceEKGvCSE2(odch,5);
end

odchylka_T_konec_pro_kont =odchylka_T_konec ;
odchylka_T_konec(67)=[] ; %odstraneni nerozmerenych signalu
odchylka_T_konec(69)=[] ;%odstraneni nerozmerenych signalu 1
% kontrola spravnych detekci pro vyrazeni z pocitani
% v pripade ze odchylka od referencni pozice je prilis vysoka, tak je
% signal zarazen do chybne detekce, nezapocitava se do smerodatne odchylky,
% ale snizi hodnotu senzitivity
chybna_detekce_konec=0;
odchylka_T_konec_upr =[];
casova_hranice=30;

for i = 1: length(odchylka_T_konec)
if abs(odchylka_T_konec(i))<= casova_hranice
    odchylka_T_konec_up=odchylka_T_konec(i);
     odchylka_T_konec_upr=[odchylka_T_konec_upr,odchylka_T_konec_up];
else abs(odchylka_T_konec(i))> casova_hranice;
    chybna_detekce_konec = chybna_detekce_konec+1;
end
end
%senzitivita
Senzitivita_konec_T=length(odchylka_T_konec_upr)/123
prumerna_odchylka = mean(odchylka_T_konec_upr);
casova_prumerna_odchylka = prumerna_odchylka *(1/fvz);
smerodatna_konec = std( odchylka_T_konec_upr);
casova_odch_konec=smerodatna_konec *(1/fvz);
disp(['smerodatna odchylka T konec ' num2str(casova_odch_konec*1000) ' ms'])
disp(['prumerna odchylka T konec ' num2str(casova_prumerna_odchylka*1000) ' ms'])



