From 428ce439cb0267173cae80fe6a33ea46ffe42f87 Mon Sep 17 00:00:00 2001
From: marwan <>
Date: Tue, 1 Jul 2008 13:09:27 +0000
Subject: [PATCH] initial import
---
twinsurr.m | 157 +++++++++++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 157 insertions(+)
create mode 100644 twinsurr.m
diff --git a/twinsurr.m b/twinsurr.m
new file mode 100644
index 0000000..7b52ca5
--- /dev/null
+++ b/twinsurr.m
@@ -0,0 +1,157 @@
+function y = twinsurr(varargin)
+% TWINSURR Creates twin surrogates for statistical tests.
+% Y=TWINSURR(X) creates twin surrogates Y based on the vector X
+% using recurrences. The matrix Y contains 100 columns of
+% 100 twin surrogates. If X is a PxQ matrix, the resulting
+% surrogate matrix is Px100xQ.
+%
+% Y=TWINSURR(X,M,T,E,...) creates twin surrogates using
+% embedding dimension M, delay T, recurrence threshold E. The
+% input arguments are similar to those of the command CRP.
+%
+% Y=TWINSURR(X,M,T,E,...,N) creates N surrogates (default is 100).
+%
+% Example: x = rand(3,1);
+% a = [.8 .3 -.25 .9]';
+% for i = 4:1000,
+% x(i) = sum(a(1:3) .* x(i-1:-1:i-3)) + a(end) * randn;
+% end
+% xs = twinsurr(x,1,1,.1,'euc',10);
+%
+% See also CRP, RECONS.
+%
+% References:
+% Thiel, M., Romano, M. C., Kurths, J., Rolfs, M., Kliegl, R.:
+% Twin Surrogates to Test for Complex Synchronisation,
+% Europhys. Lett., 75, 2006.
+
+% Copyright (c) 2008 by AMRON
+% Norbert Marwan, Potsdam University, Germany
+% http://www.agnld.uni-potsdam.de
+%
+% $Date$
+% $Revision$
+%
+% $Log$
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% programme properties
+
+global errcode props
+
+init_properties
+nsur_init = 100;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% check the input
+
+error(nargchk(1,7,nargin));
+if nargout>1, error('Too many output arguments'), end
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% read the input
+
+% 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
+
+
+varargin{8}=[];
+i_double=find(cellfun('isclass',varargin,'double'));
+i_char=find(cellfun('isclass',varargin,'char'));
+check_meth={'ma','eu','mi','nr','rr','fa','in','om','op','di'}; % maxnorm, euclidean, nrmnorm, fan, distance
+check_norm={'non','nor'}; % nonormalize, normalize
+
+if isnumeric(varargin{1}) % read commandline input
+ % check the text input parameters for method, gui
+ temp_meth=0;
+ temp_norm=0;
+ if ~isempty(i_char)
+ for i=1:length(i_char),
+ varargin{i_char(i)}(4)='0';
+ temp_norm=temp_norm+strcmpi(varargin{i_char(i)}(1:3),check_norm');
+ temp_meth=temp_meth+strcmpi(varargin{i_char(i)}(1:2),check_meth');
+ end
+ method_n=min(find(temp_meth));
+ nonorm=min(find(temp_norm))-1;
+ for i=1:length(i_char); temp2(i,:)=varargin{i_char(i)}(1:3); end
+
+ if isempty(nonorm), nonorm=1; end
+ if nonorm>1, nonorm=1; end
+ if isempty(method_n), method_n=1; end
+ if method_n>length(check_meth), method0=length(check_meth); end
+ method=check_meth{method_n};
+ norm_str = check_norm{nonorm+1};
+ else
+ method = 'max'; norm_str = 'non';
+ end
+
+ % get the parameters for creating RP
+ if max(size(varargin{1}))<=3
+ disp('Error using ==> twinsurr')
+ disp('To less values in data X.')
+ return
+ end
+ x=double(varargin{1});
+
+ if ~isempty(varargin{i_double(2)}), m=varargin{i_double(2)}(1); else m=1; end
+ if ~isempty(varargin{i_double(3)}), t=varargin{i_double(3)}(1); else t=1; end
+ if ~isempty(varargin{i_double(4)}), e=varargin{i_double(4)}(1); else e=.1; end
+ if ~isempty(varargin{i_double(5)}), nsur=varargin{i_double(5)}(1); else nsur=nsur_init; end
+else
+ disp('Error using ==> twinsurr')
+ disp('No valid arguments.')
+ return
+end
+
+
+N = length(x);
+if size(x,1) < size(x,2), x = x'; end
+m = size(x,2);
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% make surrogates
+
+X = crp(x,m,t,e,method,norm_str,'sil');
+NX = length(X);
+
+%% find twins
+h = waitbar(0,'Searching Twins');
+for i = 1:NX, waitbar(i/NX)
+ A = repmat(X(:,i),1,NX);
+ S{i} = find(all(X == A));
+end
+delete(h)
+
+for k = 1:nsur
+ %% chose randomly first point
+ i = ceil(NX * rand);
+ xs = x(i,:);
+
+ %% jump randomly between twins until length of surrogate is reached
+ while length(xs) < N
+ j = S{i};
+ j_i = ceil(length(j) * rand);
+ i = j(j_i);
+ i = i + 1;
+ if i > NX
+ i = ceil(NX * rand);
+ continue
+ end
+ xs = [xs; x(i,:)];
+ end
+ y(:,k,1:m) = reshape(xs(:),N,1,m);
+end
+
--
GitLab