function varargout=corrint(varargin) % % [y1,y2,y3]=corrint(x,embeddedDim,timeLag,timeStep,distanceThreshold,neighboorSize,estimationMode,findScaling) % % Correlation integral analysis of a time series. Based on: % % [1] Kaplan, Daniel, and Leon Glass. Understanding nonlinear dynamics. Vol. 19. Springer, 1995. % [2] Kantz, Holger, and Thomas Schreiber. Nonlinear time series analysis. Cambridge university press, 2004. % % Required input parameter: % x % Nx1 matrix (doubles) of time series to be analyzed. % % Optional Parameters are: % % embeddedDim % 1x1 Integer specifying the embedded dimension size to use (default % =2). % % timeLag % 1x1 Integer specifying the minimum time lag distance (in samples) of the point to % be estimated. Default is 2. If timeLag=-1 the timeLag is estimated % from the first zero-crossing point of the autocorrelation of x. % % timeStep % 1x1 Integer specifying time lag distance (in samples) within % each point used in the embeddedDimm vector. For example, if embeddedDim % is 3 and timeStep =2, then the embedded dimension vector will consists of % 3 samples separated by 2 samples each, covering a window of size of 7 samples. % % distanceThreshold % 1x1 double specifying the distance threshold between embedded % points. The points who's distance is less than distanceThreshold are considered in % the same neighborhood and used for either prediction, recurrence, or the % estimation of the embedded dimension. % % neighboorSize % 1x1 Integer specifying the number of neighbors to be used for % prediction and smoothing (see 'estimationMode' parameter). % % estimationMode % String specifying what analysis type to be done in the time series. % Options are: % 'recurrence' -Calculates recurrence data to be used % in for recurrence plots (default). % 'dimension' -Generates statistics for the estimation of the correlation dimension % of the time series and it's scaling % regions. % 'prediction' -Predicts second half of the time % series using the first half as a model % and neighboorSize nearest points. % 'smooth' -Predicts all point of the times % series using all other points as a % model and neighboorSize nearest % points. % % findScaling % 1x1 Boolean flag to be passed when using 'dimension' mode. If set to % true, the scaling region will be searched automatically, using % r1=std(x)/4 and r2 -> C(r1)/C(r2) ~ 5. Default value is false. % % The output returned by CORRINT is dependent on the 'estimationMode' % parameter, so that the description of the output below is broken down into the % different possible options for the 'estimationMode' parameter. % % Output Parameters - 'recurrence' mode % y1 % Lx1 Vector of integers for state i. % y2 % Lx1 Vector of integers for state state j that is a neighbor of state i (first column). % % Output Parameters - 'dimension' mode % y1 % Lx1 Vector of doubles of log(distanceThreshold). % y2 % Lx1 Vector of doubles for log(neighborhood size) given the distanceThreshold used in column 1. % y3 % 1x1 double. Optional, estimated slope of y1 and y2 % % Output Parameters - 'prediction' mode % y1 % Lx1 Vector of doubles of estimated second half of the time series. % y2 % Lx1 Vector of doubles for original second half of the time series. % y3 % 1x1 double. Optional, variance of the prediction error divided by variance of the second half of the time series. % % % Output Parameters - 'smooth' mode % y1 % Lx1 Vector of doubles of smoothed the time series. % y2 % Lx1 Vector of doubles for original time series. % y3 % 1x1 double. Optional, variance of the prediction error divided by variance of the time series. % % % %%% Beging Example %%% % % N=500; %Number of points for each process % model_names={'linearModel','nonlinearModel'}; % % %Linear Auto Regressive model with measurement noise % linearModel=zeros(N,1); % x=77; % linearModel(1)=x; % for n=2:N % x=4 + 0.95*x; % linearModel(n)= x + randn(1)*2; % end % % %Non-linear model of dimension ~ 3.9 % nonlinearModel=zeros(N,1); % x=0.2;y=0.2;z=0.2;v=0.2;model_five(1)=x; % for n=2:N % m=0.4 - 6/(1+ x^2 + y^2); % xold=x;yold=y;zold=z;vold=v; % x= 1 + 0.7*(xold*cos(m)-yold*sin(m)) + 0.2*zold; % y=0.7*(xold*sin(m) + yold*cos(m)); % z=1.4 + 0.3*vold - zold^2; % v=zold; % nonlinearModel(n)= x + 0.3*z + randn(1)*0.05; % end % % %Plot time series % figure(1) % for i=1:2 % subplot(2,1,i) % eval(['plot(' model_names{i} ');legend(''' model_names{i} ''')']) % title('Time Plot');xlabel('time') % end % % %Plot cross correlation % figure(2) % for i=1:2 % subplot(2,1,i) % eval(['x=' model_names{i} ';']) % R=xcorr(x-mean(x),'coeff'); % plot(R(round(N):end)) % eval(['legend(''' model_names{i} ''')']) % title('Autocorelation'); xlabel('lag') % end % % %Plot Phase Plots % figure(3) % for i=1:2 % subplot(2,1,i) % eval(['x=' model_names{i} ';']) % scatter(x(1:end-1),x(2:end)) % eval(['legend(''' model_names{i} ''')']) % title('Phase Plot');xlabel('x(t)');ylabel('x(t+1)') % end % % %Plot prediction errors vs surrogate % timeLag=1; % timeStep=1; % distanceThreshold=[]; % embeddedDim=4; % estimationMode='smooth'; % figure(4) % K=[1:20 25 30 50 70 100]; % D=length(K); % surrN=10; % for i=1:2 % eval(['x=' model_names{i} ';']) % err=zeros(D,1)+NaN; % surr_data=zeros(D,surrN); % SURR=surrogate(x,surrN); % for d=1:D; % neighboorSize=K(d); % [y1,y2,y3]=corrint(x,embeddedDim,timeLag,timeStep,distanceThreshold,neighboorSize,estimationMode); % err(d)=y3; % for s=1:surrN % [y1,y2,y3]=corrint(SURR(:,s),embeddedDim,timeLag,timeStep,distanceThreshold,neighboorSize,estimationMode); % surr_data(d,s)=y3; % end % end % subplot(2,1,i) % plot(K,err,'o-');hold on % errorbar(K,mean(surr_data,2),var(surr_data,[],2)./sqrt(10),'r') % eval(['legend(''' model_names{i} ''',''surrogate'')']) % xlabel('Embedded Dimension') % ylabel('err/var') % % end % % %%% End Example %%% % % Written by Ikaro Silva, 20134 % Last Modified: November 23, 2014 % Version 1.0 % % Since 0.9.8 % % % See also SURROGATE, DFA, MSENTROPY %endOfHelp persistent javaWfdbExec config if(isempty(javaWfdbExec)) [javaWfdbExec,config]=getWfdbClass('corrint'); end %Set default pararamter values inputs={'x','embeddedDim','timeLag','timeStep','distanceThreshold','neighboorSize','estimationMode','findScaling'}; outputs={'y1','y2','y3'}; embeddedDim=[]; timeLag=[]; timeStep=[]; distanceThreshold=[]; neighboorSize=[]; estimationMode='recurrence'; findScaling=0; wfdb_argument={}; y1=[]; y2=[]; y3=[]; for n=1:nargin if(~isempty(varargin{n})) eval([inputs{n} '=varargin{n};']) end end if(~isempty(embeddedDim)) wfdb_argument{end+1}='-d'; wfdb_argument{end+1}=num2str(embeddedDim); end if(~isempty(timeLag)) wfdb_argument{end+1}='-t'; wfdb_argument{end+1}=num2str(timeLag); end if(~isempty(timeStep)) wfdb_argument{end+1}='-s'; wfdb_argument{end+1}=num2str(timeStep); end if(~isempty(distanceThreshold)) wfdb_argument{end+1}='-r'; wfdb_argument{end+1}=num2str(distanceThreshold); end if(~isempty(neighboorSize)) wfdb_argument{end+1}='-n'; wfdb_argument{end+1}=num2str(neighboorSize); end switch estimationMode case 'recurrence' wfdb_argument{end+1}='-p'; case 'dimension' wfdb_argument{end+1}='-v'; y3=' '; if(findScaling) wfdb_argument{end+1}='-a'; end case 'prediction' wfdb_argument{end+1}='-P'; case 'smooth' wfdb_argument{end+1}='-S'; y3=' '; otherwise error(['Unkown estimation mode: ' estimationMode]) end javaWfdbExec.setArguments(wfdb_argument); if(config.inOctave) x=cellstr(num2str(x)); x=java2mat(javaWfdbExec.execWithStandardInput(x)); Nx=x.size; out=cell(Nx,1); for n=1:Nx out{n}=x.get(n-1); end else out=cell(javaWfdbExec.execWithStandardInput(x).toArray); end M=length(out); if(~isempty(y3)) y3=out{end}; out(end)=[]; M=M-1; if(strcmp(estimationMode,'smooth')) tmp=y3; sep=regexp(tmp,'\s'); y3=str2num(tmp(sep(end):end)); end end if(~isempty(strfind(out{1},'Possibly'))) warning(out{1}) out(1)=[]; M=M-1; end y1=zeros(M,1)+NaN; y2=zeros(M,1)+NaN; for m=1:M str=out{m}; sep=regexp(str,'\s'); y1(m)=str2num(str(1:sep)); y2(m)=str2num(str(sep(1):end)); end for n=1:nargout eval(['varargout{n}=y' num2str(n) ';']) end