Newer
Older
function [a_out, b_out]=pss(varargin)
%PSS Computes phase space size.
% [Y Z]=PSS(X [,M,T]) computes the maximal (Y) and the
% averaged (Z) phase space diameter of embedded data X,
% with the embedding parameters dimension M and lag T.
%
% [Y Z]=PSS(..., NORM) computes the maximal and averaged
% phase space diameter by using a specified norm, where
% NORM can be
% maxnorm - Maximum norm.
% euclidean - Euclidean norm (default).
% minnorm - Minimum norm.
%
% See also PHASESPACE, CRP, CRQA.
% Copyright (c) 2008-2009
% Norbert Marwan, Potsdam Institute for Climate Impact Research, Germany
% http://www.pik-potsdam.de
%
% Copyright (c) 1998-2008
% Norbert Marwan, Potsdam University, Germany
% http://www.agnld.uni-potsdam.de
%
% $Date$
% $Revision$
%
% $Log$
% Revision 2.1 2004/11/10 07:10:06 marwan
% initial import
%
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
%
% 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.
error(nargchk(1,4,nargin));
if nargout>2, error('Too many output arguments'), end
i_double=find(cellfun('isclass',varargin,'double'));
i_char=find(cellfun('isclass',varargin,'char'));
check_meth={'ma','eu','mi'}; % maxnorm, euclidean, minnorm
temp_meth=0;
if ~isempty(i_char)
for i=1:length(i_char),
varargin{i_char(i)}(3)='0';
temp_meth=strcmpi(varargin{i_char(i)}(1:2),check_meth');
end
method=min(find(temp_meth));
if isempty(method), method=2; end
if method>3, method=3; end
else
method=2;
end
x=double(varargin{i_double(1)});
if length(i_double)>1 & max(size(varargin{i_double(2)}))==1
m=varargin{i_double(2)};
else
m=1;
end
if length(i_double)>2 & max(size(varargin{i_double(3)}))==1
t=varargin{i_double(3)};
else
t=1;
end
Nx=length(x);NX=Nx-t*(m-1);
i=(1:NX)';j=0:t:0+(m-1)*t;
i=reshape(i(:,ones(1,m))+j(ones(NX,1),:),m*NX,1);
x1=x(i);
x2=reshape(x1,NX,m);
[NX, mx] = size(x2);
px = permute(x2, [ 1 3 2 ]);
py = permute(x2, [ 3 1 2 ]);
s1 = px(:,ones(1,NX),:) - py(ones(1,NX),:,:);
switch(method)
case{1}
max_dist = max(x(1:end-(m-1)*t))-min(x(1:end-(m-1)*t));
s = max(abs(s1),[],3);
case{2}
s = sum(s1.^2, 3);
max_dist = sqrt(max(s(:)));
case{3}
s = sum(abs(s1), 3);
max_dist = max(sum(x2,2))-min(sum(x2,2));
end
mean_dist = sqrt(mean(s(:)));
if nargout==2
b_out=mean_dist;
end
if nargout>0
a_out=max_dist;
else
max_dist
end