clc
clear variables
%% Nactení dat
% Nacitani dat je provedeno pro soubory dat a mat. V pripade dat souboru se
% pouzije nacitani funkci rddata. V pripade jiz rozmerenych dat ulozenych v
% mat souboru se preskoci akce filtrace a rozmereni signalu.

[filename, pathname] = uigetfile({'*.*'},'Vyber soubor ktery se ma nacist'); % open dialog pro vyber souboru
jmeno=filename(1:length(filename)-4);   %odříznutí přípony
pripona=filename(length(filename)-3:length(filename));
if strcmp(pripona,'.mat')==1 % nactení souboru
    display('Hodnoty nacitany');
    load([pathname filename]) % nacítani mat souboru
else
    b1=exist([pathname jmeno '.atr'],'file');
    b2=exist([pathname jmeno '.hea'],'file');
    b3=exist([pathname jmeno '.dat'],'file');
    if b1==2 && b2==2 && b3==2
        [ekgn,atr]=rddata(jmeno,pathname); % nacitaní binárního souboru z physionet.org
    else
        error('Potrebne soubory nenalezeny. Zkontrolujte dostupnost souboru typu atr,dat a hea');
    end
end
clc

%% Vypocty
a1=exist('prumer','var');
a2=exist('Qpozc','var');
a3=exist('Rpozc','var');
a4=exist('Spozc','var');
a5=exist('Jpozc','var');
a6=exist('ekgn','var');
a7=exist('atr','var');
if a1==1 && a2==1 && a3==1 && a4==1 && a5==1 && a7==1
    display('Hodnoty nacteny');
    
else
    if a6==1
      ekg=ekgfil(ekgn(:,2),atr);
      posunuti=mean(ekg(1:10000));
      ekg=ekg-posunuti;
      [Qpoz,Qhod,Rpoz,Rhod,Spoz,Shod,Jpoz,stredy,Tpoz]=rdet(ekg,atr);
      n=1;
      pocetR=length(Rpoz);
      prumer=zeros(floor(pocetR/10),500);
      pocet_prumeru=10*floor(pocetR/10);
        for i=1:10:(pocet_prumeru) % pro vyber intervalu vypoctu
          if (stredy(i+1)-stredy(i))>450 
              podil(n)=round((stredy(i+1)-stredy(i))/250);
              delka=round((stredy(i+1)-stredy(i))/podil(n));  
          else
            delka=stredy(i+1)-stredy(i);
          end
          if delka>450
            podil(n)=round((stredy(i+1)-stredy(i))/345);
            delka=(stredy(i+1)-stredy(i))/podil(n);
          end
            for j=1:5 % prumer z peti po sobe jdoucích ekg po kazdych 10 cyklech 
              vypis(j,:)=ekg(Rpoz(j+i)-0.5*delka:Rpoz(j+i)+0.5*delka);  
            end
%             plot(1:length(vypis),vypis(1,:),1:length(vypis),vypis(2,:),1:length(vypis),vypis(3,:),1:length(vypis),vypis(4,:),1:length(vypis),vypis(5,:))
%             pause
%             close all
          prumer(n,1:length(vypis))=(sum(vypis(1:5,:),'double'))/5;
          prumer(n,500)=length(vypis(1,:));%Uloží na konec matice délku aktuálního řádku bez nul
          prumer(n,499)=(sum(prumer(n,1:prumer(n,500))))/prumer(n,500);
          clear vypis;
          akt_delka=prumer(n,500);  
%% Výber T a P vlny - nejsou v programu momentalne potrebne 
%          plot(prumer(n,1:akt_delka-1))
%          title(['Vyber Zacatek a konec ' num2str(n) '. P vlny z ' num2str(pocet_prumeru/50)]);
%          [Ppoz,Phod]=ginput(2);
%           plot(prumer(n,1:akt_delka))
%           title(['Vyber Zacatek a konec ' num2str(n) '. T vlny z ' num2str(pocet_prumeru/50)]);
%           [Tpoz,Thod]=ginput(2);
%           Pzacc(:,n)=Ppoz(1,1);
%           Pkonc(:,n)=Ppoz(2,1);
%           Tzacc(:,n)=Tpoz(1,1);
%           Tkonc(:,n)=Tpoz(2,1);
%% 
          as=1;
          pr=[prumer(n,1:akt_delka) prumer(n,1:akt_delka)];
          if length(pr)<501
             pr(end:505)=zeros; 
          end
          [Qpozb,Qhodb,Rpozb,Rhodb,Spozb,Shodb,Jpozb,stredyb,Tpozb]=rdet(pr,as);
          Qpozc(:,n)=Qpozb(1);
          Rpozc(:,n)=Rpozb(1);
          Spozc(:,n)=Spozb(1);
          Jpozc(:,n)=Jpozb(1);
          Tpozc(:,n)=Tpozb(1);
%         n
%          Vykresleni prumeru s oznacenim klicovych bodu
%           plot(pr)
%           hold on
%           line([Spozc(:,n) Spozc(:,n)],[2.5 3],'Color','r');
%           line([Rpozc(:,n) Rpozc(:,n)],[2.5 3],'Color','g');
%           line([Tpozc(:,n) Tpozc(:,n)],[2.5 3],'Color','b');
%           pause 
%           close all
          n=n+1;
          %clc 
       end   
    else
        error('Potrebné data nenalezeny, soubor neobsahuje EKG signál (vektor ekgn)');
    end %if (23)
end     % if (20) 
%% výpocet sumacniho trendu S-T

for i=1:length(prumer(:,1))
    prumer(i,1:prumer(i,500))=prumer(i,1:prumer(i,500))-mean(prumer(i,499));    
end
delka=(prumer(1:length(prumer(:,1)),500))';
f=linspace(0.01,100,100);

for i=1:length(prumer(:,1))    
      %vypocet spektra
      c = abs(cwt(prumer(i,1:delka(1,i)),f,'morl'));%morl,meyr,coif1
      
%% trendy      
        for k=1:100 % sumace radku
            pomST(1,k)=(sum(c(k,Spozc(1,i):Tpozc(1,i))));
        end
      ST_A(:,i)=sum(pomST(1:100))/((Tpozc(1,i)-Spozc(1,i))*100);    % trend vývoje ST
      ST_B(:,i)=sum(pomST(80:100))/((Tpozc(1,i)-Spozc(1,i))*21);  
      clc
      trend80P(1,i)=(sum(prumer(i,(Jpozc(1,i)+18):(Jpozc(1,i)+22)))/5);%-prumer(i,Jpozc(1,i)); % průměrovaná metoda 80 ms (5 vzorků)
%% Vypocty distribude maxima 
      m_d(:,i)=max(max(c(1:40,Spozc(1,i):Tpozc(1,i))));
      m_h(:,i)=max(max(c(41:100,Spozc(1,i):Tpozc(1,i)-20)));
      [m_dm(i,1) m_dm(i,2)]=find(c(1:40,Spozc(1,i):Tpozc(1,i))==m_d(:,i));
      [m_hm(i,1) m_hm(i,2)]=find(c(40:100,Spozc(1,i):Tpozc(1,i)-20)==m_h(:,i));
      m_dm(i,2)=m_dm(i,2)+Spozc(1,i)-1;
      m_hm(i,1)=m_hm(i,1)+40-1;
      m_hm(i,2)=m_hm(i,2)+Spozc(1,i)-1;

%% vypocet trendu veberem okruhu dat 
%maxima v horní a dolní casti, amplituda a pozice 
      vyber=c(1:100,Spozc(i):Tpozc(i));
      for radek=1:100 % nulovani hodnot pod 0.5
         for  sloupec=1:length(vyber(1,:))
          if vyber(radek,sloupec)>=0.5
            vyber(radek,sloupec)=vyber(radek,sloupec);
          else
            vyber(radek,sloupec)=0;
          end
         end
      end
      pomSTC(i)=sum(sum(vyber(1:100,1:length(vyber(1,:)))));      
      pomSTD(i)=sum(sum(vyber(80:100,1:length(vyber(1,:)))));  
      ST_C(i)=pomSTC(i)/length(find(vyber(1:100,1:length(vyber(1,:))))>0);%stredni hodnoty na jeden bod
      ST_D(i)=pomSTD(i)/length(find(vyber(80:100,1:length(vyber(1,:)))>0));
%     soucty po v intervalech po deseti
%     a1(i)=sum(sum(vyber(1:10,1:length(vyber(1,:))))); 
%     a2(i)=sum(sum(vyber(10:20,1:length(vyber(1,:))))); 
%     a3(i)=sum(sum(vyber(20:30,1:length(vyber(1,:))))); 
%     a4(i)=sum(sum(vyber(30:40,1:length(vyber(1,:))))); 
%     a5(i)=sum(sum(vyber(40:50,1:length(vyber(1,:))))); 
%     a6(i)=sum(sum(vyber(50:60,1:length(vyber(1,:))))); 
%     a7(i)=sum(sum(vyber(60:70,1:length(vyber(1,:))))); 
%     a8(i)=sum(sum(vyber(70:80,1:length(vyber(1,:))))); 
%     a9(i)=sum(sum(vyber(80:90,1:length(vyber(1,:))))); 
%     a10(i)=sum(sum(vyber(90:100,1:length(vyber(1,:))))); 
end % konec cyklu vypoctu
%% Vypocet distribuce
dis1=zeros(40,max(m_dm(:,2)));
for i=1:length(m_dm)
   dis1(m_dm(i,1),m_dm(i,2))=dis1(m_dm(i,1),m_dm(i,2))+1;
end

dis2=zeros(100,max(m_hm(:,2)));
for i=1:length(m_hm)
   dis2(m_hm(i,1),m_hm(i,2))=dis2(m_hm(i,1),m_hm(i,2))+1;
end

xx=ST_A(2:end)-ST_A(1:end-1);
yy=ST_A(3:end)-ST_A(2:end-1);
xx=xx(1:end-1);

%% Vystupy
    osa=Rpoz(1:10:(pocet_prumeru));
    STp=find(atr(:,1)==18);
    Tp=find(atr(:,1)==19);
    STchange=atr(STp,3);
    Tchange=atr(Tp,3);

%----------------------------------------------------------------------
        figure(1)
        plot(osa,(trend80P-min(trend80P)))
        title('Vyvojove trendy ziskane odectenim hodnoty J+80 ms');
        xlabel('n[-]')
        ylabel('U[V]')
        hold on
        line([STchange STchange],[max(abs(trend80P)) 1.1*max(abs(trend80P))],'Color','r');
        line([Tchange Tchange],[max(abs(trend80P)) 1.1*max(abs(trend80P))],'Color','g');    
%----------------------------------------------------------------------        
        figure(2)
        plot(osa,ST_A)
        title('ST_A metoda (rozsah-1:100,bez prahování)');
        xlabel('n[-]')
        ylabel('Amplituda')
        hold on
        line([STchange STchange],[max(ST_A) 1.1*max(ST_A)],'Color','r');
        line([Tchange Tchange],[max(ST_A) 1.1*max(ST_A)],'Color','g');        
 %----------------------------------------------------------------------     
        figure(3);
        plot(osa,ST_B)
        title('ST_B metoda (rozsah-80:100,bez prahování)');
        xlabel('n[-]')
        ylabel('Amplituda')
        hold on
        line([STchange STchange],[max(ST_B) 1.1*max(ST_B)],'Color','r');
        line([Tchange Tchange],[max(ST_B) 1.1*max(ST_B)],'Color','g'); 
 %----------------------------------------------------------------------  
        figure(4);
        plot(osa,ST_C)
        title('ST_C metoda (rozsah-1:100,s prahováním)');
        xlabel('n[-]')
        ylabel('Amplituda')
        hold on
        line([STchange STchange],[max(ST_C) 1.1*max(ST_C)],'Color','r');
        line([Tchange Tchange],[max(ST_C) 1.1*max(ST_C)],'Color','g');
%----------------------------------------------------------------------        
        figure(5);
        plot(osa,ST_D)
        title('ST_D metoda (rozsah-80:100,s prahováním)');
        xlabel('n[-]')
        ylabel('Amplituda')
        hold on
        line([STchange STchange],[max(ST_D) 1.1*max(ST_D)],'Color','r');
        line([Tchange Tchange],[max(ST_D) 1.1*max(ST_D)],'Color','g');
%----------------------------------------------------------------------
        figure(6);
        imagesc(dis1)
        title('Distribuce maxim 1, rozsah 1:40');
        xlabel('n[-]')
        ylabel('SCALE')
%----------------------------------------------------------------------     
        figure(7);
        imagesc(dis2)
        title('Distribuce maxim 2, rozsah 41:100');
        xlabel('n[-]')
        ylabel('SCALE')
        axis([1 length(dis2(1,:)) 40 100])
%----------------------------------------------------------------------
        figure(8);
        scatter(xx, yy);
        title('Poincareho mapa pro metodu STB');       
        xlabel('x_n_+_1-x_n')
        ylabel('x_n_+_2-x_n_+_1')

%% Zdrojový kód pro vykreslení hranic do prumerovaneho ekg cyklu
% %% Zdrojový kód pro vykreslení hranic do spektra 
% for i=1:length(prumer(:,1)) 
%    c = abs(cwt(prumer(i,1:delka(1,i)),f,'meyr'));%morl,meyr,coif1
%    figure(4);
%    imagesc (c); figure(gcf)   
% end
% figure(6)
% d=length(vyber);
% for i=1:length(S1_poz)
%    line([1000*S1_poz(i) 1000*S1_poz(i)],[-1 1],'Color','r');
%    hold on 
%    line([1000*S2_poz(1,i) 1000*S2_poz(1,i)],[-0.7 0.7],'Color','b');
%    line([1000*S2_poz(2,i) 1000*S2_poz(2,i)],[-0.7 0.7],'Color','g');
% end
% plot(vyber)
% title('Ukázka detekce - vykreslení nadetekovaných pozic v prvotním signálu') 

% Vykresleni hranic do spektra jednotlivé
%    imagesc (c);
%    hold on 
%    line([Spozc(i) Spozc(i)],[1 50],'Color','r');
%    line([Tzacc(1,i) Tzacc(1,i)],[1 50],'Color','r');
