*** empty log message ***

erqa.m

@@ -42,21 +42,21 @@ function xout=erqa(varargin)
 %    Parameters not needed to specify.
 %
 %    Output:
-%      Y(:,1)=K1
-%      Y(:,2)=K2
-%      Y(:,3)=D2
-%      Y(:,4)=I2
+%      Y(:,1) = K1
+%      Y(:,2) = K2
+%      Y(:,3) = D2
+%      Y(:,4) = I2
 %      
-%    Examples: a=randn(300,1);
-%              crqa(a,1,1,.2,40,2,'max')
+%    Examples: a = randn(300,1);
+%              erqa(a,1,1,.2,40,2,'max')
 %
-%              N=300; w=40; ws=2;
-%              a=3.4:.6/(N-1):4;
-%              b=.5; for i=2:N, b(i)=a(i)*b(i-1)*(1-b(i-1));end
-%              y=crqa(b,3,2,.1,w,ws,'nogui');
+%              N = 300; w = 40; ws = 2;
+%              a = 3.4:.6/(N-1):4;
+%              b = .5; for i = 2:N, b(i) = a(i)*b(i-1)*(1-b(i-1)); end
+%              y = erqa(b,3,2,.1,w,ws,'nogui');
 %              subplot(2,1,1), plot(a,b,'.','markersize',.1)
 %              title('logistic map'), axis([3.4 4 0 1])
-%              subplot(2,1,2), plot(a(1:ws:N-w),y(1:ws:N-w,1))
+%              subplot(2,1,2), plot(a(1:ws:N-w), y(1:ws:N-w,1))
 %              ylabel('recurrence rate'), axis([3.4 4 0 1])
 %
 %      
@@ -66,6 +66,11 @@ function xout=erqa(varargin)
 %    Thiel, M., Romano, M. C., Kurths, J.: 
 %    Analytical description of Recurrence Plots of stochastic
 %    and chaotic processes, to be publ. 2002
+%
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% % Please note, that this routine is not yet finished! %
+% % Don't use the results of this script!               %
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 % Copyright (c) 1998-2003 by AMRON
 % Norbert Marwan, Potsdam University, Germany
@@ -75,6 +80,9 @@ function xout=erqa(varargin)
 % $Revision$
 %
 % $Log$
+% Revision 1.9  2004/11/10 07:06:03  marwan
+% initial import
+%
 %
 % This program is free software; you can redistribute it and/or
 % modify it under the terms of the GNU General Public License
@@ -198,7 +206,7 @@ end
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% compute
 
-try 
+%try 
 
 de=.1;
 if nogui~=2, h=waitbar(0,'1');h1=get(h,'chil');h1=get(h1,'title'); end
@@ -209,8 +217,17 @@ set(h1,'str',txt); drawnow, end
 
 try
   X1=crp(x(i:i+w-1,:),y(i:i+w-1,:),m,t,e,varargin{i_char},'silent');
+  if nogui~=2, txt=['CRP1a - ',num2str(i),'/',num2str(Nx-w)]; set(h1,'str',txt); drawnow, end
+  X1a=crp(x(i:i+w-1,:),y(i:i+w-1,:),m,t,e,varargin{i_char},'silent');
+  X2 = X1 .* X1a; X1 = X2;
   if nogui~=2, txt=['CRP2 - ',num2str(i),'/',num2str(Nx-w)]; set(h1,'str',txt); drawnow, end
-  X2=crp(x(i:i+w-1,:),y(i:i+w-1,:),m,t,e+de,varargin{i_char},'silent');
+ X2=crp(x(i:i+w-1,:),y(i:i+w-1,:),m,t,e+de,varargin{i_char},'silent');
+ X1a=crp(x(i:i+w-1,:),y(i:i+w-1,:),m,t,e+de,varargin{i_char},'silent');
+  X2a = X2 .* X1a; X2 = X2a;
+  if nogui~=2, txt=['CRP3 - ',num2str(i),'/',num2str(Nx-w)]; set(h1,'str',txt); drawnow, end
+  X3=crp(x(i:i+w-1,:),y(i:i+w-1,:),m+2,t,e,varargin{i_char},'silent');
+  X1a=crp(x(i:i+w-1,:),y(i:i+w-1,:),m+2,t,e,varargin{i_char},'silent');
+  X2a = X3 .* X1a; X3 = X2a;
 catch
   if nogui~=2, close(h), end 
   error(lasterr)
@@ -219,6 +236,7 @@ end
 txt=['ERQA - ',num2str(i),'/',num2str(Nx-w)];
 if nogui~=2, set(h1,'str',txt); drawnow, end
 warning off 
+
 N=length(X1);
 [ml NL]=dl(X1);
 pl=hist(NL,[0:N]);
@@ -233,8 +251,8 @@ i2=[round(.3*Np):round(.85*Np)]+1;
 %i1=2:round(.12*Np); 
 if length(i1)==1, i1=[]; end
 if length(i2)==1, i2=[]; end
-if ~isempty(i1),K1=regress(log(plc(i1))',[i1',ones(length(i1),1)]);else K1=NaN; end
-if ~isempty(i2),K2=regress(log(plc(i2))',[i2',ones(length(i2),1)]);else K2=NaN; end
+if ~isempty(i1),K1=-regress(log(plc(i1))',[i1',ones(length(i1),1)]);else K1=NaN; end
+if ~isempty(i2),K2=-regress(log(plc(i2))',[i2',ones(length(i2),1)]);else K2=NaN; end
 %K2=regress(log(plc(i2))',[i2',ones(length(i1),1)]);
 %K2=log(plc(i2));
 
@@ -257,6 +275,16 @@ end
 
 I2=2*H2-H2t;
 
+% ApEn
+B=(sum(X1,1)-1)/(length(X1)+1); A=(sum(X3,1)-1)/(length(X3)+1);
+
+ApEn = sum(log(B(find(B))))/(length(X1)+1) - sum(log(A(find(A))))/(length(X3)+1);
+
+%SampEn
+B=(sum(X1,1)-1)/(length(X1)-1); A=(sum(X3,1)-1)/(length(X1)-1);
+SampEn = -log(mean(A)/mean(B));
+
+
 warning on
 
 Y(i,1)=K1(1); 
@@ -265,6 +293,9 @@ Y(i,3)=D2;
 Y(i,4)=I2(1);
 Y(i,5)=I2(2);
 Y(i,6)=I2(3);
+Y(i,7)=ApEn;
+Y(i,8)=SampEn;
+
 end
 
 if nogui~=2
@@ -281,7 +312,8 @@ end
 
 if nargout, xout=Y; end
 
-catch
+if 0
+%catch
   if ishandle(h),close(h),end
   error(lasterr)
 end