Code covered by the BSD License

# Fast permutation entropy

### Valentina Unakafova (view profile)

Efficient computing of empirical permutation entropy

PE(x, Tau, d, WS)
% efficienct computing of empirical permutation entropy

% Ref: Unakafova, V. A., & Keller, K. (2013).
% Efficiently Measuring Complexity on the Basis of Real-World Data.
% Entropy, 15(10), 4392-4415.

% x - time series, Tau - delay $\tau$, d - order of ordinal patterns $d$,
% WS - size $M$ of a sliding window
% ePE -  values of the empirical permutation entropy
function ePE = PE(x, Tau, d, WS)
load(['table' num2str(d) '.mat']);% the precomputed table
pTbl = eval(['table' num2str(d)]);
Length = numel(x);                % length of the time series
d1 = d+1;
dTau = d*Tau;
nPat = factorial(d1);             % amount of ordinal patterns of order d
opd = zeros(1, nPat);             % distribution of ordinal patterns
ePE = zeros(1, Length);           % empirical permutation entropy
op = zeros(1, d);                 % ordinal pattern $(i_1,i_2,...,i_d)$
prevOP = zeros(1, Tau);           % previous ordinal patterns for $1:\tau$
opW = zeros(1, WS);               % ordinal patterns in the window
ancNum = nPat./factorial(2:d1);   % ancillary numbers
peTbl(1:WS) = -(1:WS).*log(1:WS); % table of values $g(j)$
peTbl(2:WS) = (peTbl(2:WS)-peTbl(1:WS-1))./WS;
for iTau = 1:Tau
cnt = iTau;
op(1) = (x(dTau+iTau-Tau) >= x(dTau+iTau));
for j = 2:d
op(j) = sum(x((d-j)*Tau+iTau) >= x((d1-j)*Tau+iTau:Tau:dTau+iTau));
end
opW(cnt) = sum(op.*ancNum);              % the first ordinal pattern
opd(opW(cnt)+1) = opd(opW(cnt)+1)+1;
for j = dTau+Tau+iTau:Tau:WS+dTau        % loop for the first window
cnt = cnt+Tau;
posL = 1;                 % the position $l$ of the next point
for i = j-dTau:Tau:j-Tau
if(x(i) >= x(j))
posL = posL+1;
end
end
opW(cnt) = pTbl(opW(cnt-Tau)*d1+posL);
opd(opW(cnt)+1) = opd(opW(cnt)+1)+1;
end
prevOP(iTau) = opW(cnt);
end
ordDistNorm = opd/WS;
ePE(WS+Tau*d) = -nansum(ordDistNorm(1:nPat).*log(ordDistNorm(1:nPat)));

iTau = 1;                   % current shift $1:\tau$
iPat = 1;                   % position of the current pattern in the window
for t = WS+Tau*d+1:Length   % loop over all points
posL = 1;                 % the position $l$ of the next point
for j = t-dTau:Tau:t-Tau
if(x(j) >= x(t))
posL = posL+1;
end
end
nNew = pTbl(prevOP(iTau)*d1+posL); % "incoming" ordinal pattern
nOut = opW(iPat);                  % "outcoming" ordinal pattern
prevOP(iTau) = nNew;
opW(iPat) = nNew;
nNew = nNew+1;
nOut = nOut+1;
if nNew ~= nOut            % update the distribution of ordinal patterns
opd(nNew) = opd(nNew)+1; % "incoming" ordinal pattern
opd(nOut) = opd(nOut)-1; % "outcoming" ordinal pattern
ePE(t) = ePE(t-1)+(peTbl(opd(nNew))-peTbl(opd(nOut)+1));
else
ePE(t) = ePE(t-1);
end
iTau = iTau+1;
iPat = iPat+1;
if(iTau > Tau) iTau = 1; end
if(iPat > WS) iPat = 1; end
end
ePE = ePE(WS+Tau*d:end);