diff --git a/rpde.m b/rpde.m
new file mode 100644
index 0000000000000000000000000000000000000000..765b73e33a90c9bc8e9e994a415e7cc5415681be
--- /dev/null
+++ b/rpde.m
@@ -0,0 +1,197 @@
+function out=rpde(varargin)
+% RPDE Computes the recurrence period density entropy.
+% Y=RPDE(X [,Y] [,param1,param2,...])
+% Calculates the normalised entropy Y of the
+% recurrence time distribution of time series X.
+%
+% Examples: a = sin(0:.1:80);
+% b = sin(0:.1:80) + .1 * randn(1,801);
+% rpde(a,3,15,.1)
+% rpde(b,3,15,.1)
+%
+% See also CRQAD, TT.
+%
+% References:
+% Little, M., McSharry, P., Roberts, S., Costello, D., Moroz, I.:
+% Exploiting Nonlinear Recurrence and Fractal Scaling Properties
+% for Voice Disorder Detection, Biomed. Eng. Online, 6, 2007.
+
+% Copyright (c) 2010-
+% Norbert Marwan, Potsdam Institute for Climate Impact Research, Germany
+% http://www.pik-potsdam.de
+%
+% $Date$
+% $Revision$
+%
+% $Log$
+%
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% programme properties
+
+global props
+
+init_properties
+
+lmin=1;
+w=[]; method='max'; method_n=1; t=1; m=1; e=.1;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% check the input
+
+error(nargchk(1,10,nargin));
+if nargout>1, error('Too many output arguments'), end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% check and read the input
+
+ varargin{11}=[];
+ % transform any int to double
+ intclasses = {'uint8';'uint16';'uint32';'uint64';'int8';'int16';'int32';'int64'};
+ flagClass = [];
+ for i = 1:length(intclasses)
+ i_int=find(cellfun('isclass',varargin,intclasses{i}));
+ if ~isempty(i_int)
+ for j = 1:length(i_int)
+ varargin{i_int(j)} = double(varargin{i_int(j)});
+ end
+ flagClass = [flagClass; i_int(:)];
+ end
+ end
+ if ~isempty(flagClass)
+ disp(['Warning: Input arguments at position [',num2str(flagClass'),'] contain integer values']);
+ disp(['(now converted to double).'])
+ end
+ i_double=find(cellfun('isclass',varargin,'double'));
+ i_char=find(cellfun('isclass',varargin,'char'));
+ nogui=0;
+
+ if nargin & isnumeric(varargin{1})
+
+ % check the text input parameters for method, gui
+ check_meth={'ma','eu','mi','nr','rr','fa','in','om','op','di'}; % maxnorm, euclidean, nrmnorm, fan, distance
+ check_gui={'gui','nog','sil'}; % gui, nogui, silent
+ temp_meth=0;
+ temp_gui=0;
+ if ~isempty(i_char)
+ for i=1:length(i_char),
+ varargin{i_char(i)}(4)='0';
+ temp_gui=temp_gui+strcmpi(varargin{i_char(i)}(1:3),check_gui');
+ temp_meth=temp_meth+strcmpi(varargin{i_char(i)}(1:2),check_meth');
+ end
+ method_n=min(find(temp_meth));
+ nogui=min(find(temp_gui))-1;
+ for i=1:length(i_char); temp2(i,:)=varargin{i_char(i)}(1:3); end
+ i_char(strmatch(check_gui(find(temp_gui)),temp2))=[];
+ if isempty(nogui), nogui=0; end
+ if isempty(method_n), method_n=1; end
+ if nogui>2, nogui=1; end
+ if method_n>length(check_meth), method0=length(check_meth); end
+ method=check_meth{method_n};
+ else
+ nogui=0;
+ if nargout
+ nogui=1;
+ action='compute';
+ end
+ end
+ if nogui==0
+ action='init';
+ else
+ action='compute';
+ 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 sum(double(diff(x(:,1))<=0)), embed_flag=0; 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
+ if ~isempty(varargin{i_double(5)}), w=varargin{i_double(5)}(1); else w=varargin{i_double(5)}; end
+% if ~isempty(varargin{i_double(6)}), wstep=varargin{i_double(6)}(1); else wstep=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
+ if ~isempty(varargin{i_double(6)}), w=varargin{i_double(6)}(1); else w=varargin{i_double(6)}; end
+% if ~isempty(varargin{i_double(7)}), wstep=varargin{i_double(7)}(1); else wstep=1; end
+ end
+ else
+ error('No valid arguments.')
+ end
+
+
+ Nx=length(x); Ny=length(y);
+ if size(x,1)<size(x,2), x=x'; end
+ if size(y,1)<size(y,2), y=y'; end
+
+ if size(x,2)>=2
+ xscale=x(:,1);
+ if ~isempty(find(diff(xscale)<0)), embed_flag=0;end
+ else
+ xscale=(1:length(x))';
+ end
+ if size(y,2)>=2
+ yscale=y(:,1);
+ if ~isempty(find(diff(yscale)<0)), embed_flag=0;end
+ else
+ yscale=(1:length(y))';
+ end
+
+ if max(size(x))~=max(size(y)),
+ if ~nogui, errordlg('Data must have the same length.','Check Data'), else error('Data must have the same length.'), end
+ end
+ if e<0,
+ e=1;
+ if ~nogui
+ warndlg('The threshold size E can not be negative and is now set to 1.','Check Data')
+ h=findobj('Tag','crqa_eps');
+ set(h(1),'String',str2num(e))
+ else
+ disp('The threshold size E can not be negative and is now set to 1.'),
+ end
+ end
+ if t<1,
+ t=1;
+ if ~nogui
+ warndlg('The delay T can not be smaller than one and is now set to 1.','Check Data')
+ h=findobj('Tag','crqa_maxLag');
+ set(h(1),'String',str2num(t))
+ else
+ disp('The delay T can not be smaller than one and is now set to 1.')
+ end
+ end
+ if isempty(w), w=.5*Nx; wstep=1; end
+% if w<2,
+% w=2;
+% if ~nogui, warndlg('The window size W exceeds the valid range.','Check Data')
+% else, disp('The window size W exceeds the valid range.'), end
+% end
+ if w>Nx,
+ w=Nx; wstep=1;;
+ if ~nogui, warndlg('The window size W exceeds the valid range.','Check Data')
+ else, disp('The window size W exceeds the valid range.'), end
+ end
+ t=round(t); m=round(m); w=round(w);% wstep=round(wstep);
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% compute
+
+% calculate recurrence plot (RP)
+X = crp(x,m,t,e,method,'sil');
+
+% get recurrence times (in terms of vertical white lines in the RP)
+[dummy1 dummy2 w] = tt(X);
+w(w<=lmin) = [];
+
+% get histogram (normalisation will be done in ENTROPY)
+h = histc(w,[0:max(w)]+.5);
+
+% calculate entropy
+out = entropy(h) / log(max(w));
+