--- a
+++ b/Supporting Functions/lyapunov.m
@@ -0,0 +1,52 @@
+function lam = lyapunov(x,Fs)
+% calculate lyapunov coefficient of time series
+
+%x = sin(1:dt:10)'; % this is 68.83
+%x = (1:dt:10)';    % this is 0
+%x = randn(1001,1); % this is ~300
+%x = squeeze(EEG.data(1,1:1001,1))'; %this is 241.
+dt = 1/Fs;
+[ndata nvars]=size(x);
+steps_forward = 5;
+N2 = floor(ndata/2);
+N4 = floor(ndata/4);
+TOL = 1.0e-6;
+
+exponent = zeros(N4+1,1);
+tic
+for i=N4:N2  % second quartile of data should be sufficiently evolved
+   
+   %get all points but this one.
+   js = setdiff(1:ndata-steps_forward,i);
+   %find the index of the nearest neighbor
+   [idx] = knnsearch(x(js,:),x(i,:));
+   indx = js(idx);
+   
+%   tic
+%    dist = norm(x(i+1,:)-x(i,:));
+%    indx = i+1;
+%    for j=i:ndata-5
+%        if (i ~= j) && norm(x(i,:)-x(j,:))<dist
+%            dist = norm(x(i,:)-x(j,:));
+%            indx = j; % closest point!
+%        end
+%    end
+%    toc
+   expn = 0.0; % estimate local rate of expansion (i.e. largest eigenvalue)
+   for k=1:steps_forward
+       if norm(x(i+k,:)-x(indx+k,:))>TOL && norm(x(i,:)-x(indx,:))>TOL
+           expn = expn + (log(norm(x(i+k,:)-x(indx+k,:)))-log(norm(x(i,:)-x(indx,:))))/k;
+       end
+   end
+   exponent(i-N4+1)=expn/steps_forward;   
+end
+toc
+
+sum=0;  % now, calculate the overal average over N4 data points ...
+for i=1:N4+1
+    sum = sum+exponent(i);
+end
+
+lam=sum/((N4+1)*dt);  
+% return the average value
+% if lam > 0, then system is chaotic