add new function to calculate microstates

Contents.m

 % CRP Toolbox
-% Version 5.28 (R37) 05-Sep-2023
+% Version 5.28 (R37b) 08-Mar-2024
 %
 %    ace            -    Finds optimal transformation and maximal correlation.
 %    adjust         -    Adjusts two two-column vectors.
@@ -69,4 +69,4 @@
 % along with this program; if not, write to the Free Software
 % Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 
-% Modified at 05-Sep-2023 16:18:11 by MAKEINSTALL
+% Modified at 08-Mar-2024 10:35:33 by MAKEINSTALL

crpclean.m

@@ -16,7 +16,7 @@ function crpclean
 % Norbert Marwan, Potsdam University, Germany
 % http://www.agnld.uni-potsdam.de
 %
-% Generation date: 05-Sep-2023 16:18:11
+% Generation date: 08-Mar-2024 10:35:33
 % $Date$
 % $Revision$
 

microstates.m

+function [M, P] = microstates(varargin)
+% MICROSTATES   Set of microstates from a recurrence plot.
+%    M=MICROSTATES(X) finds matrix M of all microstates from 
+%    recurrence matrix X. Microstates are numbered from 1 to 
+%    N, where N is the max. number of available microstates.
+%    Per default microstates have size K=2, thus N = 16.
+%
+%    [M,P]=MICROSTATES(X) further provides the set of microstate
+%    patterns in P as a (2x2xN) matrix. The microstate index in M
+%    corresponds to the pattern index in the 3rd dimension.
+%
+%    ...=MICROSTATES(X,K) uses microstates of size KxK. The 
+%    total number of microstates is N = K^(2^K).
+%
+%    Examples: a = sin(linspace(0,5*2*pi,1050)); % sine
+%              b = rand(1000,1); % noise
+%              Xa = crp(a,2,50,.2,'nonorm','nogui'); % recurrence plot for sine
+%              Xb = crp(b,1,1,.2,'nonorm','nogui'); % recurrence plot for noise
+%              K = 3; % size of microstates
+%              [Ma, P] = microstates(Xa, K); % microstates for sine
+%              Mb = microstates(Xb, K); % microstates for noise
+%              Ha = hist(Ma(:), 1:2^(K^2)); % histogram of microstates for sine
+%              Hb = hist(Mb(:), 1:2^(K^2)); % histogram of microstates for noise
+%
+%              % show microstates histograms
+%              subplot(2,1,1)
+%              bar(log10(Ha))
+%              ylabel('Frequency')
+%              title('Histogram microstates sine')
+%
+%              subplot(2,1,2)
+%              bar(log10(Hb))
+%              ylabel('Frequency'), xlabel('Microstate index')
+%              title('Histogram microstates noise')
+%
+%              % show most frequent microstates
+%              [~, idx] = sort(Ha, 'desc'); % sorted histogram
+%              clf
+%              for i = 1:16
+%                 subplot(4,4,i)
+%                 imagesc(P(:,:,idx(i))), caxis([0 1])
+%                 title(sprintf('Percentage: %2.2f', 100*Ha(idx(i))/sum(Ha)))
+%              end
+%              colormap([1 1 1; 0 0 0])
+%              sgtitle('Most frequent microstates')
+
+%
+%    See also CRQA, DL, TT.
+%
+%    References: 
+%    Corso, F., et al.:
+%    Quantifying entropy using recurrence matrix microstates, Chaos, 28, 2018.
+
+% Copyright (c) 2024-
+% Norbert Marwan, Potsdam Institute for Climate Impact Research, Germany
+% https://tocsy.pik-potsdam.de
+
+% This program is free software; you can redistribute it and/or
+% modify it under the terms of the GNU General Public License
+% as published by the Free Software Foundation; either version 2
+% of the License, or any later version.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% check and read the input
+narginchk(1,2)
+nargoutchk(0,2)
+
+% default values
+K = 2;
+
+warning off
+X = logical(varargin{1});
+if nargin == 2
+   K = varargin{2};
+end
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% find microstate patterns
+
+% create microstate reference pattern
+K2 = K.^2;
+
+% find microstate patterns
+patternTemplate = reshape(2.^(0:(K2-1)), K, K);
+M = conv2(X,patternTemplate, 'same') + 1;
+
+
+% create the corresponding microstate matrices for visualisation
+seq = dec2bin(0:(2^K2-1)); % translate integers to binary sequences
+patternSeq = zeros(K2, length(seq)); % sequential matrix of final patterns
+pattern = zeros(K, K, length(seq)); % matrix of final patterns
+for i = 1:length(seq)
+   for j = 1:K2
+      patternSeq(j,i) = str2num(seq(i,j));
+   end
+   pattern(:,:,i) = reshape(patternSeq(:,i), K, K);
+end
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% output results
+if nargout==2
+     P=pattern;
+end
+
+
+warning on
+

phasesynchro.m

-function CPRout = phasesynchro(varargin);
+function [CPRout RR] = phasesynchro(varargin);
 % PS   Indicator of phase synchronisation by means of recurrences.
 %    CPR=PHASESYNCHRO(X,Y [,param1,param2,...]) calculates the 
 %    index of phase synchronisation based on recurrences.
@@ -76,7 +76,7 @@ lmin=2;
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% check the input
 
 narginchk(1,9)
-nargoutchk(0,1)
+nargoutchk(0,2)
 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% read the input
@@ -287,3 +287,7 @@ end
 if nargout | nogui == 2
     CPRout = CPR;
 end
+
+if nargout == 2 | nogui == 2
+    RR = [rr1(:) rr2(:)];
+end