High dimensional chaos (HDS)


Welcome to my web site concerning the development of the algorithms estimating correlation dimension d of the continuous high-dimensional chaotic signals.

The page makes it posibble to download free the MATLAB HDS-toolkit that make possible to estimate the Dimensional Complexity  (d) of high dimensional signals.

The HDS-toolkit can be free used and distributed both for academic and comercial use. The only request is to cite 3 papers below, describing the method and algorithm:

The details of the algorithm are presented in my papers:

1.      K. Michalak, (2011) Modifications of the Takens-Ellner algorithm for medium- and high-dimensional signals, Physical Review E, Vol. 83(2):026206,

2.     K.P. Michalak. How to estimate the correlation dimension of high-dimensional signals?  AIP Chaos, 24(2):033118, 2014.

3.  K.P. Michalak. Estimating correlation dimension of high-dimensional signals - quick algorithm. in preparation, sent to a journal

The notation being used in the functions dhds.m and dhdsall.m  corresponds to the notation being used in these articles.



This is the beta version of the toolkit.

If something does not go, please write me the email with the description of your problem. I try help.


 There are 2 functions dedicated for estimating correlation dimension d

dhds.m - the short version of algorithm is implemented


dhdsall.m - the full version of algorithm is implemented


The general call of dhds.m function is:

[ dcopt, Wopt, dcW, dcM, dcMc, WV, Pmax, other ] =

dhds (series, quick precise, MaxIterationTime, Da,  DeflErr,  norm,  NW Wmax WV,  Pmin,  dcW,  fitting,  figures)


However, you can use the simplest call:   d = dhds ( series )


Use empty string:  ''    if you want use the default values of subsequent input parameters.

   Use the function surrog2.m to create the surrogates.


Let us analyze the input parameters:
  (1)  series - the signal being analyzed

(2) quick precise: this parameter makes it possible to set just one of 3 options of speed/precision:

  1 or 'quick' - quick, low precision,

  2 or 'medium' - medium (default),

  3 or 'precise' - slow, high precision.

Use 1 for preliminary estimation and 3 for nal estimation. Change other parameters if the results are not satisfactory.

(3) MaxIterationTime: this parameter makes it possible to schedule the calculations with respect to its duration. It denotes the maximal duration (in seconds) of every W-iteration. There is more di cult to nd distances that are shorter than assumed max when the estimated complexity for the given W increases and the calculations become signi cantly longer. The solution is to increase the accepted de ection error e which determines higher max. The estimation of d becomes possible, however, it is performed with a bit lower precision due to higher e. After MaxIterationTime is reached without nding Da distances, the algorithm takes decision: (a) if more than 50% of Da is reached until now, the calculation follows up to reaching Da, however, e is increased for next iteration; if less than 50% distances is reached, W-iteration is repeated with higher e. e increases by 1% by default and cannot excess 30%. The default MaxIterationTime is 60 s, 200 s or 1000 s depending on quick precise value. It is possible, also, to change the HEADER's variable 'UserAction'. If 'UserAction'=0 than the calculations are repeated using increased De Err. If 'UserAction'=1, the algorithm stops and asks for the user's decision: 'Follow' using current e or 'Break and repeat using increased e'.

(4) Da: assumed nbr of distances being shorter than the estimated distance dmax. dmax is estimated based on assumed deflection error e and used embedding dimension m. When omitted, default value of Da depends on the used quick_precise parameter: Da = 2000 or 5000 or 10000.

(5) DeflErr(1:2): the vector with 2 elements. The rst element describes the assumed deflection error e for the first iterations of calculations to estimate current
dmax. The default value is 3%. It is increased automatically when the complexity for the given W increases and it is difficult to find assumed number of distances Da in the assumed MaxIterationTime. Increasing e reduces the precision of calculations. It must be pointed, however, that to small e may shift the results into `noise` region of d=fn(Pmax) relation. You can use higher values initial values e.g. 5-10% if the expected complexity d is very high and/or the signal is short and/or you want obtain robust results in very quick time. Use lower values if the complexity of the signal is rather low, the signal is noise-free and you would like reduce the error connected with e-correction. High values accelerate strongly the calculations, especially if dw is high, however, 'e-correction error' occurs.
The second element of this input parameter De Err(2) = ecut represents the deflection error for which the
dcut is estimated. It must be equal or greater than e. dcut determines the cut-off distance and thus, the length of the actual dmax vector in the memory space and the duration of its sorting, after Da is reached. The distances shorter than cut and longer than max are the reserve if increase in e is necessary for the given W-iteration. The iteration has not to be repeated. It influences also the graphs being printed as a output of the function. If you use higher value, the relation d=fn(Pmax) is broader. The acceptable values lies between 0 and 30. The unit is [%]. By default, it is 4%. Both these parameters are automatically increased if the assumed number of distances Da shorter than max is not reached in the assumed iteration time.
The calculation for the given W-iteration is repeated using DeflErr increased by 1 making the estimation of dw possible, however, with a bit lower precision.

(6) normalization : preliminary normalisation of the signal. Acceptable values are 0 and 1. If 1 (default) - the series is converted to possess a normal distribution (0,1). Normalisation can reduce the calculation time, because the higher ratio of distances is shorter than assumed
dmax and Pmax is higher for the assumed e. Be careful if the histogram of the signal is strongly non-gaussian, and especially, if the histogram possesses the maximum/a on its borders. Normalisation makes it also possible to use preliminary estimated histogram slopes (stored in attached deltamaxM_PmaxM_normal.mat file). They have not to be estimated individually for the given signal.

(7) NW_Wmax_WV: Here, you can set the vector WV with your Window Width values used for calculations. You can just set the vector WV with the assumed W values. E.g. you want to use the same W values for different similar signals.

If the number of elements is less than 4 the algorithm finds default values depending on the number of elements:
- 1-element number : it is interpreted as NW
- 2-element vector it is interpreted as: [NW, Wmax]
- 3-element vector: it is interpreted as: [NW, Wmax, Wmaxdensity]

NW is the number of elements in WV vector (default is 40), Wmax is the maximal value in WV (default is 10 (autocorelation time of the signal))

Wmax - maximal value of the vector WV, by default 10t.
  Wmaxdensity denotes the value of W for which W values are chosen more dense.

If you expect the range of plateau in dw=fn(W) then set Wmaxdensity to be the center of the plateau in order to increase a bit the precision of polynomial fitting (default is 2.5 ).
- if [0,Wmax] - dw's are calculated for increasing W's calculated using a rule:
Wnext = Wprevious * Wincr... up to Wmax value starting from W=5samples.
Use this option as the most universal and independent on . Wincr is equal to 1.1-1.3 depending on quick precise value.
- if 'fntau' or 'tau' or ": (DEFAULT) sets 40 values W from 0.1 to 10 due to estimated autocorrelation time of the signal, unequally distributed with the maximum density for about 2.5 .
The distribution of default W-values can be modi ed in this method by changing the default parameters in the HEADER of the function: NW=40,  Wmax_over_tau=10 and defaultWdense_over_tau=2.5.

(8) Pmin: This parameter determines the method of the choice of Pmin parameter. This parameter does not influence strongly the result of final d estimation if Pmax is high enough. As a rule, it is necessary to exclude the initial part of the vector with the very short distances which do not create a regular plateau. In the case of noise contaminated signal, Pmin represents the range of the strong influence of noise.

The interpretation of this parameter depends on its value:

If 0 < Pmin < 0.1 - than it is treated as Pmin. You can set this value if you have observed the begin of the in uence of noise on the graph d=fn(Pmax,W) in the preliminary run of the algorithm. 
  If Pmin > 0.1 - it determines the deflection error (in %) for which Pmin should be estimated.
  If Pmin=0  - than no Pmin-cut-off is used and d is estimated using all the distances shorter than
dmax (Imin=1). Use this value if you want see on the graph all the relation d=fn(Pmax,W) with the noise contamination region.
  If Pmin=-1 - The default value for Pmin. It determines the automatic approximate determination of Pmin, separately for every W iteration. Pmin is chosen to be the smallest Pmax for which the estimated d lies in the range 66-150% of the preliminary/roughly estimated d. Look for the variable PminCut=1.5 in the HEADER of the function to modify this range.

(9) dcWmethod : This parameter determines the method of the estimation dw for the current W iteration from the relation d=fn(Pmax).
The question arises, which fragment of this relation should be the base for the estimation of dw?

It is observed that the relation d=fn(Pmax) possesses lacunarities in the plateau region. Thus, in order to increase the precision, a number of points from this relation should be taken to estimate dw as a mean over these points. The parameter dcW determines the the choice of Imaxdown Imaxup

range for these points.
1 - range is maximally broad: Imaxdown = Ioptmin+200,   Imaxup = Icut,
2 - (default) range is broad: Imaxdown = 0.05*Imaxopt ,   Imaxup = 1.5*Imaxopt
3 - range is medium broad: Imaxdown = 0.12*Imaxopt ,   Imaxup = 1.2*Imaxopt
4- range is narrow:               Imaxdown = 0.3*Imaxopt ,     Imaxup =    1 *Imaxopt
5- use 1 point only: Imax = Imaxopt

(10) fitting_method: This parameter determines the method of polynomial fitting applied to the relation d=fn(W).
if 0:   look for optimal fitting over different polynomials and self-defined functions in the function 'optimal_fitting.m`.

Here, you can define other functions than polynomials to find optimal fitting. The procedure works on symbolic expressions and solving differential equations, thus be careful with defining variables being optimized.
if [4,5,6] (or higher): de ne the degrees of the polynomials to be fitted.
if [4]: there is only 1 possible minimum in the 3-degree 1st derivative. You reach
unambiguous result: 0 or 1 plateau and dopt estimated.
if [5] or [6]: there are 0,1 or 2 possible minima in the 4th- and 5th-degree 1st derivatives.
if [7] or [8]: there are 0,1,2 or 3 possible minima in the 6th- and 7th-degree 1st derivatives.
if [4,5] [4 6] [4 5 9] [4 6 8 10]: you can use every combination of the degrees to be used. The graphs and results will be printed for every polynomial used. Increase NW when using high degree polynomials to increase fitting precision.
DEFAULT is [4 6 8].

(10) figures :
0 - turns off plotting figures
1 - figures are plotted once after the last W-iteration.
2 - turns on plotting intermediate figures after every W-iteration.


(12) signal_numer : important only for figure names. This number is added to figure names to distinguish figure names while calling to function in the loop for subsequent signals. Other case - figure are over-saved for subsequent calls of the function.


(13) fighandleset: vector with fighandles to use while printing figures.

when ommitted or 0 or -1 or NaN - new figures are created.

(1:3) figures refreshed after every iteration

(4:5) figures plotted while using 'optimal_fitting'procedure

(6:7) first polynomial degree, (8,9) - second degree,..





I. dopt {i,j}(k) : final estimators of dimensional complexity
i=1 - results without e-correction, i=2 - results with e-correction
j - the number of the fitted polynomial, index into ` fitting` vector
k - the number of plateaus for the given polynomial fitting. In case of higher degree polynomials, the number of the plateaus may be higher than 1.
II. Wopt {i,j}(k) : W's corresponding to dopt.   i, j, k as above.
III. dcW(1; :) - vector with dw's estimated for subsequent W's without e-correction.
      dcW(2; :) - with e-correction.
IV. dcM - matrix (sizePmaxV, sizeWV) with d estimated without e-correction for the individual pairs Pmax, W.

Used for printing d=fn(Pmax,W) relation.
V. dcM_c - the same but with e-correction
VI. WV - vector with W's for which dw were estimated to build the relation dw = fn(W)
VII. PmaxV - defined in the HEADER of the function as: PmaxV_default = 10(8:0:01:0:01). Necessary for printing graphs.
VIII. other - structure keeping many other intermediate parameters. For details - see file dhds.m, lines 1081-1109.


Exemplary use:

Analyzed signal - is the sum of 2 Lorenz signals possessing different time scales ratio 1:5

blue and green - components,

red - sum of the blue and green, analyzed signal (extract  - first 3000 samples)

The short relation d=fn(Pmax,W) after 40 W-iterations:



The relation d=fn(Pmax) for the last W-iteration


The upper relation reflects the lower one being increased by about 10% due to the deflection errors e=10% corresponding to such very small Pmax's.


Polynomial fitting applied to d=fn(W) relation:


Fitted 8-degree polynomial:

Comparison with the surrogates:

Difference between dsurr and dorig  - it points to the chaotic structure being visible especially in the range of W 400-800 px.








If you would like to be informed about the new events on this web site, please write the email to the author :  kmichalak@amu.edu.pl