%----------------------------------------------------------
%     DETEKTOR QRS A VÝPOČET TEPOVÉ FREKVENCE Z EKG SIGNÁLU
%----------------------------------------------------------
% Diplomová práce: Bezkontaktní detekce fyziologických parametrů z obrazových sekvencí.
% Tereza Bršlicová, ID 133958

clear all
close all
clc

load x.mat          % načtení dat vybraného dobrovolníka, za x dosadíme číslo dobrovolníka (1 až 10)
x = x(:,4);         % volba 4. kanálu, představujícího EKG
fvz = 1000;         % vzorkovací frekvence 

% Vykreslení původního signálu
% figure;
% plot(x);
% title('Původní EKG signál','FontWeight','bold');
% xlabel('Počet vzorků');
% ylabel('[mV]');

%% Realizace filtru pásmová propust

pp = fir1(120, [8 22]/(fvz/2));     % nastavení parametrů filtru pásmová prospust (PP)
f_pp = filter( pp, 1, x );          % filtrace PP

%% Umocnění signálu

y = f_pp.^2;          

%% Vyhlazeni umocněneho signálu filtrem dolní propust 

dp = fir1(70, 20/(fvz/2),'low');    % nastavení parametrů filtru dolní propust (DP)
yy = filter( dp, 1, y );            % filtrace DP

%% Detekce komplexu QRS 

% inicializace proměnných
index = 1;                      % priprava indexu pro while cyklus
N = length(yy);                 % délka signálu
polohy = [];                    % proměnná pro ukládání poloh R kmitů
hodnoty = [];                   % proměná pro ukládání hodnot R kmitů
okno = 100;                     % delka okna, ve kterem se budou vyhledávat maxima
pr_vzd = (1*fvz);               % odhad průměrné vzdálenosti
prah = (max(yy(1:fvz)));        % počáteční nalezení maxima pro určení prahu
pocetQRS = 0;                   % promena pro ukládání počtu QRS komplexu



while index < N-okno
    
     if yy(index) > prah*0.4     % podmínka pro rozhodování, jestli prvek větší než práh
                                 
            okno1 = yy(index:index + okno); % je-li větší, použije se okno 
            
            if pocetQRS == 0                        % první ORS komplex
                [hodnota,poloha] = max(okno1);      % hledání maxima v okně
                polohy = [polohy, index+poloha];    % uložení polohy maxima
                hodnoty = [hodnoty, hodnota];       % uložení hodnoty maxima
                index = index + okno;               % zvýšení indexu o délku okna pro další vyhledávání 
                pocetQRS = pocetQRS + 1;            % QRS komplex je započítán
                prah = hodnota;                     % prah je nastaven na novou hodnotu
            else 
                [hodnota_m,poloha_m] = max(okno1);  

                if ((poloha_m+index)-(polohy(length(polohy))))>(0.15*pr_vzd)  % minimální vzdálenost sousedních komplexů QRS je větší než 15% průměrné vzdálenosti
                    polohy = [polohy, poloha_m+index];                        % uložení polohy maxima
                    hodnoty = [hodnoty, hodnota_m];                           % uložení hodnoty maxima
                    index = index + round(pr_vzd/2);       % velikost okna
                    pocetQRS = pocetQRS + 1;               % QRS komplex je započítán
                    
                    if hodnota_m>0.5*prah       % podmínka, která ošetřuje snižování průměrné velikosti R kmitu, aby nedocházelo k detekci menšího kmitu neodpovídajícímu R kmitu
                        prah = (prah+hodnota_m)/2;          % průměrná velikost R kmitu
                    end
                    
                    if pocetQRS==2
                        pr_vzd = polohy(2)-polohy(1);       % výpočet průměrné vzdálenosti prvních dvou R kmitů
                    end
                    
                    pr_vzd = (pr_vzd+(polohy(pocetQRS)-polohy(pocetQRS-1)))/2;    % QRS komplex započítán

                else
                    index = index + 1; % není-li splněna podmínka, nastane posun na další vzorek
                end
            end
            
     else                           % není-li splněna podmínka, že bude vzorek vyšší než práh
         if pocetQRS==0             % není-li detekován žádný komplex, zvýší se index
             index = index + 1;    
         else
             if (index-(polohy(pocetQRS)))>(1.66*pr_vzd)   % maximální vzdálenost sousedních komplexů - menší než 166% průměrné vzdálenosti
                 index=polohy(pocetQRS)+round(pr_vzd/2);   % detekce se sníženým prahem
             else
                index = index + 1;
             end 
             
         end
     end
end
                        

% Ošetření zpoždění
zp1=(length(pp)-1)/2;
zp2=(length(dp)-1)/2;

zpoz=zp1+zp2;

%% Vykresleni kroků detekce QRS komplexů

figure;
subplot(4,1,1);
plot(x,'k');
xlabel('Počet vzorků');
ylabel('[mV]');
hold on
stem( polohy-zpoz, 0.5*ones(1,length(polohy)), 'b' );
title ('Původní EKG signál s detekovanými QRS komplexy','FontWeight','bold');
subplot(4,1,2); 
plot(f_pp,'k');
xlabel('Počet vzorků');
ylabel('[mV]');
title ('Signál EKG po filtraci PP','FontWeight','bold');
subplot(4,1,3);
plot(y,'k');
xlabel('Počet vzorků');
ylabel('[mV]');
title ('Signál po umocnění','FontWeight','bold');
subplot(4,1,4); 
plot(yy,'k');
xlabel('Počet vzorků');
ylabel('[mV]');
hold on
plot(polohy,hodnoty,'b*');
title ('Umocněný a filtrovaný signál s detekovanými QRS komplexy','FontWeight','bold');


%% vypocet tepove frekvence

% vytvoření vektoru vzdáleností
vekt_vzd=zeros(1,length(polohy)-1);

% výpočet vektoru vzdáleností jednotlivých R kmitů
for i=2:length(polohy)
   
    vekt_vzd(i-1)=polohy(i) - polohy(i-1); 
end

% výpočet tepové frekvence podle poloh R kmitů v signálu
tep = 60/(mean(vekt_vzd)/1000);

%% Zobrazení počtu komplexů a výsledné hodnoty tepové frekvence
disp(['Počet nalezených komplexů: '  num2str(pocetQRS)])
disp(['Tepová frekvence = ' num2str(tep) '/min'])

