function [locs, TF] = getR(EKG, EKG_fil, predR, zaR, fvz)
% Funkce vraci pozice vhodnych R vln EKG s fvz 512 Hz. Vhodnymi se mysli
% takove, ktere jsou soucasti nezkreslenych tepu a tyto tepy jsou tedy
% vhodne k prumerovani. Zaroven jsou vynechany R vlny, ktere jsou priliz
% blizko okraje signalu (tj. pred R vlnou musi byt minimalne predR vzorku
% signalu a za R vlnou minimalne zaR vzorku signalu). 
%
% Vhodne tepy jsou urceny jednoduchym pravidlem - divame se na medianovou
% hodnotu R vln v originalnim i filtrovanem signalu a pokud se dany tep
% nejak zvlast vymyka, radsi ho odstranime.
%
% EKG ... originalni EKG, vektor o delce 10240 vzorku (20 s pri fvz = 512 Hz)
% EKG_fil ... EKG filtrovane PP 1-40 Hz, vektor o delce 10240 vzorku (20 s pri fvz = 512 Hz)
% predR, zaR ... skalary udavajici kolik vzorku chceme pred/za R
%
% Detektor R vln je inspirovan predmetem AABS.

%% Detektor R vln
rad_filtru = round(fvz * 1101/512);
[b,a] = fir1(rad_filtru, [11/(fvz/2) 21/(fvz/2)], 'bandpass'); % vyber QRS frekvenci
Rzvyraznen = filtfilt(b,a,EKG_fil);
Rzvyraznen = Rzvyraznen.^2; 
rad_filtru = round(fvz * 60/512);
[b,a] = fir1(rad_filtru,3/(fvz/2),'low'); % vyber obalky
Rzvyraznen = filtfilt(b,a,Rzvyraznen);
prah=quantile(Rzvyraznen,0.8);
vzdalenost = round(fvz * 300/512);
[~,locs] = findpeaks(Rzvyraznen, 'MinPeakHeight',prah,'MinPeakDistance',vzdalenost); % vyber peaku

% Zpresneni polohy R vlny
okoli = round(fvz * 50/512); % okoli kolem detekovane pozice R, budu se divat po maximu EKG v danem okoli

for n = 1:length(locs)
    if locs(n)-predR<1 || locs(n)+zaR>length(EKG_fil) % R je blizko kraji, nemame cely tep - vynechame
        locs(n) = 0;
        continue
    else
        vektor = locs(n)-okoli:locs(n)+okoli; 
        [~,k] = max(EKG_fil(vektor)); % nalezeni opravdove R vlny (extrem v danem okoli)
        % pri okoli = 50 bude promenna vektor delky 101 a puvodne nalezena pozice R vlny bude na 51. miste
        if k==okoli+1
            continue % zadna zmena
        elseif k>okoli+1
            locs(n) = (k-(okoli+1))+locs(n); % (okoli+1) je pozice puvodne nalezene R vlny v promenne vektor
        else
            locs(n) = locs(n)-((okoli+1)-k);
        end
    end
end
locs(locs==0) = [];

% zjisteni TF
vzdalenost_R = median(diff(locs)); % vzdalenost R ve vzorcich, reprezentujici hodnota median
perioda = vzdalenost_R/fvz;
frekvence = 1/perioda; % 1/s
TF = 60*frekvence;

% kontrola tepu blizko kraju po uprave pozic R vln
if locs(1)-predR<1  % R je blizko kraji, nemame cely tep - vynechame
    locs(1) = [];
end
if locs(end)+zaR>length(EKG_fil)
    locs(end) = [];
end

% figure(5)
% subplot(211)
% plot(EKG)
% title('Original signal','FontSize',15)
% hold on
% stem(locs,EKG(locs),'LineWidth',0.8,'Color','red')
% subplot(212)
% plot(EKG_fil)
% title('Filtered signal','FontSize',15)
% hold on
% stem(locs,EKG_fil(locs),'LineWidth',0.8,'Color','red')

%% Vyber vhodnych R vln
med_R_orig = median(EKG(locs));
med_R_fil = median(EKG_fil(locs));

pom_locs = locs;
for n = 1:length(locs)
    if EKG(locs(n))>med_R_orig+3*med_R_orig || EKG(locs(n))<med_R_orig-3*med_R_orig ... % vyrazna zmena v originalnim signalu
            || EKG_fil(locs(n))>med_R_fil+0.3*med_R_fil || EKG_fil(locs(n))<med_R_fil-0.3*med_R_fil  % ci mensi zmena ve filtrovanem signalu
        locs(n) = 0;    % odstraneni pochybnych tepu z dalsiho zpracovavani
    end
end
locs(locs==0) = [];
prah_orig = 3;
prah_filt = 0.3;
while isempty(locs) % zabraneni NaN v priznakove matici
    locs = pom_locs;
    prah_orig = prah_orig + 0.5;
    prah_filt = prah_filt + 0.05;
    if EKG(locs(n))>med_R_orig+prah_orig*med_R_orig || EKG(locs(n))<med_R_orig-prah_orig*med_R_orig ... % vyrazna zmena v originalnim signalu
            || EKG_fil(locs(n))>med_R_fil+prah_filt*med_R_fil || EKG_fil(locs(n))<med_R_fil-prah_filt*med_R_fil  % ci mensi zmena ve filtrovanem signalu
        locs(n) = 0;    % odstraneni pochybnych tepu z dalsiho zpracovavani
    end
    locs(locs==0) = [];
end

% figure(5)
% x = [0 10240];
% subplot(211)
% y = [med_R_orig med_R_orig];
% line(x,y,'Color','#77AC30','LineWidth',0.7)
% y = [(med_R_orig+3*med_R_orig) (med_R_orig+3*med_R_orig)];
% line(x,y,'Color','#77AC30','LineStyle','--','LineWidth',0.7)
% y = [(med_R_orig-3*med_R_orig) (med_R_orig-3*med_R_orig)];
% line(x,y,'Color','#77AC30','LineStyle','--','LineWidth',0.7)
% hold off
% subplot(212)
% y = [med_R_fil med_R_fil];
% line(x,y,'Color','#77AC30','LineWidth',0.7)
% y = [(med_R_fil+0.3*med_R_fil) (med_R_fil+0.3*med_R_fil)];
% line(x,y,'Color','#77AC30','LineStyle','--','LineWidth',0.7)
% y = [(med_R_fil-0.3*med_R_fil) (med_R_fil-0.3*med_R_fil)];
% line(x,y,'Color','#77AC30','LineStyle','--','LineWidth',0.7)
% hold off
% 
% figure(6)
% subplot(211)
% plot(EKG_fil)
% title('Locations of R waves which are usable for beat cumulation','FontSize',15)
% hold on
% stem(locs,EKG_fil(locs),'LineWidth',0.8,'Color','red')
% hold off

end