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