Code covered by the BSD License

# Fast permutation entropy

### Valentina Unakafova (view profile)

Efficient computing of empirical permutation entropy

PEeq(x, Tau, d, WS)
% efficienct computing of empirical permutation entropy
% with modified ordinal patterns

% 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 $d$ of modified ordinal patterns,
% WS - size $M$ of a sliding window
% ePE - values of empirical permutation entropy
function ePE = PEeq(x, Tau, d, WS)
load(['tableEq' num2str(d) '.mat']); % the precomputed table
opTbl = eval(['tableEq' num2str(d)]);% of successive ordinal patterns
L = numel(x);                        % length of time series
dTau = d*Tau;
nPat = 1;
for i = 3:2:2*d+1
nPat = nPat*i;
end
opd = zeros(1, nPat);      % distribution of the modified ordinal patterns
ePE = zeros(1, L);         % empirical permutation entropy
b = zeros(Tau, d);         % indicator of equality $(b_1,b_2,\ldots,b_d)$
prevOP = zeros(1, Tau);    % previous modified ordinal patterns for $1:\tau$
opW = zeros(1, WS);        % modified ordinal patterns in the window
ancNum = ones(1, d);       % ancillary numbers
for j = 2:d
ancNum(j) = ancNum(j-1)*(2*j-1);
end
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           % all shifts
cnt = iTau;
mOP = zeros(1, d);
t = dTau+iTau;           % current time $t$ of the last point in mOP
for j = 1:d              % determining modified ordinal patterns
for i = j-1:-1:0
if(i == 0 || b(iTau, i) == 0)
if(x(t-j*Tau) > x(t-i*Tau))
mOP(j) = mOP(j)+2;
elseif(x(t-j*Tau) == x(t-i*Tau))
b(iTau, j) = 1;
end
end
end
end
mOP(1:d) = mOP(1:d)+b(iTau, 1:d);   % add equality indicator
opW(cnt) = sum(mOP.*ancNum);
opd(opW(cnt)+1) = opd(opW(cnt)+1)+1;
cnt = cnt+Tau;
for t = iTau+Tau*(d+1):Tau:WS+Tau*d % loop for the first window
b(iTau, 2:d) = b(iTau, 1:d-1);    % renew $(b_1,b_2,\ldots,b_d)$
b(iTau, 1) = 0;
posL = 1;                         % position $L$ of the next point
eqFlag = 0;                       % indicator of equality $B$
for i = 1:d;                      % determining the position $L$
if(b(iTau, i) == 0)
if(x(t-i*Tau) > x(t))
posL = posL+2;
elseif(x(t) == x(t-i*Tau))
eqFlag = 1;
b(iTau, i) = 1;
end
end
end
posL = posL+eqFlag;  % position $L$ of the next point
opW(cnt) = opTbl(opW(cnt-Tau)*(2*d+1)+posL);
opd(opW(cnt)+1) = opd(opW(cnt)+1)+1;
cnt = cnt+Tau;
end
prevOP(iTau) = opW(t-dTau);
end
OPDnorm = opd/WS;        % normalization of the ordinal distribution
ePE(WS+Tau*d) = -nansum(OPDnorm(1:nPat).*log(OPDnorm(1:nPat)));

iTau = 1;                % current shift $1:\tau$
iOP = 1;                 % position of the current pattern in the window
for t = WS+Tau*d+1:L     % loop for all points in a time series
b(iTau, 2:d) = b(iTau, 1:d-1);
b(iTau, 1) = 0;
posL = 1;
eqFlag = 0;            % x(j)==x(i)?
for i = 1:d;           % determining the position $L$
if(b(iTau, i) == 0)
if(x(t-i*Tau) > x(t))
posL = posL+2;
elseif(x(t) == x(t-i*Tau))
eqFlag = 1;
b(iTau, i) = 1;
end
end
end
posL = posL+eqFlag;                      % position $L$ of the next point
nNew = opTbl(prevOP(iTau)*(2*d+1)+posL); % "incoming" ordinal pattern
nOut = opW(iOP);                         % "outcoming" ordinal pattern
prevOP(iTau) = nNew;
opW(iOP) = nNew;
nNew = nNew+1;
nOut = nOut+1;
if nNew ~= nOut                     % if nNew == nOut, ePE does not change
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;
iOP = iOP+1;
if(iTau > Tau) iTau = 1; end
if( iOP > WS) iOP = 1; end
end
ePE = ePE(WS+Tau*d:end);