Code covered by the BSD License  

Highlights from
Fast permutation entropy

Fast permutation entropy

by

 

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);

Contact us