Code covered by the BSD License  

Highlights from
Fast permutation entropy

Fast permutation entropy

by

 

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

Contact us