clc;
clearvars;
close all;

N=10;
TS = 0.01;
NP = 1;
NZ = 0;
IODelay = NaN;
%%

DelayVect = zeros(1,N);

% 1/p^2
%     stringPath = 'Merane_data\BP/Tuhost/PRBS/PRBS_T_1-p^2_N';
     stringPath = 'Merane_data\BP\Tuhost\Sinus\SIN_T_1-p^2_N';
%     stringPath = 'Merane_data\BP\BezTuhosti\PRBS\PRBS_BT_1-p^2_N';
%     stringPath = 'Merane_data\BP\BezTuhosti\sinus\SIN_BT_1-p^2_N';

% K
    %stringPath = 'Merane_data\NaZpracovanie\K\K10__prbs2_';
% K/p
    %stringPath = 'Merane_data\NaZpracovanie\K--p\K5--p__prbs3_';
    %stringPath = 'Merane_data\NaZpracovanie\sinus\K5--p__amp0.2_freq1_2_3rad__N';

% K/p+1
    %stringPath = 'Merane_data\NaZpracovanie\K--p+1\K10--p+1__prbs3_';

for i=1:N
% 1/p^2
    %load(['Merane_data\BP\BezTuhosti\PRBS\PRBS_BT_1-p^2_N',num2str(i),'.mat']);
    %load(['Merane_data\NaZpracovanie\K--p^2\K1--p^2__prbs10_',num2str(i),'.mat']);
    load([stringPath, num2str(i),'.mat']);
    %load(['Merane_data\BP\BezTuhosti\sinus\SIN_BT_1-p^2_N',num2str(i),'.mat']);



    data = iddata(out.action.data(1:end-1),diff(out.error.data)/TS, TS);
    del = 0.1:0.05:0.8;
    fit = zeros(1,length(del));
    fInd  = 1;
    for d = del
        
        sys = tfest(data, NP, NZ, d);
        fit(fInd) = sys.Report.Fit.FitPercent;
        fInd = fInd+1;

    end
    figure(i)
    plot(del,fit);
    bestDelay = del(fit == max(fit));
    DelayVect(i) = bestDelay(1);
end