function out = ridici_v3(filename)
%% nacteni dat ze souboru
% vytvorit cell a indexovat, aby probehlo pro vsechny
Data = readmatrix (filename);    
 
% Data = readmatrix (filename);
Time = Data(1000:end-50,1)* 1e-3;
SplineDistance = Data(1000:end-50,4);
WheelAngle=Data(1000:end-50,6)*50;

%% 
% ulozi se data do matice, ale bez textoveho nadpisu, s absolutní hodnotou

data = abs(SplineDistance);
data_wheel = abs(WheelAngle);
data_wheel_filt  = filtfilt(ones(250,1)/250,1,data_wheel);
Ts=mean(Time(3:end)-Time(2:end-1));
f_s=1/Ts;
thr = 2.5;
thrW = 0.4;
%thr = quantile(SplineDistance,0.75); 3.quartil pro prahovou hodnotu

% thr = 0.8 * max(SplineDistance);

peaks = data > thr;
peaks_w = data_wheel_filt > thrW;

%zacatek a konec obdelniku, kladna a zaporna je konec obdelniku
idx_peaks = diff([peaks(1);peaks]);
idx_peaks_w = diff([peaks_w(1);peaks_w]);
% vypíše indexy špiček odchylky
idx_start = find(idx_peaks > 0.5);
idx_stop =  find(idx_peaks < -0.5);
% vypíše indexy špiček úhlu volantu
idx_startWheel = find(idx_peaks_w > 0.5);
idx_stopWheel =  find(idx_peaks_w < -0.5);

% vyšetření špiček pro odchylku, opravuje a smaže vadné prvky
if length(idx_start) > length(idx_stop)
    idx_start(end) = [];
elseif    length(idx_start) < length(idx_stop)
    idx_stop(1) = [];
end

 if idx_start(1) > idx_stop(1)
     idx_start(end) = [];
     idx_stop(1)=[];
 end

% vyšetření špiček pro úhly volantu
if length(idx_startWheel) > length(idx_stopWheel)
    idx_startWheel(end) = [];
elseif    length(idx_startWheel) < length(idx_stopWheel)
    idx_stopWheel(1) = [];
end

if idx_startWheel(1) > idx_stopWheel(1)
    idx_startWheel(end) = [];
    idx_stopWheel(1)=[];
end

% naleznutí maxima špiček
idx_peaks_time = idx_start;
idx_peakW_time =idx_startWheel;

ss_peaks = [idx_start idx_stop];
ss_peaksW =[idx_startWheel idx_stopWheel];
idx_peaksWheel = zeros(size(idx_startWheel));

idx_peaks = zeros(size(idx_start));

if length(idx_start) ~= length(idx_peaksWheel)
    
err_size = zeros(length(idx_peaks),length(idx_peaksWheel));
    if length(idx_peaks) > length(idx_peaksWheel)
        for k = 1:length(idx_peaks)
            
            err_size(k,:) = (idx_startWheel-idx_start(k)).^2;
            
        end
        
        err_size = min(err_size,[],2);
        [~,err_size] = maxk(err_size,-length(idx_peaksWheel)+length(idx_peaks)) ;
       idx_peaks_time(err_size) = []; 
       
       ss_peaks(err_size,:) =[];
        idx_peaks(err_size) =[];
    else 
       for k = 1:length(idx_peaksWheel)
            
            err_size(:,k) = ((idx_startWheel(k)-idx_start).^2)';
            
        end 
        err_size = min(err_size,[],1);
        [~,err_size] = maxk(err_size,-length(idx_peaks)+length(idx_peaksWheel)) ;
       idx_peakW_time(err_size) = []; 
       
       
       ss_peaksW(err_size,:) =[];
        
    end
   
    
    
end

for k = 1:length(ss_peaks)
    [~, idxmW] = max(data(ss_peaks(k,1):ss_peaks(k,2)));
    idx_peaks_time(k) = ss_peaks(k,1)+idxmW-1;
    
end

%naleznuti referencnih indexu pro troj. detektor
idx = zeros(length(ss_peaks),2); 
N = 1000; % podle periody vzorkovani Ts
for k = 1:length(ss_peaks)  
    
%     start
if ss_peaks(k,1) - N < 1
    
    idx(k,1) = 1;
else
    %[~,idx_k] = mink(data(idx_start(k)-N:idx_start(k)),1);
    idx(k,1) = ss_peaks(k,1) - N;
    
end
    
    
    
%     stop
   if ss_peaks(k,2) + N > length(data)
    
    idx(k,2) = length(data);
else
    %[~,idx_k] = mink(data(idx_stop(k):idx_stop(k)+N),1);
    idx(k,2) = N + ss_peaks(k,2);
   end
   
end

%% start/stop celeho peaku
% do detektoru se posle segment dat, radek 70
[start, stop] = deal(zeros(size(ss_peaks(:,1))));

for k = 1:length(ss_peaks)
    seg = data(idx(k,1):idx_peaks_time(k));
    [start(k),~] = triangle_detector(seg);
    start(k) = start(k) + idx(k,1) - 1;
    
    seg = data(idx_peaks_time(k):idx(k,2));
    [stop(k),~] = triangle_detector(seg);
    stop(k) = stop(k) + idx_peaks_time(k) - 1;
end
%% reakční čas
idx_timeW=zeros(1,length(idx_peakW_time));
for k=1:length(idx_peakW_time)
    
    plus_peak = SplineDistance(idx_peaks_time(k)) > 0;
    
    if plus_peak
        [myMax,idx_max] = max(WheelAngle(ss_peaksW(k,1):ss_peaksW(k,2)));
        
        myMax = find(WheelAngle(ss_peaksW(k,1):ss_peaksW(k,2)+idx_max)>0.25*myMax,1,'first');   
    
        idx_timeW(k) = ss_peaksW(k,1)+myMax-1;
        
    else
        
         [myMax,idx_max] = min(WheelAngle(ss_peaksW(k,1):ss_peaksW(k,2)));
        
        myMax = find(WheelAngle(ss_peaksW(k,1):ss_peaksW(k,2)+idx_max)<0.15*myMax,1,'first');   
    
        idx_timeW(k) = ss_peaksW(k,1)+myMax-1;
        
    end
    
end


reactTime = abs(Time(idx_timeW) - Time(ss_peaks(:,1)));

%% doba náběhu

 idx_time_rise=zeros(1,length(idx_peaks_time));
 for k=1:length(idx_peaks_time)
     [myReact,idx_react] = max(data(ss_peaks(k,1):ss_peaks(k,2)));
     second_peak = find(data(idx_react+ss_peaks(k,1):stop(k))>0.10*myReact,1,'last');
     idx_time_rise(k) = ss_peaks(k,1)+second_peak + idx_react-1;
 end
riseTime= Time(idx_time_rise) - Time(idx_peaks_time);

%% doba překmitu
prekmit_bod=zeros(length(stop),1);
%např. použití znova troj. detektor
prekmit_stop=stop;
prekmit_smazat = [];
for k = 1:length(ss_peaks)
    segment= data(stop(k):stop(k)+500);
    idx_min = detect_overVibrate(segment);
    if idx_min
        
        prekmit_bod(k) =idx_min+stop(k); 
    else
        
        prekmit_smazat = [prekmit_smazat k];
    end
end
prekmit_bod(prekmit_smazat)=[];
prekmit_stop(prekmit_smazat)=[];
OverShoot_time = Time(prekmit_bod) - Time(prekmit_stop);
% ještě přidat překmit v amplitudě

%% plocha peaku

S_1 = zeros(size(ss_peaks,1),1); %obsah
S_2 = S_1;
S_3 = S_1;
clear ident_tool
for k = 1:length(S_1)
    S_1(k) = sum(data(start(k):stop(k))) * Ts; %absolutni kriterium
    S_2(k) = sum(data(start(k):stop(k)))^(2) * Ts; % kvadraticke kriterium
    S_3(k) = sum(data(start(k):stop(k)) .* Time(start(k):stop(k)))  * Ts^(2); % ITAE kriterium   
%     ident_hodnoty= iddata(SplineDistance(start(k):stop(k)),WheelAngle(start(k):stop(k)),0.0099);
%      ident_tool = tfest(ident_hodnoty,6,6);
%      [n p k]=zpkdata(ident_tool,'v');
end
%% vytvoření struktur, statistika a histogramy

out.int_krit.abs = S_1;
out.int_krit.kvadr = S_2;
out.int_krit.ITAE = S_3;
out.reactTime = reactTime;
out.riseTime = riseTime;
out.OverShootTime=OverShoot_time;
out.stat.mean = [mean(S_1) mean(S_2) mean(S_3)];
out.stat.median = [median(S_1) median(S_2)  median(S_3)];
out.stat.quart = [quantile(S_1,[0.25 0.5 0.75]') quantile(S_2,[0.25 0.5 0.75]')...
quantile(S_3,[0.25 0.5 0.75]') ]; %ctvrtina polovina a 3/4 dat,  

out.stat.std = [std(S_1) std(S_2) std(S_3)];   %smerodatna odchylka
% předělat na jeden velký vektor nebo matici
out.stat.rt = [mean(reactTime) median(reactTime) std(reactTime) ...
                quantile(reactTime, [0.25 0.5 0.75])];

% pro strojove učení
% ident_hodnoty= iddata(SplineDistance(start(k):stop(k)),WheelAngle(start(k):stop(k)),0.0099);
% ident_tool = tfest(ident_hodnoty,2,0);
% T=sqrt(ident_tool.Denominator(1)/(ident_tool.Denominator(3)));
% ksi=ident_tool.denominator(2)/((2*T)*ident_tool.Denominator(3));
% Gain = ident_tool.Numerator/(ident_tool.Denominator(3));
% out.vektor = [Gain; T; ksi; reactTime; riseTime; OverShoot_time; S_1; S_2; S_3];
for k=1:length(S_1)
ident_hodnoty= iddata(SplineDistance(start(k):stop(k)),WheelAngle(start(k):stop(k)),0.0099);
ident_tool = tfest(ident_hodnoty,2,0);
T(k)=sqrt(ident_tool.Denominator(1)/(ident_tool.Denominator(3)));
ksi(k)=ident_tool.denominator(2)/((2*T(k))*ident_tool.Denominator(3));
Gain(k) = ident_tool.Numerator/(ident_tool.Denominator(3));
end
out.Gain = Gain;
out.T = T;
out.ksi = ksi;



end