clc
clear all
close all
%%%%% input data feeding
load ('Sz1.mat')
x1=EEG;

load('Sz2.mat')
x2=EEG;

load('Sz3.mat')
x3=EEG;

load('Sz4.mat')
x4=EEG;

load('Sz5.mat')
x5=EEG;





Fs=512; %%% sampling frequency
WW=1; %%% Window length
num_channel=size(x1,2);
num_channel_select=23;





%%%%% Separating the input data of each channel

for i=1:num_channel
    g1{i}=x1(:,i);
    g2{i}=x2(:,i);
    g3{i}=x3(:,i);
    g4{i}=x4(:,i);
    g5{i}=x5(:,i);
   
    
end



for i=1:num_channel
    g_tot{i}=[g1{i};g2{i};g3{i};g4{i};g5{i}];
end


%%% Picking up the selected channels (it must be uncommented after
%%% channel_selelection)

% for i=1:num_channel_select
%     g_tot_new{i}=g_tot{channel_rank_index_ascend};
% end


%%%%% Defining the length of signal and labeling

num1_window=length(g1{1})/512;
num2_window=length(g2{1})/512;
num3_window=length(g3{1})/512;
num4_window=length(g4{1})/512;
num5_window=length(g5{1})/512;



num1_seizure=num1_window-360;
num2_seizure=num2_window-360;
num3_seizure=num3_window-360;
num4_seizure=num4_window-360;
num5_seizure=num5_window-360;


num_sample_tot=size(g_tot{1},1);
num_window_tot=num_sample_tot/(Fs*WW);

num_seizure_tot=num1_seizure+num2_seizure+num3_seizure+num4_seizure+num5_seizure;

sz1_start=181;
sz1_end=num1_window-180;

sz2_start=num1_window+181;
sz2_end=num1_window+num2_window-180;

sz3_start=num1_window+num2_window+181;
sz3_end=num1_window+num2_window+num3_window-180;

sz4_start=num1_window+num2_window+num3_window+181;
sz4_end=num1_window+num2_window+num3_window+num4_window-180;

sz5_start=num1_window+num2_window+num3_window+num4_window+181;
sz5_end=num1_window+num2_window+num3_window+num4_window+num5_window-180;



basicVector=[zeros(1,180), ones(1,num1_seizure), zeros(1,360), ones(1,num2_seizure), zeros(1,360),ones(1,num3_seizure), zeros(1,360),ones(1,num4_seizure), zeros(1,360),ones(1,num5_seizure),zeros(1,180)];



%%% Average the value of all channels

x_concat=cat(3,g_tot{:});
input_data_tot=mean(x_concat,3);


%%% Average the value of selected channels (uncommented after channel
%%% selection)

% x_concat=cat(3,g_tot_new{:});
% input_data_tot=mean(x_concat,3);


%%%% Windowing the input signal with channel selection

k=1;
for nn=1:Fs*WW:num_sample_tot-Fs*WW+1

  y_final{k}=input_data_tot(nn:nn+Fs*WW-1);
  k=k+1;
end



t_window_tot=1:num_window_tot;

t1_seizure=180+WW:WW:num1_window-180/WW;
t2_seizure=num1_window+180+WW:WW:num1_window+num2_window-180/WW;
t3_seizure=num1_window+num2_window+180+WW:WW:num1_window+num2_window+num3_window-180/WW;
t4_seizure=num1_window+num2_window+num3_window+180+WW:WW:num1_window+num2_window+num3_window+num4_window-180/WW;
t5_seizure=num1_window+num2_window+num3_window+num4_window+180+WW:WW:num1_window+num2_window+num3_window+num4_window+num5_window-180/WW;


 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Coastline feature extraction (feature 1) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:num_window_tot
    for o=2:Fs*WW
        y_subtract{i}(o)=abs(y_final{i}(o)-y_final{i}(o-1));
        y_subtract{i}(1)=0;
    end 
end

for i=1:num_window_tot
    F_coastline(i)=sum(y_subtract{i},'all');
end


F_coastline_baseline=[F_coastline(1,1:sz1_start-1) , F_coastline(1,sz1_end+1:sz2_start-1),F_coastline(1,sz2_end+1:sz3_start-1) , F_coastline(1,sz3_end+1:sz4_start-1),F_coastline(1,sz4_end+1:sz5_start-1) , F_coastline(1,sz5_end+1:num_window_tot)];

    
F_coastline_seizure=[F_coastline(1,sz1_start:sz1_end) , F_coastline(1,sz2_start:sz2_end) , F_coastline(1,sz3_start:sz3_end) , F_coastline(1,sz4_start:sz4_end),F_coastline(1,sz5_start:sz5_end)];


F1_coastline_seizure=F_coastline_seizure(1,1:num1_seizure);
F2_coastline_seizure=F_coastline_seizure(1,num1_seizure+1:num1_seizure+num2_seizure);
F3_coastline_seizure=F_coastline_seizure(1,num1_seizure+num2_seizure+1:num1_seizure+num2_seizure+num3_seizure);
F4_coastline_seizure=F_coastline_seizure(1,num1_seizure+num2_seizure+num3_seizure+1:num1_seizure+num2_seizure+num3_seizure+num4_seizure);
F5_coastline_seizure=F_coastline_seizure(1,num1_seizure+num2_seizure+num3_seizure+num4_seizure+1:num1_seizure+num2_seizure+num3_seizure+num4_seizure+num5_seizure);


%%%% defining threshold of coastline feature

threshold_coastline_baseline= mean(F_coastline_baseline);
threshold_coastline_seizure= mean(F_coastline_seizure);
threshold_coastline=(4*threshold_coastline_baseline + 1*threshold_coastline_seizure)/5;

%%% Definnig feature discriminative parameter (FDP)

FDP_coastline=threshold_coastline_seizure/threshold_coastline_baseline;


%%% plotting the results of Coastline feature

figure
plot(t_window_tot,F_coastline,'linewidth',3)
title(' Coastline feature extraction')
xlabel('Time (sec)')
ylabel('Coastline Amplitude')


hold on

plot(t1_seizure,F1_coastline_seizure,t2_seizure,F2_coastline_seizure,t3_seizure,F3_coastline_seizure,t4_seizure,F4_coastline_seizure,t5_seizure,F5_coastline_seizure,'linewidth',3,'color','r')

hold off

%%%% Coastline detection

for i=1:num_window_tot
    if F_coastline(i)>=threshold_coastline
        det_coastline(i)=1;
    else
        det_coastline(i)=0;
    end
end

coastline_num_seizure=sum(det_coastline==1);
coastline_num_baseline=sum(det_coastline==0);

FP_coastline=sum(basicVector-det_coastline==-1);
FN_coastline=sum(basicVector-det_coastline==1);

TP_coastline=coastline_num_seizure-FP_coastline;
TN_coastline=coastline_num_baseline-FN_coastline;

sensitivity_coastline= TP_coastline/(TP_coastline+FN_coastline)
specificity_coastline=TN_coastline/(TN_coastline+FP_coastline)
accuracy_coastline=(sensitivity_coastline+specificity_coastline)/2


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Energy feature extraction (feature2) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


for i=1:num_window_tot
    y_square{i}=(y_final{i}).^2;
    F_energy(i)=mean(y_square{i});
end


%%% Determining values of baseline and seizure states

F_energy_baseline=[F_energy(1,1:sz1_start-1) , F_energy(1,sz1_end+1:sz2_start-1),F_energy(1,sz2_end+1:sz3_start-1) , F_energy(1,sz3_end+1: sz4_start-1), F_energy(1,sz4_end+1:sz5_start-1) , F_energy(1,sz5_end+1:num_window_tot)];

F_energy_seizure=[F_energy(1,sz1_start:sz1_end) , F_energy(1,sz2_start:sz2_end) , F_energy(1,sz3_start:sz3_end) , F_energy(1,sz4_start:sz4_end),F_energy(1,sz5_start:sz5_end)];

%%% Seperating the energy feature for each seizure occurence                  
                 
F1_energy_seizure=F_energy_seizure(1,1:num1_seizure);
F2_energy_seizure=F_energy_seizure(1,num1_seizure+1:num1_seizure+num2_seizure);
F3_energy_seizure=F_energy_seizure(1,num1_seizure+num2_seizure+1:num1_seizure+num2_seizure+num3_seizure);
F4_energy_seizure=F_energy_seizure(1,num1_seizure+num2_seizure+num3_seizure+1:num1_seizure+num2_seizure+num3_seizure+num4_seizure);
F5_energy_seizure=F_energy_seizure(1,num1_seizure+num2_seizure+num3_seizure+num4_seizure+1:num1_seizure+num2_seizure+num3_seizure+num4_seizure+num5_seizure);

%%%% defining threshold of enregy feature

threshold_energy_baseline= mean(F_energy_baseline);
threshold_energy_seizure= mean(F_energy_seizure);
threshold_energy= (4*threshold_energy_baseline + 1*threshold_energy_seizure)/5;

%%% Definnig feature discriminative parameter (FDP)

FDP_energy=sqrt(threshold_energy_seizure/threshold_energy_baseline);

%%% plotting results of Energy feature

figure
plot(t_window_tot,F_energy,'linewidth',3)
title(' Energy feature extraction')
xlabel('Time (sec)')
ylabel('Eenrgy Feature result')


hold on

plot(t1_seizure,F1_energy_seizure,t2_seizure,F2_energy_seizure,t3_seizure,F3_energy_seizure,t4_seizure,F4_energy_seizure,t5_seizure,F5_energy_seizure,'linewidth',3,'color','r')

hold off

%%%% energy detection

for i=1:num_window_tot
    if F_energy(i)>=threshold_energy
        det_energy(i)=1;
    else
        det_energy(i)=0;
    end
end

energy_num_seizure=sum(det_energy==1);
energy_num_baseline=sum(det_energy==0);

FP_energy=sum(basicVector-det_energy==-1);
FN_energy=sum(basicVector-det_energy==1);

TP_energy=energy_num_seizure-FP_energy;
TN_energy=energy_num_baseline-FN_energy;

sensitivity_energy= TP_energy/(TP_energy+FN_energy)
specificity_energy=TN_energy/(TN_energy+FP_energy)
accuracy_energy=(sensitivity_energy+specificity_energy)/2


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Nonlinear Energy feature extraction (feature 3) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:num_window_tot
    for o=2:Fs*WW-1
    
    NE{i}(o)=(y_final{i}(o).^2)-y_final{i}(o-1)*y_final{i}(o+1);
    NE{i}(1)=0;
    NE{i}(Fs*WW)=0;
    end
end
 
for i=1:num_window_tot
    F_NE(i)=mean(NE{i});
end


%%%% Determining baseline and seizure states for nonlinear Energy feature

F_NE_baseline=[F_NE(1,1:sz1_start-1) , F_NE(1,sz1_end+1:sz2_start-1),F_NE(1,sz2_end+1:sz3_start-1) , F_NE(1,sz3_end+1: sz4_start-1), F_NE(1,sz4_end+1:sz5_start-1) , F_NE(1,sz5_end+1:num_window_tot)];

F_NE_seizure=[F_NE(1,sz1_start:sz1_end) , F_NE(1,sz2_start:sz2_end) , F_NE(1,sz3_start:sz3_end) , F_NE(1,sz4_start:sz4_end),F_NE(1,sz5_start:sz5_end)];

%%%%%%%%%%%%% Determining baseline and seizure states for nonlinear Energy

F1_NE_seizure=F_NE_seizure(1,1:num1_seizure);
F2_NE_seizure=F_NE_seizure(1,num1_seizure+1:num1_seizure+num2_seizure);
F3_NE_seizure=F_NE_seizure(1,num1_seizure+num2_seizure+1:num1_seizure+num2_seizure+num3_seizure);
F4_NE_seizure=F_NE_seizure(1,num1_seizure+num2_seizure+num3_seizure+1:num1_seizure+num2_seizure+num3_seizure+num4_seizure);
F5_NE_seizure=F_NE_seizure(1,num1_seizure+num2_seizure+num3_seizure+num4_seizure+1:num1_seizure+num2_seizure+num3_seizure+num4_seizure+num5_seizure);

%%%%%% Defining threshold of Nonlinear Energy feature

threshold_NE_baseline= mean(abs(F_NE_baseline));
threshold_NE_seizure= mean(abs(F_NE_seizure));
threshold_NE=(1*threshold_NE_baseline + 1*threshold_NE_seizure)/2;
%%% Definnig feature discriminative parameter (FDP)

FDP_NE=sqrt(threshold_NE_seizure/threshold_NE_baseline);

%%% plotting results of nonlinear Energy feature

figure
plot(t_window_tot,F_NE,'linewidth',3)
title(' Nonlinear feature extraction')
xlabel('Time (sec)')
ylabel('Nonlinear Feature result')


hold on

plot(t1_seizure,F1_NE_seizure,t2_seizure,F2_NE_seizure,t3_seizure,F3_NE_seizure,t4_seizure,F4_NE_seizure,t5_seizure,F5_NE_seizure,'linewidth',3,'color','r')

hold off

%%%% Nonlinear energy detection

for i=1:num_window_tot
    if abs(F_NE(i))>=threshold_NE
        det_NE(i)=1;
    else
        det_NE(i)=0;
    end
end

NE_num_seizure=sum(det_NE==1);
NE_num_baseline=sum(det_NE==0);

FP_NE=sum(basicVector-det_NE==-1);
FN_NE=sum(basicVector-det_NE==1);

TP_NE=NE_num_seizure-FP_NE;
TN_NE=NE_num_baseline-FN_NE;

sensitivity_NE= TP_NE/(TP_NE+FN_NE)
specificity_NE=TN_NE/(TN_NE+FP_NE)
accuracy_NE=(sensitivity_NE+specificity_NE)/2

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Variance feature extraction (feature 4) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:num_window_tot
    
    F_var(i)=var(y_final{i});
end

%%%% Determining baseline and seizure states of Variance feature

F_var_baseline=[F_var(1,1:sz1_start-1) , F_var(1,sz1_end+1:sz2_start-1),F_var(1,sz2_end+1:sz3_start-1) , F_var(1,sz3_end+1: sz4_start-1), F_var(1,sz4_end+1:sz5_start-1) , F_var(1,sz5_end+1:num_window_tot)];

F_var_seizure=[F_var(1,sz1_start:sz1_end) , F_var(1,sz2_start:sz2_end) , F_var(1,sz3_start:sz3_end) , F_var(1,sz4_start:sz4_end),F_var(1,sz5_start:sz5_end)];

%%%%%%%%%%%%% Determining baseline and seizure states for Variance Feature

F1_var_seizure=F_var_seizure(1,1:num1_seizure);
F2_var_seizure=F_var_seizure(1,num1_seizure+1:num1_seizure+num2_seizure);
F3_var_seizure=F_var_seizure(1,num1_seizure+num2_seizure+1:num1_seizure+num2_seizure+num3_seizure);
F4_var_seizure=F_var_seizure(1,num1_seizure+num2_seizure+num3_seizure+1:num1_seizure+num2_seizure+num3_seizure+num4_seizure);
F5_var_seizure=F_var_seizure(1,num1_seizure+num2_seizure+num3_seizure+num4_seizure+1:num1_seizure+num2_seizure+num3_seizure+num4_seizure+num5_seizure);


%%% Defining threshold of the Variance feature

threshold_var_baseline= mean(F_var_baseline);
threshold_var_seizure= mean(F_var_seizure);
threshold_var=(1*threshold_var_baseline + 1*threshold_var_seizure)/2;

%%% Definnig feature discriminative parameter (FDP)

FDP_var=sqrt(threshold_var_seizure/threshold_var_baseline);

%%% plotting results of Variance feature

figure
plot(t_window_tot,F_var,'linewidth',3)
title(' Variance feature extraction')
xlabel('Time (sec)')
ylabel('Variance Feature result')

hold on

plot(t1_seizure,F1_var_seizure,t2_seizure,F2_var_seizure,t3_seizure,F3_var_seizure,t4_seizure,F4_var_seizure,t5_seizure,F5_var_seizure,'linewidth',3,'color','r')

hold off

%%%% Variance detection

for i=1:num_window_tot
    if F_var(i)>=threshold_var
        det_var(i)=1;
    else
        det_var(i)=0;
    end
end

var_num_seizure=sum(det_var==1);
var_num_baseline=sum(det_var==0);

FP_var=sum(basicVector-det_var==-1);
FN_var=sum(basicVector-det_var==1);

TP_var=var_num_seizure-FP_var;
TN_var=var_num_baseline-FN_var;

sensitivity_var= TP_var/(TP_var+FN_var)
specificity_var=TN_var/(TN_var+FP_var)
accuracy_var=(sensitivity_var+specificity_var)/2



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Minimum feature extraction (feature 5) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:num_window_tot
    F_min(i)=min(y_final{i});
    
end

%%%% Determining baseline and seizure states of Minimum feature

F_min_baseline=[F_min(1,1:sz1_start-1) , F_min(1,sz1_end+1:sz2_start-1),F_min(1,sz2_end+1:sz3_start-1) , F_min(1,sz3_end+1: sz4_start-1), F_min(1,sz4_end+1:sz5_start-1) , F_min(1,sz5_end+1:num_window_tot)];

F_min_seizure=[F_min(1,sz1_start:sz1_end) , F_min(1,sz2_start:sz2_end) , F_min(1,sz3_start:sz3_end) , F_min(1,sz4_start:sz4_end),F_min(1,sz5_start:sz5_end)];

%%%%%%%%%%%%% Determining baseline and seizure states for Minimum feature

F1_min_seizure=F_min_seizure(1,1:num1_seizure);
F2_min_seizure=F_min_seizure(1,num1_seizure+1:num1_seizure+num2_seizure);
F3_min_seizure=F_min_seizure(1,num1_seizure+num2_seizure+1:num1_seizure+num2_seizure+num3_seizure);
F4_min_seizure=F_min_seizure(1,num1_seizure+num2_seizure+num3_seizure+1:num1_seizure+num2_seizure+num3_seizure+num4_seizure);
F5_min_seizure=F_min_seizure(1,num1_seizure+num2_seizure+num3_seizure+num4_seizure+1:num1_seizure+num2_seizure+num3_seizure+num4_seizure+num5_seizure);



%%% Defining threshold of the Minimum feature

threshold_min_baseline= mean(abs(F_min_baseline));
threshold_min_seizure= mean(abs(F_min_seizure));
threshold_min=(1*threshold_min_baseline + 3*threshold_min_seizure)/4;

%%% Defining feature discriminative parameter (FDP)

FDP_min=threshold_min_seizure/threshold_min_baseline;

%%% plotting results of Minimum feature

figure
plot(t_window_tot,F_min,'linewidth',3)
title(' Minimum feature extraction')
xlabel('Time (sec)')
ylabel('Minimum Feature result')

hold on

plot(t1_seizure,F1_min_seizure,t2_seizure,F2_min_seizure,t3_seizure,F3_min_seizure,t4_seizure,F4_min_seizure,t5_seizure,F5_min_seizure,'linewidth',3,'color','r')

hold off

%%%% Minimum detection

for i=1:num_window_tot
    if abs(F_min(i))>=threshold_min
        det_min(i)=1;
    else
        det_min(i)=0;
    end
end

min_num_seizure=sum(det_min==1);
min_num_baseline=sum(det_min==0);

FP_min=sum(basicVector-det_min==-1);
FN_min=sum(basicVector-det_min==1);

TP_min=min_num_seizure-FP_min;
TN_min=min_num_baseline-FN_min;

sensitivity_min= TP_min/(TP_min+FN_min)
specificity_min=TN_min/(TN_min+FP_min)
accuracy_min=(sensitivity_min+specificity_min)/2



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Maximum feature extraction (feature 6) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:num_window_tot
    
    F_max(i)=max(abs(y_final{i}));
    
end

%%%% Determining baseline and seizure states of Maximum feature

F_max_baseline=[F_max(1,1:sz1_start-1) , F_max(1,sz1_end+1:sz2_start-1),F_max(1,sz2_end+1:sz3_start-1) , F_max(1,sz3_end+1: sz4_start-1), F_max(1,sz4_end+1:sz5_start-1) , F_max(1,sz5_end+1:num_window_tot)];

F_max_seizure=[F_max(1,sz1_start:sz1_end) , F_max(1,sz2_start:sz2_end) , F_max(1,sz3_start:sz3_end) , F_max(1,sz4_start:sz4_end),F_max(1,sz5_start:sz5_end)];

%%%%%%%%%%%%% Determining baseline and seizure states for maximum feature

F1_max_seizure=F_max_seizure(1,1:num1_seizure);
F2_max_seizure=F_max_seizure(1,num1_seizure+1:num1_seizure+num2_seizure);
F3_max_seizure=F_max_seizure(1,num1_seizure+num2_seizure+1:num1_seizure+num2_seizure+num3_seizure);
F4_max_seizure=F_max_seizure(1,num1_seizure+num2_seizure+num3_seizure+1:num1_seizure+num2_seizure+num3_seizure+num4_seizure);
F5_max_seizure=F_max_seizure(1,num1_seizure+num2_seizure+num3_seizure+num4_seizure+1:num1_seizure+num2_seizure+num3_seizure+num4_seizure+num5_seizure);

%%% Defining threshold of the Maximum feature

threshold_max_baseline= mean(abs((F_max_baseline)));
threshold_max_seizure= mean(abs((F_max_seizure)));
threshold_max=(1*threshold_max_baseline + 3*threshold_max_seizure)/4;

%%% Defining feature discriminative parameter (FDP)

FDP_max=mean(F_max_seizure)/mean(F_max_baseline);

%%% plotting results of Maximum feature

figure
plot(t_window_tot,F_max,'linewidth',3)
title(' Maximum feature result')
xlabel('Time (sec)')
ylabel('Minimum Feature magnitude')

hold on

plot(t1_seizure,F1_max_seizure,t2_seizure,F2_max_seizure,t3_seizure,F3_max_seizure,t4_seizure,F4_max_seizure,t5_seizure,F5_max_seizure,'linewidth',3,'color','r')

hold off

%%%% Max detection

for i=1:num_window_tot
    if abs(F_max(i))>=threshold_max
        det_max(i)=1;
    else
        det_max(i)=0;
    end
end

max_num_seizure=sum(det_max==1);
max_num_baseline=sum(det_max==0);

FP_max=sum(basicVector-det_max==-1);
FN_max=sum(basicVector-det_max==1);

TP_max=max_num_seizure-FP_max;
TN_max=max_num_baseline-FN_max;

sensitivity_max= TP_max/(TP_max+FN_max)
specificity_max=TN_max/(TN_max+FP_max)
accuracy_max=(sensitivity_max+specificity_max)/2



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Range feature extraction (feature 7) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:num_window_tot
    
    F_range(i)=F_max(i)-F_min(i);
    
end

%%%% Determining baseline and seizure states of Range feature

F_range_baseline=[F_range(1,1:sz1_start-1) , F_range(1,sz1_end+1:sz2_start-1),F_range(1,sz2_end+1:sz3_start-1) , F_range(1,sz3_end+1: sz4_start-1), F_range(1,sz4_end+1:sz5_start-1) , F_range(1,sz5_end+1:num_window_tot)];

F_range_seizure=[F_range(1,sz1_start:sz1_end) , F_range(1,sz2_start:sz2_end) , F_range(1,sz3_start:sz3_end) , F_range(1,sz4_start:sz4_end),F_range(1,sz5_start:sz5_end)];

%%%%%%%%%%%%% Determining baseline and seizure states for the Range feature

F1_range_seizure=F_range_seizure(1,1:num1_seizure);
F2_range_seizure=F_range_seizure(1,num1_seizure+1:num1_seizure+num2_seizure);
F3_range_seizure=F_range_seizure(1,num1_seizure+num2_seizure+1:num1_seizure+num2_seizure+num3_seizure);
F4_range_seizure=F_range_seizure(1,num1_seizure+num2_seizure+num3_seizure+1:num1_seizure+num2_seizure+num3_seizure+num4_seizure);
F5_range_seizure=F_range_seizure(1,num1_seizure+num2_seizure+num3_seizure+num4_seizure+1:num1_seizure+num2_seizure+num3_seizure+num4_seizure+num5_seizure);



%%% Defining threshold of the Range feature

threshold_range_baseline= mean(abs(F_range_baseline));
threshold_range_seizure= mean(abs(F_range_seizure));
threshold_range=(1*threshold_range_baseline + 3*threshold_range_seizure)/4;

%%% Definnig feature discriminative parameter (FDP)

FDP_range=threshold_range_seizure/threshold_range_baseline;

%%% plotting results of Range feature

figure
plot(t_window_tot,F_range,'linewidth',3)
title(' Range feature result')
xlabel('Time (sec)')
ylabel('Range magnitude')

hold on

plot(t1_seizure,F1_range_seizure,t2_seizure,F2_range_seizure,t3_seizure,F3_range_seizure,t4_seizure,F4_range_seizure,t5_seizure,F5_range_seizure,'linewidth',3,'color','r')

hold off

%%%% range detection

for i=1:num_window_tot
    if abs(F_range(i))>=threshold_range
        det_range(i)=1;
    else
        det_range(i)=0;
    end
end

range_num_seizure=sum(det_range==1);
range_num_baseline=sum(det_range==0);

FP_range=sum(basicVector-det_range==-1);
FN_range=sum(basicVector-det_range==1);

TP_range=range_num_seizure-FP_range;
TN_range=range_num_baseline-FN_range;

sensitivity_range= TP_range/(TP_range+FN_range)
specificity_range=TN_range/(TN_range+FP_range)
accuracy_range=(sensitivity_range+specificity_range)/2


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Average feature extraction (feature 8) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:num_window_tot
    
    F_mean(i)=mean(y_final{i});
    
end



%%%% Determining baseline and seizure states of Average feature

F_mean_baseline=[F_mean(1,1:sz1_start-1) , F_mean(1,sz1_end+1:sz2_start-1),F_mean(1,sz2_end+1:sz3_start-1) , F_mean(1,sz3_end+1: sz4_start-1), F_mean(1,sz4_end+1:sz5_start-1) , F_mean(1,sz5_end+1:num_window_tot)];

F_mean_seizure=[F_mean(1,sz1_start:sz1_end) , F_mean(1,sz2_start:sz2_end) , F_mean(1,sz3_start:sz3_end) , F_mean(1,sz4_start:sz4_end),F_mean(1,sz5_start:sz5_end)];

%%%%%%%%%%%%% Determining baseline and seizure states for the mean feature

F1_mean_seizure=F_mean_seizure(1,1:num1_seizure);
F2_mean_seizure=F_mean_seizure(1,num1_seizure+1:num1_seizure+num2_seizure);
F3_mean_seizure=F_mean_seizure(1,num1_seizure+num2_seizure+1:num1_seizure+num2_seizure+num3_seizure);
F4_mean_seizure=F_mean_seizure(1,num1_seizure+num2_seizure+num3_seizure+1:num1_seizure+num2_seizure+num3_seizure+num4_seizure);
F5_mean_seizure=F_mean_seizure(1,num1_seizure+num2_seizure+num3_seizure+num4_seizure+1:num1_seizure+num2_seizure+num3_seizure+num4_seizure+num5_seizure);



%%% Defining threshold of the Average feature

threshold_mean_baseline= mean(abs(F_mean_baseline));
threshold_mean_seizure= mean(abs(F_mean_seizure));
threshold_mean=(1*threshold_mean_baseline + 3*threshold_mean_seizure)/4;

%%% Definnig feature discriminative parameter (FDP)

FDP_mean=threshold_mean_seizure/threshold_mean_baseline;

%%% plotting results of Average feature

figure
plot(t_window_tot,F_mean,'linewidth',3)
title(' Avg feature result')
xlabel('Time (sec)')
ylabel('Avg magnitude')

hold on

plot(t1_seizure,F1_mean_seizure,t2_seizure,F2_mean_seizure,t3_seizure,F3_mean_seizure,t4_seizure,F4_mean_seizure,t5_seizure,F5_mean_seizure,'linewidth',3,'color','r')

hold off

for i=1:num_window_tot
    if abs(F_mean(i))>=abs(threshold_mean)
        det_mean(i)=1;
    else
        det_mean(i)=0;
    end
end

mean_num_seizure=sum(det_mean==1);
mean_num_baseline=sum(det_mean==0);

FP_mean=sum(basicVector-det_mean==-1);
FN_mean=sum(basicVector-det_mean==1);

TP_mean=mean_num_seizure-FP_mean;
TN_mean=mean_num_baseline-FN_mean;

sensitivity_mean= TP_mean/(TP_mean+FN_mean)
specificity_mean=TN_mean/(TN_mean+FP_mean)
accuracy_mean=(sensitivity_mean+specificity_mean)/2


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Correlation feature extraction (feature 9) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for  i=3:num_window_tot
    
    HV_cor(i)= min([F_max(i),F_max(i-1),F_max(i-2)]);
    HV_cor(1)=F_max(1);
    HV_cor(2)=F_max(2);
    
    LV_cor(i)=max([F_min(i) , F_min(i-1) , F_min(i-2)]);
    LV_cor(1)=min(F_min(1)); 
    LV_cor(2)=min(F_min(2));
    
    F_cor(i)=HV_cor(i)-LV_cor(i);
    
end


%%%% Determining baseline and seizure states of Correlation feature

F_cor_baseline=[F_cor(1,1:sz1_start-1) , F_cor(1,sz1_end+1:sz2_start-1),F_cor(1,sz2_end+1:sz3_start-1) , F_cor(1,sz3_end+1: sz4_start-1), F_cor(1,sz4_end+1:sz5_start-1) , F_cor(1,sz5_end+1:num_window_tot)];

F_cor_seizure=[F_cor(1,sz1_start:sz1_end) , F_cor(1,sz2_start:sz2_end) , F_cor(1,sz3_start:sz3_end) , F_cor(1,sz4_start:sz4_end),F_cor(1,sz5_start:sz5_end)];

%%%%%%%%%%%%% Determining baseline and seizure states for the cor feature

F1_cor_seizure=F_cor_seizure(1,1:num1_seizure);
F2_cor_seizure=F_cor_seizure(1,num1_seizure+1:num1_seizure+num2_seizure);
F3_cor_seizure=F_cor_seizure(1,num1_seizure+num2_seizure+1:num1_seizure+num2_seizure+num3_seizure);
F4_cor_seizure=F_cor_seizure(1,num1_seizure+num2_seizure+num3_seizure+1:num1_seizure+num2_seizure+num3_seizure+num4_seizure);
F5_cor_seizure=F_cor_seizure(1,num1_seizure+num2_seizure+num3_seizure+num4_seizure+1:num1_seizure+num2_seizure+num3_seizure+num4_seizure+num5_seizure);




%%% Defining threshold of the Correlation feature

threshold_cor_baseline= mean(abs(F_cor_baseline));
threshold_cor_seizure= mean(abs(F_cor_seizure));
threshold_cor=(1*threshold_cor_baseline + 1*threshold_cor_seizure)/2;

%%% Definnig feature discriminative parameter (FDP)

FDP_cor=threshold_cor_seizure/threshold_cor_baseline;

%%% plotting results of Correlation feature

figure
plot(t_window_tot,F_cor,'linewidth',3)
title(' Correlation feature result')
xlabel('Time (sec)')
ylabel('Correlation magnitude')

hold on

plot(t1_seizure,F1_cor_seizure,t2_seizure,F2_cor_seizure,t3_seizure,F3_cor_seizure,t4_seizure,F4_cor_seizure,t5_seizure,F5_cor_seizure,'linewidth',3,'color','r')

hold off

%%% Correlation detection

for i=1:num_window_tot
    if abs(F_cor(i))>=threshold_cor
        det_cor(i)=1;
    else
        det_cor(i)=0;
    end
end

cor_num_seizure=sum(det_cor==1);
cor_num_baseline=sum(det_cor==0);

FP_cor=sum(basicVector-det_cor==-1);
FN_cor=sum(basicVector-det_cor==1);

TP_cor=cor_num_seizure-FP_cor;
TN_cor=cor_num_baseline-FN_cor;

sensitivity_cor= TP_cor/(TP_cor+FN_cor)
specificity_cor=TN_cor/(TN_cor+FP_cor)
accuracy_cor=(sensitivity_cor+specificity_cor)/2


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Feature ranking %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%% Ranking based on FDP

FDP_tot=[FDP_coastline, sqrt(FDP_energy), sqrt(FDP_NE), sqrt(FDP_var), FDP_min, FDP_max, FDP_range, FDP_mean, FDP_cor];
[FDP_tot_sort,FDP_tot_sort_index]=sort(FDP_tot,'descend')



%%%%% Selecting top 4 rank FDP features

 rank_FDP_tot=FDP_tot_sort_index(1:4)
 
 
 
 
 
 
 
 
 