
% FILTRACE

h=fir1(51,[0.1 1]/100); %filtrace pasmovou propusti 14-25Hz 
PPGf=filtfilt(h,1,ppg); %filtrace DP, mezni frekvence 5 Hz, vyhlazeni signalu 

% ODSTRANĚNÍ DRIFTU 

prodlouzenyPPGf=[PPGf(1:350) PPGf PPGf(end-349:end)];   % prodlouzeni EDA z obou stran o 2500 vzorku, opakovanim zacatku a konce, prodlouzeni kvuli vypoctu ve for cyklu
for i=351:length(prodlouzenyPPGf)-350                   % projeti prodlouzenyEDAf od vzroku 2500 do vzorku end-2500, odpovida to tedy normalnimu EDA 
    prumer(i-350)=mean(prodlouzenyPPGf(i-350:i+349));   % vypocet driftu signalu, vezme se okno delky 5000 vzorku a projizdi se signalem, z tech 5000 hodnot se vezme minimalni hodnota, ta je urcena jako drift pro prislusny vzorek   
end

PPG0=PPGf-prumer;      % odecteni driftu od EDA filtrovany

% DETEKCE PÍKŮ

pruchod=0; j=0;    % pomocne promenne pro vypocet piku
for i=1:N
    if PPG0(i)>0 && pruchod==0   % detekuji se pruchody signalu nulou, nejdriv pruchod ze zapornych do kladnych 
        pruchod=1;
        p1=i;                   % oznaceni vzorku na kterem doslo k prvnimu pruchodu 
    elseif PPG0(i)<0 && pruchod==1  % detekce z kladnych do zapornych
        pruchod=2;
        p2=i;                   % oznaceni vzorku na kterem doslo k druhemu pruchodu 
    end
    
    if pruchod==2                       % pokud jsou detekovany oba pruchody, tak se urci maximalni hodnota (pik) mezi pruchody a jeji pozice
        [cislo pozice]=max(PPG0(p1:p2)); 
        pruchod=0;                      % zruseni detekovaneho prechodu, nastaveni pro novou detekci pruchodu          
        j=j+1;
        ppgx(j)=pozice+p1;   %zapsani pozice piku upravenou pro cely signal, posun o p1
        ppgy(j)=cislo;       %zapsani hodnoty piku
    end
end
%Kód výše je použitím a úpravou: 
%BUBNÍK, K.; KOLÁŘOVÁ, J., Software pro zpracování a filtraci PPG signálu, Brno: FEKT VUT v Brně 2012. Dostupné z WWW: http://www.dbme.feec.vutbr.cz/vyzkum-a-vyvoj/produkty

 

%  VÝPOČET HRV
polohy=ppgx; %vektor poloh detekovaných píků
cas=t(polohy); %časy detekce
i=1; %iterační proměnná
while i<(length(cas))
    tep(i)=60/(cas(i+1)-cas(i)); %výpočet tepové frekvence odečtením sousedních hodnot R vln
    i=i+1;
    if i>5 % ošetření chybné detekce
        pom=mean(tep); %průměr TF předcházejících tepů pro ošetření chybné detekce
    
        if tep(i-1)<pom*1.2 && tep(i-1)>pom*0.8 %v případě, že vzdálenost mezi detekcí...
                                                %...je větší než 1.2x průměr nebo menší než 0.8x průměr detekce je chybná
            tep(i-1)=tep(i-1); %zachová se předcházející hodnota
            i=i;
        else %pokud ne je detecke správná
            tep(i-1)=pom; %na pozici s chybnou detekcí se uloží průměr předchozích hodnot        
            i=i;
        end
    end
    t_cas(i-1)=cas(i); %získání hodnot pro osu x při vykreslení tepové frekvence
     
end
rozdil=zeros(1,length(tep)-1); %proměnná pro rozdíly tepových frekvencí
for i=2:length(tep) %procházení vektoru TF
    rozdil(i-1)=tep(i)-tep(i-1); %výpočet rozdílů tepových frekvencí pro názornější vykreslení
end
% PROKLÁDÁNÍ SPLAJNEM, VÝPOČET VLIVU SYMPATIKU A PARASYMPATIKU

x =t(polohy(1:end-2)) ;
y =rozdil; %upřesnění prokládaných bodů
xx =t(1:polohy(end-2));% počet bodů křivky
PPGs = spline(x,y,xx); %proložení hodnot rozdílů kubickým splajnem


[lmaxima,indices] = localmax(PPGs);     %nalezení lokálních maxim PPGs
parasympaticus= PPGs(indices)>0;        %hledání maxim (>0)
sympaticus= PPGs(indices)<0;            %hledání minim (<0)
parasympaticus=indices(parasympaticus); %polohy maxim
sympaticus=indices(sympaticus);         %polohy minim

o =t(parasympaticus);
p =PPGs(parasympaticus);
oo =t(1:polohy(end-2));
parasymp = spline(o,p,oo);  %proložení maxim kubickým splajnem

k =t(sympaticus) ;
l =PPGs(sympaticus);
kk =t(1:polohy(end-2));
symp = spline(k,l,kk);      %proložení minim kubickým splajnem

vyslANS=parasymp+symp;      %odečtení vlivu smypatiku a parasympatiku
vyslANS=[xx;vyslANS];       %vytvoření matice s osou a výsledným rozdílem
