Skip to content
Snippets Groups Projects
Commit 932cad0c authored by marwan's avatar marwan
Browse files

initial import

parent f5de7e1e
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
phasesynchro.m 0 → 100644
+ 289
− 0
View file @ 932cad0c
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
recons.m 0 → 100644
+ 228
− 0
View file @ 932cad0c
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
taucrp.m 0 → 100644
+ 274
− 0
View file @ 932cad0c
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)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment