From 3ec33ec9c0cdad9ac69a8f6db4c429a78a62085b Mon Sep 17 00:00:00 2001
From: marwan <>
Date: Tue, 1 Mar 2016 16:26:03 +0000
Subject: [PATCH] added support for multi-column phase space vectors
---
twinsurr.m | 60 +++++++++++++++++++++++++++++++-----------------------
1 file changed, 35 insertions(+), 25 deletions(-)
diff --git a/twinsurr.m b/twinsurr.m
index ec94fc0..8369108 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);
--
GitLab