1 line
2.4 KiB
Matlab
1 line
2.4 KiB
Matlab
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);
|
|
|
|
|