% METODA TEPOVE FREKVENCE, TVAROVEHO FAKTORU A BINARNIHO KODOVANI
%vsechny signaly jsou na CD 
clear all
clc
close all
H=zeros(9,46); % po skonceni cyklu for v teto matice budou vysledky, radky popsany dale rady 227-235
HR_vsechny_PVC=zeros(1,2);
HR_vsechny_ostatni=zeros(1,2186);
FF_vsechny_PVC=zeros(1,2);
FF_vsechny_ostatni=zeros(1,2186);
save HR_vsechny_PVC HR_vsechny_PVC
save HR_vsechny_ostatni HR_vsechny_ostatni
save FF_vsechny_PVC FF_vsechny_PVC
save FF_vsechny_ostatni FF_vsechny_ostatni
save H H;
%4:49 vsechny signaly
% 1)Pro prohlizeni hodnot parametru na box plotech vsech 36 signalu
% nechat projit for cyklus pres vsechny signaly a pote odkomentovat rady 
% skriptu 355-380
% 2)pro prohlizeni box plotu jednotliveho signalu zakomentovat rady
% 333-335 clear all clc close all , vybrat si signal treba q=7, a pote
% odkomentovat rady 355-380 
% 3)pro prohlizeni celkove stredni senzitity a specificity ze
% vsech 36 signalu odkomentovat rady 345-347
% 4) pro prohlizeni detekce v jednotlivych signalech, vybrat si signal,
% zakomentovat 333-335 a odkomentovat  %%ukazka vysledku rada 242
% 5)pokud si vyberete signal, ktery neobsahuje PVC, tak nic nezobrazi 
for q=4:49
files = dir;
directoryNames = {files([files.isdir]).name};
directoryNames = directoryNames(~ismember(directoryNames,{'.','..'}));
k=struct(files);
sign=files(q).name;
o=char(sign)
load(o)
load('H.mat')
load('FF_vsechny_PVC.mat')
load('FF_vsechny_ostatni.mat')
load('HR_vsechny_PVC.mat')
load('HR_vsechny_ostatni.mat')
time = 1/360:1/360:tm(end,1);  %casova osa signalu v sekundach
[~,y1]=size(time);
% figure(1)
% plot(time,x(1:y1,1)) % zobrazeni 1. signalu
% xlabel('t[s]')
% ylabel('U[mV]')

%pozice a pocty ruznych patologii
[pozice_V]=find(type=='V');%pvc
%% pokud signal bude obsahovat 0 pocet PVC --> konec
if length(pozice_V)>0 
[pozice_N]=find(type=='N');%normalni
[pozice_x]=find(type=='x');
[pozice_R]=find(type=='R');%Right brunch bundle boundaries
[pozice_L]=find(type=='L');%Left brunch bundle boundaries
[pozice_J]=find(type=='J');%Junctional premature	
[pozice_j]=find(type=='j');%Junctional escape	
[pozice_F]=find(type=='F');%fusion pvc
[pozice_A]=find(type=='A');%APC
[pozice_plus]=find(type=='+'); %sum a dalsi
[pozice_vlnka]=find(type=='~');%sum a dalsi

pozice_V=(pozice_V)';
% if length(pozice_V)>0
[y,~]=size(ann);
A=zeros(1,y);
Bpo=zeros(1,y);
Cpred=zeros(1,y);
Detekce=zeros(1,y);

for i=1:y 
A(i)=time(1,ann(i,1));
end

for j=3:length(ann)-2
Bpo(1,j)=A(j+1)-A(j); % deskriptor RR po aktualni pozice QRS
Cpred(1,j)=A(j)-A(j-1); % deskriptor RR pred aktualni pozice QRS
end 

for i=1:y
Rvse(i)=ann(i,1)   ;
% plot(time(1,Rvse(i)),x(Rvse(i),1),'*b');
end

%vymena PVC na 1 a Ostatni na 2
for i=1:length(ann)-2
    if type(i,1)=='V'
type_1(i,1)=1;
else
    type_1(i,1)=2;
    end
end

%% Detekce zacatku QRS, 73 vzorku = 200ms , vypocet nejvetsi diference
for j=3:length(Rvse)-2
txR=time(1,Rvse(1,j));
tyR=x(Rvse(1,j),1);
tx=time(1,Rvse(1,j)-72);
ty=x(Rvse(1,j)-72,1);

xuyrok(1,1:73)=x(Rvse(1,j)-72:Rvse(1,j),1);
krok=(abs(xuyrok(1,1))+xuyrok(1,end))/73;
kroknasobeni(1,:)=[1:1:150];

for i=2:74
xrovnomerne(1,1)=x(Rvse(1,j)-72,1);
xrovnomerne(1,i)=xrovnomerne(1,1)+krok*kroknasobeni(1,i-1);
end

for i=1:73
diference(1,i)=abs(xrovnomerne(1,i)-(xuyrok(1,i)));
end
m=max(diference);
[~,ymaxl]=find(diference==m);
[~,maxl]=size(ymaxl);
if maxl>1
    ymaxl=ymaxl(1);
end
ymax=ymaxl;
dvesektime=time(1,Rvse(1,j)-(73-ymaxl));
% xx=[dvesektime,dvesektime];
% yy=[xuyrok(1,ymax),2];
% plot(xx,yy,'r')
% line([dvesektime,dvesektime],[xuyrok(1,ymax),2])
padsekbodtime=time(1,Rvse(1,j)-(73-ymax)-18);
padsekbodampl=x(Rvse(1,j)-(73-ymax)-18,1);
% xx=[padsekbodtime,dvesektime];
% yy=[padsekbodampl,xuyrok(1,ymax)];
% plot(xx,yy,'r')
% line([padsekbodtime,dvesektime ],[padsekbodampl,xuyrok(1,ymax)])
%% 19 vzorku = 50 ms, vypocet nejvetsi diference -> zacatek detekovan
krokpades=(abs(x(Rvse(1,j)-(73-ymax),1)-padsekbodampl))/18;
for i=2:19
rovnopades(1,1)=x(Rvse(1,j)-(73-ymax)-18,1);
rovnopades(1,i)=rovnopades(1,1)-krokpades*kroknasobeni(1,i-1);
% plot(time(Rvse(1,3)-72:Rvse(1,3)),xrovnomerne(1,i-1),'*r')
end
xpadesxyrok(1,1:19)=x(Rvse(1,j)-(73-ymax)-18:Rvse(1,j)-(73-ymax),1);
for i=1:19
diferencepades(1,i)=abs(rovnopades(1,i)-(xpadesxyrok(1,i)));
end
mpad=max(diferencepades);
[~,ymaxx]=find(diferencepades==mpad);
[~,maxx]=size(ymaxx);
if maxx>1
ymaxx=ymaxx(1);
end
ymaxpad=ymaxx;
onset(1,j)=time(1,Rvse(1,j)-(73-ymax)-(19-ymaxpad));
pozice_a1(1,j)=Rvse(1,j)-(73-ymax)-(19-ymaxpad);
[~,tim]=find(time==onset(1,j));
timm(1,j)=tim;
posun(1,j)=abs(x(timm(1,j),1));
%% vykreskeni cervenou hvezdickou zacatku QRS, nedorucuju spoustet to pres vsechna QRS, bude to trvat moc dlouho
% hold on
% plot(onset(1,j),x(Rvse(1,j)-(73-ymax)-(19-ymaxpad),1),'*r')

%% Detekce konce QRS, 90 vzorku = 248 ms a vypocet nejvetsi diference
a3=diff(x(Rvse(1,j):Rvse(1,j)+90,1));
mdvapet=max(a3);
[ymaxdvapet]=find(a3==mdvapet);
if length(ymaxdvapet)>1
    ymaxdvapet=ymaxdvapet(1);
end
dvapetktime=time(1,Rvse(1,j)+ymaxdvapet);
offsettime(1,j)=dvapetktime;
pozice_a2(1,j)=Rvse(1,j)+ymaxdvapet;
end

%% tvarovy faktor
for i=3:length(Rvse)-2
a0=x(pozice_a1(i):pozice_a2(i),1);
a0=(a0)';
a=diff(x(pozice_a1(i):pozice_a2(i),1));
a=(a)'; % 1 derivace
[~,delkaa]=size(a(1,:));
ti=time(1,pozice_a1(i):pozice_a1(i)+delkaa-1);
%% vypocet a vykresleni derivaci
% hold on
% plot(ti(1,:),a(1,:),'green')
a2=diff(diff(x(pozice_a1(i):pozice_a2(i),1)));
a2=(a2)'; % 2 derivace
[~,delkaa2]=size(a2(1,:));
ti1=time(1,pozice_a1(i):pozice_a1(i)+delkaa2-1);
% plot(ti1(1,:),a2(1,:),'red')
rozptyl_0=var(a0); % rozptyl signalu bez derivace
rozptyl_1=var(a);  % rozptyl signalu po 1 derivaci
rozptyl_2=var(a2); % rozptyl signalu po 2 derivaci
FF(1,i)=(rozptyl_2/rozptyl_1)/(rozptyl_2/rozptyl_0); % vypocet tvaroveho faktoru
end
%%  FF 
for i=1:length(FF)
if FF(1,i)>12 
    poz_kom(1,i)=1;
else
    poz_kom(1,i)=0;
end
end
%% vypocet tepove frekvence
for i=1:length(FF)
HR(1,i)=60/Cpred(1,i) ;

if HR(1,i)>90
poz_HR(1,i)=1;
else
poz_HR(1,i)=0;            
end
end

for i=1:length(FF)
if type_1(i)==1 && poz_kom(i)==1 && poz_HR(i)==1 
    vysl_spr_det(1,i)=1;
    spr_det_extr(1,i)=time(1,Rvse(1,i));
else
    spr_det_extr(1,i)=0;  
    vysl_spr_det(1,i)=0;
end 
end

for i=1:length(FF)
if type_1(i)==2 && poz_kom(i)==0 | type_1(i)==2 && poz_HR(i)==0 
    vysl_spr_nedet(1,i)=1;
    spr_nedet_extr(1,i)=time(1,Rvse(1,i));
else
    spr_nedet_extr(1,i)=0; 
    vysl_spr_nedet(1,i)=0;
end 
end

senzetivita=sum(vysl_spr_det==1)/(sum(vysl_spr_det==1)+(sum(type_1==1)-sum(vysl_spr_det==1)));
specificita=sum(vysl_spr_nedet==1)/(sum(vysl_spr_nedet==1)+(sum(type_1==2)-sum(vysl_spr_nedet==1)));

H(1,q-3)=senzetivita;
H(2,q-3)=specificita;
H(3,q-3)=sum(type_1==1); % pocet_PVC
H(4,q-3)=sum(vysl_spr_det==1);  % pocet detekovanych PVC
H(5,q-3)=sum(type_1==1)-sum(vysl_spr_det==1);% pocet nedetekovanych PVC
H(6,q-3)=sum(type_1==2); % pocet ostatnich
H(7,q-3)=sum(vysl_spr_nedet==1); % detekovane ostatni
H(8,q-3)=sum(type_1==2)-sum(vysl_spr_nedet==1); % nedetekovane ostatni
H(9,:)=[102,103,104,105,106,107,108,109,111,112,113,114,115,116,117,118,119,121,122,123,124,200,201,202,203,205,207,208,209,210,212,213,214,215,217,219,220,221,222,223,228,230,231,232,233,234];
% 9 radek obsahuje cisla signalu 
%% ukazka vysledku 
% if sum(vysl_spr_det==1)>0
% figure(1)
% plot(time,x(1:y1,1)) % zobrazeni 1. signalu
% xlabel('t[s]')
% ylabel('U[mV]')
% title('Vysledky detekce')
% hold on
% plot(spr_det_extr(spr_det_extr~=0),2.5,'V','Color',[0 1 0])
% % plot(A(pozice_V),2,'x','Color',[1 0 0]) % spatne zarazene, jenom pro ukazku
% plot(spr_nedet_extr(spr_nedet_extr~=0),3,'o','Color',[1 0 0])
% % nedetekovane PVC nebo Ostatni nebudou na miste R piku vubec zobrazeny zadnym symbolem.
% else
%     %pokud nenadetekuje zadne PVC
%  figure(1)
% plot(time,x(1:y1,1)) % zobrazeni 1. signalu
% xlabel('t[s]')
% ylabel('U[mV]')
% title('Vysledky detekce')
% hold on
% plot(spr_nedet_extr(spr_nedet_extr~=0),3,'o','Color',[1 0 0])   
% end
%% pouzit pro vykresleni boxplotu parametru FF 36 signalu
for i=1:length(pozice_V)-2
FF_u_PVC(1,i)=FF(pozice_V(1,i));
FF(pozice_V(1,i))=0;
end
j=1;
for i=1:length(FF)
    if FF(1,i)~=0
        FF_ostatni(1,j)=FF(1,i); 
        j=j+1;
    else
    end
end
if length(pozice_V)>2
if q==4
FF_vsechny_PVC=FF_u_PVC;
FF_vsechny_ostatni=FF_ostatni;
else
kpn=length(FF_vsechny_PVC); 
hipn=length(FF_u_PVC);
FF_vsechny_PVC(1,kpn+1:kpn+hipn)=FF_u_PVC;
kpn2=length(FF_vsechny_ostatni) ;
hipn2=length(FF_ostatni);
FF_vsechny_ostatni(1,kpn2+1:kpn2+hipn2)=FF_ostatni;
end
end
%% pouzit pro vykresleni boxplotu parametru HR 36 signalu 
for i=1:length(FF)
    HR(1,i)=60/Cpred(1,i) ;
end
for i=1:length(pozice_V)-2
HR_u_PVC(1,i)=HR(pozice_V(1,i));
HR(pozice_V(1,i))=0;
end
j=1;
for i=1:length(HR)
    if HR(1,i)~=0
        HR_ostatni(1,j)=HR(1,i); 
        j=j+1;
    else
    end
end
if length(pozice_V)>2
if q==4
HR_vsechny_PVC=HR_u_PVC;
HR_vsechny_ostatni=HR_ostatni;
else
kp=length(HR_vsechny_PVC); 
his=length(HR_u_PVC);
HR_vsechny_PVC(1,kp+1:kp+his)=HR_u_PVC;
kp2=length(HR_vsechny_ostatni) ;
hi2=length(HR_ostatni);
HR_vsechny_ostatni(1,kp2+1:kp2+hi2)=HR_ostatni;
end
end
%%
save H H
save HR_vsechny_PVC HR_vsechny_PVC
save HR_vsechny_ostatni HR_vsechny_ostatni
save FF_vsechny_PVC FF_vsechny_PVC
save FF_vsechny_ostatni FF_vsechny_ostatni
clear all
clc
close all
else
%       sprintf('spuste,prosim,dalsi signal,v tomto signalu PVC nejsou')
    H(:,q-3)=0; % v signalu neobsahujicim extrasystoly budou vsude nuly ve sloupci matice H
    save H H
end
end
load('H.mat')
mean_senz=sum(H(4,:))/(sum(H(4,:))+sum(H(5,:)))% stredni hodnota senzetivity
mean_spec=sum(H(7,:))/(sum(H(7,:))+sum(H(8,:))) % stredni hodnota specificity

% mean_senz =
% 
%     0.9022
% 
% 
% mean_spec =
% 
%     0.6972
%%
% load('HR_vsechny_PVC')
% load('HR_vsechny_ostatni')
% load('FF_vsechny_PVC')
% load('FF_vsechny_ostatni')
% %  boxplot FF
% FF_vsechny_PVC(FF_vsechny_PVC==0)=[];
% FF_vsechny_ostatni(FF_vsechny_ostatni==0)=[];
% figure(2)
% p11 = (FF_vsechny_PVC)';
% p22 = (FF_vsechny_ostatni)';
% xxp = [p11;p22];
% gip = [ones(size(p11)); 2*ones(size(p22))];
% boxplot(xxp,gip,'labels',{'Komorove extrasystoly' 'Ostatni'})
% title('Tvarovy faktor')
% %  boxplot HR
% HR_vsechny_PVC(HR_vsechny_PVC==0)=[];
% HR_vsechny_ostatni(HR_vsechny_ostatni==0)=[];
% figure(3)
% p1 = (HR_vsechny_PVC)';
% p2 = (HR_vsechny_ostatni)';
% xx = [p1;p2];
% gi = [ones(size(p1)); 2*ones(size(p2))];
% boxplot(xx,gi,'labels',{'Komorove extrasystoly' 'Ostatni'})
% ylabel('[tepu/min]')
% title('Tepova frekvence')
%%




