Code: Select all
%% Sample code to apply a windowing recurrence analysis with CRP Toolbox for MATLAB
% Allows to change the embedding and recurrence parameters in the
% running windows.
% CRP Toolbox for MATLAB required: http://tocsy.pik-potsdam.de/CRPtoolbox/
%% Simple system with changing dynamics (AR(1) process)
% AR-coefficient
a = [.9 * ones(1000,1); .97 * ones(300,1); .8 * ones(1000,1)];
% length of time series
N = length(a);
% calculate AR(1) time series
x = zeros(N,1);
x(1) = randn;
for i = 2:N
x(i) = a(i) * x(i-1) + randn;
end
% plot time series
plot(x)
%% (1) windowing RQA - simple approach
w = 200; % window length
ws = 10; % moving step
cnt = 1; % counter (iterates with each window)
clear Y % in this variable the RQA results will be stored
for i = 1:ws:N-w % run through each window
xs = x(i:i+w-1); % subsequence of the time series (= window)
Y(cnt,:) = crqa(xs,1,1,.1,'euc','sil'); % perform RQA for the current window
cnt = cnt + 1; % iterate counter
end
%%
% plot results
% (note that here the time point of the window is the starting of the window)
plot(1:ws:N-w,Y(:,2))
grid on
%% (2) windowing RQA - preset windows in separate variable
w = 200; % window length
ws = 50; % moving step
windows = 1+floor(w/2):ws:N-floor(w/2); % mid points of the windows
%%
% in this variable the RQA results will be stored;
% the preallocation with zeros can speed up MATLAB
Y = zeros(length(windows),13);
cnt = 1;
for i = windows
xs = x(i-floor(w/2):i+floor(w/2)); % subsequence of the time series (= window)
Y(cnt,:) = crqa(xs,3,1,.2,'euc','sil'); % perform RQA for the current window
cnt = cnt + 1;
end
%%
% plot results
% (note that here the time point of the window is the center of the window)
plot(windows,Y(:,2))
grid on
%% (3) windowing RQA with changing parameters
% in this example, the embedding dimension is found by FNN for
% each window separately
w = 200;
ws = 50;
windows = 1+floor(w/2):ws:N-floor(w/2);
Y = zeros(length(windows),13);
cnt = 1;
for i = windows
xs = x(i-floor(w/2):i+floor(w/2)); % subsequence of the time series (= window)
mAll = fnn(xs,'sil'); % get the false nearest neighbours as a function of the embedding dimension
m(cnt) = find(mAll > 0, 1, 'last') + 1; % find the embedding dimension (where FNN vanishes)
Y(cnt,:) = crqa(xs,m(cnt),1,.2*sqrt(m),'rr','sil'); % apply RQA with changing embedding dimension
cnt = cnt + 1;
end
plot(windows,Y(:,2))
grid on