crp2.m

function xout=crp2(varargin)
%CRP2   Creates a cross recurrence plot/ recurrence plot and the LOS.
%    CRP(X [,Y [,param1,param2,...]) creates a cross recurrence  
%    plot/ recurrence plot from the embedding vectors X and Y. 
%    Results can be stored into the workspace. Further it
%    is possible to estimate the line of synchronization (LOS)
%    in order to get the nonparametric time-relationship between
%    the two considered systems.
%
%    R=CRP(X,M2,T,E) uses the additionally dimension M2, delay T 
%    and the size of neighbourhood E and creates a recurrence 
%    plot of X.
%    
%    R=CRP(X,Y,'distance','nonormalize') creates a 
%    distance coded matrix plot without normalization
%    of the data.
%
%    Allows to change the parameters interactively by 
%    using a GUI.
%
%    The embedding dimension M is given by the size of the
%    N x M matrix X and Y; if the matrix Y is not specified, 
%    a simple recurrence plot is created.
%
%    Parameters: additionally dimension M2, delay T and the 
%    size of neighbourhood E are the first three numbers 
%    after the data series; further parameters can be used
%    to switch between various methods of finding the
%    neighbours of the phasespace trajectory, to suppress
%    the normalization of the data and to suppress the 
%    GUI (useful in order to use this programme by other 
%    programmes).
%
%    Methods of finding the neighbours.
%      maxnorm     - Maximum norm.
%      euclidean   - Euclidean norm.
%      minnorm     - Minimum norm.
%      nrmnorm     - Euclidean norm between normalized vectors
%                    (all vectors have the length one).
%      fan         - Fixed amount of nearest neighbours.
%      omatrix     - Order matrix (disabled).
%      opattern    - Order patterns recurrence plot.
%      distance    - Distance coded matrix (global CRP, Euclidean norm).
%
%    Normalization of the data series.
%      normalize   - Normalization of the data.
%      nonormalize - No normalization of the data.
%
%    Suppressing the GUI.
%      gui         - Creates the GUI and the output plot.
%      nogui       - Suppresses the GUI and the output plot.
%      silent      - Suppresses all output.
%
%    Parameters not needed to be specified.
%
%    Current limitation: for higher speed in
%    output the whole matrix of the recurrence
%    plot is in the work space - this limits
%    the application of long data series.
%
%    Examples: a=sin((1:1000)*2*pi/200);  % pendulum's location vector
%              b=cos((1:1000)*2*pi/200);  % pendulum's velocity vector
%              plot(a,b,'.')
%              crp2(a(1:500),b(1:500),'nonorm','euclidean')
%
%              b=sin(.01*([1:1000]*2*pi/67).^2);
%              crp2(b(1:500),a(1:700),3,10,.06,'fan')
%
%    See also CRP, CRP_BIG and TRACKPLOT.
%
%    References: 
%    Marwan, N., Thiel, M., Nowaczyk, N.: 
%    Cross Recurrence Plot Based Synchronization of Time Series,
%    Nonlin. Proc. Geophys., 9, 2002.

% Copyright (c) 1998-2004 by AMRON
% Norbert Marwan, Potsdam University, Germany
% http://www.agnld.uni-potsdam.de
%
% $Date$
% $Revision$
%
% $Log$
% Revision 4.7  2004/11/10 07:04:40  marwan
% initial import
%
%
% This program is part of the new generation XXII series.
%
% 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.



warning off
global errcode props

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% programme properties

errcode=0;
init_properties
hCRP=[];hCtrl=[];nogui=[];obj=[];mflag=[];
set(0,'ShowHidden','On')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% check and read the input
error(nargchk(1,9,nargin));
if nargout>1, error('Too many output arguments'), end

check_meth={'ma','eu','mi','nr','fa','in','om','op','di'}; 	% maxnorm, euclidean, nrmnorm,  fan, distance
check_norm={'non','nor'};				% nonormalize, normalize
check_gui={'gui','nog','sil'};				% gui, nogui, silent

if isnumeric(varargin{1}) 		% read commandline input
   varargin{9}=[];
   i_double=find(cellfun('isclass',varargin,'double'));
   i_char=find(cellfun('isclass',varargin,'char'));

   % check the text input parameters for method, gui and normalization
   temp_meth=0;
   temp_norm=0;
   temp_gui=0;
   if ~isempty(i_char)
      for i=1:length(i_char), 
         varargin{i_char(i)}(4)='0';
         temp_meth=temp_meth+strcmpi(varargin{i_char(i)}(1:2),check_meth'); 
         temp_norm=temp_norm+strcmpi(varargin{i_char(i)}(1:3),check_norm'); 
         temp_gui=temp_gui+strcmpi(varargin{i_char(i)}(1:3),check_gui'); 
      end
      method=min(find(temp_meth));
      nonorm=min(find(temp_norm))-1;
      nogui=min(find(temp_gui))-1;
      if isempty(method), method=1; end
      if isempty(nonorm), nonorm=1; end
      if isempty(nogui), nogui=0; end
      if method>length(check_meth), method=length(check_meth); end
      if nonorm>1, nonorm=1; end
      if nogui>2, nogui=2; end
   else
      method=1; nonorm=1; nogui=0;
   end
   if nogui==0 & nargout>0, nogui=1; end

   % get the parameters for creating RP
     if max(size(varargin{1}))<=3
        error('To less values in data X.')
     end
     x=double(varargin{1});
     if isempty(varargin{2}) | ~isnumeric(varargin{2}), y=x; else
     y=double(varargin{2}); end

     if (isnumeric(varargin{2}) & max(size(varargin{2}))==1) | ~isnumeric(varargin{2})
       y=x;
       if ~isempty(varargin{i_double(2)}), m0=varargin{i_double(2)}(1); else m0=1; end
       if ~isempty(varargin{i_double(3)}), t=varargin{i_double(3)}(1); else t=1; end
       if ~isempty(varargin{i_double(4)}), e=varargin{i_double(4)}(1); else e=.1; end
     else
       if ~isempty(varargin{i_double(3)}), m0=varargin{i_double(3)}(1); else m0=1; end
       if ~isempty(varargin{i_double(4)}), t=varargin{i_double(4)}(1); else t=1; end
       if ~isempty(varargin{i_double(5)}), e=varargin{i_double(5)}(1); else e=.1; end
     end
     t=round(t); m0=round(m0); mflag=method;
     if e<0, e=1; disp('Warning: The threshold size E cannot be negative and is now set to 1.'), end
     if t<1, t=1; disp('Warning: The delay T cannot be smaller than one and is now set to 1.'), end
     action='init';

    if size(x,1)==1, x=x'; end, if size(y,1)==1, y=y'; end 
    m=max([size(x,2) size(y,2)]);

  if ~isempty(find(isnan(x)))
     disp('NaN detected (in first variable) - will be cleared.')
     for k=1:size(x,2),  x(find(isnan(x(:,k))),:)=[]; end
  end
  if ~isempty(find(isnan(y)))
     disp('NaN detected (in second variable) - will be cleared.')
     for k=1:size(y,2),  y(find(isnan(y(:,k))),:)=[]; end
  end
  if size(x,1) < t*(m-1)+1 | size(y,1) < t*(m-1)+1
     error(['Too less data',10,...
            'Either too much NaN or the number of columns in the vectors do not match.'])
  end

    Nx=size(x,1); Ny=size(y,1);
    NX=Nx-t*(m-1);NY=Ny-t*(m-1);  
    x0=zeros(Nx,m);y0=zeros(Ny,m);
    x0(1:size(x,1),1:size(x,2))=x; 
    y0(1:size(y,1),1:size(y,2))=y; 

    if nonorm==1, 
	 x=(x0-repmat(mean(x0),Nx,1))./repmat(std(x0),Nx,1);
	 y=(y0-repmat(mean(y0),Ny,1))./repmat(std(y0),Ny,1);
    end

  if ~isempty(find(isnan(x))), for k=1:size(x,2),  x(find(isnan(x(:,k))),:)=[]; end, end
  if ~isempty(find(isnan(y))), for k=1:size(y,2),  y(find(isnan(y(:,k))),:)=[]; end, end
  if size(x,1) < t*(m-1)+1 | size(y,1) < t*(m-1)+1
     error(['Too less data',10,...
            'Either too much NaN or the number of columns in the vectors do not match.'])
  end