A Robust Algorithm for Optimisation and Customisation of Fractal Dimensions of Time Series Modified by Nonlinearly Scaling Their Time Derivatives: Mathematical Theory and Practical Applications

Abstract

Standard methods for computing the fractal dimensions of time series are usually tested with continuous nowhere differentiable functions, but not benchmarked with actual signals. Therefore they can produce opposite results in extreme signals. These methods also use different scaling methods, that is, different amplitude multipliers, which makes it difficult to compare fractal dimensions obtained from different methods. The purpose of this research was to develop an optimisation method that computes the fractal dimension of a normalised (dimensionless) and modified time series signal with a robust algorithm and a running average method, and that maximises the difference between two fractal dimensions, for example, a minimum and a maximum one. The signal is modified by transforming its amplitude by a multiplier, which has a non-linear effect on the signal’s time derivative. The optimisation method identifies the optimal multiplier of the normalised amplitude for targeted decision making based on fractal dimensions. The optimisation method provides an additional filter effect and makes the fractal dimensions less noisy. The method is exemplified by, and explained with, different signals, such as human movement, EEG, and acoustic signals.

1. Introduction

The Hausdorff-Besicovitch dimension, 𝐷_𝐻, is defined by most efficient covering [1] of irregular curves and surface profiles, usually approximated by the box-counting, circle-counting, or yardstick methods. In these methods, an irregular curve is covered by 𝑁 number of boxes of size 𝑟, circles of radius 𝑟, or rulers of length 𝑟. The Hausdorff dimension, 𝐷_𝐻, is then determined from 𝐷_𝐻=lim_𝑟→0log𝑁().log1/𝑟 (1)

This method dates back to Felix Hausdorff who coined the term “fractal dimension” (“gebrochene Dimension,” [2]) by extending Carathéodory’s [3] 𝑝-dimensional measure to noninteger values of 𝑝. It was reinvented by Richardson [4], for investigating the complexity and ruggedness of coastlines with yardstick methods, who empirically found the following equation: Σ𝑙∝𝑙^−𝛼, (2)

where 𝑙 is the length of the yardstick (e.g., in km), Σ𝑙 is the sum of all yardsticks covering the total length, and the exponent 𝛼 characterises the irregularity of the coastline or the frontier. It was Mandelbrot [5] who recognised that this exponent corresponds to a fractal dimension, to be calculated from the gradient of log𝑁 against log1/𝑟.

Rewriting Richardson’s equation (2) in the form of 𝐿_𝑟=𝜆𝑟^1−𝐷, (3)

where 𝐿 is the length of a curve (a coastline or frontier in Richardson’s cases), 𝜆 is a constant, and 𝐷=𝐷_𝐻. Considering that 𝐿=𝑁𝑟, (3) implies that 1𝑁=𝜆𝑟^𝐷,1log𝑁=log𝜆+𝐷log𝑟, (4)(5)

where log𝜆 is the intercept and 𝐷 is the gradient when plotting log𝑁 against log1/𝑟, or −log𝑁log(𝑟)=−log𝜆log(𝑟)+𝐷, (6)

which shows that plotting −log𝑁/log𝑟 against log𝑟 delivers a reciprocal function of log 𝑟, with the asymptotic value of 𝐷, comparable to finite element convergence tests. In fact, increasing the number of elements 𝑁 results in 𝑟→0.

The Hausdorff dimension 𝐷_𝐻 is therefore conveniently explained as the asymptotic value of the ratio of log 𝑁 to log1/𝑟 as 𝑟→0. The fact that 𝐷 is the exponent in (4) links 𝐷_𝐻 to the Hurst exponent 𝐻, considering that 𝐷_𝐻=2−𝐻.

In addition to the standard box-counting, circle-counting, or yardstick methods, more effective methods were developed, such as the methods by Katz [6], Higuchi [7], Sevcik [8], and Raghavendra and Dutt (multiresolution box-counting MRBC and multiresolution length-based MRL methods [9]). The choice of mathematical methods gained fresh momentum when Raghavendra and Dutt [9, 10] compared existing methods and found Katz’ method [6] to be highly inaccurate. They also demonstrated the bad correlation of a hypnogram with the corresponding sleep EEG’s fractal dimensions calculated with Katz’ method, whereas Higuchi’s method provided a good correlation [10]. Based on Raghavendra and Dutt’s [10] results, Castiglioni [11] investigated Katz’ method, identified flaws, and argued that Katz’ fractal dimension is “strongly influenced by amplitude, duration and units of measure of the waveform, resulting practically useless for any real biomedical application.” As Katz [6] calculates the Euclidean distance between consecutive points, Castiglioni [11] criticises the fact that this method “sums together terms with different units” (e.g., the unit of the time scale such as seconds and the unit of the signal amplitude such as millivolt). Castiglioni [11] suggests an alternative method, based on Mandelbrot’s [12] approach that considers a mono-dimensional space, that is, the signal amplitude only, instead of a two-dimensional approach (signal amplitude and time scale). Higuchi’s [7] method is also calculated in a monodimensional space. According to Blaszczyk and Klonowski [13], “…it is important that scaling of the signal amplitude, has no influence on the results, since it causes only parallel shifting of the regression line along the ln(𝐿(𝑘)) axis [annotation: ln(𝐿(𝑘))=log𝑁 in (5)], without changing its angular coefficient [annotation: angular coefficient = gradient 𝐷 of (5)].” This refers to the principle of a mono-dimensional space, where only the signal amplitude is considered for calculating the fractal dimension. In contrast to that, scaling of the signal amplitude when using the Euclidean distance affects the fractal dimension (Castiglioni [11]).

The problem of calculating the length of a signal through the Euclidean distance of amplitude and time dimensions is that quantities of different units, for example, milliVolt (mV) and seconds (s), are squared and summed up, and finally the sum is square rooted. This does not only apply to Katz’ [6] method but also to Raghavendra and Dutt’s [9] multiresolution length-based MRL method. Moreover, Raghavendra and Dutt’s [9] multiresolution box-counting MRBC method calculates the ratio of quantities of different units: the number of boxes required to cover the distance ℎ between two consecutive data points is ceil(|ℎ|/𝑑𝑡), where ceil is the ceiling function and 𝑑𝑡 is the reciprocal value of the sampling frequency. Subsequently, the amplitude of the next data point 𝑦_(𝑖+1) is updated by summing up the amplitude of the previous datum 𝑦_(𝑖), the distance ℎ, and 𝑑𝑡: 𝑦_(𝑖+1)=𝑦_(𝑖)±|ℎ|∓𝑑𝑡, thereby summing up different units (e.g., mV and s).

These problems can be overcome when using Sevcik’s method [14], by normalising the amplitude of a signal to its range (𝑦_max−𝑦_min) and the time axis to the duration of the signal (𝑡_max), thereby transforming the signal into a unit square. This procedure solves the problem of confounding different units but renders signals with different data ranges and periods incomparable. If two signals recorded at 100 Hz for periods of 5 s and 10 s are compared, then Δ𝑥 of the unit square is 0.002 and 0.001, respectively. Assuming the signals have the same 𝑦_max and 𝑦_min, then signal 1 appears to be stretched by a factor of 2 in 𝑥-direction (former time scale) compared to signal 2, that is, the 𝑥-scale of signal 1 is multiplied by two. The same principle applies to equal period but different data range signals: if the data range of signal 2 is twice the one of signal 1, then Δ𝑦 of the unit square of signal 1 is twice as large as Δ𝑦 of signal 2, and the amplitude of signal 1 is multiplied by two (multiplier 𝑚=2) compared to signal 2. As already mentioned by Castiglioni [11], scaling of the signal’s amplitude, that is, applying different multipliers to 𝑥- and 𝑦-axes, results in different fractal dimensions when using the Euclidean distance for calculating the fractal dimensions. Whereas the unit square method appears to be a uniform method across all signals, Sevcik’s method in fact applies variable multipliers 𝑚 to signal amplitudes 𝑦, namely, 𝑚=𝑡_max/(𝑦_max−𝑦_min). However, in order to compare the fractal dimensions across signals, the multiplier 𝑚 must be the same, when using the Euclidean distance for calculating the fractal dimensions. In contrast to that, when applying a mono-dimensional space, that is, Higuchi’s method, to determining the fractal dimension of a signal, then the multiplier 𝑚 has no effect [13].

The comparability problem when transforming a signal into a unit square could be overcome when using a normalisation factor that is common to all signals that are to be compared in terms of their fractal dimensions. Such normalisation factors could be as follows: (a) the resolution of the recording device (𝑦-axis); (b) standard deviation of the signal; and (c) the window width (𝑥-axis) used for a running average method. The disadvantage when using the resolution of the recording device is that different research teams might use recording devices with different resolution, or a company releases the next generation of recording devices with better resolution. In these cases, the fractal dimensions of signals are no longer comparable. The disadvantage when using the standard deviation of the signal is that the standard deviation is different in similar signals and even changes within the signal’s time series. The disadvantage when using the window width of a running average method is that the optimal window width has to be determined beforehand, for example, from sensitivity analysis, and fractal dimensions of a signal cannot be compared any more when introducing a second factor, the multiplier, together with, and in addition to, the window width. Doubling the window width corresponds to multiplying the 𝑥-axis by 0.5. The only way of avoiding scaling problems when normalising the signal as explained above and still obtaining a dimensionless number is to normalise the signal to unit amplitude and unit time. In the end, it does not matter which parameter the time series signal is normalised to, be it amplitude range, maximum, or resolution; time period or window width; unit amplitude and unit time; as the whole exercise serves only to transform the signal data into dimensionless values in both amplitude and time dimension. This means that a time series signal, measured in milliVolts (amplitude axis) and seconds (time axis), has its amplitude and time normalised to 1 mV and 1 s, respectively.

The aim of this paper is to test Castiglioni’s [11] hypothesis, namely, that scaling of a signal, that is, applying different multipliers to normalised and dimensionless 𝑦- or 𝑥-axes, results in different fractal dimensions. Furthermore, it will be shown how scaling of a signal’s amplitude affects the Euclidean distance between data points and how this principle modifies the signal when viewed and assessed in a monodimensional space for calculating the modified signal’s fractal dimension with Higuchi’s method [7]. Finally, it will be explored whether scaling of normalised and dimensionless 𝑦- or 𝑥-axes provides an optimisation method for calculating fractal dimensions of the modified signal, which enables better decision making, classification, and quantification of events occurring throughout a time series. If decisions hinge on the fractal dimension of a biomedical signal (e.g., automatic detection of ventricular tachycardia), then the difference in fractal dimension between physiological and pathological signals (or more generally, between signals of different classes) should be as clear as possible. Thus, the question arises, whether we can “design” the fractal dimensions by scaling the amplitude of related signal classes such that the difference in fractal dimensions of the modified signal increases and becomes clearer and more pronounced? It has to be clearly stated at that point that the targeted exploration of the optimisation method for fractal dimensions refers to the fractal dimensions of a modified signal and no longer to the fractal dimensions of the original signal.

2. Statement of the Problem

Figure 1 shows an acceleration signal and the corresponding 𝐷_𝐻 calculated with three different methods: Higuchi’s method [7], Raghavendra and Dutt’s MRBC [9], and the method and algorithm developed in this study. In some parts of the signal, the 𝐷_𝐻 of all three methods are practically identical in magnitude and shape; in other parts, however, they diverge significantly. Figure 1 serves only the purpose of demonstrating these differences, without making a statement of which method is the best. In fact, as will be shown in Section 7 of this paper, all three 𝐷_𝐻 are of the same quality and therefore valid and correct (with the exception of 𝐷_𝐻<1). The point here is not the identification of the best method but rather that method which serves best for one’s own means, such as for strategic decision making and improved signal classification.

Figure 1 

Acceleration of a wheelchair recorded at 60 Hz during a wheelchair rugby match (data from [15]) against time (bottom subfigure); fractal dimension 𝐷𝐻 (running average with a window width of 151 data points over 2.5 s) of the acceleration signal calculated with three different methods (middle subfigure): Higuchi’s [7] method, Raghavendra and Dutt’s [9] MRBC, and the method developed in this study (cf. Section 4); top subfigure: absolute difference in 𝐷𝐻 between the three methods.

The principle of the MRBC method is to derive the number of boxes covering the signal amplitude change Δℎ per time step Δ𝑡 by making the ratio of Δℎ to Δ𝑡 integer with a ceiling function and summing up the number of boxes. The principle of Higuchi’s method is to sum up the change in amplitude Δℎ normalised to the time step Δ𝑡. This principle is still comparable to box-counting methods, however, by using “boxes” with a “noninteger” height or side length. If a string of data consists only of 0 and 1 (positive and negative), with a device resolution of 1, that is, zero signal with a slight noise, then the number of boxes is zero if two consecutive data are equal. In the MRBC method, relatively small changes in amplitude (with respect to Δ𝑡) always deliver a single box.

The problem for both methods is that a constant signal still has a length—in time direction, but not in amplitude direction. Nevertheless, both methods return zero Δℎ in such cases. This accounts for zero number of “boxes” (integer or non-integer ones) despite the apparent length of the signal. This problem is not only imminent in longer segments with consecutively identical data, but also influences the 𝐷_𝐻 if only two consecutive data points are of equal value.

Summing up the number of boxes within a data window at frequencies of 𝑓, 𝑓/2, and 𝑓/4, the gradient of log𝑁 against log𝑓 (Figure 2(a)) can be smaller than 1 in the MRBC method and even 2 in Higuchi’s method, even at high correlations (𝑟>0.998). It has to be noted that in time series, 1≤𝐷_𝐻≤2. From 369.6 s to 369.95 s in Figure 2, the average difference between Higuchi’s method and the MRBC method was 0.969 (Higuchi: 1.986, MRBC: 1.017) with a maximum difference of 0.99 (Higuchi: 2, MRBC: 1.01). It is surprising that two well-performing algorithms deliver opposite dimensions in extreme signals.

(a)

(b)

Figure 2 

(a) Log𝑁 against log𝑓; 𝑁=sumofamplitudes (Higuchi’s method) or boxes (MRBC) of the signal at different frequencies; 𝑓=60 Hz (original recording frequency), 30 Hz, and 15 Hz; gradient of 𝐷𝐻=0.79 (MRBC) was taken at 372.5 s in Figure 1; gradient of 𝐷𝐻=1.95 (Higuchi’s method) was taken at 370 s in Figure 1; (b) log𝑁 against time (window width = 151 data).

The reason for diverging 𝐷_𝐻 values of Higuchi’s method and the MRBC method is explained as follows. If the ratio of Δℎ to Δ𝑡 at the original sampling frequency 𝑓 is slightly smaller than 1 (e.g., 0.92), then the number of boxes (𝑁) in the MRBC method is 1, whereas the corresponding parameter (𝑁) of Higuchi’s method remains at 0.92. The ratio of Δℎ to Δ𝑡 at 𝑓/4 is close to 0.23(=0.92/4) in Higuchi’s method but still 1 in the MRBC method due to the ceiling function. Logarithming the sum of 𝑁 derived at 𝑓 leads to comparable data with MRBC and Higuchi’s methods (Figure 3(b)). However, logarithming the sum of 𝑁 derived at 𝑓/4 results in far smaller data values when using Higuchi’s method than when applying the MRBC method. This is apparent in Figure 2(b) between seconds 385 and 404 where the data derived at 𝑓/4 diverge significantly. The larger the difference in log𝑁 data, the steeper the gradient is and thus the larger the fractal dimension is. This explains why Higuchi’s method can deliver a fractal dimension of 2, whereas the MRBC method produces a fractal dimension of 1 of the same part of the signal. The ceiling function of the MRBC method leads to different scaling factors at 𝑓, 𝑓/2, and 𝑓/4 if |Δℎ|/Δ𝑡 is slightly smaller than 1, which consequently affects the fractal dimension.

(a)

(b)

(c)

Figure 3 

(a) Acceleration signal of Figure 1 from 337 s-338 s; 𝑦=amplitude; (b) 𝑦-differential (time derivative of the signal) 𝐿=|Δ𝑦| (Higuchi’s method) and Euclidean distance 𝐿=𝑚2Δ𝑦2+Δ𝑥2 at different 𝑚 against time (𝐿 is normalised such that the average 𝐿 and the area under the curves are identical in all 3 datasets); (c) Euclidean distance 𝑚2Δ𝑦2+Δ𝑥2 at different 𝑚 against 𝑦-differential |Δ𝑦|.

Esteller et al. [16] and Raghavendra and Dutt [9, 10] tested different fractal dimension algorithms with continuous nowhere differentiable functions (CNDFs, [17, 18]) of known and controllable fractal dimension 𝐷_𝐻 and with specially designed waveforms (with difference in amplitude and frequency) but not with actual time series signals (biomedical, acoustic, etc.). Such CNDFs are the deterministic Weierstrass cosine, Weierstrass-Mandelbrot cosine, and the Knopp function, as well as the stochastic Brownian motion function. The best performing methods were the Higuchi’s method and Raghavendra and Dutt’s MRBC and MRL. This proves that these methods are perfectly applicable to the aforementioned CNDFs but not necessarily how well these methods behave in general time series.

3. Solution to the Problem

In order to assess the influence of scaling a signal on its fractal dimension 𝐷_𝐻, 𝐷_𝐻 has to be determined in a 2-dimensional space rather than in a mono-dimensional one. It is therefore required that the length of a signal is calculated from the Euclidean distance between data points of time series. This makes 𝐷_𝐻 variable when scaling the amplitude [11], to be shown subsequently. The main question in this case is not “what is the correct method of calculating fractal dimensions?,” but rather “how can we maximise the information to be obtained from fractal dimensions?.” The term “maximising the information” is seen from an engineering point of view, that is, using the difference in 𝐷_𝐻 of two signals (or parts thereof when using a sliding window method) for practical application, for example, for automated decision making. This would apply to classifying sleep stadia from EEG signals, or identifying arrhythmias from EKG signals, by setting off an alarm at the onset of the latter.

As 1≤𝐷_𝐻≤2, two dimensions are anyway required, of the same physical quantity and measurement unit, as a 𝐷_𝐻 of 2 does not apply to a mono-dimensional space. This is achieved when plotting signals on paper or on the screen, and then they become real two-dimensional objects. Coastlines, rivers, and the Brownian motion are 2D curves with the same units in the two orthogonal directions. Time series signals have different units (e.g., mV and s) and therefore have to be transformed into a dimensionless space by normalisation. The difficulties of normalisation were already pointed out in Section 1. Therefore, the normalisation suggested in Section 1 is applied: 𝑦=signalamplitudee.g.mV(s)unitamplitudee.g.mV,𝑥=signalperiodunittime(s), (7)

where 𝑦 and 𝑥 are dimensionless and the reciprocal value of 𝑥 equals the dimensionless frequency 𝑓₀=1/𝑥. The suggested normalisation to unit amplitude and unit time serves only for avoiding prescaling of the amplitude, which will be modified subsequently by a multiplier 𝑚. It is essential that this multiplier is the same across all signals to be compared.

After normalisation, a variable multiplier 𝑚 is applied to the signal’s amplitude for scaling purposes and for modifying the signal. As the Euclidean distance 𝐿 between two data points is Δ𝑦²+Δ𝑥² (where 𝑥 and 𝑦 are dimensionless and unitless), it becomes 𝑚²Δ𝑦²+Δ𝑥² after applying the multiplier 𝑚. The smaller 𝑚 is, the more 𝐿 is dominated by Δ𝑥, the normalised and dimensionless time of (7). The larger 𝑚 is, the more 𝐿 is dominated by Δ𝑦, the normalised and dimensionless signal amplitude of (7). If 𝑚→∞, then 𝐿→𝑚Δ𝑦, and the fractal dimension calculated from this modified signal corresponds to Higuchi’s fractal dimension of the original signal. As Higuchi’s fractal dimension is calculated in a mono-dimensional space, the data differential is |Δ𝑦|; the fractal dimension based on which is unaffected by scaling.

In order to demonstrate how the signal is modified by scaling the amplitude with the multiplier 𝑚, the amplitude differential |Δ𝑦| and the Euclidean distance 𝐿=𝑚²Δ𝑦²+Δ𝑥² of the signal shown in Figure 3(a) are plotted in Figure 3(b). The difference between amplitude differentials at 𝑚=∞ and 𝑚=1 is extremely small, resulting in Higuchi’s fractal dimension if 1≤𝑚≤∞ in this specific part of the signal. When reducing 𝑚, the normalised 𝐿 (Figure 3(b)) below the average closes in on the normalised average 𝐿 faster than the normalised 𝐿 above the average. This non-linear scaling process is also shown in Figure 3(c). The time derivative of the original signal is therefore non-linearly modified when replacing the amplitude differential by the Euclidean distance. In other words, signal modification by replacing the amplitude differential by the Euclidean distance reduces the difference between mean amplitude differential and differential values below average faster than the difference between differential values above average and mean differential.

If the amplitude multiplier 𝑚=∞, then the Euclidean distance is identical to the scaled amplitude differential |Δ𝑦|, and therefore the fractal dimensions of the modified signal are identical to Higuchi’s fractal dimensions. If the amplitude differential |Δ𝑦| is sufficiently large compared to Δ𝑥, then multipliers 𝑚 between infinity and values far smaller than infinity (e.g., 1 in the signal shown in Figure 3(a)) deliver Euclidean distances very close to |Δ𝑦|, and therefore fractal dimensions are similar to the original signal. The smaller 𝑚 is, the more the Euclidean distance approaches Δ𝑥, and the more the signal is converted to a horizontal line. The resulting fractal dimension of the modified signal is therefore 1 if 𝑚=0. Consequently, the fractal dimension 𝐷_𝐻𝑚 of the modified signal (at any 𝑚) ranges between the actual fractal dimension 𝐷_𝐻𝑜 of the original signal and unity: 1<𝐷_𝐻𝑚<𝐷_𝐻𝑜, where 𝐷_𝐻𝑚 applies to 0<𝑚<∞, and 𝐷_𝐻𝑜 applies to 𝑚=∞, as 𝐷_𝐻𝑚→𝐷_𝐻𝑜 if 𝑚→∞.

Determining and selecting the optimal multiplier 𝑚 depends on the maximal difference between the fractal dimensions of two signals (modified by applying the amplitude multiplier 𝑚) to be compared for decision making. This optimisation process is explained subsequently, after specifying the algorithm and testing its accuracy with continuous nowhere differentiable functions (CNDFs, [17, 18]).

4. Robust Algorithm

A signal (Figure 4(a)) with amplitude 𝑦 (dimensionless) is recorded at frequency 𝑓₀ (dimensionless) over a window width of 𝑤 data points.

(a)

(b)

(c)

(d)

Figure 4 

Frequency reduction of a signal; (a) original signal at 𝑓=60 Hz; (b) signal reduced to 𝑓=30 Hz with two solutions (insert): every other even point (red) and every other odd point (green); insert: ↔= 2nd half of that green segment which precedes the first green line connecting points 2 and 4; (c) signal reduced to 𝑓=15 Hz with four solutions; (d) magnified part of Figure 4(c); the signal starts at 14.5 s; red, yellow, green, and blue lines start at points 1, 2, 3, and 4, respectively; therefore, 1/4, 1/2, and 3/4 of the preceding yellow, green, and blue segments, respectively, have to be included such that each coloured signal starts at 14.5 s exactly.

For modifying the signal, the magnitude of the amplitude 𝑦 is multiplied by the multiplier 𝑚. This method will be subsequently referred to as modified amplitude fractal dimension method (MAFDM).

The Euclidean distance 𝐿 between two data points is 𝐿_𝑖=𝑚²Δ𝑦²+Δ𝑥²=𝑚²𝑦_𝑖+1−𝑦_𝑖²+1𝑓₀², (8)

where 𝑦_𝑖 is the normalised amplitude of an 𝑖th data point.

The relative length 𝑅 of distance 𝐿 normalised to Δ𝑥 is 𝑅_𝑖=𝐿_𝑖Δ𝑥=𝑓₀𝑚²𝑦_𝑖+1−𝑦_𝑖²+1𝑓₀². (9)

The relative length over window width 𝑤 starting at datum 𝑖 results from 𝑅_𝑤=𝑅_𝑖+𝑅_𝑖+1+⋯+𝑅_{𝑖+𝑤−1}=^{𝑖+𝑤−1}_𝑖=1𝑅_𝑖. (10)

Reducing the signal to half 𝑓₀ (Figure 4(b)) by taking every other datum yields 𝐿_𝑖=𝑚²𝑦_𝑖+2−𝑦_𝑖²+2𝑓₀²,𝑅_𝑖=𝑓₀2𝑚²𝑦_𝑖+2−𝑦_𝑖²+2𝑓₀². (11)

The relative length over window width 𝑤 starting at datum 𝑖 delivers 𝑅_𝑤=𝑅_𝑖+𝑅_𝑖+1+⋯+𝑅_{𝑖+𝑤−2}+𝑅_𝑖−1(1/2)+𝑅_{𝑖+𝑤−2+1}(1/2)2=∑^{𝑖+𝑤−2}_𝑖=1𝑅_𝑖+𝑅_𝑖−1(1/2)+𝑅_{𝑖+𝑤−2+1}(1/2)2. (12)

Equation (12) is divided by two as 𝑅_𝑤 is the average of two solutions (Figure 4(b)), resulting from taking every other even datum or every other odd datum (insert in Figure 4(b)). Both solutions are embedded in Σ𝑅_𝑖 in (12). The terms 0.5𝑅_𝑖−1 and 0.5𝑅_{𝑖+𝑤−2+1} are due to the fact that the first green 𝐿 in Figure 4(b) starts only at point no. 2 (insert in Figure 4(b)), and the last green 𝐿 ends at the point preceding the last one of the shown window. This means that half-𝐿 segments were missing before point 2 (↔ in Figure 4(b)/insert) and after the second last point, if they were not included in the first place.

Reducing the signal to a quarter 𝑓₀ (Figure 4(c)) yields 𝐿_𝑖=𝑚²𝑦_𝑖+4−𝑦_𝑖²+4𝑓₀²,𝑅_𝑖=𝑓₀4𝑚²𝑦_𝑖+4−𝑦_𝑖²+4𝑓₀². (13)

The relative length over window width 𝑤 starting at datum 𝑖 delivers 𝑅_𝑤=𝑅_𝑖+𝑅_𝑖+1+⋯+𝑅_{𝑖+𝑤−4}+𝑅_𝑖−134+𝑅_𝑖−224+𝑅_𝑖−314+𝑅_{𝑖+𝑤−4+1}34+𝑅_{𝑖+𝑤−4+2}24+𝑅_{𝑖+𝑤−4+3}14×(4)⁻¹,𝑅_𝑤=^{𝑖+𝑤−4}_𝑖=1𝑅_𝑖+𝑅_𝑖−134+𝑅_𝑖−224+𝑅_𝑖−314+𝑅_{𝑖+𝑤−4+1}34+𝑅_{𝑖+𝑤−4+2}24+𝑅_{𝑖+𝑤−4+3}14×(4)⁻¹. (14)(15)

Equation (15) is divided by four as 𝑅_𝑤 is the average of four solutions (Figure 4(c)). The terms 0.75𝑅_𝑖−1, 0.75𝑅_{𝑖+𝑤−4+1}, 0.5𝑅_𝑖−2, 0.5𝑅_{𝑖+𝑤−4+2}, 0.25𝑅_𝑖−3, and 0.25𝑅_{𝑖+𝑤−4+3} are required because three-quarters of the preceding blue 𝐿 have to be added to the first blue 𝐿 (which starts at point 4; Figure 4(d)), half of the preceding green 𝐿 has to be added to the first green 𝐿 (which starts at point 3), and a quarter of the preceding yellow 𝐿 has to be added to the first yellow 𝐿 (which starts at point 2, whereas the first red 𝐿 starts at point 2).

Generalising the procedure, with 𝑤 = window width, 𝑗 = start datum of window, and 𝑛=multiplier of Δ𝑥 (if 𝑛=4, then 4Δ𝑥=4/𝑓₀), 𝑖 and 𝑙 are counters. Note that 𝑅_{𝑗,𝑤,𝑛}=^{𝑗+𝑤−𝑛}_𝑖=𝑗𝑅_𝑖+^𝑛_𝑙=1𝑅_𝑗−𝑙𝑛−𝑙𝑛+^𝑛_𝑙=1𝑅_{𝑗+𝑤−𝑛+𝑙}𝑛−𝑙𝑛×(𝑛)⁻¹=^{𝑗+𝑤−𝑛}_𝑖=𝑗𝑅_𝑖+^𝑛_𝑙=1𝑛−𝑙𝑛𝑅_𝑗−𝑙+𝑅_{𝑗+𝑤−𝑛+𝑙}(4)⁻¹, (16)

where 𝑅_𝑖=𝑓₀𝐿_𝑖𝑛,𝐿_𝑖=𝑚²𝑦_𝑖+𝑛−𝑦_𝑖²+𝑛𝑓₀². (17)

If 𝑛=1, (16) reduces to 𝑅_{𝑗,𝑤,𝑛}=∑^{𝑗+𝑤−𝑛}_𝑖=𝑗𝑅_𝑖, that is, to (10).

It is important to consider equal distribution of 𝑛 on a logarithmic scale, that is, 𝑛=1,2,4,8…, or 𝑛=1,3,9…, and not 𝑛=1,2,3,4,5, as in the latter case 𝑛 is weighted towards smaller frequencies.

The fractal dimension of the modified signal, 𝐷_𝐻𝑚, results from log𝑅_{𝑗,𝑤,𝑛}=𝐷_𝐻𝑚𝑓log⁰𝑛+𝐼, (18)

where 𝐷_𝐻𝑚 is the gradient and 𝐼 is the intercept.

In summary, the new method is characterised by the following:(1) multiplier 𝑚 applied to modifying the signal’s amplitude 𝑦;(2) Euclidean distance 𝐿 between data points;(3) different time resolutions and frequencies (𝑓₀,Δ𝑥=original,2Δ𝑥=simulatedhalffrequency,…,Δ𝑥=1/𝑓₀);(4) averaging of sum of distances for smaller frequencies (half frequency results in 2 datasets of distances); (5) running the average of 𝐷_𝐻𝑚 versus 𝑥;(6)𝐷_𝐻𝑚 is gradient of −log(Σ𝐿_Δ𝑥) versus log(Δ𝑥) or log(Σ𝐿_Δ𝑥) versus log(𝑓₀) checked by the goodness of fit (correlation coefficient), considering that there are no higher sampling frequencies available other than the original 𝑓₀;(7) equidistant lower frequencies such as 1,2,4,8,16…(not1,2,3,4,5…) or 1,3,9,27… on a log scale.

5. Influence of the Signal Amplitude on 𝐷_𝐻𝑚

In a triangular signal (Figure 5) of constant amplitude ℎ (unity) and Δ𝑥=1/𝑓₀, the amplitude ℎ is varied by the multiplier 𝑚.

Figure 5 

Triangular signal with amplitude ℎ against normalised time 𝑥 and normalised time resolution 𝑘, where 𝑘=𝑛Δ𝑥.

The time resolution is odd, that is, 𝑘=1Δ𝑥,3Δ𝑥,5Δ𝑥,7Δ𝑥,9Δ𝑥…𝑛Δ𝑥.

The length of diagonal at 𝑛=3,5,7,9… is 𝐿_𝑘>1Δ𝑥=(𝑚ℎ)²+𝑛𝑓₀². (19)

The corresponding length of the original signal (𝑘=1Δ𝑥) up to 𝑛 is 𝐿_𝑘=1Δ𝑥=𝑛(𝑚ℎ)²+1𝑓₀². (20)

The square root is multiplied by 𝑛 as the total signal length corresponds to a time period the length of which is 𝑛Δ𝑥. For example, at 𝑛=11 (Figure 5), the original triangular signal (𝑘=1Δ𝑥) consists of 11 segments, whereas the green diagonal (𝑘=1Δ𝑥, Figure 5) consists of a single line.

The relative length 𝑅 is defined as follows: 𝑅_𝑘>1Δ𝑥=𝑓₀𝑛(𝑚ℎ)²+𝑛𝑓₀²,𝑅_𝑘=1Δ𝑥=𝑓₀𝑛(𝑚ℎ)²+1𝑓₀². (21)

The slope 𝐷_𝐻𝑚 of a line connecting two points, with the first point defined by relative length 𝑅_𝑘=1Δ𝑥 and frequency 𝑓₀=1/(1Δ𝑥) and the second point defined by 𝑅_𝑘>1Δ𝑥 and 𝑓₀=1/(𝑛Δ𝑥) with 𝑛=3,5,7…, is calculated from 𝐷_𝐻𝑚=log𝑓₀𝑛(𝑚ℎ)²+1𝑓₀²𝑓−log⁰𝑛(𝑚ℎ)²+𝑛𝑓₀²×log𝑓₀𝑓−log⁰𝑛⁻¹,𝐴_1,2=lim_{𝑚→0𝑚→∞}log𝑓₀𝑛(𝑚ℎ)²+1𝑓₀²𝑓−log⁰𝑛(𝑚ℎ)²+𝑛𝑓₀²×log𝑓₀𝑓−log⁰𝑛⁻¹, (22)(23)

where 𝐴_1,2 are the asymptotic values at 𝑚→0 and 𝑚→∞, respectively. Note that 𝐴_1,2=lim_{𝑚→0𝑚→∞}𝑓log⁰𝑛(𝑚ℎ)²+1/𝑓₀²𝑓₀/𝑛(𝑚ℎ)²+𝑛/𝑓₀²×𝑓log⁰𝑓₀/𝑛⁻¹. (24)

Cancelling out 𝑓₀ and calculating each 𝐴 individually, 𝐴₁=lim_𝑚→0log𝑛²(𝑚ℎ)²+1/𝑓₀²(𝑚ℎ)²+𝑛/𝑓₀²×log(𝑛)⁻¹,𝐴₁=log𝑛²1/𝑓₀²/𝑛/𝑓₀²(𝑛)𝑛log=log²1/𝑓₀𝑛/𝑓₀×log(𝑛)⁻¹=log𝑛²𝑓₀/𝑓₀𝑛(𝑛)=(𝑛)logloglog(𝑛)=1. (25)

This result corresponds to the fractal dimension of a straight line (as 𝑚, and thus the amplitude, is 0), which has a 𝐷_𝐻𝑚 of 1. Note that 𝐴₂=lim_𝑚→∞log𝑛²(𝑚ℎ)²+1/𝑓₀²(𝑚ℎ)²+𝑛/𝑓₀²×log(𝑛)⁻¹,𝐴₂=log𝑛²𝐶(𝑛),log (26)

where 𝐶=lim_𝑚→∞(𝑚ℎ)²+1/𝑓₀²(𝑚ℎ)²+𝑛/𝑓₀²=1. (27)

Thus, 𝐴₂=log𝑛²(𝑛)=log2log𝑛log𝑛=2. (28)

This result corresponds to a densely filled area as Δ𝑥 becomes relatively small with respect to amplitude of 𝑚→∞; the 𝐷_𝐻𝑚 of such an area equals 2.

If 𝑚→0 and 𝑚→∞, then 𝐷_𝐻𝑚 asymptotes to 1 and 2, respectively, in this specific triangular signal (Figure 5), the 𝐷_𝐻𝑜 of which calculated with Higuchi’s method equals 2. Scaling the amplitude via 𝑚 corresponds to stretching or compressing the signal between an area and a straight horizontal line, respectively. This proves that the multiplier 𝑚 of the amplitude decisively influences the signal’s fractal dimension which varies between 1 and 2 (Figure 6) in this triangular signal and spans, in general, the range between 𝐷_𝐻𝑜 and 1 in different parts of a time series. Figure 6 shows the spectrum of 𝐷_𝐻𝑚, calculated from (22), against the multiplier 𝑚. This spectrum (with a positive gradient) must not be confused with the fractal spectrum (with a negative gradient), calculated from generalised Rényi’s entropy (cf. Kulish et al. [19]). The 𝐷_𝐻𝑚 spectrum (Figure 6) is decisive for the new (modified amplitude fractal dimension method) MAFDM as it enables the selection of 𝑚 for optimisation and customisation purposes, to be explained further later.

Figure 6 

Amplitude spectrum of the fractal dimension 𝐷𝐻𝑚 of the triangular signal of Figure 5 at different time resolutions and different multipliers 𝑚 of the signal amplitude.

6. Accuracy of the Proposed Algorithm

The accuracy of the algorithm for the MAFDM, (16)–(18) at 𝑚=1, was tested with four CNDFs [17, 18] in order to benchmark it against the performances of Higuchi’s and MRBC methods, as reported by Raghavendra and Dutt [9, 10].

(1) Knopp Function 𝐾(𝑡): 𝐾(𝑡)=^∞_𝑘=1𝑎^𝑘𝑏^𝑘𝑡=^∞_𝑘=1𝑎^𝑘𝑏^𝑘𝑡−𝑏^𝑘𝑡+0.5, (29)

where ⟨⟨𝑥⟩⟩ denotes a triangular wave function, that is, the distance from 𝑥 to the nearest integer, and ⌊𝑥⌋ denotes the floor function; 0.5≤𝑎≤1, 𝑏>1, and 𝑎𝑏≥1.

If 0≤𝑡≤1, then the Minkowski-Bouligand dimension 𝐷_MB is 𝐷_MB=log4𝑎.log𝑏 (30)

Keeping 𝑏 constant and calculating 𝑎 as a function of 𝐷_MB𝑏𝑎=^𝐷MB4. (31)

Thus, 𝐾𝑡,𝐷_MB=^∞_𝑘=1𝑏^𝐷MB4^𝑘𝑏^𝑘𝑡. (32)

𝐷_MB=𝐷_𝐻 if the fractal is strictly self-similar. 𝑏 was set to 2, 𝑘 was limited to 50, and 𝐾(𝑡,𝐷_MB) was calculated for one second at a frequency of 1 kHz. The results are shown in Figure 7(a). For the proposed algorithm of the MAFDM, 𝐷_𝐻𝑚 is over/underestimated at small (close to 1) and large (close to 2) theoretical 𝐷_𝐻. The MAFDM is more accurate than Higuchi’s and MRBC methods, except for 𝐷_𝐻≥1.9.

(a)

(b)

(c)

(d)

Figure 7 

Accuracy of the new amplitude spectrum fractal dimension method (MAFDM; red curves; 𝑚 = 1) compared to Higuchi’s method and the MRBC method, tested with the following continuous nowhere differentiable functions (CNDFs): Knopp/Takagi function (a); Weierstrass cosine function (b); Weierstrass-Mandelbrot function (c); Brownian motion (d).

(2) Weierstrass Cosine Function 𝑊(𝑡,𝐻): 𝑊(𝑡,𝐻)=^∞_𝑘=0𝛾^−𝑘𝐻cos2𝜋𝛾^𝑘𝑡^𝑘, (33)

where 𝛾>1 and 𝐻 is the Hurst exponent (=2−𝐷_𝐻). 𝛾 was set to 5, 𝑘 was limited to 50, and 𝑊(𝑡,𝐻) was calculated for one second at a frequency of 1 kHz. The results are shown in Figure 7(b). At medium 𝐷_𝐻, the MAFDM is marginally less accurate (maximally by 2.5% of the theoretical 𝐷_𝐻) than Higuchi’s and MRBC.

(3) Weierstrass-Mandelbrot Function 𝑊𝑚(𝑡,𝐻): 𝑊𝑚(𝑡,𝐻)=^∞_𝑘=−∞1−cos𝑏^𝑘𝑡𝑏^𝐻𝑘. (34)

For comparative reasons, 𝑏 was set to 1.5, 𝑘 was limited to 50, and 𝑊𝑚(𝑡,𝐻) was calculated for one second at a frequency of 1 kHz. The results are shown in Figure 7(c). The MAFDM has the same accuracy as the MRBC method and performs better than Higuchi’s method.

(4)  Brownian Motion Function. The time series was calculated in Matlab R2010b (by MathWorks, Natick, MA, USA) with the command wbfm(𝐻,𝑛), where 𝐻 is the Hurst exponent (=2−𝐷_𝐻) and 𝑛 is the number of data generated. As the Brownian motion function is stochastic, the shapes of the curves are not identical, and therefore ten different fractal curves were tested for each 𝐷_𝐻 (from 1 to 2, in 0.1 increments). The results are shown in Figure 7(d). The MAFDM has roughly the same performance as the MRBC method and Higuchi’s method.

In summary (Figure 8), the MAFDM did not show any disadvantage compared to the MRBC method and Higuchi’s method, delivered comparable results, and is therefore considered sufficiently accurate. The multiplier 𝑚 was set to 1 in order to compare 𝐷_𝐻𝑚 of the unscaled signal to 𝐷_𝐻𝑜 of Higuchi’s method. 𝐷_𝐻𝑚, however, is 𝑚 dependent (Figure 9(a)), and the 𝐷_𝐻𝑚 values at 𝑚=1 are at the beginning of the right-hand asymptotic segment. The sampling frequency of CNDFs influences the accuracy as well (Figure 9(b)); the higher the frequency is, the closer the theoretical 𝐷_𝐻 to the estimated one is. In contrast to that, 𝐷_𝐻𝑚 is relatively insensitive to window widths (Figure 9(b)). The 𝐷_𝐻𝑚 data obtained from the 10 Hz signal are generally higher than the ones of the 100 Hz signal. The reason for this is that the deterministic CNDFs achieve a higher 𝐷_𝐻𝑚 by increasing the amplitude of their higher harmonics. The smaller the sampling frequency, the more irregular (“chaotic”) is the signal and the higher is the estimated 𝐷_𝐻𝑚.

Figure 8 

Summary of the accuracy of the modified amplitude fractal dimension method (MAFDM; 𝑚=1) tested with the four CNDFs (nowhere differentiable functions).

(a)

(b)

Figure 9 

(a) 𝐷𝐻𝑚 of Knopp function against log 𝑚; (b) 𝐷𝐻 of Knopp function against sampling frequency at different number of samples [window width 𝑤𝑤]; 10 Hz∶𝑤𝑤=100,1000,10000 data; 100 Hz∶𝑤𝑤=100,1000 data; 1000 Hz∶𝑤𝑤=1000,10000 data; 10000 Hz∶𝑤𝑤=10000 data.

7. Optimisation and Customisation Method

The optimisation method serves for improved decision making. For example, for automated distinction between a physiological signal and a pathological one (possibly due to a life-threatening condition), their 𝐷_𝐻𝑚 must not overlap, that is, must not result in false positive or false negative diagnoses. Therefore, the 𝐷_𝐻𝑚 of the physiological signal must be as small as possible and the one of the pathological on as large as possible (or vice versa). This is achieved by optimising the amplitude multiplier 𝑚. The optimal multiplier 𝑚 is neither selected arbitrarily nor subjectively but rather follows engineering optimisation methods, namely,(i) maximising the difference between maximal and minimal 𝐷_𝐻𝑚 of signals or parts thereof; (ii) maximising the difference between maximal and average 𝐷_𝐻𝑚;(iii) maximising the ratio of maximal to average 𝐷_𝐻𝑚.

Maximally separating the 𝐷_𝐻𝑚 of physiological and pathological signals is achieved by plotting the 𝑚-dependent 𝐷_𝐻𝑚 spectra of both signals, and by identifying the maximal 𝐷_𝐻𝑚 differential, Δ𝐷_𝐻𝑚=𝐷_𝐻𝑚max−𝐷_𝐻𝑚min, (35)

where 𝐷_𝐻𝑚min and 𝐷_𝐻𝑚max are the average 𝐷_𝐻𝑚 of the physiological and the pathological signals, respectively, if the pathological signal is expected to have the higher 𝐷_𝐻𝑚, which anyway results from the spectra. For identification of events within a signal, it depends on the number of events and whether they are relatively evenly distributed. If there are more than two different events and none of them is extremely rare, then those parts of the signal that are expected to have the smallest and largest 𝐷_𝐻𝑚 are identified and their Δ𝐷_𝐻𝑚 is maximised according to (35). This is exemplified in Figure 10, in the same signal shown in Figure 2.

(a)

(b)

(c)

Figure 10 

Acceleration of a wheelchair recorded at 60 Hz during a wheelchair rugby match (data from [15]) against time (c); the two arrows refer to the time data with the lowest and highest 𝐷𝐻 used for determining the amplitude spectrum of Figure 11; (b) fractal dimension 𝐷𝐻𝑚 (running average with a window width of 151 data points over 2.5 s) of the acceleration signal calculated for different amplitude multipliers 𝑚 (cf. Figure 11); note that the fractal dimension at 𝑚→∞ corresponds to the one calculated with Higuchi’s method in Figure 2; (a) optimised 𝐷𝐻𝑚 at 𝑚=0.8.

Figure 10 shows the spectrum across different 𝑚 values against the time, as well as the events that were considered to have minimal and maximal 𝐷_𝐻𝑚. The maximal Δ𝐷_𝐻𝑚 is identified from Figure 11, where minimal and maximal 𝐷_𝐻𝑚 were plotted against the logarithm of 𝑚 and the optimised value of 𝑚 was determined from maximal Δ𝐷_𝐻𝑚. At 𝑚→0 and 𝑚→∞, both the highest and smallest 𝐷_𝐻 asymptote to one and two, respectively. The maximal 𝐷_𝐻𝑚 range amounts to Δ𝐷_𝐻𝑚 = 0.78. In this specific case, it turned out that the optimal 𝑚 is close to 1, that is, 𝑚=0.8, which means that the original signal is coincidentally close to the optimised one. The unit of the unnormalised acceleration signal is 𝑔 (9.81 m/s²). If the signal was recorded in m/s², then 𝑚 would have been 0.08, and Δ𝐷_𝐻𝑚 at 𝑚=0.8 would have been reduced to one-third of Δ𝐷_𝐻𝑚max. When transforming the signal into a unit square, that is, using Sevcik’s unit square method [14], the multiplier of the normalised 𝑦-amplitude would have been 𝑚=43.3 (log𝑚=1.64, Figure 11, amplitude range 4.62𝑔, time period 200 s, 𝑚=200/4.62), and Δ𝐷_𝐻𝑚 would have been reduced to 6.5% of Δ𝐷_𝐻𝑚max. When using Higuchi’s method [7], the multiplier would have been 𝑚>1000, and Δ𝐷_𝐻𝑜 would have been reduced to zero (Figures 10 and 11). The multiplier 𝑚 of Higuchi’s method results from not calculating the Euclidean distance and thereby reducing the normalised time dimension to zero, which in turn causes 𝑚→∞. The multiplier 𝑚 of Sevcik’s method results from transforming the signal to a unit square, thereby stretching the normalised amplitude axis with respect to the normalised time axis.

Figure 11 

Fractal dimension 𝐷𝐻𝑚 against the amplitude multiplier 𝑚; 𝐷𝐻𝑚max=𝐷𝐻𝑚 with the highest value (at 371.1 s ±75 data of Figure 10, window width of 2.5s=151data); 𝐷𝐻𝑚min=𝐷𝐻𝑚 with the smallest value (at 423.5s±75 data of Figure 10); 𝐷𝐻𝑚max−𝐷𝐻𝑚min=Δ𝐷𝐻𝑚, the 𝐷𝐻𝑚 differential; optimal amplitude multiplier 𝑚 is defined at the peak of the 𝐷𝐻𝑚 differential (maximal 𝐷𝐻𝑚 distance between 𝐷𝐻𝑚max and 𝐷𝐻𝑚min).

Plotting the 𝐷_𝐻𝑚 values at different 𝑚 against time reveals the influence of 𝑚 on the magnitude of 𝐷_𝐻𝑚. In Figure 10, 𝐷_𝐻𝑚 coloured in lime corresponds to the maximal 𝐷_𝐻𝑚 range. For 𝐷_𝐻𝑚 coloured in red (maximal 𝐷_𝐻𝑚 across the signal) or purple (minimal 𝐷_𝐻𝑚), the range approaches zero and is therefore not ideal for distinguishing different events. This is due to the fact that the highest and smallest 𝐷_𝐻𝑚 of this specific signal asymptote to one and two, respectively. At this point, it has to be noted that, whereas 𝐷_𝐻𝑚 values at 𝑚→0 always asymptote to one, 𝐷_𝐻𝑚 values at 𝑚→∞ asymptote to 𝐷_𝐻𝑜, the fractal dimension of the original unmodified signal. Therefore, the asymptotic value of 𝐷_𝐻𝑚max at 𝑚→∞ is not necessarily larger than the one of 𝐷_𝐻𝑚min (as will be shown in Section 8). Before reaching their asymptotic values, 𝐷_𝐻𝑚max>𝐷_𝐻𝑚min, but at larger 𝑚, they can converge and even switch their position. This complicates the usage of Higuchi’s method for obtaining clearly separated data of different signals (or parts thereof), as shown in Figures 2 and 10.

If there are only two different events (e.g., physiological and pathological) embedded in a signal, with the latter being rare, then the average 𝐷_𝐻𝑚 of the entire signal, 𝐷_𝐻𝑚avg, is compared to 𝐷_𝐻𝑚 of the rare event (𝐷_𝐻𝑚max if the rare event produces a higher 𝐷_𝐻𝑚 than the rest of the signal). 𝐷_𝐻𝑚 of the rare event should not significantly influence the average 𝐷_𝐻𝑚. Δ𝐷_𝐻𝑚 then results from Δ𝐷_𝐻𝑚=𝐷_𝐻𝑚max−𝐷_𝐻𝑚avg. (36)

In order to suppress 𝐷_𝐻𝑚avg further, thereby accentuating 𝐷_𝐻𝑚max, 𝐷_𝐻𝑚 can be optimised as to the 𝐷_𝐻𝑚 ratio as follows: 𝐷_𝐻𝑚𝐷ratio=^𝐻𝑚max−1𝐷_𝐻𝑚.avg−1 (37)

The three different optimisation methods, (35)–(37) facilitate that the 𝐷_𝐻𝑚 of the different signals are clearly separated by applying a zooming effect to the signals. It will be shown in Section 8, that this zooming effect also provides a filter effect.

The basic principle of optimising the amplitude of signals is to find the largest possible difference between 𝐷_𝐻𝑚 of different signals or parts thereof. The optimisation methods are exemplified by, and explained in, four different cases in Section 8. The 𝐷_𝐻𝑚 is calculated with the MAFDM in all four cases.

8. Results

Case 1 (assessment of emotional reactions with EEG). Fuss and Kulish (unpublished data) recorded the EEG of test persons during watching short movies with unexpected scary events (the so-called prank videos or screamers) in order to measure the intensity of the emotional pressure. By quickly plotting the 𝐷_𝐻𝑚 at different 𝑚, it becomes evident that the 𝐷_𝐻𝑚 correlates with the emotional pressure during the aftermath of startling. It is therefore advisable to suppress the magnitude of the small 𝐷_𝐻𝑚, by keeping the one of the maximal 𝐷_𝐻𝑚. The aim is therefore to maximise the 𝐷_𝐻𝑚 range between the highest 𝐷_𝐻𝑚 and the average 𝐷_𝐻𝑚 (Figure 12). The average 𝐷_𝐻𝑚 is calculated across the entire signal and is, in this example, smaller than the average of the highest and smallest 𝐷_𝐻𝑚. This can be achieved by maximising the 𝐷_𝐻𝑚 differential of highest 𝐷_𝐻𝑚, and average 𝐷_𝐻𝑚 whereas keeping the 𝐷_𝐻𝑚 differential of average 𝐷_𝐻𝑚 and smallest 𝐷_𝐻𝑚 as small as possible. Alternatively, the 𝐷_𝐻𝑚 ratio of highest 𝐷_𝐻𝑚 minus 1 to average 𝐷_𝐻𝑚 minus 1 can be calculated according to (37). Figure 13 shows the 𝐷_𝐻𝑚 at the original signal amplitude, the range optimised 𝐷_𝐻𝑚 and the 𝐷_𝐻𝑚 optimised to the maximal range between the highest 𝐷_𝐻𝑚 and the average 𝐷_𝐻𝑚. The effect of the latter is that the magnitude of the highest peak is kept, whereas the magnitude of smaller 𝐷_𝐻𝑚 is reduced. A further effect is that the noise level of the 𝐷_𝐻𝑚 decreases, specifically at small 𝐷_𝐻𝑚. The optimisation method therefore provides an additional filter effect. The 𝐷_𝐻𝑚 ratio increases as 𝑚 decreases and finally asymptotes at small 𝑚. The optimal 𝑚 is located at the beginning of the asymptotic segment: 𝑚=0.001. The maximal 𝐷_𝐻𝑚 at this multiplier is no longer located at 24.2 s (Figure 13) but rather at 21.9 s, where the EEG signal has its maximal amplitude. Also the 𝐷_𝐻𝑚 noise before and after startling is reduced to a minimum.

Figure 12 

Fractal dimension 𝐷𝐻𝑚 against the amplitude multiplier 𝑚; 𝐷𝐻𝑚max=𝐷𝐻𝑚 with the highest value (at 24.2s±127 data of Figure 13, window width of 1s=255data); 𝐷𝐻𝑚min=𝐷𝐻𝑚 with the smallest value (at 5.5s±127 data of Figure 13); 𝐷𝐻𝑚 average (avg)=𝐷𝐻 averaged over 60 s; max−min, max−avg, and avg−min=𝐷𝐻𝑚 differentials; (max−1)/(avg−1)=𝐷𝐻𝑚 ratio.

Figure 13 

EEG signal recorded at 256 Hz against time (bottom subfigure); the two arrows refer to the time data with the lowest and highest 𝐷𝐻 used for determining the amplitude spectrum of Figure 12; top subfigure: fractal dimension 𝐷𝐻𝑚 (running average with a window width of 255 data points over 1 s) of the EEG signal calculated for different amplitude multipliers 𝑚 (cf. Figure 12).

Case 2 (comparison of the same signal recorded with different sensors). Fuss and Chua (unpublished data) investigated the sensors influence on 𝐷_𝐻𝑚. They recorded the acceleration signal of a rugby wheelchair at different activities (collisions, pushing, coast down, and zero activity) at 100 Hz with three different sensors mounted on the frame of the chair (3G Apple iPhone, 4G Apple iPod Touch, and Minimax by Catapult). For the smart phones, it is evident that the generation number (3G and 4G) was decisive and not whether the device was an iPhone or an iPod Touch. Subsequently, the signal was reduced to 50 Hz by taking every other point. Both signals were range optimised (Figures 14 and 15), the results of which are shown in Table 1.

Table 1 

Optimal multipliers 𝑚 of the normalised acceleration amplitude of the three different devices tested.

(a)

(b)

Figure 14 

Fractal dimension 𝐷𝐻𝑚 against the amplitude multiplier 𝑚; 𝐷𝐻𝑚differential=𝐷𝐻𝑚max−𝐷𝐻𝑚min; left subfigure: original signal at 100 Hz; right subfigure: frequency reduced to 50 Hz.

(a)

(b)

(c)

(d)

(e)

Figure 15 

Wheelchair acceleration signal recorded with three different devices at 100 Hz against time (a); the two arrows refer to the time data with the lowest and highest 𝐷𝐻𝑚 values used for determining the amplitude spectrum of Figure 14; (b): unoptimised fractal dimension 𝐷𝐻𝑚 of the 100 Hz signal (𝑚=1, running average with a window width of 251 data points over 2.5 s); (c): optimised fractal dimension 𝐷𝐻𝑚 of the 100 Hz signal; (d): fractal dimension 𝐷𝐻𝑜 of the original 100 Hz signal at 𝑚=1000; (e): optimised fractal dimension 𝐷𝐻𝑚 of the signal reduced to 50 Hz (running average with a window width of 125 data points over 2.5 s) calculated for different amplitude multipliers 𝑚 (cf. Figure 14 and Table 1).

The unoptimised 𝐷_𝐻𝑚 of the 100 Hz signals (Figure 15(b)) showed pronounced differences between the 𝐷_𝐻𝑚 values. This could be due to the different device resolutions: 3G iPhone: 0.018112𝑔, Minimax: 0.006𝑔, and 4G iPod Touch: 0.000015𝑔. The high resolution of the iPod Touch could be a result of data averaging or filtering. The 𝐷_𝐻𝑚 of the 3G iPhone signal is higher, on average, compared to the two other signals. After optimisation, the 𝐷_𝐻𝑚 of the 3G iPhone signal dropped and fell within the range of the 𝐷_𝐻𝑚 of other two signals. Reducing the frequency to one half merely increases 𝐷_𝐻𝑚 but does not markedly change the trend.

Furthermore, the 𝐷_𝐻𝑚 optimisation diagram (Figure 14) shows an interesting result: the two 𝑚-spectrum curves intersect, as the 𝐷_𝐻𝑚 of minimal activity at maximal 𝑚 is higher than the one of maximal activity at maximal 𝑚. After the intersection point, the 𝐷_𝐻𝑚 difference is negative. This does not only explain the high values of Higuchi’s method (Figure 15(d)) at near-zero values of the acceleration signal, but also does this justify the optimisation procedure: select those parts of the signal which should provide different 𝐷_𝐻𝑚 and increase the 𝐷_𝐻𝑚 differential to a maximum through optimisation. Higuchi’s method would have been useless in this specific case (Figure 15(d); 𝑚=1000), as collisions and low activity would have produced 𝐷_𝐻𝑚 in the same range, and the noise would have also been more pronounced. In fact, it reflects the result of Higuchi’s method shown in Figures 2 and 10, namely, that Higuchi’s method produces maximal 𝐷_𝐻𝑚 in acceleration signals of low-level activity.

Case 3 (acoustic signature of a ball impact). Fuss [20] recorded the impact sound of golf balls at 11.025 kHz and correlated it to the hardness of the balls. In addition to FFT (for determining the power spectrum and the frequencies of the impact sound), the fractal dimension of the impact sound can be calculated. Figure 16 shows the range optimisation diagram with the optimal 𝑚 at 0.2. Reducing the original amplitude to one-fifth results in a smoother signal, spanning a larger 𝐷_𝐻𝑚 range (Figure 17). Again, the filter effect is apparent.

Figure 16 

Fractal dimension 𝐷𝐻𝑚 against the amplitude multiplier 𝑚.

Figure 17 

Acoustic signal recorded at 11.025 kHz against time (bottom subfigure); top subfigure: fractal dimension 𝐷𝐻𝑚 (running average with a window width of 55 data points over 5 ms) of the acoustic signal calculated for different amplitude multipliers 𝑚 (cf. Figure 16).

Case 4 (quality of noise). Fard, Subic, Lo and Fuss (unpublished data, submitted) recorded the acoustic sound pressure levels of a car seat at different excitation frequencies (27, 31, and 35 Hz). The 𝐷_𝐻𝑚 of the acoustic signal of the rattling seat suggested that the frequency of 31 Hz produces the highest 𝐷_𝐻𝑚, whereas the 𝐷_𝐻𝑚 values of the signal after 27 and 35 Hz excitation were almost identical. The difference in 𝐷_𝐻𝑚 between the three signals was optimized, and the average 𝑚-value of the maximal range between the 31 Hz excitation and 27&35 Hz excitations was selected (Figure 18). For comparative reasons, 𝑚 must be the same for all three signals. Even before optimisation, the 𝐷_𝐻𝑚 values of the three excitation frequencies were significantly different at 𝑃<0.05 (due to the high recording frequency of 10 kHz).
Figure 19 shows the 𝐷_𝐻𝑚 of the original signal (𝑚=1) and the optimised signal amplitude (𝑚=0.0135). Signal amplitude optimisation not only separates the 𝐷_𝐻𝑚 levels of the three different signals but also clearly exhibits a filter effect. 𝐷_𝐻𝑚 at the original signal amplitude suggests a periodic change of 𝐷_𝐻𝑚 only in the 31 Hz excitation signal, whereas the 27 and 35 Hz excitations signals are far from periodic and rather noisy and chaotic in nature. After applying the optimised 𝑚 value, all signals are clearly periodic.

Figure 18 

Fractal dimension 𝐷𝐻𝑚 (calculated from the mean of the running average) of the three signals against the amplitude multiplier 𝑚.

(a)

(b)

Figure 19 

Acoustic signals and their corresponding fractal dimensions 𝐷𝐻𝑚 (running average with a window width of 201 data points over 20 ms) recorded at 10 kHz against time; SPL: sound pressure level in Pascal; mSPL: SPL multiplied by amplitude multiplier 𝑚; (a) fractal dimensions 𝐷𝐻𝑚 of the acoustic signal calculated from original signals (amplitude multiplier 𝑚=1; 𝐷𝐻𝑚≈ 𝐷𝐻𝑜) with mean 𝐷𝐻𝑚 values of 27 Hz: 1.377, 35 Hz: 1.394, and 31 Hz: 1.441; (b) fractal dimensions 𝐷𝐻𝑚 calculated from amplitude multiplier 𝑚=0.0135 (cf. Figure 18) with mean 𝐷𝐻𝑚 values of 27 Hz: 1.057, 35 Hz: 1.109, and 31 Hz: 1.348.

9. Discussion

This paper introduces a method for optimising the fractal dimension 𝐷_𝐻𝑚 of a signal for improved decision making by scaling the signal amplitude after normalising amplitude and time axes to dimensionless values, such that the Euclidean distance between consecutive data points can be calculated. Scaling time series signals for calculating 𝐷_𝐻𝑚 is not new.

(1) Scaling was applied by Sevcik [14] when normalising the signal to the amplitude range and maximal time period, thereby transforming the signal into a unit square, such that amplitude and time axes have the same length of 1. This scaling method is chosen arbitrarily and is therefore a subjective method. The fractal dimension 𝐷_𝐻 of this transformed signal is only applicable to its unit square, but 𝐷_𝐻 of signals with different amplitude ranges and/or different time periods cannot be compared any more due to different amplitude multipliers (m = time period/amplitude range).

(2) Scaling is applied by mono-dimensional methods such as Higuchi’s method [7]. As the time dimension of the signal is not considered, the signal is transformed into a vertically oriented rectangle with horizontal side length of zero, thereby transforming the signal into a vertical line. The same effect is achieved by multiplying the amplitude axis of the signal by a very large number (𝑚→∞) such that the time axis becomes comparatively infinitesimally small and approaches zero.

As there is no rule for scaling—otherwise Sevcik’s and Higuchi’s methods would stand in contrast to each other—any scaling method can be applied, provided that it is neither arbitrary nor subjective. The method introduced in this paper follows an engineering optimisation approach, by maximising the difference between the smallest and largest 𝐷_𝐻𝑚 values across signals to be compared, or across parts of a signal. This is achieved by applying an optimised multiplier 𝑚 to the normalised amplitude of the signals, and this multiplier must be the same in all signals (or parts thereof) to be compared. In most cases, Higuchi’s method narrows down the difference between the smallest and largest 𝐷_𝐻𝑚 values and therefore is not suitable for decision making. In most cases, the optimal multiplier 𝑚 of the normalised amplitude 𝑦 is smaller than 1 (depending on the original unit if the recorded data); the smaller 𝑚 is, the less noisy the 𝐷_𝐻𝑚 dataset is (obtained from a sliding window method) and the more regular the 𝐷_𝐻𝑚 of periodic signals is due to the filter effect of small 𝑚-values. The proposed optimisation technique enables the researcher to customise the normalised signal amplitude such that maximal information is obtained, specifically for improved decision making and this with an additional filtering option.

It is suggested that the normalised amplitude 𝑦, and not the normalised time 𝑥, is scaled by 𝑚. Increasing the normalised frequency 𝑓₀ by a factor of 𝑚(=reducingΔ𝑥/𝑚) delivers the same 𝐷_𝐻𝑚 as increasing the amplitude 𝑦 by a factor of 𝑚. This explains the effect of Higuchi’s method, which does not calculate the Euclidean distance between data points and is therefore unaffected by Δ𝑥. Reducing Δ𝑥 to zero means that 𝑚→∞. Therefore, Higuchi’s method delivers the asymptotic value on the right side of the range spectrum, at maximal 𝑚 (Figures 11 and 14, with the effect seen in Figures 2, 10, and 15(d)). In fact, 𝐷_𝐻𝑜 calculated from Higuchi’s method is identical to the maximal 𝐷_𝐻𝑚 in Figure 10. Scaling the frequency 𝑓₀ via 𝑚 instead of the amplitude, however, is not advisable, as 𝑓₀ is connected to the actual sampling frequency of the measurement device.

It is evident that the multiplier 𝑚 of the normalised signal amplitude must be the same in all related signals for comparative reasons. For example, this applies to parts of a signal which are compared, as well as to signal data related to a therapeutic period (of one or more patients) or to a specific pathology across a cohort of patients. This is not the case in Sevcik’s unit square method, as different multipliers result if the amplitude range and/or time periods of signals to be compared are not identical. If 𝐷_𝐻𝑚 is determined as a function of normalised time with a running average analysis, then the window width has to be taken care of. If the window width is widened, then the 𝐷_𝐻𝑚 of the signal is less noisy, the valleys become narrower (and the peaks wider) along the time axis, and the valleys become shallower (and the peaks flatter) along the amplitude axis. The optimal multiplier 𝑚 has to be determined at the same window width (e.g., Figure 18) as subsequently used for the optimised fractal dimension analysis (e.g., Figure 19). For the optimisation procedure, the 𝐷_𝐻𝑚 averaged over time or over the window width is used (e.g., Figure 18). Different recording devices are likely to deliver different 𝐷_𝐻𝑚 if their sensor specifications do not match. Optimisation of the 𝐷_𝐻𝑚 signal can mitigate the problem.

Finally, it has to be pointed out again that the fractal dimension 𝐷_𝐻𝑚, obtained after modifying the signal by multiplying its amplitude by the multiplier 𝑚, is no longer the fractal dimension 𝐷_𝐻𝑜 of the original signal but rather the fractal dimension 𝐷_𝐻𝑚 of the modified or transformed signal. 𝐷_𝐻𝑚 serves only for comparing two or more signals (or parts of the same signal) for the purpose of decision making and for classifying or quantifying the intensity of events occurring throughout the time series signal. Thus, 𝐷_𝐻𝑚 does not serve for accurate evaluation of the fractal dimension 𝐷_𝐻𝑜 of the original signal. Yet, fractal dimensions of signals are often compared in order to draw conclusions, the practical applications of which ultimately result in decision making. The latter is optimised with the method described and presented in this study.

10. Conclusions

Different, usually well-performing, methods for computing the fractal dimensions of time series can produce opposite results in extreme signals, namely, dimensions of 1 (or even <1) and 2. These results can be reproduced when calculating the Euclidean distance between data points and modifying their amplitude by a multiplier. Equally, standard methods for computing fractal dimensions use different scaling methods, that is, different amplitude multipliers, which makes it difficult to compare fractal dimensions obtained from different methods.

In order to overcome this problem, a fractal dimension optimisation method is proposed that computes the fractal dimension of a normalised (dimensionless) and modified time series signal with a robust algorithm and a running average method, and maximises the difference between two fractal dimensions, for example, a minimum and a maximum one.

The optimisation method developed in this study is essential for decision making, particularly when the decision hinges on the condition that two different signals (or parts of a signal) have different fractal dimensions. The optimisation method is applied to calculating the maximal fractal dimension differential, as a function of the amplitude multiplier. Alternatively, the fractal dimension differential between maximal dimension and average dimension can be maximised. This optimisation method enhances fractal dimensions above average and suppresses those below average. Maximising the ratio of maximal to average dimension provides a further optimisation effect.

The smaller the amplitude multiplier, the stronger is the filter effect on the fractal dimension. Large multipliers result in a very noisy fractal dimension signal. The filter effect even reveals periodic events of a signal and enhances larger fractal dimensions.

The fractal dimension 𝐷_𝐻𝑚, obtained after modifying the signal by multiplying its amplitude by the multiplier 𝑚, is no longer the fractal dimension 𝐷_𝐻𝑜 of the original signal but rather the fractal dimension 𝐷_𝐻𝑚 of the modified signal only, used solely for the purpose of decision making and for classifying or quantifying the intensity of events occurring throughout the time series signal.

Conflict of Interests

The author reports no conflict of interests.

Acknowledgment

The author thanks the two reviewers for their invaluable critical comments and suggestions, which significantly contributed to the improvement of the paper.