Computational and Mathematical Methods in Medicine

Volume 2013 (2013), Article ID 178476, 19 pages

http://dx.doi.org/10.1155/2013/178476

Research Article

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

Franz Konstantin Fuss

School of Aerospace, Mechanical and Manufacturing Engineering, RMIT University, Bundoora Campus, Plenty Road, P.O. Box 71, Bundoora, Melbourne VIC 3083, Australia

Received 15 January 2013; Revised 18 March 2013; Accepted 30 May 2013

Academic Editor: Ernst Niebur

Copyright © 2013 Franz Konstantin Fuss. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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

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:

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 against .

Rewriting Richardson’s equation (2) in the form of

where is the length of a curve (a coastline or frontier in Richardson’s cases), is a constant, and . Considering that , (3) implies that

where is the intercept and is the gradient when plotting against , or

which shows that plotting against 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 .

The Hausdorff dimension is therefore conveniently explained as the asymptotic value of the ratio of log  to as . The fact that is the exponent in (4) links to the Hurst exponent , considering that .

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 axis [annotation: 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 is updated by summing up the amplitude of the previous datum , the distance , and : , 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 and the time axis to the duration of the signal (), 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 and , 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 ) 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, 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 ). 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 , , and , the gradient of against (Figure 2(a)) can be smaller than 1 in the MRBC method and even 2 in Higuchi’s method, even at high correlations . It has to be noted that in time series, . 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) against ; (Higuchi’s method) or boxes (MRBC) of the signal at different frequencies;  Hz (original recording frequency), 30 Hz, and 15 Hz; gradient of (MRBC) was taken at 372.5 s in Figure 1; gradient of (Higuchi’s method) was taken at 370 s in Figure 1; (b) 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 is close to 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 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 diverge significantly. The larger the difference in 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 , , and 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; ; (b) -differential (time derivative of the signal) (Higuchi’s method) and Euclidean distance 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 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 , 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:

where and are dimensionless and the reciprocal value of equals the dimensionless frequency . 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 is extremely small, resulting in Higuchi’s fractal dimension if 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 . Consequently, the fractal dimension of the modified signal (at any ) ranges between the actual fractal dimension of the original signal and unity: , where applies to , 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  Hz; (b) signal reduced to  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  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, , , and 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

where is the normalised amplitude of an th data point.

The relative length of distance normalised to is

The relative length over window width starting at datum results from

Reducing the signal to half (Figure 4(b)) by taking every other datum yields

The relative length over window width starting at datum delivers

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 and 0.5 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

The relative length over window width starting at datum delivers

Equation (15) is divided by four as is the average of four solutions (Figure 4(c)). The terms , , , , , and 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 of (if , then ), and are counters. Note that

where

If , (16) reduces to , that is, to (10).

It is important to consider equal distribution of on a logarithmic scale, that is, , or , and not , as in the latter case is weighted towards smaller frequencies.

The fractal dimension of the modified signal, , results from

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 ;(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 versus or versus 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 or on a log scale.

5. Influence of the Signal Amplitude on

In a triangular signal (Figure 5) of constant amplitude (unity) and , 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, .

The length of diagonal at is

The corresponding length of the original signal up to is

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

The relative length is defined as follows:

The slope of a line connecting two points, with the first point defined by relative length and frequency and the second point defined by and with , is calculated from

where are the asymptotic values at and , respectively. Note that

Cancelling out and calculating each individually,

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

where

Thus,

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 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 , 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 :

where denotes a triangular wave function, that is, the distance from to the nearest integer, and denotes the floor function; , , and .

If , then the Minkowski-Bouligand dimension is

Keeping constant and calculating as a function of

Thus,

if the fractal is strictly self-similar. was set to 2, was limited to 50, and ) 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 .

(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 :

where and is the Hurst exponent . 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 :

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 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 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; ) 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 data; 100 Hz data; 1000 Hz data; 10000 Hz 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,

where and 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 .

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 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, , 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 would have been reduced to one-third of . 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 (, Figure 11, amplitude range 4.62, time period 200 s, ), and would have been reduced to 6.5% of . When using Higuchi’s method [7], the multiplier would have been , 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 ; with the highest value (at 371.1 s ±75 data of Figure 10, window width of ); with the smallest value (at data of Figure 10); , the differential; optimal amplitude multiplier is defined at the peak of the differential (maximal distance between and ).

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 always asymptote to one, values at asymptote to , the fractal dimension of the original unmodified signal. Therefore, the asymptotic value of at is not necessarily larger than the one of (as will be shown in Section 8). Before reaching their asymptotic values, , 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, , is compared to of the rare event ( 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

In order to suppress further, thereby accentuating , can be optimised as to the ratio as follows:

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: . 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 ; with the highest value (at data of Figure 13, window width of ); with the smallest value (at data of Figure 13); average averaged over 60 s; , , and differentials; 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.

Optimal at 100 Hz
Optimal at 50 Hz

4G iPod Touch
1.0
2.0

3G iPhone
0.4
1.0

Minimax
1.2
2.5

(a)

(b)

Figure 14: Fractal dimension against the amplitude multiplier ; ; 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 , 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 ; (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); ), 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 (due to the high recording frequency of 10 kHz).

Figure 19 shows the of the original signal () and the optimised signal amplitude (). 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 ; ≈ ) with mean values of 27 Hz: 1.377, 35 Hz: 1.394, and 31 Hz: 1.441; (b) fractal dimensions calculated from amplitude multiplier (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.

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).

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 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.

References

H. Takayasu, In the Physical Sciences, Manchester University Press, Manchester, UK, 1990.

F. Hausdorff, “Dimension und äußeres Maß,” Mathematische Annalen, vol. 79, pp. 157–179, 1919.

C. Carathéodory, “Über das lineare Maß von Punktmengen—eine Verallgemeinerung des Längenbegriffs,” Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-Physikalische Klasse, vol. 1914, pp. 404–426, 1914.

L. F. Richardson, “The problem of contiguity: an appendix of statistics of deadly quarrels,” General Systems Yearbook, vol. 6, pp. 139–187, 1961.

B. B. Mandelbrot, “How long is the coast of Britain? Statistical self-similarity and fractional dimension,” Science, vol. 156, no. 3775, pp. 636–638, 1967.

M. J. Katz, “Fractals and the analysis of waveforms,” Computers in Biology and Medicine, vol. 18, no. 3, pp. 145–156, 1988.

T. Higuchi, “Approach to an irregular time series on the basis of the fractal theory,” Physica D, vol. 31, no. 2, pp. 277–283, 1988.

C. Sevcik, “On fractal dimension of waveforms,” Chaos, Solitons and Fractals, vol. 28, no. 2, pp. 579–580, 2006.

B. S. Raghavendra and D. N. Dutt, “Computing fractal dimension of signals using multiresolution box-counting method,” International Journal of Engineering and Mathematical Sciences, vol. 6, pp. 53–68, 2010.

B. S. Raghavendra and D. Narayana Dutt, “A note on fractal dimensions of biomedical waveforms,” Computers in Biology and Medicine, vol. 39, no. 11, pp. 1006–1012, 2009.

P. Castiglioni, “What is wrong in Katz's method? Comments on: ‘a note on fractal dimensions of biomedical waveforms’,” Computers in Biology and Medicine, vol. 40, no. 11-12, pp. 950–952, 2010.

B. B. Mandelbrot, The Fractal Geometry of Nature, Freeman, New York, NY, USA, 1983.

J. W. Blaszczyk and W. Klonowski, “Postural stability and fractal dynamics,” Acta Neurobiologiae Experimentalis, vol. 61, no. 2, pp. 105–112, 2001.

C. Sevcik, “A procedure to estimate the fractal dimension of waveforms,” Complexity International, vol. 5, pp. 1–19, 1998.

F. K. Fuss, A. Subic, and J. C. C. Chua, “Analysis of wheelchair rugby accelerations with fractal dimensions,” Procedia Engineering, vol. 34, pp. 439–442, 2012.

R. Esteller, G. Vachtsevanos, J. Echauz, and B. Litt, “A Comparison of waveform fractal dimension algorithms,” IEEE Transactions on Circuits and Systems I, vol. 48, no. 2, pp. 177–183, 2001.

H.-H. Kairies, “Functional equations for peculiar functions,” Aequationes Mathematicae, vol. 53, no. 1-2, pp. 207–241, 1997.

P. Kannappan, “Functional equations and inequalities with applications,” in Springer Monographs in Mathematics, Nondifferentiable Functions, chapter 8, pp. 359–370, Springer, Heidelberg, Germany, 2009.

V. Kulish, A. Sourin, and O. Sourina, “Human electroencephalograms seen as fractal time series: mathematical analysis and visualization,” Computers in Biology and Medicine, vol. 36, no. 3, pp. 291–302, 2006.

F. K. Fuss, “Non-linear viscoelastic properties of golf balls,” in The Impact of Technology on Sport II, F. K. Fuss, A. Subic, and S. Ujihashi, Eds., pp. 207–222, Taylor and Francis Group, London, UK, 2008.