recurrence_plot.m

来自「主成分分析和偏最小二乘SquaresPrincipal成分分析（ PCA ）和偏」· M 代码 · 共 796 行 · 第 1/2 页

%    Examples: a=sqrt(100^2-(-71:71).^2); b=1:100;
%              b(101:100+length(a))=-(a)+170;
%              b(end+1:end+100)=100:-1:1;
%              crp_big(b,1,1,.1,'max')
%
%    See also CRP, CRP2 and CRQA.

% Copyright (c) 1998-2003 by AMRON
% Norbert Marwan, Potsdam University, Germany
% http://www.agnld.uni-potsdam.de
%
% Last Revision: 2004-01-23
% Version: 4.5.8
%
% 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.
%
% last modified 12.09.04 by Max Logunov

warning off

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

nogui=[];obj=[];mflag=[];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% check and read the input
error(nargchk(1,9,nargin));

if isnumeric(varargin{1})==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
   check_meth={'ma','eu','mi','nr','fa','in','di'}; 	% maxnorm, euclidean, nrmnorm,  fan, distance
   check_norm={'non','nor'};			% nonormalize, normalize
   check_gui={'gui','nog','sil'};		% gui, nogui, silent
   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>7, method=7; 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; end
%     if ~isempty(varargin{2}) & isnumeric(varargin{2}), y=double(varargin{2}); end
     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)}), m=varargin{i_double(2)}(1); else m=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)}), m=varargin{i_double(3)}(1); else m=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
     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
     t=round(t); m=round(m); mflag=method;
     action='init';

    Nx=length(x); Ny=length(y);
    NX=Nx-t*(m-1);NY=Ny-t*(m-1);  
    if (NX<1 | NY<1) & nogui==0;
         errordlg('The embedding vectors cannot be created. Dimension M and/ or delay T are to big. Please use smaller values.','Dimension/ delay to big')
         waitforbuttonpress
    end
  if size(x,1)<size(x,2), x=x'; end
  if size(y,1)<size(y,2), y=y'; end

  if ~isempty(find(isnan(x)))
     disp('Warning: 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('Warning: NaN detected (in second variable) - will be cleared.')
     for k=1:size(y,2),  y(find(isnan(y(:,k))),:)=[]; end
  end

  if size(x,2)>=2
     xscale=x(:,1); 
     if ~isempty(find(diff(xscale)<0)), error('First column of the first vector must be monotonically non-decreasing.'),end
     if nonorm==1, x=(x(:,2)-mean(x(:,2)))/std(x(:,2)); else x=x(:,2); end
  else
     if nonorm==1, x=(x-mean(x))/std(x); end
     xscale=(1:length(x))'; 
  end
  if size(y,2)>=2
     yscale=y(:,1); 
     if ~isempty(find(diff(yscale)<0)), error('First column of the second vector must be monotonically non-decreasing.'),end
     if nonorm==1, y=(y(:,2)-mean(y(:,2)))/std(y(:,2)); else y=y(:,2); end
  else
      if nonorm==1, y=(y-mean(y))/std(y); end
      yscale=(1:length(y))';
  end
  ds=eye(m);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% nogui

if nogui>0
   hCRP=9999;
   if nogui~=2 
       tx(1)={'Maximum norm'}; 
       tx(2)={'Euclidean norm'}; 
       tx(3)={'Minimum norm'}; 
       tx(4)={'Euclidean norm of normalized distance'}; 
       tx(5)={'fixed amount of nearest neighbours'};
       tx(6)={'interdependent neighbours'};
       tx(7)={'distance plot'};
       disp(['use method: ', char(tx(method))]);
       if nonorm==1, disp('normalize data'); else disp('do not normalize data'); end
   end
   action='compute';
   if (NX<1 | NY<1)
         disp('Warning: The embedding vectors cannot be created.')
	 disp('Dimension M and/ or delay T are to big. Please use smaller values.')
         action='';
   end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% switch routines

switch(action)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% initialization
  case 'init'

  errcode=1;
  xshuttle(:,1)=xscale;
  yshuttle(:,1)=yscale;
  xshuttle(:,2)=x;
  yshuttle(:,2)=y;
  ds=eye(m);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% computation
  
  case 'compute'

  errcode=11;
  txt_cross='Cross ';
  if size(x)==size(y) 
    if x==y, txt_cross=''; end
  end

  maxN=sqrt(m)*max(NX,NY);
  if maxN>800
    n=round(log(maxN/350)/log(2)); % orig war 150
    maxN=maxN/(2^n);
    maxN=round(maxN);NX0=NX;NY0=NY;smax=sqrt(m)*(max(abs(x))+max(abs(y)));
  else
    maxN=round(sqrt(m)*800);NX0=NX;NY0=NY;smax=sqrt(m)*(max(abs(x))+max(abs(y)));
  end
  s_norm=sqrt(m)*abs(max(x(:))-min(x(:)));

  X=uint8(zeros(NY+(m-1)*t,NX+(m-1)*t));
%   set(findobj('Tag','Status','Parent',findobj('Parent',hCRP,'Tag','CRPPlot')),'String','Build Embedding Vectors'),drawnow
    if m*Nx>maxN
      ax=ceil(m*Nx/maxN);
      Nx2=ceil(Nx/ax);
      x(Nx+1:Nx2*ax+t*(m-1))=x(Nx);
    else
      ax=1; Nx2=Nx;
    end
    x(Nx+1:Nx2*ax+t*(m-1))=x(Nx);
    x0=x;

    if m*Ny>maxN
      ay=ceil(m*Ny/maxN);
      Ny2=ceil(Ny/ay);
    else
      ay=1; Ny2=Ny;
    end
    y(Ny+1:Ny2*ay+t*(m-1))=y(Ny);
    y0=y;

%   set(findobj('Tag','Status','Parent',findobj('Parent',hCRP,'Tag','CRPPlot')),'String','Compute CRP'),drawnow
  k=0;
  for i=1:ax,
    
    x=x0(1+Nx2*(i-1):Nx2+Nx2*(i-1)+t*(m-1));
      for j=1:ay,
        y=y0(1+Ny2*(j-1):Ny2+Ny2*(j-1)+t*(m-1));
	k=k+1;

    Nx=length(x); Ny=length(y);
    NX=Nx-t*(m-1);NY=Ny-t*(m-1);  
%%

  ii=(1:NX)';jj=0:t:0+(m-1)*t;
  ii=reshape(ii(:,ones(1,m))+jj(ones(NX,1),:),m*NX,1);
  x1=x(ii);
  x2=reshape(x1,NX,m);

  ii=(1:NY)';jj=0:t:0+(m-1)*t;
  ii=reshape(ii(:,ones(1,m))+jj(ones(NY,1),:),m*NY,1);
  y1=y(ii);
  y2=reshape(y1,NY,m);
 
  y2=y2*ds; % switch vectors

  [NX, mx] = size(x2);
  [NY, my] = size(y2);
	
  clear jx jy

  switch(mflag)
  
  %%%%%%%%%%%%%%%%% local CRP, fixed distance
  
  case {1,2,3}

  errcode=111;
  px = permute(x2, [ 1 3 2 ]);
  py = permute(y2, [ 3 1 2 ]);
  s1 = px(:,ones(1,NY),:) - py(ones(1,NX),:,:);
  switch mflag
    case 1
  %%%%%%%%%%%%%%%%% maximum norm
    s = max(abs(s1),[],3);
      matext=[num2str(round(100*e)/100) '\sigma (fixed distance maximum norm)'];
    case 2
  %%%%%%%%%%%%%%%%% euclidean norm
      s = sqrt(sum(s1.^2, 3));
      matext=[num2str(round(100*e)/100) '\sigma (fixed distance euclidean norm)'];
    case 3
  %%%%%%%%%%%%%%%%% minimum norm
      s = sum(abs(s1), 3);
      matext=[num2str(round(100*e)/100) '\sigma (fixed distance minimum norm)'];
  end

  X1=255*s/max(s(:))<(255*e/max(s(:)));
  X0=(uint8(X1))'; clear s s1 x1 y1 px py X1
  X(1+Ny2*(j-1):Ny2+Ny2*(j-1),1+Nx2*(i-1):Nx2+Nx2*(i-1))=X0;
  X(NY0+1:end,:)=[];
  X(:,NX0+1:end)=[];
  NX=NX0; NY=NY0;
%  matext=[num2str(round(100*e)/100) '\sigma (fixed distance maximum norm)'];

  %%%%%%%%%%%%%%%%% local CRP, normalized distance euclidean norm
  
  case 4

%  set(findobj('Tag','Status','Parent',findobj('Parent',hCRP,'Tag','CRPPlot')),'String',['Normalize Embedding Vectors ',num2str(ax*(i-1)+j)]),drawnow
  Dx=sqrt(sum(((x2(:,:)).^2)'))';
  Dy=sqrt(sum(((y2(:,:)).^2)'))';
  x1=x2./repmat(Dx,1,m);x2=x1;
  y1=y2./repmat(Dy,1,m);y2=y1; clear Dx Dy y1 x1
  

  px = permute(x2, [ 1 3 2 ]);
  py = permute(y2, [ 3 1 2 ]);
  s1 = px(:,ones(1,NY),:) - py(ones(1,NX),:,:);
  s = sqrt(sum(s1.^2, 3));
  
  X0=(uint8(255*s/max(s(:)))<(255*e/max(s(:))))'; clear s s1 x1 y1 px py

  X(1+Ny2*(j-1):Ny2+Ny2*(j-1),1+Nx2*(i-1):Nx2+Nx2*(i-1))=X0;
  X(NY0+1:end,:)=[];
  X(:,NX0+1:end)=[];
  NX=NX0; NY=NY0;
  matext=[num2str(round(100*e)/100) '\sigma (normalized distance euclidean norm)'];


  %%%%%%%%%%%%%%%%% local CRP, fixed neigbours amount
  
  case 5
  
  if e>=1 
    e=round(e)/100;
    txt=['The value for fixed neigbours amount has to be smaller '...
         'than one. Continue the computation with a value of ' ...
	 num2str(e)];
    if nogui==0
      warndlg(txt,'Threshold value mismatch');
      drawnow
      waitforbuttonpress
      set(findobj('Tag','Size','Parent',gcf),'String',num2str(e))
    end
  end
  

  px = permute(x2, [ 1 3 2 ]);
  py = permute(y2, [ 3 1 2 ]);
  s1 = px(:,ones(1,NY),:) - py(ones(1,NX),:,:);
  s = sqrt(sum(s1.^2, 3));
  
  mine=round(NY*e);
  [SS, JJ]=sort(s'); JJ=JJ';

  X1(NX*NY)=uint8(0); X1(JJ(:,1:mine)+repmat([0:NY:NX*NY-1]',1,mine))=uint8(1);

  X0=reshape(X1,NY,NX); clear X1 SS JJ s px py
  X(1+Ny2*(j-1):Ny2+Ny2*(j-1),1+Nx2*(i-1):Nx2+Nx2*(i-1))=X0;
  X(NY0+1:end,:)=[];
  X(:,NX0+1:end)=[];
  matext=[num2str(round(1000*e)/10) '% (fixed neighbours amount)'];
  NX=NX0; NY=NY0;

  
  %%%%%%%%%%%%%%%%% local CRP, interdependent neigbours 
  
  case 6
  
  if e>=1 
    e=round(e)/100;
    txt=['The value for fixed neigbours amount has to be smaller '...
         'than one. Continue the computation with a value of ' ...
	 num2str(e)];
  end
  
  px = permute(x2, [ 1 3 2 ]);
  py = permute(y2, [ 1 3 2 ]);
  px2 = permute(x2, [ 3 1 2 ]);
  py2 = permute(y2, [ 3 1 2 ]);

  s1 = px(:,ones(1,NX),:) - px2(ones(1,NX),:,:);
  sx = sqrt(sum(s1.^2, 3));
  s1 = py(:,ones(1,NY),:) - py2(ones(1,NY),:,:);
  sy = sqrt(sum(s1.^2, 3));

  mine=round(min(NX,NY)*e);
  [SSx, JJx]=sort(sx);%SSx(1,:)=[]; JJx(1,:)=[];
  [SSy, JJy]=sort(sy);%SSy(1,:)=[]; JJy(1,:)=[];
  ee=mean(SSy(mine:mine+1,:));
  X0=uint8(zeros(NY,NX));
  for ii=1:min(NX,NY),
    JJx((JJx(1:mine,ii)>min(NX,NY)),ii)=1;
    X0(ii,JJx(1:mine,ii))=sy(ii,JJx(1:mine,ii))<=ee(ii);
  end

  X(1+Ny2*(j-1):Ny2+Ny2*(j-1),1+Nx2*(i-1):Nx2+Nx2*(i-1))=X0;
  X(NY0+1:end,:)=[];
  X(:,NX0+1:end)=[];
  NX=NX0; NY=NY0;
  clear X1 SS* JJ* s sx sy s1 px py 
  matext=[num2str(round(1000*e)/10) '% (interdependent neighbours)'];

  
  %%%%%%%%%%%%%%%%% global CRP
  
  case 7
  
  px = permute(x2, [ 1 3 2 ]);
  py = permute(y2, [ 3 1 2 ]);
  s1 = px(:,ones(1,NY),:) - py(ones(1,NX),:,:);
  s = sqrt(sum(s1.^2, 3))';
  X0=uint8(255*s/s_norm);  

  X(1+Ny2*(j-1):Ny2+Ny2*(j-1),1+Nx2*(i-1):Nx2+Nx2*(i-1))=X0;
  X(NY0+1:end,:)=[];
  X(:,NX0+1:end)=[];
  NX=NX0; NY=NY0;
  matext='';

  
  end
  
  
  end, end
  
  %%%%%%%%%%%%%%%%% show CRP 

 xout=X;
  
end
warning on