crqa.m

     if time_scale_flag
        x_var=var(x(i:i+w-1,2:end)); y_var=var(y(i:i+w-1,2:end)); 
        temp=cov(x(i:i+w-1,2:end),y(i:i+w-1,2:end));
     else
        x_var=var(x(i:i+w-1,1)); y_var=var(y(i:i+w-1,1)); 
        temp=cov(x(i:i+w-1,1),y(i:i+w-1,1));
     end
     xy_var=temp(1,2);

          
     do_norm = {'non';'nor'};

     % if plugin exist and method is MAX, MIN or EUC
     if plugin_exist & ( method_n < 4 | method_n == 9 ) & length(x) == length(y) 

         errcode=21;
         warning off
         while 1
             tmp_xdatafile = tempname;
             if ~exist(tmp_xdatafile), break, end
         end
         while 1
             tmp_ydatafile = tempname;
             if ~exist(tmp_ydatafile), break, end
         end
         while 1
             tmp_rqadatafile = tempname;
             if ~exist(tmp_rqadatafile), break, end
         end
         
         
         x_tmp = x(i:i+w-1,:);
         y_tmp = y(i:i+w-1,:);
         
         % save data in temporary file
         save(tmp_xdatafile,'x_tmp','-ascii','-tabs');
         save(tmp_ydatafile,'y_tmp','-ascii','-tabs');

         % call extern rp programme
         m_str = {'MAX', 'EUC', 'MIN', 'NR', 'RR', 'FAN', 'IN', 'OM', 'OP', 'EUC'};

         [status ] = system([plugin_path,filesep,plugin_name,' -m ',num2str(m), ...
                                       ' -t ',num2str(t), ...
                                       ' -e ',num2str(e), ...
                                       ' -n ',m_str{method_n}, ...
                                       ' -w ',num2str(theiler_window), ...
                                       ' -l ',num2str(lmin), ...
                                       ' -v ',num2str(vmin), ...
                                       ' -i ',tmp_xdatafile, ...
                                       ' -j ',tmp_ydatafile, ...
                                       ' -o ',tmp_rqadatafile, ...
                                       ' -s']);
         errcode=22;
         % import RQA
         rqa_in = [];
         try
             fid = fopen(tmp_rqadatafile,'r'); % open RQA data file
             while 1
                 dataStr = fgetl(fid); % read RQA data
                 if ~ischar(dataStr), break, end % leave the loop if end of file
                 if isempty(findstr(dataStr, '#')) % neglect comments line (e.g. header)
                     rqa_in = str2num(dataStr); % import RQA measures into local variable rqa
                 end
             end
             fclose(fid); % close RQA data file

             RR = rqa_in(1);
             DET = rqa_in(2);
             LAM = rqa_in(4);
             Lmax = rqa_in(6);
             L = rqa_in(7);
             ENTR = rqa_in(8);
             Vmax = rqa_in(10);
             TT = rqa_in(11);
             t1 = rqa_in(13);
             t2 = rqa_in(14);
             RT = rqa_in(16);
             RTmax = rqa_in(15);
             RF = rqa_in(19);
             ENTW = rqa_in(17);
             warning off
             delete(tmp_rqadatafile);
             delete(tmp_xdatafile);
             delete(tmp_ydatafile);
             warning on
         catch 
             warning off
             delete(tmp_rqadatafile);
             delete(tmp_xdatafile);
             delete(tmp_ydatafile);
             warning on
         end

         if nogui~=2 & ishandle(h1), set(h1,'str',[num2str(i),'/',num2str(Nx-w)]); waitbar(i/(Nx-w)); drawnow, end



     % use builtin implementation
     else
     
         errcode=25;
         try
         %  X=crp_big(x(i:i+w-1,:),y(i:i+w-1,:),m,t,e,'fan','silent');
             if time_scale_flag
               if length(x(i:i+w-1,:)) > 2000
                  X=crp_big(x(i:i+w-1,:),y(i:i+w-1,:),m,t,e,method,'nonorm','silent');
               else
                  X=crp(x(i:i+w-1,:),y(i:i+w-1,:),m,t,e,method,'nonorm','silent');
               end
             else
               X=crp2(x(i:i+w-1,:),y(i:i+w-1,:),m,t,e,method,'nonorm','silent');
             end
           %  X=crp(x(i:i+w-1,:),y(i:i+w-1,:),m,t,e,varargin{i_char},'silent');
           warning off 
           if nogui~=2 & ishandle(h1), set(h1,'str',[num2str(i),'/',num2str(Nx-w)]); waitbar(i/(Nx-w)); drawnow, end

         %if 0
         catch
             error(lasterr)
             if nogui~=2 & ishandle(hw), close(hw), end
         end

         N=size(X);
         if theiler_window > 0
             X_theiler=double(triu(X,theiler_window))+double(tril(X,-theiler_window));
         else
             X_theiler=X;
         end
         errcode=26;

         % compute recurrence times of 1st and 2nd type
         if size(X_theiler,2) > 1000
             t1=[];t2=[];
             rps2=find(diff(double(X_theiler(:)))==1);
             rps=find(X_theiler(:));
             t1=diff(rps);
             t2=diff(rps2);
         else
             t1 = []; t2 = [];
             for i2=1:size(X_theiler,2)
                 if nogui~=2 & (Nx-w < 2), waitbar(i2/size(X_theiler,2)), end
                 rps2=find(diff(double(X_theiler(:,i2)))==1);
                 rps=find(X_theiler(:,i2));
                 t1=[t1;diff(rps)];
                 t2=[t2;diff(rps2)];
             end
         end
         t1=mean(t1);
         t2=mean(t2);

         errcode=27;
         [a b]=dl(X_theiler);

         warning off 
         errcode=271;
         b(find(b<lmin))=[];
         [c d dw]=tt(X_theiler);
         warning off 
         errcode=272;
         d(find(d<vmin))=[];

         errcode=273;
         RR=sum(X_theiler(:))/(N(1)*N(2));
         %b(find(b>=max(N)-lmin))=[]; if isempty(b), b=0; end
         if isempty(b), b=0; end
         errcode=274;
         if sum(sum(X_theiler)) > 0
           DET=sum(b)/sum(sum(X_theiler));
         else
           DET=NaN;
         end
         errcode=275;
         L=mean(b);
         histL=hist(b(:),[1:min(N)]);
         ENTR=entropy(histL(:));
         errcode=276;
         if sum(d)>0
           LAM=sum(d)/sum(sum(X_theiler));
         else
           LAM=NaN;
         end
       
         % recurrence times
         RT = mean(dw);
         [dwh dwi] = hist(dw,[1:max(dw)]);
         if dwh
             [dws dwsi] = sort(dwh);
             RTp = dwi(dwsi(end)); % most probable recurrence time
         else
             RTp = 0;
         end
         RTmax = max(dw); % maximal recurrence time
         RF = 1/RTmax; % minimal recurrence frequency
         ENTW = entropy(dwh(:));

         errcode=277;
         TT=mean(d);
         b=[b;0]; Lmax=max(b);
         d=[d;0]; Vmax=max(d);
         
     end % end plugin
  
     warning on

     errcode=28;
     Y(i,1)=RR; 
     Y(i,2)=DET;
     Y(i,3)=L;
     Y(i,4)=Lmax;
     Y(i,5)=ENTR;
     Y(i,6)=LAM;
     Y(i,7)=TT;
     Y(i,8)=Vmax;
     Y(i,9)=t1;
     Y(i,10)=t2;
     Y(i,11)=RT;
     Y(i,12)=RTmax;
     Y(i,13)=RF;
     Y(i,14)=ENTW;
     Y(i,15)=x_var;
     Y(i,16)=y_var;
     Y(i,17)=xy_var;
end % end window loop
if ishandle(hw), waitbar(1); drawnow; close(hw), end % close waitbar

  if ~nogui
    h=findobj('tag','crqa_button_apply');
    set(h(1),'ToolTip','Starts the computation.','String','Apply','Callback','crqa compute')
    for j=1:length(obj); 
        h=findobj('Tag',obj{j},'Parent',h_fig(1)); 
        if ~isempty(h)
          set(h,'Enable','On')
        end
    end
  end

  if ~nogui, 
    h=findobj('Tag','crqa_Fig'); 
    if ~isempty(h), set(0,'CurrentFigure',h(1)); end
    setptr(gcf,'arrow'), 
  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% plot
if ~nogui
    set(0,'showhidden','on')
    errcode=30;
    h=findobj('Tag','crqa_Fig'); if ~isempty(h), set(0,'CurrentFigure',h(1)); end
    tx={'RR';'DET';'L';'ENTR';'LAM';'TT';'T_1';'T_2';'Variance';'Covariance'};
    index=[1,2,3,5,6,7,9,10,11,12,13];
    tags={'crqa_axes_RR','crqa_axes_DET','crqa_axes_L','crqa_axes_ENTR','crqa_axes_LAM','crqa_axes_TT','crqa_axes_T1','crqa_axes_T2','crqa_axes_Var','crqa_axes_CoVar'};
    h=findobj('Tag','crqa_axes_RR','Parent',gcf); h_axes.h(1)=h(1);
    h=findobj('Tag','crqa_axes_DET','Parent',gcf); h_axes.h(2)=h(1);
    h=findobj('Tag','crqa_axes_L','Parent',gcf); h_axes.h(3)=h(1);
    h=findobj('Tag','crqa_axes_ENTR','Parent',gcf); h_axes.h(4)=h(1);
    h=findobj('Tag','crqa_axes_LAM','Parent',gcf); h_axes.h(5)=h(1);
    h=findobj('Tag','crqa_axes_TT','Parent',gcf); h_axes.h(6)=h(1);
    h=findobj('Tag','crqa_axes_T1','Parent',gcf); h_axes.h(7)=h(1);
    h=findobj('Tag','crqa_axes_T2','Parent',gcf); h_axes.h(8)=h(1);
    h=findobj('Tag','crqa_axes_Var','Parent',gcf); h_axes.h(9)=h(1);
    if ~all(x(:)==y(:)), h=findobj('Tag','crqa_axes_CoVar','Parent',gcf); h_axes.h(10)=h(1); end
   for i=1:9,
     set(gcf,'CurrentAxes',h_axes.h(i))
     if size(Y,1)==1
       cla
       text(0.5,0.5,sprintf('%6.4f',Y(index(i))),'FontWeight','bold','HorizontalAlign','Center')
%        if ~all(x(:)==y(:)) & i==9
%          cla 
%              text(0.5,0.5,sprintf('%6.5f, %6.5f',Y(index(i)),Y(index(i+1))),'FontWeight','bold','HorizontalAlign','Center')
%              set(gcf,'CurrentAxes',h_axes.h(i+1));cla
%          set(h_axes.h(i+1),'visible','off');
%        end
     else
       if i==9
         cla
             h2=stairs(xscale(1:wstep:length(Y)),Y(1:wstep:end,index(i)));
             set(h2,'color',props.line.Color)
         if ~all(x(:)==y(:))
           set(gca,'Color','none')
           hold on; h2=stairs(xscale(1:wstep:length(Y)),Y(1:wstep:end,index(10)),'r');
           set(gcf,'CurrentAxes',h_axes.h(i+1));cla
           set(h_axes.h(10),'visible','on');
           
           h3=stairs(xscale(1:wstep:length(Y)),Y(1:wstep:end,index(11))); set(h3,'color',[0 .4 0]);
           ylabel(tx(i+1));
           set(h_axes.h(i+1),'Tag',tags{i+1},'Units','Norm','Color','none','YAxisLocation','right','YColor',[0 .4 0])
         end
       else
         plot(xscale(1:wstep:length(Y)),Y(1:wstep:end,index(i)),'color',props.line.Color)
       end
     end
     set(gcf,'CurrentAxes',h_axes.h(i));
     ylabel(tx(i));
     set(gca,'Tag',tags{i},'color',props.axes.Color,'Units','Norm')
   end
   h=findobj('Tag','crqa_button_store');
   set(h(1),'UserData',Y,'Enable','On')
   xlim = get(h_axes.h(9),'xlim');
   if length(Y(:,1)) > 1, set(h_crqa_axes_Data,'xlim',xlim), end
   if ~all(x(:)==y(:)) & length(Y(:,1)) == 1, set(h_axes.h(10),'visible','off'); delete(get(h_axes.h(10),'children')), end
else
  if nargout, xout=Y(:,1:14); end
end

if nargout, xout=Y(:,1:14); end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% the end 

end




%if 0
catch
  if nogui~=2, if ishandle(hw), close(hw), end , end
      
  z=whos;x=lasterr;y=lastwarn;in=varargin{1};
  print_error('crqa',z,x,y,in,method,action)
  try, if ~nogui
    h=findobj('tag','crqa_button_apply');
    set(h(1),'ToolTip','Starts the computation.','String','Apply','Callback','crqa compute')
    for j=1:length(obj); 
        h=findobj('Tag',obj{j},'Parent',h_fig(1)); 
        if ~isempty(h)
          set(h,'Enable','On')
        end
    end, end
  end
  set(0,'showhidden','on')
    h=findobj('Tag','crqa_Fig'); if ~isempty(h), setptr(h(1),'arrow'); end
  set(0,'showhidden','off')
end


try, set(0,props.root), end



try
  set(0,'ShowHidden','off')
end