Sample code to apply a windowing recurrence analysis with CR

Everything about quantification of recurrence plots.

Sample code to apply a windowing recurrence analysis with CR

Postby Norbert » Tue Jul 7, 2015 21:27

I post here a MATLAB coding example that was developed within the hands-on session of the 6th RP symposium in Grenoble. It allows a windowed recurrence analysis with changing parameters for each window.

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
User avatar
Norbert
Expert
 
Posts: 173
Joined: Wed Jan 4, 2006 11:03
Location: Potsdam Institute for Climate Impact Research, Germany

Return to Quantification

Who is online

Users browsing this forum: No registered users and 1 guest

cron