function [QRS etime]=obalka_multi(signaly)

%%%VEVERKA VOJTĚCH 2015
global QRS etime
stime=tic;
%%%PŘEDZPRACOVÁNÍ%%%
%nastavení filtru
fvz=500; %vzorkovací frekvence
N=251; %řád filtru
svodu = size(signaly,1);

pp=fir1(N,[15 30]/(fvz));%natavení parametrů filtru
k = fvz/5;   %ošetření začátku a konce
for i=1:svodu  
signal=signaly(i,:);
zacatek(1:k) = signal(1);
konec(1:k) = signal(end);
signal = [zacatek signal konec];
signals(i,:)=signal;
end
%filtrace pásmovou propusti
for i=1:svodu  
signal=signals(i,:);
y(i,:) = filter( pp, 1, signal );
end

%eliminace zpoždění způsobeného zpožďovacími členy filtru FIR
delay = mean(grpdelay(pp));
for i=1:svodu  
pomzp=y(i,:);
pomzp(1:delay) = [];
pomzp=pomzp(100:5075);
fin(i,:)=pomzp;
end


 %výpočet obálky filtrovaného signálu
 
for i=1:svodu  
 yy=fin(i,:);
 yy = hilbert(yy);
 yy = abs(yy);
 yy=yy.^2;
 allobal(i,:)=yy;
end


%%%DETEKCE%%%
for i=1:svodu
poloha=detekuj(allobal(i,:));
polohy(1,i) = {poloha};
end
    QRS = [];
    shluky = shlukovka(polohy);
    for i = 1:length(shluky)
        if length(shluky{1,i})>svodu/2
           QRS = [QRS median(shluky{1,i})];
        end
    end
etime=toc(stime);



% signal=signaly(2,:);
% signal=signal(1:5000);
% y=allobal(2,:);
% % y=y(1:5000);
% QRSref=[255 708 1147 1394 1988 2392 2639 3188 3536 3778 4337 4765]
% QRS2 = QRS*0.002;
% QRSref=QRSref*0.002;
% t=0:1/fvz:length(signal)/fvz-1/fvz;
% t1=0:1/fvz:length(y)/fvz-1/fvz;
% figure(1)
% subplot(2,1,1)
% plot(t,signal)
% hold on
% line ([QRSref ;QRSref],[ylim],'color','red')
% subplot(2,1,2)
% plot(t1,y)
% hold on
% line ([QRS2 ;QRS2],[ylim],'color','red')
%     xlabel('Čas [s]');
%     ylabel('Napětí [uV]');

end


    


function [polohy]=detekuj(signal)
fvz=500;    
prah = mean(signal);  %vypocteni prahove hodnoty
d = length(signal);   %delka signalu  
krok = 1;             %krok pro while cyklus  
polohy = [];          %promena pro ulozeni poloh QRS  
okno = 120;           %delka okna pro vyhledávání maximální hodnoty
ref = 0.2 *fvz;       %refrakterni doba
signal=[signal zeros(1,okno)];    % rzsireni signalu o okno

while krok<d;
    
    if signal(krok)>prah; % kontrola nadprahovych hodnot
        
        kus = signal(krok:krok+okno); % vyber useku pro vyber maxima
        [hodnota,poloha] = max( kus ); 
        polohy = [polohy, krok+poloha-1]; % ulozeni polohy
        krok=polohy(end)+ref;
       if length(polohy)>=2 
           RR=polohy(end)-polohy(end-1); 
           RR1=RR*1.66;
       end          
       
  %%% Proběhne kontrola rekrakterni faze, pokud nepřijde QRS po refrakterni
  %%% fazi sníží se hodnota prahu
        if length(polohy)>=2
            if  krok==(polohy(end)+RR1)
                prah=0.5*prah;
                krok=polohy(end)+ref;
            end
        end
    end
    krok=krok+1;
end
end

function shluky = shlukovka(polohy)
mezery=30;
polohy2 = sort(cell2mat(polohy));
if length(polohy2)>=2
    vzdalenosti = pdist(polohy2'); 
    shlukovani = linkage(vzdalenosti); 
    shluk = cluster(shlukovani,'cutoff', mezery,'criterion','distance')'; % nastaveni shlukove analyzy
    [poradi,pocet] = unique(shluk,'first'); 
    [~,poradi2] = sort(pocet); 
    shluky = cell(1,length(poradi)); % ulozeni shluku
    for i = poradi
        shluky{i} = sort(polohy2(shluk==poradi2(i)));
    end
else
    shluky = {polohy2};
end
end

