1 line
2.4 KiB
Mathematica
1 line
2.4 KiB
Mathematica
|
function divergence=div_calc_shorttimeseries(state,ws_sec,fs,period_sec,progress)
% calculate local dynamic stability (max lyapunov exponent).
% Input:
% state: appropriate state space
% ws: window size over which divergence should be calculated(in seconds)
% fs: sample frequency
% period: ..... (can be dominant period in the signal(in seconds))
% plotje: show a graph.
%
% Note that in this version ws should be larger than 4*period, as the long
% term divergence is calculated from 4*period to ws*fs.
%
% Output:
% divergence: the divergence curve
% lds: the 2 estimates of the maximum lyapunov exponents (short term and long term)
%
% Earlier versions by KvS en SMB. Made the routine faster, 14/01/2011, Sietse Rispens
% Use less memory to prevent memory error, 01/02/2011, Sietse Rispens
% Take mean of the log of divergence instead of log of the mean of divergence, 26/05/2011, Sietse Rispens
% Allow non-integer ws_sec and period_sec and change fitting-range, 01/08/2011, Sietse Rispens
% Do not try to find neighbours that need to be followed beyond end of time series, 01/11/2012, Sietse Rispens
% Exclude neighbours closer than Period in time (instead of period/2), 01/11/2012, Sietse Rispens
[m,n]=size(state);
ws=round(ws_sec*fs);
period=round(period_sec*fs);
mcompletewindow = m-ws+1;
statecw = state(1:mcompletewindow,:);
divergence_sum=zeros(1,ws);
divergence_count=zeros(1,ws);
diff_state = statecw*0;
diff_state_sqr = diff_state;
diff_total = zeros(size(diff_state,1),1);
for i = 1:mcompletewindow
if ~isnan(state(i,:))
start=round(max([1,i-period+1]));
stop=round(min([mcompletewindow,i+period-1]));
for j=1:n,
diff_state(:,j) = statecw(:,j)-statecw(i,j);
diff_state_sqr(:,j)=diff_state(:,j).^2;
end
diff_total(:,1)=sum(diff_state_sqr,2);
diff_total(start:stop,1)=NaN;
[mini,index]=min(diff_total);
if i+ws>m || index+ws>m
% do not use these data
else
divergence_sum = divergence_sum + log(sqrt(sum((state(i:i+ws-1,:)-state(index:index+ws-1,:)).^2,2)))';
divergence_count = divergence_count + 1;
end
end
if progress > 0
if mod(i,round(progress))==0
if i>round(progress)
fprintf('\b\b\b\b\b\b\b\b\b\b');
end
fprintf('i=%8d', i);
end
end
end
if progress > 0
fprintf('\n');
end
divergence= (divergence_sum./divergence_count);
|