initial import

twinsurr.m

+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
+