From 932cad0c25f85e0f196450df012426a2614151e1 Mon Sep 17 00:00:00 2001
From: marwan <>
Date: Fri, 25 Jan 2008 12:47:26 +0000
Subject: [PATCH] initial import
---
phasesynchro.m | 289 +++++++++++++++++++++++++++++++++++++++++++++++++
recons.m | 228 ++++++++++++++++++++++++++++++++++++++
taucrp.m | 274 ++++++++++++++++++++++++++++++++++++++++++++++
3 files changed, 791 insertions(+)
create mode 100644 phasesynchro.m
create mode 100644 recons.m
create mode 100644 taucrp.m
diff --git a/phasesynchro.m b/phasesynchro.m
new file mode 100644
index 0000000..d9eaafd
--- /dev/null
+++ b/phasesynchro.m
@@ -0,0 +1,289 @@
+function CPRout = 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.
+%
+% CPR=PHASESYNCHRO(X,Y,M,T,E,W) uses the dimension M, delay T,
+% the size of neighbourhood E and the range W of past and future
+% time steps.
+%
+% If X and Y are multi-column vectors then they will be
+% considered as phase space vectors (TAUCRP can be used
+% for real phase space vectors without embedding).
+%
+% The call of PHASESYNCHRO without output arguments plots the
+% tau-recurrence rate and the CPR value in the current figure.
+%
+% Parameters: dimension M, delay T, the size of neighbourhood
+% E and the range W are the first four numbers after the data
+% series; further parameters can be used to switch between
+% various methods of finding the neighbours of the phasespace
+% trajectory and to suppress the normalization of the data.
+%
+% Methods of finding the neighbours/ of plot.
+% maxnorm - Maximum norm.
+% euclidean - Euclidean norm.
+% minnorm - Minimum norm.
+% maxnorm - Maximum norm, fixed recurrence rate.
+% fan - Fixed amount of nearest neighbours.
+%
+% Normalization of the data series.
+% normalize - Normalization of the data.
+% nonormalize - No normalization of the data.
+%
+% Suppressing the plot of the results.
+% silent - Suppresses the plot.
+%
+% Parameters not needed to be specified.
+%
+%
+% Examples: a = sin((1:1000) * 2 * pi/67);
+% b = sin((1:1000) * 2 * pi/67) + randn(1,1000);
+% phasesynchro(a,b,2,17,'nonorm','euclidean');
+%
+% See also CRP, CRP2, CRP_BIG, JRP, CRQA.
+%
+% References:
+% Marwan, N., Romano, M. C., Thiel, M., Kurths, J.:
+% Recurrence Plots for the Analysis of Complex Systems, Physics
+% Reports, 438(5-6), 2007.
+%
+% Romano, M. C., Thiel, M., Kurths, J., Kiss, I. Z., Hudson, J.:
+% Detection of synchronization for non-phase-coherent and
+% non-stationary data, Europhysics Letters, 71(3), 2005.
+
+% Copyright (c) 2008 by AMRON
+% Norbert Marwan, Potsdam University, Germany
+% http://www.agnld.uni-potsdam.de
+%
+% $Date$
+% $Revision$
+%
+% $Log$
+%
+%
+% Examples: a = sin((1:1000) * 2 * pi/67);
+% b = sin((1:1000) * 2 * pi/67) + rand(1,1000);
+% X = phasesynchro(a,b,2,17,'nonorm','euclidean');
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% programme properties
+
+global errcode props
+
+init_properties
+m_init = 1;
+tau_init = 1;
+eps_init = 0.1;
+w_init = 100;
+method_str={'Maximum Norm','Euclidean Norm','Minimum Norm','Maximum Norm, fixed RR','FAN'};
+norm_str = {'nor','non'};
+lmin=2;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% check the input
+
+error(nargchk(1,9,nargin));
+if nargout>1, error('Too many output arguments'), end
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% read the input
+
+
+check_meth={'ma','eu','mi','rr','fa'}; % maxnorm, euclidean, nrmnorm, fan
+check_norm={'non','nor'}; % nonormalize, normalize
+check_gui={'gui','nog','sil'}; % gui, nogui, silent
+
+if isnumeric(varargin{1}) % read commandline input
+ varargin{9}=[];
+ % transform any int to double
+ intclasses = {'uint8';'uint16';'uint32';'uint64';'int8';'int16';'int32';'int64'};
+ flagClass = [];
+ for i = 1:length(intclasses)
+ i_int=find(cellfun('isclass',varargin,intclasses{i}));
+ if ~isempty(i_int)
+ for j = 1:length(i_int)
+ varargin{i_int(j)} = double(varargin{i_int(j)});
+ end
+ flagClass = [flagClass; i_int(:)];
+ end
+ end
+ if ~isempty(flagClass)
+ disp(['Warning: Input arguments at position [',num2str(flagClass'),'] contain integer values']);
+ disp(['(now converted to double).'])
+ end
+ i_double=find(cellfun('isclass',varargin,'double'));
+ i_char=find(cellfun('isclass',varargin,'char'));
+
+ % check the text input parameters for method, gui and normalization
+ temp_meth=0;
+ temp_norm=0;
+ temp_gui=0;
+ if ~isempty(i_char)
+ for i=1:length(i_char),
+ varargin{i_char(i)}(4)='0';
+ temp_meth=temp_meth+strcmpi(varargin{i_char(i)}(1:2),check_meth');
+ temp_norm=temp_norm+strcmpi(varargin{i_char(i)}(1:3),check_norm');
+ temp_gui=temp_gui+strcmpi(varargin{i_char(i)}(1:3),check_gui');
+ end
+ method=min(find(temp_meth));
+ nonorm=min(find(temp_norm))-1;
+ nogui=min(find(temp_gui))-1;
+ if isempty(method), method=1; end
+ if isempty(nonorm), nonorm=1; end
+ if isempty(nogui), nogui=0; end
+ if method>length(check_meth), method=length(check_meth); end
+ if nonorm>1, nonorm=1; end
+ if nogui>2, nogui=2; end
+ else
+ method=1; nonorm=1; nogui=0;
+ end
+ if nogui==0 & nargout>0, nogui=1; end
+
+ % get the parameters for creating RP
+ if max(size(varargin{1}))<=3
+ error('To less values in data X.')
+ end
+ x=double(varargin{1});
+ if isempty(varargin{2}) | ~isnumeric(varargin{2}), y=x; else
+ y=double(varargin{2}); end
+
+ if (isnumeric(varargin{2}) & max(size(varargin{2}))==1) | ~isnumeric(varargin{2})
+ y=x;
+ if ~isempty(varargin{i_double(2)}), m0=varargin{i_double(2)}(1); else m0=m_init; end
+ if ~isempty(varargin{i_double(3)}), t=varargin{i_double(3)}(1); else t=tau_init; end
+ if ~isempty(varargin{i_double(4)}), e=varargin{i_double(4)}(1); else e=eps_init; end
+ if ~isempty(varargin{i_double(5)}), w=varargin{i_double(5)}(1); else w=w_init; end
+ else
+ if ~isempty(varargin{i_double(3)}), m0=varargin{i_double(3)}(1); else m0=m_init; end
+ if ~isempty(varargin{i_double(4)}), t=varargin{i_double(4)}(1); else t=tau_init; end
+ if ~isempty(varargin{i_double(5)}), e=varargin{i_double(5)}(1); else e=eps_init; end
+ if ~isempty(varargin{i_double(6)}), w=varargin{i_double(6)}(1); else w=w_init; end
+ end
+ t=round(t); m0=round(m0); mflag=method;
+ if e<0, e=1; disp('Warning: The threshold size E cannot be negative and is now set to 1.'), end
+ if t<1, t=1; disp('Warning: The delay T cannot be smaller than one and is now set to 1.'), end
+ if m0 < 1, m0 = 1; end
+ if t < 1, t = 1; end
+ if size(x,1)==1, x=x'; end, if size(y,1)==1, y=y'; end
+ m=max([size(x,2) size(y,2)]);
+ if w > size(x,1); w = size(x,1); end
+
+ if method==8 & (m*m0) > 1,
+ m0=1;
+ error(['The neighbourhood criterion ''Oder matrix''',10,'is not implemented - use crp or crp_big instead.'])
+ end
+ if method==9 & (m*m0) == 1,
+ m0=2;
+ disp(['Warning: For order patterns recurrence plots the dimension must',10,...
+ 'be larger than one. ',...
+ 'Embedding dimension is set to ',num2str(m0),'.'])
+ end
+ action='init';
+
+ if ~isempty(find(isnan(x)))
+ disp('NaN detected (in first variable) - will be cleared.')
+ for k=1:size(x,2), x(find(isnan(x(:,k))),:)=[]; end
+ end
+ if ~isempty(find(isnan(y)))
+ disp('NaN detected (in second variable) - will be cleared.')
+ for k=1:size(y,2), y(find(isnan(y(:,k))),:)=[]; end
+ end
+ if size(x,1) < t*(m-1)+1 | size(y,1) < t*(m-1)+1
+ error(['Too less data',10,...
+ 'Either too much NaN or the number of columns in the vectors do not match.'])
+ end
+
+ Nx=size(x,1); Ny=size(y,1);
+ NX=Nx-t*(m0-1);NY=Ny-t*(m0-1);
+ x0=zeros(Nx,m);y0=zeros(Ny,m);
+ x0(1:size(x,1),1:size(x,2))=x;
+ y0(1:size(y,1),1:size(y,2))=y;
+
+
+ if ~isempty(find(isnan(x))), for k=1:size(x,2), x(find(isnan(x(:,k))),:)=[]; end, end
+ if ~isempty(find(isnan(y))), for k=1:size(y,2), y(find(isnan(y(:,k))),:)=[]; end, end
+ if size(x,1) < t*(m0-1)+1 | size(y,1) < t*(m0-1)+1
+ error(['Too less data',10,...
+ 'Either too much NaN or the number of columns in the vectors do not match.'])
+ end
+
+else
+ error('No valid input given!')
+end
+
+
+
+
+% compute RR-tau
+X = taucrp(x, m, t, e, w, check_meth(method),check_norm(nonorm+1));
+X_ = taucrp(y, m, t, e, w, check_meth(method),check_norm(nonorm+1));
+
+
+%
+X = [zeros(2*w+1,1) X zeros(2*w+1,1)];
+X_ = [zeros(2*w+1,1) X_ zeros(2*w+1,1)];
+x_diff = diff(X,[],2); % mark start and end points of line structures with -1 and +1
+x_diff_ = diff(X_,[],2); % mark start and end points of line structures with -1 and +1
+
+[line_start col1] = find(x_diff' == 1); % index of start points
+[line_end col2] = find(x_diff' == -1); % index of end points
+[line_start_ col1_] = find(x_diff_' == 1); % index of start points
+[line_end_ col2_] = find(x_diff_' == -1); % index of end points
+
+
+l = line_end - line_start; % length of all lines
+l_corr = l; l_corr(l < lmin) = 0; % only lines longer than lmin
+tau = unique(col1); % lag at which the lines where found
+l_ = line_end_ - line_start_; % length of all lines
+l_corr_ = l_; l_corr_(l_ < lmin) = 0; % only lines longer than lmin
+tau_ = unique(col1_); % lag at which the lines where found
+
+
+RR = zeros(2*w+1,1); RR_ = zeros(2*w+1,1);
+for i = tau'; % waitbar(i/max(tau))
+ j = col1==i;
+ N_recpoints = sum(l(j)); % number of recurrence points at a given lag
+ N_pointsline = sum(l_corr(j)); % number of recurrence points forming lines
+ DET(i) = N_pointsline/N_recpoints;
+ L(i) = mean(l_corr(j));
+ RR(i) = sum(l(j))/(NX-abs(w-i));
+end
+for i = tau_';
+ j = col1_==i;
+ N_recpoints = sum(l_(j)); % number of recurrence points at a given lag
+ N_pointsline = sum(l_corr_(j)); % number of recurrence points forming lines
+ DET_(i) = N_pointsline/N_recpoints;
+ L_(i) = mean(l_corr_(j));
+ RR_(i) = sum(l_(j))/(NX-abs(w-i));
+end
+RR(find(isnan(RR)))=0; RR_(find(isnan(RR)))=0;
+
+
+
+rr1 = normalize(RR(w+1:end));
+rr2 = normalize(RR_(w+1:end));
+
+
+
+% ACF of RR-tau
+a1 = xcov(rr1,'unbiased'); a1(1:w) = [];
+i1 = find(a1 < 1/exp(1),1);
+a2 = xcov(rr2,'unbiased'); a2(1:w) = [];
+i2 = find(a2 < 1/exp(1),1);
+i3 = max([i1 i2]);
+
+
+% correlation coefficient
+C = corrcoef(rr1(i3:end), rr2(i3:end));
+CPR = C(1,2);
+
+% plot
+if nogui ~= 2
+ plot(0:i3,rr1(1:i3+1),':',0:i3,rr2(1:i3+1),':'), hold on
+ plot(i3:w,rr1(i3+1:end),i3:w,rr2(i3+1:end)), hold off
+ title(['CPR: ',sprintf('%2.2f',CPR)], 'fontw','bold')
+ xlabel('Lag'), ylabel('RR_{\tau}')
+end
+if nargout | nogui == 2
+ CPRout = CPR;
+end
diff --git a/recons.m b/recons.m
new file mode 100644
index 0000000..3ce53b6
--- /dev/null
+++ b/recons.m
@@ -0,0 +1,228 @@
+function xout=recons(varargin)
+%RECONS Reconstruct a time series from a recurrence plot.
+% Y = RECONS(X) reconstructs a time series Y from the
+% recurrence plot in the matrix X.
+%
+% Y = RECONS(X,NAME) reconstructs the time series using
+% the named cumulative distribution function, which can
+% be 'norm' or 'Normal' (defeault), 'unif' or 'Uniform'.
+%
+% Y = RECONS(X,P) reconstructs the time series using
+% the cumulative distribution function given by vector P.
+%
+% See also CRP, CRP2, JRP.
+%
+% References:
+% Thiel, M., Romano, M. C., Kurths, J.:
+% How much information is contained in a recurrence plot?,
+% Phys. Lett. A, 330, 2004.
+
+% Copyright (c) 2005-2006 by AMRON
+% Norbert Marwan, Potsdam University, Germany
+% http://www.agnld.uni-potsdam.de
+%
+% $Date$
+% $Revision$
+%
+% $Log$
+%
+%
+%
+% This program is part of the new generation XXII series.
+%
+% 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.
+
+
+errcode = 0;
+error(nargchk(1,2,nargin));
+if nargout>1, error('Too many output arguments'), end
+
+if isnumeric(varargin{1})==1 % read RP
+ X = double(varargin{1});
+else
+ error('Input must be numeric.')
+end
+
+
+
+N = size(X); maxN = prod(N)+1;
+if N(1) ~= N(2)
+ error('Input must be square matrix.')
+end
+
+try
+% final distribution
+errcode = 1;
+i=(-5:10/(N(1)-1):5)';
+fs=cumsum(10*(1/(sqrt(2*pi))).*exp(-(i).^2./2));
+P=interp1(fs,i,[1:length(fs)],'nearest')';
+
+if nargin > 1 & isnumeric(varargin{2})==1
+ P = varargin{2};
+elseif nargin > 1 & ischar(varargin{2})==1
+ if strcmpi(varargin{2}(1:3), 'uni')
+ P = (0:1/(N(1)-1):1)';
+ disp('uni')
+ end
+end
+P(isnan(P)) = []; % remove NaN from distribution
+
+
+h_waitbar = waitbar(0,'Remove double columns'); drawnow
+
+% remove double rows
+errcode = 2;
+[X2 uniI J]=unique(X,'rows');
+uniI = sort(uniI);
+
+% this are the removed rows
+rem = setdiff(1:length(X2),uniI);
+%X3(rem,:) = X(rem,:);
+
+% this are the corresponding columns
+for i = 1:length(rem), waitbar(i/length(rem))
+% [dummy, dblI(i), a2] = intersect(X, X(rem(i),:), 'rows');
+ X2 = X;
+ X2(rem,:) = -1;
+ dblI(i) = find(all((X2 == repmat(X(rem(i),:), size(X2,2), 1))'));
+end
+
+% number of non-neighbours
+errcode = 3;
+
+waitbar(0, h_waitbar, 'Search neighbours'); drawnow
+NN = maxN*ones(N);
+%N = zeros(size(X1));
+for i = 1:size(X,1), waitbar(i/size(X,1))
+ if ~ismember(i,rem)
+ % indices of RPs at i
+ I = find(X(i,:));
+ I(ismember(I,rem)) = [];
+
+ % number of non-neighbours
+ n = sum((X(I,:) - repmat(X(i,:),length(I),1)) == 1,2);
+ NN(i,I) = n';
+ end
+
+end
+errcode = 4;
+
+% remove entries on the main diagonal in N
+i(1:size(NN,1)) = 1; NN(find(diag(i))) = maxN;
+
+% determine the both indices where N == 0 for all i
+iNaN = NN == maxN; i = find(iNaN); NN2 = NN; NN2(i) = zeros(size(i)); % remove NaNs
+n = sum(NN2,1);
+borderI=find(n==0);
+borderI(ismember(borderI,rem)) = [];
+
+k = borderI(1);
+r = k; % start point in the rank order vector
+
+errcode = 5;
+waitbar(0, h_waitbar, 'Reconstructing time series'); drawnow
+while min(NN(:)) < maxN
+ % look for the minimal N(k,:)
+ waitbar(length(r)/N(1))
+ [dummy kneu] = min(NN(k,:));
+ if ~isempty(dummy) kneu = find(NN(k,:) == dummy); else break, end
+ if length(kneu) > 1 % if a set of minimal N exist
+ [dummy kneu_index] = min(NN(kneu,k));
+ kneu=kneu(kneu_index);
+ end
+ NN(:,k) = maxN;
+ kold = k;
+ k=kneu;
+ r = [r; k]; % add the new found index to the rank order vector
+end
+
+errcode = 6;
+% adjust length of the distribution
+P = interp1((1:length(P))/length(P), P, (1:length(r))/length(r),'cubic');
+
+% assign the distribution to the rank order series
+waitbar(0, h_waitbar, 'Assigning distribution'); drawnow
+clear xneu
+for i = 1:length(r);
+ xneu(r(i)) = P(i);
+ waitbar(i/length(r))
+end
+
+errcode = 7;
+
+% add the removed (doubled) elements
+for i =1:length(rem)
+ xneu(rem(i)) = xneu(dblI(i));
+end
+delete(h_waitbar)
+
+
+errcode = 8;
+
+if nargout
+ xout = xneu;
+else
+ xneu
+end
+
+
+
+
+%%%%%%% error handling
+
+%if 0
+catch
+ z=whos;x=lasterr;y=lastwarn;in=varargin{1};
+ if ~isempty(findobj('Tag','TMWWaitbar')), delete(findobj('Tag','TMWWaitbar')), end
+ if ~strcmpi(lasterr,'Interrupt')
+ fid=fopen('error.log','w');
+ err=fprintf(fid,'%s\n','Please send us the following error report. Provide a brief');
+ err=fprintf(fid,'%s\n','description of what you were doing when this problem occurred.');
+ err=fprintf(fid,'%s\n','E-mail or FAX this information to us at:');
+ err=fprintf(fid,'%s\n',' E-mail: marwan@agnld.uni-potsdam.de');
+ err=fprintf(fid,'%s\n',' Fax: ++49 +331 977 1142');
+ err=fprintf(fid,'%s\n\n\n','Thank you for your assistance.');
+ err=fprintf(fid,'%s\n',repmat('-',50,1));
+ err=fprintf(fid,'%s\n',datestr(now,0));
+ err=fprintf(fid,'%s\n',['Matlab ',char(version),' on ',computer]);
+ err=fprintf(fid,'%s\n',repmat('-',50,1));
+ err=fprintf(fid,'%s\n',x);
+ err=fprintf(fid,'%s\n',y);
+ err=fprintf(fid,'%s\n',[' during ==> recons']);
+ if ~isempty(errcode)
+ err=fprintf(fid,'%s\n',[' errorcode ==> ',num2str(errcode)]);
+ else
+ err=fprintf(fid,'%s\n',[' errorcode ==> no errorcode available']);
+ end
+ err=fclose(fid);
+ disp('----------------------------');
+ disp(' ERROR OCCURED ');
+ disp([' during executing recons']);
+ disp('----------------------------');
+ disp(x);
+ if ~isempty(errcode)
+ disp([' errorcode is ',num2str(errcode)]);
+ else
+ disp(' no errorcode available');
+ end
+ disp('----------------------------');
+ disp(' Please send us the error report. For your convenience, ')
+ disp(' this information has been recorded in: ')
+ disp([' ',fullfile(pwd,'error.log')]), disp(' ')
+ disp(' Provide a brief description of what you were doing when ')
+ disp(' this problem occurred.'), disp(' ')
+ disp(' E-mail or FAX this information to us at:')
+ disp(' E-mail: marwan@agnld.uni-potsdam.de')
+ disp(' Fax: ++49 +331 977 1142'), disp(' ')
+ disp(' Thank you for your assistance.')
+ warning('on')
+ end
+ try, set(0,props.root), end
+ set(0,'ShowHidden','Off')
+ return
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/taucrp.m b/taucrp.m
new file mode 100644
index 0000000..57c178f
--- /dev/null
+++ b/taucrp.m
@@ -0,0 +1,274 @@
+function X = taucrp(varargin)
+%TAUCRP Creates a close returns plot.
+% R=TAUCRP(X [,Y] [,param1,param2,...]) creates a cross
+% recurrence plot/ recurrence plot R for a limited range of
+% past and future states, also known as close returns plot.
+%
+% R=TAUCRP(X [,Y],M,T,E,W) uses the dimension M, delay T, the size
+% of neighbourhood E and the range W of past and future time
+% steps.
+%
+% If X and Y are multi-column vectors then they will be
+% considered as phase space vectors (TAUCRP can be used
+% for real phase space vectors without embedding).
+%
+% Parameters: dimension M, delay T, the size of neighbourhood
+% E and the range W are the first four numbers after the data
+% series; further parameters can be used to switch between
+% various methods of finding the neighbours of the phasespace
+% trajectory and to suppress the normalization of the data.
+%
+% Methods of finding the neighbours/ of plot.
+% maxnorm - Maximum norm.
+% euclidean - Euclidean norm.
+% minnorm - Minimum norm.
+% maxnorm - Maximum norm, fixed recurrence rate.
+% fan - Fixed amount of nearest neighbours.
+% distance - Distance coded matrix (global CRP, Euclidean norm).
+%
+% Normalization of the data series.
+% normalize - Normalization of the data.
+% nonormalize - No normalization of the data.
+%
+% Parameters not needed to be specified.
+%
+%
+% Examples: a = sin((1:1000) * 2 * pi/67);
+% w = 160;
+% X = taucrp(a,2,17,.2,w,'nonorm','euclidean');
+% imagesc(1:size(X,2),-w:w,X), colormap([1 1 1; 0 0 0])
+%
+% See also CRP, CRP2, CRP_BIG, JRP, CRQA.
+%
+% References:
+% Marwan, N., Romano, M. C., Thiel, M., Kurths, J.:
+% Recurrence Plots for the Analysis of Complex Systems, Physics
+% Reports, 438(5-6), 2007.
+
+% Copyright (c) 2008 by AMRON
+% Norbert Marwan, Potsdam University, Germany
+% http://www.agnld.uni-potsdam.de
+%
+% $Date$
+% $Revision$
+%
+% $Log$
+
+
+warning off
+global errcode props nonorm
+
+errcode=0;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% programme properties
+
+init_properties
+m_init = 1;
+tau_init = 1;
+eps_init = 0.1;
+w_init = 100;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% check and read the input
+
+error(nargchk(1,8,nargin));
+if nargout>1, error('Too many output arguments'), end
+
+check_meth={'ma','eu','mi','rr','fa','di'}; % maxnorm, euclidean, nrmnorm, fan, distance
+check_norm={'non','nor'}; % nonormalize, normalize
+check_gui={'gui','nog','sil'}; % gui, nogui, silent
+
+if isnumeric(varargin{1}) % read commandline input
+ varargin{9}=[];
+ % transform any int to double
+ intclasses = {'uint8';'uint16';'uint32';'uint64';'int8';'int16';'int32';'int64'};
+ flagClass = [];
+ for i = 1:length(intclasses)
+ i_int=find(cellfun('isclass',varargin,intclasses{i}));
+ if ~isempty(i_int)
+ for j = 1:length(i_int)
+ varargin{i_int(j)} = double(varargin{i_int(j)});
+ end
+ flagClass = [flagClass; i_int(:)];
+ end
+ end
+ if ~isempty(flagClass)
+ disp(['Warning: Input arguments at position [',num2str(flagClass'),'] contain integer values']);
+ disp(['(now converted to double).'])
+ end
+ i_double=find(cellfun('isclass',varargin,'double'));
+ i_char=find(cellfun('isclass',varargin,'char'));
+
+ % check the text input parameters for method, gui and normalization
+ temp_meth=0;
+ temp_norm=0;
+ temp_gui=0;
+ if ~isempty(i_char)
+ for i=1:length(i_char),
+ varargin{i_char(i)}(4)='0';
+ temp_meth=temp_meth+strcmpi(varargin{i_char(i)}(1:2),check_meth');
+ temp_norm=temp_norm+strcmpi(varargin{i_char(i)}(1:3),check_norm');
+ temp_gui=temp_gui+strcmpi(varargin{i_char(i)}(1:3),check_gui');
+ end
+ method=min(find(temp_meth));
+ nonorm=min(find(temp_norm))-1;
+ nogui=min(find(temp_gui))-1;
+ if isempty(method), method=1; end
+ if isempty(nonorm), nonorm=1; end
+ if isempty(nogui), nogui=0; end
+ if method>length(check_meth), method=length(check_meth); end
+ if nonorm>1, nonorm=1; end
+ if nogui>2, nogui=2; end
+ else
+ method=1; nonorm=1; nogui=0;
+ end
+ if nogui==0 & nargout>0, nogui=1; end
+
+ % get the parameters for creating RP
+ if max(size(varargin{1}))<=3
+ error('To less values in data X.')
+ end
+ x=double(varargin{1});
+ if isempty(varargin{2}) | ~isnumeric(varargin{2}), y=x; else
+ y=double(varargin{2}); end
+
+ if (isnumeric(varargin{2}) & max(size(varargin{2}))==1) | ~isnumeric(varargin{2})
+ y=x;
+ if ~isempty(varargin{i_double(2)}), m0=varargin{i_double(2)}(1); else m0=m_init; end
+ if ~isempty(varargin{i_double(3)}), t=varargin{i_double(3)}(1); else t=tau_init; end
+ if ~isempty(varargin{i_double(4)}), e=varargin{i_double(4)}(1); else e=eps_init; end
+ if ~isempty(varargin{i_double(5)}), w=varargin{i_double(5)}(1); else w=w_init; end
+ else
+ if ~isempty(varargin{i_double(3)}), m0=varargin{i_double(3)}(1); else m0=m_init; end
+ if ~isempty(varargin{i_double(4)}), t=varargin{i_double(4)}(1); else t=tau_init; end
+ if ~isempty(varargin{i_double(5)}), e=varargin{i_double(5)}(1); else e=eps_init; end
+ if ~isempty(varargin{i_double(6)}), w=varargin{i_double(6)}(1); else w=w_init; end
+ end
+ t=round(t); m0=round(m0); mflag=method;
+ if e<0, e=1; disp('Warning: The threshold size E cannot be negative and is now set to 1.'), end
+ if t<1, t=1; disp('Warning: The delay T cannot be smaller than one and is now set to 1.'), end
+ if m0 < 1, m0 = 1; end
+ if t < 1, t = 1; end
+ if size(x,1)==1, x=x'; end, if size(y,1)==1, y=y'; end
+ m=max([size(x,2) size(y,2)]);
+ if w > size(x,1); w = size(x,1); end
+
+ if method==8 & (m*m0) > 1,
+ m0=1;
+ error(['The neighbourhood criterion ''Oder matrix''',10,'is not implemented - use crp or crp_big instead.'])
+ end
+ if method==9 & (m*m0) == 1,
+ m0=2;
+ disp(['Warning: For order patterns recurrence plots the dimension must',10,...
+ 'be larger than one. ',...
+ 'Embedding dimension is set to ',num2str(m0),'.'])
+ end
+ action='init';
+
+ if ~isempty(find(isnan(x)))
+ disp('NaN detected (in first variable) - will be cleared.')
+ for k=1:size(x,2), x(find(isnan(x(:,k))),:)=[]; end
+ end
+ if ~isempty(find(isnan(y)))
+ disp('NaN detected (in second variable) - will be cleared.')
+ for k=1:size(y,2), y(find(isnan(y(:,k))),:)=[]; end
+ end
+ if size(x,1) < t*(m-1)+1 | size(y,1) < t*(m-1)+1
+ error(['Too less data',10,...
+ 'Either too much NaN or the number of columns in the vectors do not match.'])
+ end
+
+ Nx=size(x,1); Ny=size(y,1);
+ NX=Nx-t*(m0-1);NY=Ny-t*(m0-1);
+ x0=zeros(Nx,m);y0=zeros(Ny,m);
+ x0(1:size(x,1),1:size(x,2))=x;
+ y0(1:size(y,1),1:size(y,2))=y;
+
+ if nonorm==1,
+ x=(x0-repmat(mean(x0),Nx,1))./repmat(std(x0),Nx,1);
+ y=(y0-repmat(mean(y0),Ny,1))./repmat(std(y0),Ny,1);
+ end
+
+ if ~isempty(find(isnan(x))), for k=1:size(x,2), x(find(isnan(x(:,k))),:)=[]; end, end
+ if ~isempty(find(isnan(y))), for k=1:size(y,2), y(find(isnan(y(:,k))),:)=[]; end, end
+ if size(x,1) < t*(m0-1)+1 | size(y,1) < t*(m0-1)+1
+ error(['Too less data',10,...
+ 'Either too much NaN or the number of columns in the vectors do not match.'])
+ end
+
+else
+ error('No valid input given!')
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if m0>1
+ x2=x(1:end-t*(m0-1),:);
+ y2=y(1:end-t*(m0-1),:);
+ for i=1:m0-1,
+ x2(:,m*i+1:m*(i+1))=x(1+t*i:end-t*(m0-i-1),:);
+ y2(:,m*i+1:m*(i+1))=y(1+t*i:end-t*(m0-i-1),:);
+ end
+ x=x2; y=y2; Nx=size(x,1); Ny=size(y,1);
+ m=m0*m; clear x2 y2
+end
+
+x1=repmat(x,1,Ny);
+for mi=1:m, x2(:,mi)=reshape(rot90(x1(:,0+mi:m:Ny*m+mi-m)),Nx*Ny,1); end
+y1=repmat(y,Nx,1); x1=x2; clear x2
+
+
+% create embedding vectors
+NX = Nx - t*(m0-1); NY = Ny - t*(m0-1);
+
+[idx add] = meshgrid(linspace(0,(m0-1)*t,m0), 1:Nx-(m0-1)*t);
+idx = idx + add;
+x1 = []; y1 = [];
+for i = 1:m
+ x1 = [x1 reshape(x(idx,i),size(idx,1),size(idx,2))];
+ y1 = [y1 reshape(y(idx,i),size(idx,1),size(idx,2))];
+end
+
+
+% indices for which the recurrences should be computed
+[i1 i2] = meshgrid(1:NX,-w:w);
+i2 = i2 + i1;
+
+i_remove = i2 > length(y1) | i2 < 1;
+i2(i_remove) = 1;
+
+% compute distance matrix
+D = (x1(i1,:) - y1(i2,:));
+
+switch method
+ case 1
+ D2 = max(abs(D),[],2);
+ X = double(D2 < e);
+ case 2
+ D2 = sqrt(sum(D.^2, 2));
+ X = double(D2 < e);
+ case 3
+ D2 = sum(abs(D), 2);
+ X = double(D2 < e);
+ case 4
+ D2 = max(abs(D),[],2);
+ SS = sort(D2(:));
+ idx = ceil(e * length(SS));
+ e_ = SS(idx);
+ X = double(D2 < e_);
+ case 5
+ D2 = sqrt(sum(D.^2, 2));
+ D3 = reshape(D2,size(i1,1),size(i1,2));
+ [SS, JJ] = sort(D3,1); JJ = JJ';
+ mine = round((2*w+1)*e);
+ X1(NX*(2*w+1)) = 0;
+ X1(JJ(:,1:mine)+repmat([0:(2*w+1):NX*(2*w+1)-1]',1,mine)) = 1;
+ X = reshape(X1,(2*w+1),NX);
+ case 6
+ D2 = sqrt(sum(D.^2, 2));
+ X = double(D2);
+
+ end
+
+clear X1 SS JJ s px D D_ D2 D3
+
+X = reshape(X,size(i1,1),size(i1,2)); X(i_remove) = 0;
+%imagesc(X)
--
GitLab