close all;clear all;clc;
%% Inicializace
s = serial('COM9');
set(s,'BaudRate', 57600);
% s = Bluetooth('EKG',1);
load look_up_table; 
fopen(s);
fvz = 1000;
kalibrace = 1023; % inicializace kalibracni hodnoty 
pocet_vzorku_x = 2000; % Tyto vzorky slouží jako okno pro výpočet TF
zoom = 1000;
omezeni = 1;    % omezovac chyby baterie
baterka = 0;    % ovladani baterie
rr_int = 0;     % inicializace intervalu
rr_int2 = 0;    % inicializace intevalu
pom = 1;        % pomocne počitadlo   
u = 1;          % rr_ int index
z = 1;          % Zoom index
v = 1;          % pocitadlo
interval = 1;   %Interval pro výpočet TF
rr_int_tf = 500; %Inicializace proměnné výpočtu tepové frekvence
prumer_osy_y = 0;
zesileni = 1;
%%
pro_tep = 1;
tepova_frekvence = NaN; % pro ucebni fazi tepove frekvence
%%
max_y = 0;
hl_pocitadlo = 1;
zpetne_hledani = 'zatim_nic_nehledam';
%% Čekání na Wait
while (1)
    a = sscanf((fscanf(s)),'%s');
    switch a
        case 'Wait'
       fwrite(s,1,'int8');
        break
    end
end
%% Akvizice dat  
while(1)
    %% ZOOM
    while (1)
        while z <=zoom      % Pro natavení zoomu 
            volty = (3/1023) * (str2num(fscanf(s))-511.17);
            prepocet = (volty/zesileni)*1000;
            zobrazeni_h{z} = prepocet;
            y_h{pom} = zobrazeni_h{z};
            z = z + 1;      % Počítadlo zobrazení
            pom = pom + 1;  % Pomocné počítadlo pro počet vzorků
            if u > 9
            pamet_tep(pro_tep) = tepova_frekvence;
            pro_tep = pro_tep + 1;
            end
        end
        if u > 9
      figure(1)
      plot(pamet_tep)
      histogram = gca;
      histogram.XLabel.String = 'Čas [ms]';
      histogram.YLabel.String = 'BPM';
      histogram.Title.String =  'Tepová frekvence';
        
        drawnow
        end

        zobrazeni = cell2mat(zobrazeni_h);
        graf=max(zobrazeni);
        y = cell2mat(y_h);
        zaznam{hl_pocitadlo} = y;
        %% Filtrace 50 Hz
%         dolni_f = 49.5;     % Filtrace pouze pro zobrazení/nahrávání
%         horni_f = 50.5;
%         rad_f = 2;
%         Wn = [dolni_f/(fvz/2) horni_f/(fvz/2)];
%         [b,a] = butter(rad_f,Wn,'stop');
%         zobrazeni = filtfilt(b,a,zobrazeni);
%         %
figure(2)
      plot(zobrazeni)
      graf = gca;
      osa_x = graf.XTick;
      osa_y = graf.YTick;
      graf.XTickLabel = osa_x/1000;
      graf.YTickLabel = osa_y; % posun offsetu 
      graf.XLabel.String = 'Čas [ms]';
      graf.YLabel.String = 'Napětí [mV]';
      graf.Title.String =  'EKG';
      grid on 
      grid minor
        drawnow
        z = 1;zobrazeni = {};
        if pom >= pocet_vzorku_x         % Kontrola dosažení 4s záznamu
            pom = 1;
            break
        end
    end

%% Detekce R vlny
dolni_f = 5;
horni_f = 15;
rad_f = 2;
Wn = [dolni_f/(fvz/2) horni_f/(fvz/2)];
[b,a] = butter(rad_f,Wn,'bandpass');
signal = filtfilt(b,a,y);
% plot(signal)
derivace = diff(signal);
mocnina = derivace.^2;
delka_okna = 110;    % udava se pocet vzorku integracniho okna
    for int = 1:length(mocnina)-delka_okna
        integrace(int) = sum(mocnina(int:int+delka_okna));
    end
p = 1;  % Pocitadlo R_vln
i = 1;  % Inkrement hlavni smycky
pomocne_pocitadlo = 1;

    Prah_signal = 0.125*max(integrace);    % Prah signalu
    Prah = Prah_signal;                    % Prah signalu
    Prah_sum = 0.125*mean(integrace);   % Prah sumu

while i<length(integrace)- delka_okna
    if integrace(i)> Prah
        okno = integrace(i:i+delka_okna);
        Max_okno = max(okno);
        Sum_okno = mean(okno);
        poloha_peaku_pom = find(integrace==Max_okno);
            if poloha_peaku_pom == 1 || poloha_peaku_pom == 2
                 poloha_peaku_pom = 3;         
            else
                dolni_mez = integrace(poloha_peaku_pom - 1);  
                horni_mez = integrace(poloha_peaku_pom + 1);   
            if dolni_mez <= Max_okno && horni_mez <= Max_okno   
                poloha_peaku(p) = poloha_peaku_pom;
            if p >= 2
                rr_int(u) = abs(diff(poloha_peaku(end-1:end)));     % zapis rr_intervalu do pameti
                rr_int_help = abs(diff(poloha_peaku(end-1:end)));   % stávající interval
                u = u+1;
          
            end    
                hodnota_peaku(p) = integrace(poloha_peaku(p)); 
                Prah_signal = hodnota_peaku(end);    % Odhad amplitudy dalsiho peaku
                switch zpetne_hledani         % Zpetne hledani
                     case    'Hledam'
                        Prah_signal = 0.25 * Max_okno + 0.75 * Prah_signal;  % Odhad peaku pro zpentne hledani
                        zpetne_hledani = 'nasel';
                     otherwise
%% Rovnice pro výpočet prahových hodnot
                        Prah_signal = 0.125 * Max_okno + 0.8 * Prah_signal;    
                        Prah_sum = 0.125 * Sum_okno + 0.875 * Prah_sum;
                        Prah = Prah_sum + 0.25 * (Prah_signal - Prah_sum);
                end
                p = p+1;
                pomocne_pocitadlo = 1;              
           end
        end
        i = i+0.25*fvz;     % Refrakterni perioda 
    end
    if u < 9 % učební fáze
        RR_min = 0.92 * mean(rr_int);
        RR_max = 1.16 * mean(rr_int);
        RR_AVE1 = mean(rr_int);
    end
    
    if u > 9    % systém naučen
        if RR_min <= rr_int_help &&  rr_int_help <= RR_max % výpočet validní RR 
            if p > 2
            rr_int2(v) = rr_int_help;   % zápis validní rr do vektoru
            rr_int_tf(interval) = rr_int_help
            RR_AVE1 = mean(rr_int2);    % výpočet střední hodnoty intervalu
            RR_min = 0.92 * mean(rr_int2)   % minimum
            RR_max = 1.16 * mean(rr_int2)   % maximum
            rr_int_help = 0;
            v = v + 1;
            interval = interval + 1;
            end
        end
    
        if pomocne_pocitadlo > RR_AVE1*1.26
            if length(poloha_peaku) > 0
                i = poloha_peaku(end);
            end

        Prah = Prah * 0.5;
        zpetne_hledani = 'Hledam';
        pomocne_pocitadlo = 1;
            end
     
    end

i = i+1;
pomocne_pocitadlo = pomocne_pocitadlo + 1;
if interval == 3
    interval = 1;
end
if v == 8
    v = 1;
end
end
if u > 9
     rozdil = mean(rr_int_tf);
     prepocet  = rozdil/fvz;
     tepova_frekvence = round(60/mean(prepocet),0)
end
%% Nulovani promennych
    poloha_peaku = [];
    hodnota_peaku = [];
    rozdil_vektor = [];
    prepocet = [];
%     integrace = [];
ampituda=peak2peak(y)
    y = [];
    y_h = {};
%% Zapis tepove frekvence do vektoru
hl_pocitadlo = hl_pocitadlo + 1;
drawnow
%% Detekce stavu baterie
switch baterka
    case 1
        pomoc = 'Chyba baterie'
    otherwise
fprintf(s,'%s','B');    % Odeslání příkazu B
             while (1)
             a = sscanf(fscanf(s),'%s');    % čekání na návěstí BATERRY
                switch a 
                    case 'BATERRY'
                    b = str2num(fscanf(s));   % Druhý vzorek je stav Baterie
                    osa = str2num(fscanf(s))  % Třetí vzorek je vypocet osy y 
                      zesileni = 0.2222*osa^2 - 429.7*osa + 208665
                    break
                end
                if omezeni == 1000   % v případě nezaznamenání příkazu baterka po dobu 1s je cyklus čtení stavu baterie ukončen
                    baterka = 1;
                    break
                end
                omezeni = omezeni + 1;
             end
b=b(1);
napeti = round( b * (3/1023),2);
index_napeti = find(krivka(:,1)==napeti);
procenta = krivka(index_napeti,2);
end
omezeni = 1;
end
fclose(s) 
clear ('s')
delete('s')    