Larges Lyapunov Exponent from time series

dist.m

function D = dist(X)
	s = length(X(1,:));
	n = length(X(:,1));
	D = zeros(s,s);
	for i=1:s
		for j=1:i
			if i == j
				D(i,j) = 0;
			else
				D(i,j) = sqrt(sum(power(X(:,i)-X(:,j),2)));
			endif
		endfor
	endfor
	D = D + D' + eye(s,s)*1e6;
endfunction

findLyapunov.m

%
% findLyapunov()
% funkcja znajduje najwiekszy wykladnik Lapunowa dla trajektorii zanurzonej w 3 wymiarach
% x, y, z - wspolrzedne punktow odleglych o jednakowy odstep czasu
% autor: Mateusz Galazyn
% oparte na algorytmie z: 
% WOLF, SWIFT, SWINNEY, AND VASTANO "Determining lyapunov exponents from a time series."
% Physica D: Nonlinear Phenomena 16, 3 (July 1985), 285–317.
%

function L = findLyapunov(x,y,z)
	s = length(x);
	D = dist([x' ; y' ; z' ]);
	L = 0;
	counter = 0; % licznik punktow dla ktorych zostala obliczona ekspansja
	for i=1:floor(s-1)
		[v, iv] = min(D(i,:)); % najblizszy sasiad
		if iv < s
			distancePrime = dist( [x(i+1) x(iv+1) ; y(i+1) y(iv+1) ; z(i+1) z(iv+1) ] )(1,2); % L prime
			distance = dist( [x(i) x(iv) ; y(i) y(iv) ; z(i) z(iv) ] )(1,2); % L
			if distance > 0 && distancePrime > 0
				++counter;
				L += log2( distancePrime / distance	);
			endif
		endif
	endfor

	if counter == 0
		L = 0; % jezeli orbita okresowa, zaden punkt nie zostal uzyty do obliczen
	else
		L /= counter;
	endif
endfunction