Time-Series Feature Extraction by Return Map Analysis and Its Application to Bearing-Fault Detection

Abstract

Developing new features for time-series characterization is a current challenge in data science and machine learning. In this paper, we propose a new metric based on a simple and efficient algorithm, namely, the return map. The return map analysis is well established in the field of non-linear dynamics, in particular, for fitting the parameters of a chaotic system from a waveform, or to attack a chaotic communication channel. We show that our metrics work well for both non-linear dynamics and time-series feature extraction problems in the field of machine learning. In an experiment aiming to classify vibration signals of normal and damaged bearings, we compare our method with two other methods that reported to have excellent accuracy, based on entropy and statistical feature distribution, respectively. We show that our method achieves higher accuracy with almost the lowest time costs, which was confirmed in experiments with two different datasets containing three main classes of bearings: normal, with inner race faults, and with outer race faults, having different damage origins and recorded in various conditions. In particular, for the dataset supplied by Case Western Reserve University, our method reached an accuracy of 100% at signals of 5000 sample points length, with a total time of 0.4 s required for feature estimation, while the entropy-based method reached an accuracy of 95% with a time of 100 s, and a statistical feature distribution method reached an accuracy of 93% with a total time of 1.9 s. Results show that the developed method is better suited to real-time bearing condition monitoring applications than most of the methods reported to date.

1. Introduction

Extraction of features from data, such as time series and images, is a stage of data processing in many fields of study and is of particular importance in machine learning. Time series, or one-dimensional signals, play a major role in many domains and represent, for example, the output of acceleration or pressure sensors in technical systems, heart and muscle electrical activity in biological objects, signals in communication lines, seismic activity in geological applications, evolution of economic indices, and many other types of data. Processing of time series includes their characterization, which means that long signals are transformed into a compact feature space representing only important properties of a particular signal. The resulting feature sets are used for measurements, classification, and forecasting [1]. One of the possible applications of such feature sets is the characterization of non-linear and chaotic dynamics in natural, technical, or simulated systems. In engineering, extracting features from chaotic time series is useful in designing sensors [2] and chaotic communication systems [3]. Fault detection is another field where automatic diagnostics based on feature extraction are widely used [4,5].

A popular diagnostic task for machine learning is the detection of bearing faults. About 40% of all failures in electric motors are associated with bearings [6]; 70% of gearboxes and generators downtime is caused by bearing faults [7]. For providing bearing diagnostics, different methods of signal acquisition are utilized. More often, vibration sensors or accelerometers are used [4,8]. Alternatively, acoustic signals [9,10], phase currents [11,12], and bearing temperature [13,14] are acquired using appropriate sensors. Analysis of vibration signals shows high efficiency, although it is not applicable in the case of aggressive environments [15]. Bearing damages may develop slowly, and the bearing thus does not fail immediately. This allows preventive maintenance strategies based on condition monitoring to be implemented. Therefore, it is preferable to develop diagnostic methods that are sensitive not only to the fact of the presence or absence of a fault but also reveal its type and extent of damage. Bearing diagnostics are a challenging task because the shaft rotation speed of the machine in an industrial environment can vary, which makes frequency-domain methods require special adaptation [12], so avoiding them is a fine strategy, though frequency-domain features are very efficient. Also, torque loads and radial forces change characteristics of vibration in the machine, which should be considered when selecting features and training the classifier. All aforementioned forces research was intended to develop new algorithms with higher accuracy and computational performance. The technique recommended by the Society For Machinery Fail Prevention Technology (MFPT) implements a signal decomposition using FIR filters and subsequent calculation of approximate entropy (ApEn) as a metric for feature extraction, which results in up to 100% accuracy [16]. Discrete wavelet transform (DWT) decomposition with the subsequent characterization of components qualitatively [17] or by statistical features [18] is another fruitful approach, allowing classification accuracy of more than 95%. High accuracy and robustness are provided by the use of the generalized Gaussian distribution, which performs precise statistical correspondence to the coefficients of the wavelet sub-bands, having a heavy-tailed bell-shaped distribution differentiated for various classes of faults. Some works use entropy estimations together with other features, such as energetic and time-domain estimates [19].

For clarity, we presented some of the recent frequency-independent methods in Table 1, pointing out their features, strengths, and weaknesses.

The methods presented in Table 1, as well as many others found in the literature, consist of four standard stages. First, the signal is filtered, where FIR filters, empirical modes decomposition, or wavelet transform are used. Second, features, such as amplitude, statistical, or entropy estimates are calculated from the filtered sub-bands. Third, the dimensionality of the obtained feature vectors is reduced, and a classifier is applied in the fourth stage. Usually 1-NN, SVM, or neural networks are used as classifiers, but the particular classifier seems to have little impact on the resulting accuracy [18]. Along with the aforementioned, the methods presented in Table 1 can be considered as typical state-of-the-art methods for analyzing bearing vibration signals. From the literature overview, it follows that researchers focus primarily on the accuracy of the methods, but do not provide data on the computational time, though this is critical for real-time applications. Among the time-consuming operations, it is worth mentioning the need to compute multiple features with subsequent dimensionality reduction, as well as the computation of any entropy estimates. For this reason, it is of interest both to develop new methods that utilize a limited set of rapidly computable features and to compare them with the presented algorithms.

In this study, we adapt the return map technique, well known in non-linear dynamics, to produce signal features, similar to entropy, which was first used in non-linear dynamics rather than in signal analysis. In [3], we have presented a novel method of chaotic time-series comparison based on return maps [20], called quantified return map analysis (QRMA), although other chaotic time-series comparison methods utilizing return maps have been developed [21,22]. Return map is a mathematical function that maps peaks and valleys of a time series onto the 2-dimensional plane, taking into account the amplitude or amplitude and phase properties of a time series, depending on the particular algorithm realization. The aforementioned works consider return maps as a tool to estimate the similarity of two chaotic waveforms. In particular, our previous work [3] demonstrates that QRMA can effectively distinguish modes of operation of a chaotic system, including the case of additive noise, which allows this algorithm to be used to attack chaotic communication systems. Some studies [21,22] focus on estimating parameters of a model system using observed chaotic time series. We introduce another approach, which would characterize the time series by a scalar or a vector obtained by a numerical algorithm applied to the return map built from the time series. Thus, the measure of signal similarity would be performed by the comparison of obtained numerical values. This opens up the possibility of not only using return maps to estimate the parameters of the chaotic oscillators but also extends the QRMA algorithm’s scope of application to a wide class of machine-learning problems, which forms the aim of this work. The contribution and novelty of this work are as follows:

• A new algorithm for time-series characterization was developed using two types of return maps—amplitude and amplitude-phase—and matching time-series to feature vector. The method was called distance-based QRMA, or, shortly, dRMA.
• The developed algorithm was applied to a set of test problems. We examined the similarity of dRMA and the largest Lyapunov exponent (LLE) of the chaotic time series, the possibility of estimating parameters of chaotic dynamics in noisy conditions, and the classification of bearing faults.
• We compare the classification accuracy and performance of the proposed bearing-fault classification algorithm with competing methods that have high scores in the literature, namely, the statistical DWT-based method and the entropy-based method. We show that the developed method performs more accurately and faster than the competitors.

Table 1. Comparison of some recent methods for bearing-fault classification which do not use frequency-based features.

Method with Reference	Features Description	Strengths	Weaknesses
Empirical mode decomposition and statistical feature extraction [23].	Metrics applied to intrinsic mode functions, such as kurtosis, skewness, etc.	High performance for partitioning a dataset into multiple classes.	Feature dimensionality reduction is required to prevent high computational complexity and accuracy degradation.
Wavelet decomposition and statistical density modeling [18].	Parameters of generalized Gaussian distribution, which approximates the density of decomposition coefficients.	High accuracy and robustness. Choice of the particular wavelet filters, decomposition levels, and classifiers is not essential.	Computational cost of the method is not reported. Convergence issues are possible.
Entropy estimation of raw signals [24].	Sample entropy with parameter values $m=0..2$.	No pre-processing for accelerometer signal is required. The method is accurate and simple to implement.	Calculating entropy is computationally expensive.
Complementary ensemble empirical mode decomposition (CEEMD) and weighted multi-scale entropy estimation [25].	Weighted multi-scale fuzzy, permutation, and dispersion entropies with scale factors 1..10 after CEEMD filtration.	High performance for partitioning a dataset into multiple classes, including different fault severity.	Method involves multiple computationally expensive stages, also feature dimensionality reduction may be required.

Table 1. Comparison of some recent methods for bearing-fault classification which do not use frequency-based features.

Method with Reference	Features Description	Strengths	Weaknesses
Empirical mode decomposition and statistical feature extraction [23].	Metrics applied to intrinsic mode functions, such as kurtosis, skewness, etc.	High performance for partitioning a dataset into multiple classes.	Feature dimensionality reduction is required to prevent high computational complexity and accuracy degradation.
Wavelet decomposition and statistical density modeling [18].	Parameters of generalized Gaussian distribution, which approximates the density of decomposition coefficients.	High accuracy and robustness. Choice of the particular wavelet filters, decomposition levels, and classifiers is not essential.	Computational cost of the method is not reported. Convergence issues are possible.
Entropy estimation of raw signals [24].	Sample entropy with parameter values $m=0..2$.	No pre-processing for accelerometer signal is required. The method is accurate and simple to implement.	Calculating entropy is computationally expensive.
Complementary ensemble empirical mode decomposition (CEEMD) and weighted multi-scale entropy estimation [25].	Weighted multi-scale fuzzy, permutation, and dispersion entropies with scale factors 1..10 after CEEMD filtration.	High performance for partitioning a dataset into multiple classes, including different fault severity.	Method involves multiple computationally expensive stages, also feature dimensionality reduction may be required.

The further article is organized as follows. In Section 2, we describe the methods, as well as the measured and artificial data. Section 3 demonstrates the experimental results obtained, including an evaluation of the classification accuracy. Section 4 concludes the paper.

2. Materials and Methods

2.1. Return Map Analysis

2.1.1. Return Maps

Return map analysis (RMA) is a known technique for evaluating the properties of chaotic time series. It was introduced by Perez et al. in 1995 [20] for the extraction of messages incorporated in the chaotic waveform. The return map is built using peaks and valleys of a time series. Figure 1 shows the signal and its peaks and valleys. The amplitude return map (ARM) algorithm takes into account the peak and valley amplitudes of a signal, while the amplitude-phase return map (APRM) algorithm calculates both the peak amplitudes and the inter-peak and inter-valley distances of a signal.

The formulae for the ARM algorithm are as follows:

$$\begin{array}{c}\hfill A={\displaystyle \frac{X_m+Y_m}{2}};B=X_m-Y_m;\\ \hfill C=-{\displaystyle \frac{X_{m+1}+Y_m}{2}};D=X_{m+1}-Y_m.\end{array}$$

(1)

where $X_m$ is an array of peak values (amplitudes) without the last element, and $X_{m+1}$ is an array of peak amplitudes without the first element. The definitions of $Y_m$ and $Y_{m+1}$ are similar to the definitions of $X_m$ and $X_{m+1}$ except that they are the arrays of valleys of a signal.

In the APRM algorithm, the A, B, C, and D arrays are calculated as follows:

$$\begin{array}{c}\hfill A=\Delta T_{peaks+1}-\Delta T_{peaks};B=X_m;\\ \hfill C=\Delta T_{valleys+1}-\Delta T_{valleys};D=-Y_m.\end{array}$$

(2)

where $\Delta T_{peaks+1}$ is an array of peak times without the first element, $\Delta T_{peaks}$ is an array of peak times without the last element, $X_m$ is an array of peak amplitudes, and $Y_m$ is an array of valley amplitudes. The definitions of $\Delta T_{valleys+1}$ and $\Delta T_{valleys}$ are similar to the definitions of $\Delta T_{peaks+1}$ and $\Delta T_{peaks}$ with the difference being that these are the arrays of valleys of a signal.

2.1.2. QRMA and dRMA

For years, the return maps method was applied to attack chaos-based communication systems [26], outperforming the spectral methods in this task. The disadvantage of this approach was qualitative rather than quantitative results, but more recent work has suggested different estimates of return maps [3,21,22].

The QRMA algorithm proposed by Rybin et al. [3] allows the estimation of return maps. It consists of the following stages: for two compared signals, return maps are found, then a 2-dimensional plane of return maps is presented as 2-dimensional histograms with a certain bin count (usually $10\times 10$), histograms are subtracted, and the obtained difference is quantified as a percentage of the number of mismatched areas of the return map to its total area. This method is sensitive and noise-resistive but has a narrow dynamic range, which makes a comparison of two-time series with notably different return maps impossible.

In this paper, we propose a new method for time-series estimation called distance-based QRMA, or dRMA for short. The steps of the distance-RMA algorithm are as follows. First, use RMA to obtain the four arrays $(A,B,C,D)$. Then, concatenate them in pairs $[A,C]$, $[B,D]$, but for the ARM, use $-C$ instead of C. Second, compute the Euclidean distances between the two arrays obtained in the previous step. As the last step, calculate the mean value of these distances.

Figure 2 shows the execution of the second step of the dRMA algorithm. The Euclidean distance is calculated for each of the concatenated points, resulting in RMA arrays. For example, a distance between point 1 and point 4 is $d_{14}$.

The pseudocode for the dRMA algorithm is listed below, see Listing 1.

Listing 1. Distance-based QRMA algorithm.

A, B, C, D = RMA (signal) X = concatenate (A, “−”C) Y = concatenate (B, D) distances = distance (X, Y) m_dist = mean (distances)

2.2. Time-Series Decomposition and Analysis Techniques

In our study, we use two decompositions. Empirical mode decomposition (EMD) is used with our proposed dRMA algorithm, and discrete wavelet transform (DWT) is used in the competing algorithm [18] used for comparison.

2.2.1. EMD and DWT

Empirical mode decomposition (EMD) is a data-driven method used for analyzing non-stationary and non-linear time series data. It decomposes a signal into a finite number of oscillatory components called intrinsic mode functions (IMFs). Each IMF represents a different frequency component of the signal. The decomposition process is performed iteratively by extracting the oscillatory component with the highest frequency as an IMF and then subtracting it from the original signal. This process continues until the remaining signal is a monotonic function, which is considered the final IMF. EMD is widely applied in signal processing and time-frequency analysis.

The EMD consists of sequential calculation of empirical modes $c_j$ and residuals $r_j=r_{j-1}-c_j$. The result is a signal decomposition [27]:

$$X\left(t\right)=\sum _{j=1}^n{c}_j+r_n,$$

(3)

where n is the number of empirical modes.

Discrete wavelet transform (DWT) is a technique used in signal processing and image compression to analyze and decompose signals or images into different frequency components. It involves transforming a signal into a set of wavelet coefficients, which represent the signal’s intensity at different scales and resolutions. DWT is widely used in various applications, such as image and video compression, denoising, and feature extraction [18].

2.2.2. Statistics Estimation of DWT-Decomposed Signal

A high-performance bearing-fault detection algorithm using wavelet transform and generalized Gaussian distribution was proposed in [18]. It consists of applying discrete wavelet decomposition to a vibration signal, constructing a histogram for each sub-band, approximating the histogram with a generalized Gaussian distribution, and extracting statistical features from it.

The distribution function of the GGD is defined as

$$p(x;\alpha ,\beta ,\mu )={\displaystyle \frac{\beta }{2\alpha \mathrm{\Gamma }(1/\beta )}}{e}^{-{\left(\right|x-\mu |/\alpha )}^{\beta }},$$

(4)

where $\mathrm{\Gamma }$ means the Gamma function:

$$\mathrm{\Gamma }\left(z\right)={\int }_0^{\infty }{t}^{z-1}{e}^{-t}dt,z>0.$$

(5)

2.2.3. Entropy Estimation

The core benefit of entropy estimation is that it can be applied to the raw time series. Wong et al. [24,28] achieved up to 99.9% accuracy of bearing-fault classification with only sample entropy features. Seera et al. [29] examined a combination of the power spectrum and sample entropy applied also to raw vibrational signals to bearings, and achieved an accuracy 99.84% with three entropy features. In contrast, Ma et al. [19] applied ensemble intrinsic time-scale decomposition (ITD) with adaptive noise to the raw signals of rolling bearings to compute entropy at different scales. However, decomposition was required due to the lower quality of the input data for classification compared to previous papers.

The sample entropy (SampEn), or an approximation of the Kolmogorov entropy, is an indicator of data unpredictability. When new and unexpected patterns appear in a dataset, its entropy increases. The sample entropy is calculated as follows [6]:

• A data matrix is created and the correlation corresponding to zero is calculated.
• The Chebyshev distances (6) are calculated to exclude cases where differences between data points occur only along one axis, and the Heaviside function is calculated from the distances found (7):

  $$dist=max\left|tempMat\right(:,i+1:N-dim)-repmat(tempMat(:,i),1,N-dim-i)$$

  (6)

  where $tempMat(:,i+1:N-dim)$ extracts all subsequent patterns of length m from position $i+1$ to $N-dim$, and $repmat\left(\right)$ creates a matrix that repeats this pattern a certain number of times. This allows the current pattern to be compared with subsequent patterns.

  $$D=(dist<r)$$

  (7)
• The proportion of vectors within region r is calculated:

  $$coun{t}_i={\displaystyle \frac{su{m}_D}{N-dim}}$$

  (8)
• The average proportion of matches over the length of the sequence is calculated for all vectors:

  $$corre{l}_{m-dim+1}={\displaystyle \frac{su{m}_{count}}{N-dim}}$$

  (9)
• Then the sample entropy is calculated:

  $$sampEn=log\left({\displaystyle \frac{corre{l}_1}{corre{l}_2}}\right)$$

  (10)

When calculating the entropy, it is important to take into account the length of the sequence, or embedding dimension, in order to maximize accuracy. In machine learning, calculation of SampEn for different embedding dimensions (usually dimensions 0–2), is used to enlarge the features vector.

Approximate entropy (ApEn) is a complexity index that measures the complexity of the signal sequence and the probability of generating new patterns. This method is also used for feature extraction in machine learning. The specific steps for approximate entropy calculation are as follows [27]:

• The sampled signal $x\left(1\right),x\left(2\right),\dots ,x\left(N\right)$ should be reconstructed as $X\left(i\right)=\left[x\right(i),x(i+1),\dots ,x(i+m-1\left)\right]$, where $i=1,2,\dots ,N+m-1$, and m is the vector dimension.
• Calculate the distance between $X\left(i\right)$ and $X\left(j\right)$:

  $$d(i,j)=\underset{k=0\to m-1}{max}|X(i+k)-X(j+k)|$$

  (11)
• Set deviation r and calculate the number of $d(i,j)<r$ for each vector $N^m\left(i\right)$.
• Calculate the ratio between $N^m\left(i\right)$ and the total vector distance, then calculate the average of $C_i^m\left(r\right)$:

  $$C_i^m\left(r\right)={\displaystyle \frac{N^m\left(i\right)}{N-m+1}}$$

  (12)

  $$C^m\left(r\right)={\displaystyle \frac{1}{N-m+1}}\sum _{i=1}^{N-m+1}ln{C}_i^m\left(r\right)$$

  (13)
• The final equation for approximate entropy is:

  $$ApEn(m,r,N)=C^m\left(r\right)-C^{m+1}\left(r\right)$$

  (14)

Both ApEn and SampEn have the algorithm complexity of $O\left(N^2\right)$, where N is the length of the time series analyzed. Such complexity makes entropy computation time-consuming, which limits the application of these methods in real-time problems.

2.3. Chaotic Systems for Artificial Time-Series Generation

In this section, we describe the chaotic systems that have been used to generate artificial time series. The convenience of using chaotic systems as generators of artificial non-linear signals lies in the possibility to flexibly control their parameters, so one can provide the necessary difference between the generated time series. For all the experiments, the simulation was performed using the fourth-order Runge–Kutta method with a time step $h=0.01$.

2.3.1. Unified System

The unified system was introduced in 2002 by J. Lü and G. Chen [30]. A unified system is a chaotic system that can exhibit the behavior of three systems, namely, Lorenz, Lü, and Chen, depending on a single parameter $\alpha$ value. It is described by the following differential equations:

$$\begin{array}{cc}\hfill & \dot{x}=(25\alpha +10)(y-x);\hfill \\ \hfill & \dot{y}=(28-35\alpha )x-xz+(29\alpha -1)y;\hfill \\ \hfill & \dot{z}=xy-(\alpha +8)z/3.\hfill \end{array}$$

(15)

The possible ranges for $\alpha$ include:

• $\alpha \in [0,0.8)$ for obtaining the Lorentz-like dynamics.
• $\alpha =0.8$ for the Lü-like dynamics.
• $\alpha \in (0.8,1]$ for the Chen-like dynamics.

Phase portraits of the Unified system at different values of parameter $\alpha$ are presented in Figure 3. The values of the bifurcation parameter are chosen in such a way that the system smoothly transitions from the chaotic regime to the periodic one.

2.3.2. Gokyildirim System

The Gokyildirim chaotic system was presented in [31] in 2023. It is a seven-term chaotic system with quadratic non-linearity and is described by the following differential equations:

$$\begin{array}{cc}\hfill & \dot{x}=y+xy-yz;\hfill \\ \hfill & \dot{y}=x+ayz;\hfill \\ \hfill & \dot{z}=b{y}^2+xy.\hfill \end{array}$$

(16)

where the standard parameter set is as follows: $a=0.1$, $b=0.7$. Standard initial conditions are $x_0=1.5$, $y_0=3$, and $z_0=6$. The presented system shows sensitive and rich dynamic behaviors in a small range of system parameters, which makes it suitable for different applications.

Phase portraits of the Gokyildirim system at different values of parameter a are presented in Figure 4. The values of the bifurcation parameter are chosen in such a way that the system smoothly transitions from the periodic regime to the chaotic.

2.4. Classification Workflow

The flowchart of the classification process is presented in Figure 5.

As the competitive methods, we implemented entropy estimation of raw signals [24,28] and wavelet decomposition followed by statistical density modeling [18]. All three feature descriptors were designed to produce four features, subsequently feeding them to the unified KNN classifier.

For the entropy-based method, the derived features are the result of applying the sample entropy to a signal with embedding dimensions of zero, one, and two, and the result of the approximate entropy of a signal with the embedding dimension of zero. As a result of using the generalized Gaussian distribution, the values of $\alpha$ and $\beta$ were obtained for the first and second discrete wavelet decomposition levels. For the proposed dRMA method, the derived features are the values of the ARM and APRM, calculated for two close intrinsic mode functions (IMFs). Particular indices of IMFs depend on the dataset.

2.5. Bearing Datasets

2.5.1. Paderborn Dataset

In this paper, the bearing dataset of Paderborn University [32] was used as a test machine-learning problem. The Paderborn bearing dataset is a benchmark dataset for condition monitoring of rolling bearings. It includes measurements of motor currents, which are obtained from frequency inverters. The dataset distinguishes between two categories of bearing damage: artificially induced and naturally occurring during accelerated lifetime tests. Artificial damage was induced manually using precision machining tools or etching. Both types of damages are present on the inner and outer rings of the ball bearing of model SKF 6203.

The dataset consists of recordings for 32 bearings, 4 modes for operating conditions, and 20 signals for each bearing under specified operating conditions. The dataset also includes data about current and other quantities such as temperature, torque, etc. The load conditions are summarized in

Table 2. Operating parameters of bearings from Paderborn dataset.

No.	Rotational Speed [rpm]	Load Torque [Nm]	Radial Force [N]	Name of Setting
0	1500	0.7	1000	N15_M07_F10
1	900	0.7	1000	N09_M07_F10
2	1500	0.1	1000	N15_M01_F10
3	1500	0.7	400	N15_M07_F04

.

Damages to the outer race (OR) and inner race (IR) of the motor bearing are presented in Figure 6. Also, an accelerometer position is depicted.

In order to ensure an even distribution of the data, four bearing vibration data points were taken for each of the three classes. The chosen bearings in this study are:

• Without any defects (referred to as “normal”): K001, K002, K003, and K006.
• With an outer ring defect (referred to as “OR”): KA04, KA05, KA22, and KA30.
• With an inner ring defect (referred to as “IR”): KI01, KI16, KI17, and KI21.

Damaged bearings used in the dataset differ in cause and character of damage. For outer race faults, bearing KA05 was damaged by an electric engraver, while other bearings were damaged in accelerated life tests. KA04 and KA22 have pitting, and KA30 has plastic deformations and indentations. Among inner race faulty bearings, KI01 contains trenches generated by electrical discharge machining, and others have single-point pitting with different extents of damage obtained in accelerated life tests. All of this complicates the task of fault classifying, as various natures of damage create different vibration signals, even if they are in the same class.

Also note that for the bearings, data on both available rotational speeds and loads were used. Figure 5 shows the algorithms used for classification in the experiment on Paderborn data. The parameters of all methods and the features to be classified were selected to give the best classification accuracy. For all recordings, the sample rate is $f_{samp}=64$ kHz, and for all the implemented algorithms, we did not use resampling.

2.5.2. CWRU Dataset

The Case Western Reserve University (CWRU) dataset consists of vibration signals recorded from the drive end and the fan end of the motor housing. In the experiments, a motor with 2 HP was used. Bearing defects can be classified into four groups:

• normal (without any defect);
• with an inner race defect;
• with a ball defect;
• with an outer race defect. The three outer race positions relative to the load zone are:

  –

  6 o’clock (centered);

  –

  3 o’clock (orthogonal);

  –

  12 o’clock (opposite).

In this study, we used three classes: normal bearings, bearings with an outer race defect, and bearings with an inner race defect. Only signals from the drive end were used. The sizes of the defects were 0.014 inches and 0.021 inches. Two loads were used: zero and two HP. Bearings with an outer race defect were located at 6 o’clock. The sample rate was $f_{samp}=12$ kHz.

3. Results

All the experiments in this section were performed using the workstation Z790 Gaming X AX with the following characteristics:

• CPU 13-th Gen Intel® Core ™ i7-13700K 3.4 GHz, 16 cores
• RAM 32 Gb
• GPU NVIDIA GeForce RTX 4070
• SSD Samsung 980 PRO 1 Tb.

3.1. dRMA as a Measure for Chaotic Dynamics Characterization

In order to evaluate the possibility of applying the dRMA method to characterize complex time series, we investigate the behavior of the dRMA features synchronously with the change in the bifurcation parameter of chaotic systems. As waveform generators, we consider the Unified (15) and Gokyildirim (16) systems.

To indicate changes in the dynamics of chaotic systems, we use two standard tools: the bifurcation diagram and the largest Lyapunov exponent (LLE). When LLE becomes positive, the system switches from periodic to chaotic mode, and vice versa. This is a significant change in the dynamics, and we expect that an adequate metric will also follow the change in LLE. With that, the metric should not repeat the LLE dynamics completely.

Figure 7 and Figure 8 present the results of the experiments. In Figure 7, one may see that both the amplitude return map estimate (ARM) and amplitude-phase estimate (APRM) follow both slow and rapid changes in LLE. When dynamics switch from chaotic to periodic, both return map estimates and LLE do fall. For the Gokyildirim system, the reaction of ARM and APRM estimates is a bit different; they may fall or rise as a reaction to the change in LLE sign. This may be due to the complexity of the dynamics of the system, which, when changing the bifurcation parameter, repeatedly transitions from periodic to chaotic mode and back.

Table 3. Spearman $\rho$ correlation coefficients between dRMA features and LLE.

System	corr (ARM, LLE)	corr (APRM, LLE)
Unified	0.6021	−0.6625
Gokyildirim	0.4729	−0.3737

numerically estimates the correlation between dRMA features and LLE for both systems. We obtained a moderate-to-high negative correlation for the Unified system, and a low-to-moderate correlation for the Gokyildirim system. From this, we may conclude that dRMA can be used as a metric to characterize the non-linear dynamics of the system, taking into account its peculiarities and differences from LLE.

3.2. Estimation of dRMA Sensitivity in Case of Noisy Time Series

Signals from the Unified and Gokyildirim systems were artificially noised using Gaussian noise. In both experiments, signals from state variable x were used. Figure 9 and Figure 10 illustrate the difference in dRMA algorithm output as a function of signal-to-noise ratio. Dynamic regimes presented in Figure 3 and Figure 4 were used. This means that at small values of $\Delta \alpha$ or $\Delta a$, the system modes were quite similar, while at larger values they were significantly different. Values for dRMA features are normalized as follows. For regimes with parameter $\Delta p_i$, values $V_i$ were obtained for every value of SNR, and compared with values $V_0$ obtained for the chaotic system output with the basic parameter value by the following formula:

$$dRM{A}_i={\displaystyle \frac{V_i-V_0}{V_i}}·100\%$$

(17)

As a result, it can be concluded that for almost all values of signal-to-noise ratio greater than 0 dB, the system modes could be distinguished from each other using dRMA features. This indicates the high noise robustness of the proposed estimate, which is of practical importance for applied problems.

3.3. Classifier Settings

The 1-nearest neighbor classifier was used with the distance metric being the Euclidean distance. The predictor data were standardized. No weighting was performed. The tie-breaking algorithm for choosing between classes consists of using the smallest index among tied groups. The nearest neighbor algorithm was set to perform an exhaustive search of all of the distances for the new point in all of the experiments conducted. During the classification, K-Fold cross-validation was applied to the data samples with five folds.

3.4. Case I: Paderborn Dataset Classification

Classification workflow using the proposed dRMA-based and competitive methods is presented in Figure 5. EMD was used to eliminate noise and to obtain useful information from signals. It has been experimentally determined that the intrinsic mode functions 3 and 4 contained the most relevant information. Figure 11 shows the original signal from bearing with code K001 with a rotation speed of 900 rpm, a load torque of 0.7 Nm, and a radial force of 1000 N (setting “N09_M07_F10”) with 3000 time points, and corresponding IMF3 and IMF4.

The distribution of dRMA features for modes 3 and 4 obtained for all analyzed dataset data comprising 3 classes is shown in Figure 12. There are just minor overlaps between points of different classes in the feature space, which allows the use of different types of classifiers. Thus, simple and efficient KNN is a relevant choice.

A generalized feature value for bearings from different classes is shown in Figure 13, using mode 3 as an example. Amplitude-phase features show greater variation, which makes them more significant for classification, but also they can vary notably within a class. This characterizes the complexity of the classification problem.

Sub-signals of 1000, 2000, 5000, 10,000, 30,000, and 50,000 time points were used for the classification of each of the vibration signals. The number of signal samples of each type of damage is 5. The first data sample has always been taken from the starting value of half a sample size. The next data samples were taken from the end of the previous data sample. All of the four operational settings were used. No data augmentation techniques were employed.

The overall classification accuracy and total feature computation time across the dataset are shown in Figure 14.

The best result of the classification was acquired using the dRMA method for 50,000 time points. The achieved accuracy result is 94.96%. Also, this method was faster in almost all experiments, except for the longest time series. Methods based on DWT and GGD achieved only 91.13% accuracy, while entropy-based methods required more than 24 h for calculation features at $n>20,000$ time points using the dedicated workstation. Note that no parallelization mechanisms have been applied to the calculations. Therefore, although the accuracy of the entropy method for 10,000 points outperformed the result of other methods, we refused to run it for a longer time series.

3.5. Case II: Bearing Data Center Dataset Classification

To validate the applicability of the proposed methods on different datasets, we performed the classification of bearing defects using the Bearing Data Center dataset from Case Western Reserve University. The same workflow, as presented in Figure 5, was used. It is a popular dataset, and several referred works, e.g., [19,25], use it for classification as a test bench.

We experimentally determined that the first and second modes contain the most relevant information. Sub-signals of 1000, 2000, 5000, and 10,000-time points were used to classify each of the drive-end accelerometer vibration signals. The ten data samples were used for defective bearings. To balance classes, we used twenty samples for the normal bearings. The samples were taken from the random index of the signal. Two of the operational settings were used: the zero and the two-horsepower motor load. The resulting accuracy and computation time of the methods are shown in Figure 15.

In this test problem, the dRMA method shows the best accuracy at a lower computational cost compared to the other considered methods. The best result for the dRMA method is 100% accuracy obtained with sub-signals of 5000 and 10,000 time points length. The same accuracy result was obtained with the entropy-based method with a sample size of 10,000 data points, but the computational cost for this method is much higher compared to the dRMA method by a factor of more than 500. A method based on DWT coefficients distribution approximation reached an accuracy only of 98%.

4. Conclusions

In this paper, a new algorithm for time-series characterization has been developed using two types of return maps: amplitude and phase-amplitude maps. The method was called distance-based QRMA, or dRMA for short.

The comparison of dRMA and the largest Lyapunov exponent on chaotic time series demonstrates the possibility of dRMA for estimating the parameters of complex dynamics. Also, the proposed method was proven to solve this task in noisy conditions, which makes it a versatile instrument for the characterization of both simulated noiseless and real-world signals.

To test the possibility of using the proposed method in machine-learning tasks, the method based on empirical mode decomposition followed by dRMA-based feature extraction was implemented and applied for the classification of bearing faults. We used vibration signals from the open Paderborn University bearing dataset and Case Western Reserve University Bearing Data Center dataset. As a classifier, we used KNN. In the Paderborn bearing dataset, recordings from both artificially damaged bearings and bearings naturally damaged in accelerated life test experiments were used. Having this, we compared the developed method with the two competing methods: an entropy-based method, and a method based on distribution approximation of wavelet coefficients obtained as a result of discrete wavelet decomposition. We found that the developed method has superior accuracy when more than 10,000 time points are used, and accuracy of classification data into three classes reached 95%. Also, our method was generally faster than competing methods, e.g., outperformed the entropy-based method in computational time by 2–3 orders. In the Case Western dataset, we used three classes: normal bearings, bearings with an outer ring defect, and bearings with an inner ring defect. Two loads were used, corresponding to different shaft rotation speeds. We have shown that the proposed dRMA method has the best accuracy, reaching 100% in classification with the sub-signals of 5000 time points in length. Also, the dRMA method outperformed other considered methods in computational time, by two times for the distribution approximation method, and by 1–2 orders for the entropy method.

With the obtained computational speed, the proposed method can be used in various real-time controller applications demanded in industry, wind turbines, engines and gearbox monitoring, oil pipelines, etc. The disadvantages of the proposed method at the current stage of development include the need for manual selection of the most informative intrinsic mode functions. The automation of this process, as well as the reduction of computational time even more will be the direction of future work.