Recurrence Plots Using Matlab

This post has been taken from : http://www.recurrence-plot.tk/rp-tutorial.php and is only meant for the benefit of students. In case you find the topic helpful, kindly acknowledge the above source and not me.

Matlab Tutorial for Creating Recurrence Plots

For systems with dimensions higher than three, it is obviously difficult to look at their phase space trajectory. However, there is a tool which enables us to study the phase space portrait also for higher dimensional systems: The recurrence plot a visualization tool, which visualizes the recurrences of states in the phase space by a 2-dimensional plot:

R_{i, j} = Θ ( ε_i − || x_i − x_j||), x_i ∈ ℜ^m, i, j=1…N,

A recurrence is then defined when the distance between two states i and j (points on the trajectory) is smaller than a threshold ε. Thus, this is a pairwise test between each state with each other state (N² tests if we have N states). The recurrence plot is a N × N matrix of black and white dots (usually a black dot corresponds with a recurrence), with two time-axes.

Recurrences are typical for dynamical systems. After some time, the state of a system will recur as close as one wishes to a former state. Recurrence plots can help to find a first characterization of the data or to find transitions and interrelations.

For a first illustration we will import a time series (cycles.dat) with containing cycles of 100, 40 and 20 periods and linearly transform it to an equidistant time scale

x0 = load('cycles.dat'); t = 0:3:996; x = interp1(x0(:,1),x0(:,2),t,'linear'); 

We start with the assumption that the phase space is only 1-dimensional. The calculation of the distances between all points of the phase space trajectory reveals the distance matrix S.

N = length(x); S = zeros(N, N);  for i = 1:N,     S(:,i) = abs( repmat( x(i), N, 1 ) - x(:) ); end 

Here we have used the vectorization feature of Matlab, which provides a faster computation than to compute each N² distances in two loops.

Now we plot the distance matrix

imagesc(t, t, S) axis square 

and the recurrence plot by applying a threshold of 0.1 to the distance matrix

imagesc(t, t, S < 1) axis square colormap([1 1 1;0 0 0]) xlabel('Time'), ylabel('Time') 

Both plots reveal periodically occuring patterns. The distances between these periodic patterns represent the cyclicities contained in the time series. The most significant periodic structures have periods of 200 and 100 time units. The period 200 is most significant because of the superposition of the cycles 100 and 40, which are common divisors of 200. Moreover, there are small substructures in the recurrence plot, which have sizes of 40 and 20 time units.

Now we calculate a recurrence plot of the Lorenz system (lorenz.dat).

x = load('lorenz.dat'); 

We chose the resampled first component of this system and reconstruct a phase space trajectory by using an embedding of m = 3 and τ = 4 (which corresponds to a delay of 0.34 sec).

t = x(1:5:3000,1); y = x(1:5:3000,2); m = 3; tau = 4;  N = length(y); N2 = N - tau * (m - 1); 

The original data series has a length of 600, but the resulting phase space trajectory has the length 596. Now we create the phase space trajectory with

clear xe for mi = 1:m;     xe(:, mi) = y([1:N2] + tau * (mi-1)); end 

We can accelerate the pair-wise test between each points on the trajectory with a fully vectorized algorithm. For that we need to transfer the trajectory vector into two test vectors, whose component-wise test will provide the pair-wise test of the trajectory vector:

x1 = repmat(xe, N2, 1); x2 = reshape(repmat(xe(:), 1, N2)', N2 * N2, m); 

Using these vectors we calculate the recurrence plot using the Euclidean norm without any loop

S = sqrt(sum( (x1 - x2) .^ 2, 2 )); S = reshape(S, N2, N2);  imagesc(t(1:N2), t(1:N2), S < 2) axis square colormap([1 1 1;0 0 0]) xlabel('Time (sec)'), ylabel('Time (sec)') 

This recurrence plot is rather homogenous with many short diagonal lines. These lines represent epochs, where the phase space trajectory runs parallel to former or later sequences of this trajectory, i.e. the states and the dynamics is similar at these times. The distances between these diagonal lines, representing the periods of the cycles, differ and are not constant (as they are in a harmonic oscillation).

Using the CRP toolbox for Matlab, recurrence plots can be computed much simpler:

X = crp(y(:,1),3,4,2,'euc','nonorm'); imagesc(t(1:size(X,1)), t(1:size(X,2)), X) axis square xlabel('Time (sec)'), ylabel('Time (sec)') 

The arguments after the data series are (in that order) embedding dimension, delay, threshold, norm and a switch (to avoid normalisation of the data). The arguments are optional. Without the output argument, the recurrence plot can be interactively modified by a GUI.

Using this toolbox it is also possible to compute cross recurrence plots, which visualise the pair-wise test of two different phase space trajectories, and joint recurrence plots, which compare the simultaneous recurrence structure of two phase space trajectories. Such bi- and multivariate extensions of recurrence plots offer nonlinear correlation tests and synchronization analysis. A detailed introduction in recurrence plot based methods can be found at the web site www.recurrence-plot.tk.

Exercise

1. Recurrence plot of the Roessler system

1.1. Import the Roessler system from the file roessler.dat and create the corresponding recurrence plot using all three components (first column in file is the time scale!).

1.2. Reconstruct the phase space by using the first component and then by the third component. Compute the corresponding recurrence plots for three different epochs in the data and comment on the results.

2. How should a data series look like that the corresponding recurrence plot contains a circle?

3. Consider the harmonic oscillations
a = sin((1:1000) * 2 * pi/67); and
b = sin(.01 * ([1:1000] * 2 * pi/67) .^ 2);.
Compute the cross recurrence plot between both oscillations by crp(a,b). Comment on the result.