diff --git a/twinsurr.m b/twinsurr.m
index ec94fc0ce524d19e6cf3e8fc5a73ebb92ed65ee7..836910850e1041cb09afd2a35fd78b60f9030bf2 100644
--- a/twinsurr.m
+++ b/twinsurr.m
@@ -45,6 +45,9 @@ function [y nTw] = twinsurr(varargin)
 % $Revision$
 %
 % $Log$
+% Revision 5.5  2016/03/01 15:56:47  marwan
+% complete recoding based on an implementation by Maik Riedl (much faster calculations)
+%
 % Revision 5.4  2009/03/24 08:33:47  marwan
 % copyright address changed
 %
@@ -148,19 +151,20 @@ end
 
 
 if size(x,1) < size(x,2), x = x'; end
-if size(x,2) > 1
-    warning('Multicolumn vectors are not supported. Only first column will be used.')
-    x = x(:,1);
+if size(x,2) > 1 & m > 1
+    warning('Embedding for multi-column phase space vectors is not supported. Continue without embedding.')
+    m = 1;
 end
 
 N = length(x); 
 
-
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% main part: create surrogates
 
-
-surr = zeros(N,nsur);
-
+if size(x,2) > 1
+   surr = zeros(N,size(x,2),nsur);
+else
+   surr = zeros(N,nsur);
+end
 
 %% find twins
 h = [];
@@ -169,6 +173,7 @@ twinList = findTwins(x,m,t,e,method,nonorm,sil,h);
 nTwins = length(twinList); % number of twins 
 
 if nTwins < 1
+    if ~sil & exist('h') & ishandle(h), close(h), end
     error(['No twins found!',10,'The derived surrogates would identical with the original data.']')
 end
 
@@ -203,8 +208,12 @@ while cnt < nsur
     end
     if length(idxSurr)==N
         cnt = cnt + 1;
-        surr(:,cnt)=x(idxSurr);
-        if ~sil & exist('h') & mod(cnt,100)==0
+        if size(x,2) > 1
+           surr(:,:,cnt)=x(idxSurr,:);
+        else
+           surr(:,cnt)=x(idxSurr,:);
+        end
+        if ~sil & ~isempty(h) & mod(cnt,100)==0
             waitbar((2*N+N*cnt/nsur)/(3*N),h,'Create surrogates')
         end
     end
@@ -213,7 +222,7 @@ end
 if ~sil & exist('h') & ishandle(h), close(h), end
 
 
-
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% output
 if nargout > 0;
    y = surr;
 else
@@ -227,12 +236,7 @@ end
 
 
 
-
-
-
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% helper function (findTwins)
 
 
 function twinList=findTwins(x,m,t,thr,method,nonorm,sil,h)
@@ -248,7 +252,6 @@ function twinList=findTwins(x,m,t,thr,method,nonorm,sil,h)
 %threshold THR.
 %Twins are points with the same neighbourhood.
 
-%% Reconstruct phase space
 if nonorm
     if ~sil, disp('Normalize data.'), end
     x = zscore(x);
@@ -256,12 +259,19 @@ end
 
 updateTime = 100;
 N = length(x);
-Nx = N - (m-1)*t;
-xEmb = zeros(Nx,m);                
-for i = 1:m
-    xEmb(:,i) = x(((i-1)*t+1):(end-t*(m-i)));
-end
 
+%% Reconstruct phase space if m is set
+if m > 1
+    Nx = N - (m-1)*t;
+    xEmb = zeros(Nx,m);                
+    for i = 1:m
+        xEmb(:,i) = x(((i-1)*t+1):(end-t*(m-i)));
+    end
+else
+    xEmb = x;
+    Nx = N;
+end
+    
 %% Find recurrence points (= neighbours) for each state
 % this is a list of neighbours for each time point/state
 idxNeighbours = [];
@@ -270,11 +280,11 @@ for i = 1:Nx, if ~sil & ~isempty(h) & mod(i,updateTime)==0, waitbar(i/(3*N)), en
     % distance between current state and all other states
     switch method
         case 'eu'
-           d = sqrt(sum((xEmb(1:end,:)-xEmb(i,:)).^2,2));
+           d = sqrt(sum((xEmb-repmat(xEmb(i,:),Nx,1)).^2,2));
         case 'mi'
-           d = sum(abs(xEmb(1:end,:)-xEmb(i,:)),2);
+           d = sum(abs(xEmb-repmat(xEmb(i,:),Nx,1)),2);
         otherwise
-           d = max(abs(xEmb(1:end,:)-xEmb(i,:)),[],2);
+           d = max(abs(xEmb-repmat(xEmb(i,:),Nx,1)),[],2);
     end
     % indices of the Neighbours
     idxNeighbours{i} = find(d<thr);