clear all;close all;clc;
%% Program pro detekci K-komplexů s využitím S-transformace-odhad obálky umocněním na druhou
% zvolená podmínka sklon a extrém
% Diplomová práce; Kristýna Hlaváčová

load databaze % načtení databáze spánkových signálů EEG
fvz = 200; % vzorkovací frekvence
x = databaze(:,1);% VYBER SIGNÁL Z DATABÁZE!!!!!!!

% Vyber hodnocení experta dle příslušného signálu, který byl vybrán výše
% pro výběr hodnocení druhého experta odkomentuj příslušné proměnné
expert1=importdata('Visual_scoring1_excerpt1.txt');% ... 1.expert (1-10)
% expert2=importdata('Visual_scoring2_excerpt1.txt');% ... 2.expert (1-5)


pom1=expert1.data(:,1); pom1=(pom1)';% výběr začátku trvaní K-komplexu
% pom2=expert2.data(:,1); pom2=(pom2)';
komplexStart1=round(pom1*fvz);% indexy počátku K-komplexu
% komplexStart2=round(pom2*fvz);

trvani1=expert1.data(:,2); trvani1=(trvani1)';% výběr začátku doby trvaní K-komplexu
% trvani2=expert2.data(:,2); trvani2=(trvani2)';
komplexDuration1=round(trvani1*fvz);% indexy doby trvaní K-komplexu
% komplexDuration2=round(trvani2*fvz);

komplexEnd1=komplexStart1+komplexDuration1;% indexy konce K-komplexu
% komplexEnd2=komplexStart2+komplexDuration2;

%% Stockwellova transformace

% nastavení parametrů pro tranformaci
step_t = 1;      	% 1 vz
step_f = .1;     	% 1, .1 Hz
f_od = .5;           % .5           
f_do = 2;         % 0   2 Hz pro EEG
fs = 200;

[ST,t,f] = S_tr_time(x,fs,f_od,f_do,step_t,step_f);% Stockwellova transformace
t=((1:length(x))./fvz); % výpočet časové osy

figure(1)% vykreslení původního signálu s vyznačenými začátky K-komplexů
plot(t,x); hold on
for k = 1:length(komplexStart1)
    line([t(komplexStart1(k)) t(komplexStart1(k))],[-100 100],'color','r')
end
% for k = 1:length(komplexStart2)
%     line([t(komplexStart2(k)) t(komplexStart2(k))],[-100 100],'color','k')
% end
hold off
title('Původní signál s vyznačenými počátky K-komplexů','FontSize', 18),xlabel('Čas [s]','FontSize', 18),ylabel('Výchylka [µV]','FontSize', 18);
legend('Signál','Expert1');
% legend('Signál','Expert2');

figure(2)% vykreslení spektrogramu S-transformace
imagesc(t, f, abs(ST)),title('Spektrogram původního signálu','FontSize', 18);
xlabel('Čas [s]','FontSize', 18),ylabel('Frekvence [Hz]','FontSize', 18);

%% Výpočet obálky umocněním realné složky S-transformace na druhou

STreal =abs(real(ST)); % absolutní hodnota reálné složky po výpočtu S-transformace
STreal =sum(STreal); % suma sloupců

STnorm =(STreal) ./(abs(max(STreal))); % 1.normalizace 
SE = STnorm.^2; % výpočet odhadu obálky umocněním na druhou
SE =(SE) ./(abs(max(SE))); % 2.normalizace 

figure(3)% vykreslení obálky
subplot 211
plot(t,SE)
title('Obálka signálu - umocnění reálné složky S-transformace na druhou','FontSize', 18),xlabel('Čas [s]','FontSize', 18),ylabel('Normalizovaná výchylka [-]','FontSize', 18);
subplot 212
plot(t,SE); hold on
for k = 1:length(komplexStart1)
    line([t(komplexStart1(k)) t(komplexStart1(k))],[-100 100],'color','r')
end
% for k = 1:length(komplexStart2)
%     line([t(komplexStart2(k)) t(komplexStart2(k))],[-100 100],'color','k')
% end
hold off
title('Obálka signálu s vyznačenými počátky K-komplexů','FontSize', 18),xlabel('Čas [s]','FontSize', 18),ylabel('Normalizovaná výchylka [-]','FontSize', 18);
legend('Signál','Expert1');
% legend('Signál','Expert2');
ylim([0 1]);

%% Použití prahu pro detekci vrcholů v obálce umocněním na druhou

okno=30*fvz;% jelikož EEG bývá zkoumáno v rámci 30 s

SmerodatOdch=[];
for j=1:okno:length(t)  % výpočet směrodatné odchylky pro práh v daném okně
         Odch=std(SE(j:(j+okno-1)));
         SmerodatOdch=[SmerodatOdch Odch];  
end
prah=SmerodatOdch*4; % násobek směrodatné odchylky jako konečný práh pro detekci 

bx=[]; by=[];
pomx=[]; pomy=[];
vyber_usek=1:okno:length(t);
for l=1:length(prah)  
         [pomy,pomx]=findpeaks(SE(vyber_usek(l):(vyber_usek(l)+okno-1)),'MinPeakHeight',prah(l),'MinPeakDistance',fvz);
          pomx=find(ismember(SE,pomy));
          bx=[bx pomx];
end

figure(4)  % vykreslení detekovaných vrcholů s jednotlivými prahy v úseku 30s
plot(t,SE); hold on
for ll = 1:length(prah)
    line(t(vyber_usek(ll):(vyber_usek(ll)+okno-1)),prah(ll),'color','g'); hold on
end
plot(t,SE,t(bx),SE(bx),'mo');hold on
for k = 1:length(komplexStart1)
    line([t(komplexStart1(k)) t(komplexStart1(k))],[-100 100],'color','r')
end
% for k = 1:length(komplexStart2)
%     line([t(komplexStart2(k)) t(komplexStart2(k))],[-100 100],'color','k')
% end
hold off
hold off
hold off
title('Obálka s vyznačenými vrcholy a postupnými prahy','FontSize', 18),xlabel('Čas [s]','FontSize', 18),ylabel('Normalizovaná výchylka [-]','FontSize', 18);
ylim([0 1]);

%% DETEKCE K-komplexů ve FILTROVANÉM SIGNÁLU

propust = fir1(1000,[0.5 2]/100); % pásmová propust 0.5 - 2 Hz
xf = filtfilt(propust,1,x); % filtrace signálu pásmovou propustí

cs=xor((xf(1:end-1)>=0),(xf(2:end)>=0));% zerocrossing - detekce průchodů nulou
dx=find(cs);

% Detekce maxim a minim ve filtrovaném signálu mezi průchody nulou
pomminX=[]; pommaxX=[];
minX=[]; maxX=[];
neg=[]; pos=[];
for nn=2:length(dx)-1
    mez1=xf(dx(nn-1)+1:dx(nn)-1);               % interval pro minimum
    if length(mez1)>5
       pomminX=min(mez1(:));
        if pomminX<0
           mez3=xf((dx(nn)+1):(dx(nn+1)-1));   % interval pro maximum
            if length(mez3)>5
               pommaxX=max(mez3(:));        
               if pommaxX>0
                 pomminX=find(ismember(xf,pomminX));
                 pommaxX=find(ismember(xf,pommaxX));
                  minX=[minX pomminX];
                  neg=[neg dx(nn-1) dx(nn)];
                  maxX=[maxX pommaxX];
                  pos=[pos dx(nn) dx(nn+1)];
               end   
            end
        end
     end
end

maxX=unique(maxX);minX=unique(minX);
pos=unique(pos);neg=unique(neg);

figure(5)
plot(t,x,'g'); hold on % původní signál
plot(t,xf,'k');hold on % filtrovaný signál pásmovou propustí
for k = 1:length(komplexStart1)
    line([t(komplexStart1(k)) t(komplexStart1(k))],[-100 100],'color','r')
end
% for k = 1:length(komplexStart2)
%     line([t(komplexStart2(k)) t(komplexStart2(k))],[-100 100],'color','k')
% end
hold off
hold off
title('Původní a filtrovaný signál pásmovou propustí','FontSize', 18),xlabel('Čas [s]','FontSize', 18),ylabel('Výchylka [µV]','FontSize', 18);
legend('Původní signál','Filtrovaný signál','Expert1');
%legend('Původní signál','Filtrovaný signál','Expert2');

figure(6) % vykreslení detekovaných maxim a minim ve filtrovaném signálu
subplot 211
plot(t, xf, '-b');hold on
plot(t(dx), xf(dx), 'rp')
hold off
grid
legend('Filtrovaný signál', 'Průchody nulou')
title('Filtrovaný signál s vyznačenými průchody nulou','FontSize', 18),xlabel('Čas [s]','FontSize', 18),ylabel('Výchylka [µV]','FontSize', 18);

subplot 212
plot(t, xf)
hold on
plot(t,xf,t(maxX),xf(maxX),'mo');hold on
plot(t,xf,t(minX),xf(minX),'go')
hold off
title('Filtrovaný signál s detekcí maxim a minim','FontSize', 18),xlabel('Čas [s]','FontSize', 18),ylabel('Výchylka [µV]','FontSize', 18);

trvanipos=[];% trvání pozitivní doby
trvaniposInd=[];
  for gp=1:2:length(pos)-1 % vytvoření intervalů mezi 2 vrcholy
     trvanipos(gp)=pos(gp+1)-pos(gp);
     trvaniposInd=[trvaniposInd pos(gp) pos(gp+1)];
  end
trvanipos=trvanipos(trvanipos~=0);

trvanineg=[];% trvání negativní doby
trvaninegInd=[];
  for gn=1:2:length(neg)-1 % vytvoření intervalů mezi 2 vrcholy
     trvanineg(gn)=neg(gn+1)-neg(gn);
     trvaninegInd=[trvaninegInd neg(gn) neg(gn+1)];
  end
trvanineg=trvanineg(trvanineg~=0);

nP=[neg pos];
trvaninP=[trvanineg trvanipos];
trvanitotal=[]; % trvání celkové doby tzn. negativní a pozitivní dohromady
trvaninPInd=[];
for gt=1:2:length(neg)-1 % vytvoření intervalů mezi 2 vrcholy
     trvanitotal(gt)=pos(gt+1)-neg(gt);
     trvaninPInd=[trvaninPInd neg(gt) pos(gt+1)];
end
  
trvanitotal=trvanitotal(trvanitotal~=0);
trvaninPInd=unique(trvaninPInd);
intervalTotal=[];
for xv=1:length(trvanitotal) % výběr těch intervalů, které splňují danou podmínku o trvání pozitivní a negativní doby
    if (trvanitotal(xv)>=0.5*fvz && trvanitotal(xv)<=1.5*fvz)&&((trvanineg(xv)>trvanitotal(xv)/4)&&(trvanipos(xv)>trvanitotal(xv)/4))
        intervalTotal=[intervalTotal, trvanitotal(xv)];
    end
end

gb=ismember(trvanitotal,intervalTotal);
indexik=find(gb==1);
trvaniInd=[];
for lp=1:length(indexik)
    trvaniInd=[trvaniInd trvaninPInd(indexik(lp)) trvaninPInd(indexik(lp)+1)];
end

vyslmaxX=[];vyslminX=[];  % nalezená maxima a minima v intervalech, která splňují výše uvedenou podmínku
for hh=1:length(indexik)
    aa=indexik(hh);
    bb=indexik(hh)+1;
    zerosa=trvaninPInd(aa);
    zeroha=trvaninPInd(bb);
    [rowa,cola]=find((maxX>zerosa & maxX<zeroha));
    vyslmaxX=[vyslmaxX maxX(cola)];
    [rowi,coli]=find(minX>zerosa & minX<zeroha);
     vyslminX=[vyslminX minX(coli)];
   
end
vyslmaxY=xf(vyslmaxX);vyslminY=xf(vyslminX);

figure(7)% vykreslení maxim a minim, která vyhovují podmínce - doba trvání total je mezi 0,5s a 1,5s;
%trvání negativního a pozitivního intervalu je větší než 1/4 trvaní total
%intervalu
subplot 211
plot(t, xf);hold on
plot(t,xf,t(vyslmaxX),xf(vyslmaxX),'mo');hold on
plot(t,xf,t(vyslminX),xf(vyslminX),'go')
hold off
title('Filtrovaný signál s detekcí maxim a minim splňující podmínku trvání','FontSize', 18),xlabel('Čas [s]','FontSize', 18),ylabel('Výchylka [µV]','FontSize', 18);

%% Výpočet podmínky o hodnotě SKLONU mezi extrémními výchylkami a výpočet podmínky o hodnotě mezi EXTRÉMY

novemaxX=[];
for u=1:length(vyslmaxX)% výpočet sklonu a extrému
        if (((abs(vyslmaxY(u)-vyslminY(u)))/((((vyslmaxX(u))/fvz)-((vyslminX(u))/fvz))))>=(75/0.5))&& ((vyslmaxY(u)+abs(vyslminY(u)))>=(75));
             novemaxX=[novemaxX vyslmaxX(u)];
        end
end 

figure(7)% vykreslení vrcholů splňujících podmínku sklonu 
subplot 212
plot(t, xf); hold on
plot(t,xf,t(novemaxX),xf(novemaxX),'mo'); hold on
for k = 1:length(komplexStart1)
    line([t(komplexStart1(k)) t(komplexStart1(k))],[-100 100],'color','r')
end
% for k = 1:length(komplexStart2)
%     line([t(komplexStart2(k)) t(komplexStart2(k))],[-100 100],'color','k')
% end
hold off
hold off
title('Filtrovaný signál s detekovanými vrcholy splňující podmínku sklonu a extrému','FontSize', 18),xlabel('Čas [s]','FontSize', 18),ylabel('Výchylka [µV]','FontSize', 18);

%% Porovnání výsledku - konečné rozhodnutí o K-komplexech

novemaxX;bx=bx';
odVrcholu=1*fvz;
konecneX=[];

for w=1:length(bx)% hledám vrcholy (z detekce Stockwell a filtrovaného signálu), které si odpovídají v okně 1s
    range=(bx(w)-odVrcholu:bx(w)+odVrcholu);
    vrchol=ismember(novemaxX,range);
    
    if vrchol(vrchol)==1
         konecneX=[konecneX bx(w)];
    end
end


figure(8)% vykreslení detekovaných vrcholů po porovnání
plot(t,SE); hold on
for k = 1:length(komplexStart1)
    line([t(komplexStart1(k)) t(komplexStart1(k))],[-100 100],'color','r')
end
% for o = 1:length(komplexEnd1)
%      line([t(komplexEnd1(o)) t(komplexEnd1(o))],[-100 100],'color','y')
% end
% for k = 1:length(komplexStart2)
%     line([t(komplexStart2(k)) t(komplexStart2(k))],[-100 100],'color','k')
% end
% for o = 1:length(komplexEnd2)
%      line([t(komplexEnd2(o)) t(komplexEnd2(o))],[-100 100],'color','c')
% end
plot(t,SE,t(konecneX),SE(konecneX),'mo')
hold off
ylim([0 1]);
title('Obálka umocněním na druhou s detekovanými vrcholy po porovnání','FontSize', 18),xlabel('Čas [s]','FontSize', 18),ylabel('Normalizovaná výchylka [-]','FontSize', 18);
