How to calculate Yaghoobi Fractal dimension

How to calculate the Yaghoobi fractal dimension for a high dimension trajectory in the high dimensional phase space.

Article title: A new approach to measure the fractal dimension of a trajectory in the high-dimensional phase space

Article DOI: https://doi.org/10.1016/j.chaos.2021.111239

[reza yaghoobi]

function main()
clc;clear all;close all
Kmax = 8;
Y = randn(55,4);
F = YaghoobiFractalDimension(Y,Kmax);
end
function YFD = YaghoobiFractalDimension(Y,Kmax)
%% Estimate Yaghoobi fractal dimension
% Y: a high dimension trajectory
% Kmax: The length of actions
% YFD: Fractal Dimension of the input signal
dim = size(Y,2);
N = size(Y,1);
for k = 1:Kmax
    clear xx L
    for i = 1:k
        a = Y(i:k:(i+fix((N-i)/k)*k),:);
        d = diff(a);
        r = sqrt(sum(d.^2,2));
        alpha = [N-1]/(k*floor([N-i]/k));
        L(i) = alpha*sum(r)/(k^(dim-1));
    end
    LM(k) = mean(L);
end
LM(LM==0)=1e-200;
lLM = log10(LM);
lk = log10(1:Kmax);
lk(LM==0) = [];
lLM(LM==0) =[];
%% LSE
[a,b] = LSE(lk,lLM);
%% Plot
plot(lk,a*lk+b,':k');hold on
plot(lk,lLM,'.','MarkerSize',15);
xlabel('Log_1_0(k)')
ylabel('Log_1_0(L_k)')
xlim([lk(1) lk(end)])
YFD = abs(a);
title(['FD = ',num2str(YFD)])
end
%% fitting
function [a,b] = LSE(x,y)
N = size(x,2);
a = [N*x*y'-sum(x)*sum(y)]/[N*sum(x.^2)-sum(x)^2];
b = [sum(x.^2)*sum(y)-x*y'*sum(x)]/[N*sum(x.^2)-sum(x)^2];
end