
clear all
close all
clc

soubor = input ('Zadej nazev souboru: ');   % výběr EKG záznamu
load (soubor)   % načtení souboru
fvz = 500;  % určění vzorkovací frekvence - pro databázi CSE 500 Hz

svod = input ('Vyber svod (1 - 12): '); % výběr svodu

switch svod
    case 1
        x = x(1,:); % svod I
    case 2
        x = x(2,:); % svod II
    case 3
        x = x(3,:); % svod III
    case 4
        x = x(4,:); % svod V1
    case 5
        x = x(5,:); % svod V2
    case 6
        x = x(6,:); % svod V3
    case 7
        x = x(7,:); % svod V4
    case 8
        x = x(8,:); % svod V5
    case 9
        x = x(9,:); % svod V6
    case 10
        x = x(10,:);    % svod aVF
    case 11
        x = x(11,:);    % svod aVR
    case 12
        x = x(12,:);    % svod aVL
end

% realizace filtru pásmová propust
b = fir1(121, [10 20]/(fvz/2));  % nastavení parametrů filtru
z = filter( b, 1, x );  % filtrace signálu

% analytický filtr - vytvoření obálky signálu
analytitcky_signal = hilbert(z);    % vytvoření analytického signálu
obalka=abs(analytitcky_signal);     % absolutní hodnota z anal. signálu - obálka
Y = obalka.^2;   % umocneni obálky

% detekovani QRS komplexu
N = length( Y ); 
index = 1;      % priprava indexu pro cyklus WHILE
polohy = [];    % sem budeme ukladat polohy zacatku
hodnoty = [];   % proměná pro ukládání velikosti R-kmitů
velikost_okna = (0.1*fvz); % delka okna velikosti predpokladane sirky QRS (100ms) 
pocitadlo = 0; % promena pro pocitani QRS komplexu
prum_vzdal = (1*fvz);   % odhad průměrné vzdálenosti
prah = max(Y(1:2*fvz)); % počáteční nalezení maxima pro určení prahu
M=0.4;

while index<N-velikost_okna
     if Y(index)>prah*M     % rozhodování zda je vzorek větší než práh...
            cast = Y(index:index+velikost_okna); % ...pokud ano aplikuje se okno 
            
            if pocitadlo==0     % pokud se jedná o první komplex
                [hodnota,poloha] = max(cast); % hledání maxima v okně
                polohy = [polohy, index+poloha]; % uložení polohy maxima
                hodnoty=[hodnoty, hodnota]; % uložení hodnoty maxima
                index = index + velikost_okna; % zvýšení indexu pro další hledání o délku okna
                pocitadlo = pocitadlo + 1;  % započtení QRS komplexu
                prah = hodnota; % nová hodnota pro nastavení prahu
            else 
                [hodnota_pom,poloha_pom] = max(cast);

                if ((poloha_pom+index)-(polohy(length(polohy))))>(0.15*prum_vzdal)  % min. vzdál sousedních komplexů - větší než 15% prům. vzdálenosti
                    polohy = [polohy, poloha_pom+index]; % uložení polohy maxima
                    hodnoty=[hodnoty, hodnota_pom]; % uložení hodnoty maxima
                    index = index + round(prum_vzdal/2); %velikost_okna;
                    pocitadlo = pocitadlo + 1;  % započtení QRS komplexu
                    if hodnota_pom>0.5*prah     % podmínka proti snižování průměrné velikosti R-kmitu v důsledku detekce malého píku neodpovídajícímu R-kmitu
                        prah = (prah+hodnota_pom)/2;    % průměrná velikost R-kmitu
                    end
                    if pocitadlo==2
                        prum_vzdal = polohy(2)-polohy(1);   % výpočet prům. vzdálenosti prvních dvou R-kmitů
                    end
                    
                    prum_vzdal = (prum_vzdal+(polohy(pocitadlo)-polohy(pocitadlo-1)))/2;    % započtení QRS komplexu

                else
                    index = index + 1; % pokud není splněna podmínka posun na další vzorek
                end
            end
            
     else   % pokud není splněna podmínka, že vzorek je vyšší než práh
         if pocitadlo==0 % pokud není detekován žádný komplex... 
             index = index + 1; %...zvyšuje se index
         else
             if (index-(polohy(pocitadlo)))>(1.66*prum_vzdal)   % max. vzdál sousedních komplexů - menší než 166% prům. vzdálenosti
                 index=polohy(pocitadlo)+round(prum_vzdal/2);   % odkud se bude detekovat se sníženým prahem
                 if M>=0.2  % snížení koeficientu určující práh
                    M = M - 0.1;
                 else 
                    M=0.4;    % nastavení základních 40% z prům velikosti
                    prum_vzdal=prum_vzdal+prum_vzdal/10; % zvýšení prům. vzdálenosti
                 end 
             else
                index = index + 1;
             end             
         end
     end
end

disp (['V EKG signalu bylo detekovano ... ' num2str(pocitadlo) ' ... QRS komplexu.'])

% vykresleni
figure(1);
subplot(3,1,1);
plot(x); 
hold on
plot(polohy-(120/2),max(x), 'r*')
title ('Originální zaruseny EKG signal s detekovanymi QRS komplexy')
subplot(3,1,2); 
plot(z); 
hold on
plot(obalka,'r');
title ('Graf signalu EKG po pruchodu filtrem PP a obalka')
subplot(3,1,3); 
plot(Y);
hold on
plot(polohy,hodnoty, 'r*');
title ('Signal po umocneni')
