% METODA TEPOVE FREKVENCE, TVAROVEHO FAKTORU A PRAVIDLA NEJBLIZSIHO SOUSEDA
%vsechny signaly jsou na CD 
clear all
clc
close all
senz_spec=zeros(10,46);
save senz_spec senz_spec % matice vysledku
%4:49 vse 36 signalu
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('senz_spec.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);
Cpred(1,j)=A(j)-A(j-1);
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
%% 19 vzorku 50 ms -> konec QRS je detekovan
dvapetktime=time(1,Rvse(1,j)+ymaxdvapet+19);
offsettime(1,j)=dvapetktime;
% pozice_a2(1,j)=Rvse(1,j)+ymaxdvapet+19;
pozice_a2(1,j)=Rvse(1,j)+ymaxdvapet;
%% 
% hold on
% plot(dvapetktime,x(Rvse(1,j)+ymaxdvapet+19,1),'*g')
end
%% kratke nazvy deskriptoru, viz tabulka 3.1 v BP
for j=3:length(Rvse)-2
pik(1,j)=x(Rvse(1,j),1);
if pik(1,j)>=0
    pozitivpik(1,j)=pik(1,j); %Pp   
    negativpik(1,j)=x(pozice_a2(1,j),1);
 else
     pozitivpik(1,j)=x(pozice_a2(1,j)+10,1);
     negativpik(1,j)=pik(1,j); %Pn
    end
end

for j=3:length(Rvse)-2
   if x(pozice_a1(1,j),1)>0  
        pozitivpik(1,j)=pozitivpik(1,j)-posun(1,j);
     negativpik(1,j)=negativpik(1,j)-posun(1,j); %Pn  
     if x(pozice_a1(1,j),1)<0 
         pozitivpik(1,j)=pozitivpik(1,j)+posun(1,j); %Pp   
    negativpik(1,j)=negativpik(1,j)+posun(1,j);
     if x(pozice_a1(1,j),1)==0 
         pozitivpik(1,j)=pozitivpik(1,j); %Pp   
    negativpik(1,j)=negativpik(1,j);
     end
   end
end         
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
for i=1:length(FF)
    HR(1,i)=60/Cpred(1,i) ;
end
else
end

if length(pozice_V)>0
%% Vytvoreni matice deskriptoru

matice_discr(:,1)=FF;
matice_discr(:,2)=HR;
%% krizova validace a klasifikace
rng(0,'twister'); % nastaveni, aby generovani zacinalo vzdy stejne
nfold=10;

paramOut=type_1;
param=matice_discr;

cv = cvpartition(paramOut,'kfold',nfold);  %paramOut je vektor, kde jsou ulozeny anotace (tj. zda je to PVC nebo ostatni-musi byt ve forme cisel - tj. PVC nahradte treba 1 a zbytek 2
indices = zeros(size(paramOut));
for i=1:nfold
    indices(cv.test(i)) = i;
end

type='euclidean'; % nastaveni typu pocitani vzdalenosti
k=1; % pocet uvazovanych nejblizsich sousedu
accInit=zeros(nfold,1); % deficnice promnene pro ukladani uspesnosti klasifikace s originalnimi hodnotami parmametru
accNorm = zeros(nfold,1); % deficnice promnene pro ukaldani uspesnosti klasifikace s normalizovanymi hodnotami parmametru

for i=1:nfold % pro kazdou iteraci
    testIdx = (indices == i); trainIdx = ~testIdx; % vyber jedne aktualni casti jako testovaci,zbytek bude trenovaci
    ptr=param(trainIdx,:); %  aktualni trenovaci cast matice matice_discr
    
    mitr=min(ptr); % minimum z treninkoveho setu
    matr=max(ptr); % maximum z treninkoveho setu
    ptrnorm=(ptr-repmat(mitr,size(ptr,1),1))./(repmat(matr,size(ptr,1),1)-repmat(mitr,size(ptr,1),1)); % normalizace treninovaci mnoziny
    
   
    pte=param(testIdx,:); % vyber aktualni testovaci mnoziny
    ptenorm=(pte-repmat(mitr,size(pte,1),1))./(repmat(matr,size(pte,1),1)-repmat(mitr,size(pte,1),1)); % normalizace testovaci mnoziny
       
    
    knnClassInit = knnclassify(pte, ptr, paramOut(trainIdx), k, type); % trenovani a klasifikace - s originalnimi hodnotami
    knnClassNorm = knnclassify(ptenorm,ptrnorm,paramOut(trainIdx),k,type); % trenovani a klasifikace - s normalizovanymi hodnotami
  
    badInit = ~strcmp(knnClassInit,paramOut(testIdx)); % vypocet spatne klasifikovanych cyklu
    badNorm = ~strcmp(knnClassNorm,paramOut(testIdx));
    N=sum(testIdx==1);
    accInit(i) = 1-sum(badInit)/N; % vypocet uspesnosti
    accNorm(i) = 1-sum(badNorm)/N;
    [aaknnInit,grpInit] = confusionmat(paramOut(testIdx),knnClassInit) % zobrazeni matice zmen
    [aaknnNorm,grpNorm] = confusionmat(paramOut(testIdx),knnClassNorm)
    usp(1:2,1:2,i)=aaknnNorm;
end
accKNNInit1 = mean(accInit) % presnost
accKNNNorm1 = mean(accNorm)
for i=1:nfold
if usp(1,1,i)==usp(2,1,i)&& usp(1,2,i)==usp(2,2,i)
usp(1,1,i)=0;
usp(1,2,i)=0;
usp(2,1,i)=0;
end
end
senzetivita=(sum(usp(1,1,:)))/(sum(usp(1,1,:))+(sum(usp(1,2,:)))); 
specificita=(sum(usp(2,2,:)))/(sum(usp(2,2,:))+(sum(usp(2,1,:))));
pocet_PVC=sum(type_1==1);
pocet_detek_PVC=sum(usp(1,1,:));
pocet_ostatnich=sum(type_1==2);
pocet_sprdetek_ost=sum(usp(2,2,:));
senz_spec(1,q-3)=senzetivita;
senz_spec(2,q-3)=specificita;
senz_spec(3,q-3)=accKNNNorm1;%presnost
senz_spec(4,q-3)=pocet_PVC;% celkovy pocet PVC v signalu
senz_spec(5,q-3)=pocet_detek_PVC; % pocet detekovanych spravne PVC
senz_spec(6,q-3)=sum(usp(1,2,:));%zarazene do ostatnich,ale meli by patrit k PVC
senz_spec(7,q-3)=pocet_ostatnich;% celkovy pocet ostatnich patologii v signalu
senz_spec(8,q-3)=pocet_sprdetek_ost;% pocet detekovanych spravne ostatnich
senz_spec(9,q-3)=sum(usp(2,1,:));% zarazene do PVC, ale meli byt zarazeny do Ostatnich
senz_spec(10,:)=[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];
% 10 radek obsahije cisla signalu
save senz_spec senz_spec
clear all
clc
close all
else
%     sprintf('spuste,prosim,dalsi signal,v tomto signalu PVC nejsou')
    senz_spec(:,q-3)=0; % signal neobsahuje PVC
    save senz_spec senz_spec
end

end
load('senz_spec.mat')
%%
mean_senz=sum(senz_spec(5,:))/(sum(senz_spec(5,:))+sum(senz_spec(6,:)))% stredni hodnota senzetivity
mean_spec=sum(senz_spec(8,:))/(sum(senz_spec(8,:))+sum(senz_spec(9,:))) % stredni hodnota specificity

% mean_senz =
% 
%     0.6088
% 
% 
% mean_spec =
% 
%     0.9646
%%
