Нужно посчитать корреляцию функции

T=5*3600; % Продолжительность реализации, с
dt=1; % Период дискретности, с
sd=19; % CKO дрейфа гироскопа, угл.с.
Td=136; % Период корреляции дрейфа гироскопа, c
sw=1
w(1)=1.744;
w(2)=0.348;
w(3)=-1.213;
w(4)=0.133;
w(5)=0.382;
N   = round(T/dt);                  % Число моментов времени в реализации
x(N)= 0;                            % Резервирование памяти
x(1)=1;
% F   = 1-dt/Td;
% G   = sd*sqrt(2*dt/Td);
F   = exp(-dt/Td);
G   = sd*sqrt(1-F*F);
x(2)=F*x(1) + G*w(2);
x(3)=F*x(2) + G*w(3);
x(4)=F*x(3) + G*w(4);
x(5)=F*x(4) + G*w(5);
% Построение КФ истинной
Tau = 4*Td;                         % Макс. значение аргумента графика
Nkf = 500;                          % Число точек на графике КФ
tau = linspace(0,Tau,Nkf);
kfa = sd*sd*exp(-tau/Td);
% Построение КФ по точкам, введённым пользователем
tap = Td*[0 0.5 1.0 1.5 2.0 3.0]';
% Построение КФ по цветному шуму
Ndt = fix(Tau/Nkf/dt);          % Число тактов ЦШ для одной точки графика
if Ndt == 0 Ndt = 1; end
v   = x;                        % Цветной шум для сдвига при построении КФ
for i = 1:fix(Tau/(Ndt*dt))
    l = length(x);
    kfx(i)  = x*v'/(l-1) - mean(x)*mean(v);
    tax(i)  = (i-1)*Ndt*dt;
    x(l-Ndt+1:l) = [];
    v(1:Ndt)= [];
end
% Проверка введённых пользователем значений КФ
for i = 1:length(tap)
    [val, Ind(i)] = min(abs((tap(i)-tax)));
end
kfc = kfx(Ind)';                % Правильные значения КФ
dkf = abs((kfp-kfc)./kfc);      % Относительные погрешности пользователя
colr= [1 0.75 0.75];            % Цвет неправильного ответа