Improved Synchronous Sampling and Its Application in High-Speed Railway Bearing Damage Detection

: Synchronous analysis is one of the most effective and practical techniques in rotating machinery diagnostics, especially in cases with variable speed operations. A modern analog-to-digital convertor (ADC) usually digitizes an analog signal to an equal time interval data series. Synchronous resampling converts the data series from an equal time interval data series to an equal shaft rotation angle interval data series. This conversion is usually achieved in the digital domain with the aid of shaft speed information, through either direct measurement or identification from a measured vibration signal, which is a time-consuming process. In order to improve the computational efficiency as well as the data processing accuracy, in this paper, a fast synchronous time-point calculation method based on an inverse function interpolation procedure is proposed. By identifying the inverse function of the instantaneous phase with respect to time, the calculation process of synchronous time points is optimized, which results in improved calculation efficiency and accuracy. These advantages are demonstrated by numerical simulations as well as experimental verifications. The numerical simulation results show that the proposed method can improve calculation speed by about five times. The synchronous analysis based on the proposed method was applied to a bearing fault detection in a high-speed rail carriage, which demonstrated the advantages of the proposed algorithm in improving the signal-to-noise ratio (SNR) for bearing damage feature extraction.


Introduction
Rotating machinery is widely used in modern industrial, civil, and military fields.Due to its complex working environment, failures may occur, which will incur increased operation and maintenance costs and the loss of productivity [1].Sometimes it may even lead to catastrophic accidents.Therefore, fault diagnosis and condition monitoring is one of the research foci in rotating machinery [2].Nowadays, damage detection and condition monitoring based on vibration analysis is still the most commonly used method [3,4].
Synchronous analysis is an effective component of rotating machinery vibration signal processing.Usually, in order to diagnose faults of rotating machinery, it is necessary to extract the relevant fault features from the vibration response through signal processing techniques [5].These fault features are commonly associated with the shaft rotation speed.Hence, if the shaft rotating speed varies, the damage features may diminish or even vanish on spectrum due to the smearing effect.Synchronous analysis is a technique developed by integrating shaft speed into vibration signal processing, namely, synchronous resampling.In the synchronous resampling technique, the signal is not discretized in forms of equal time intervals, but in forms of equal shaft rotation angles.Thus, the shaft rotation effect is significantly reduced or eliminated, and the smearing effect caused by variations in shaft speed is reduced or eliminated.Techniques related to synchronous analysis have demonstrated their effectiveness and practicality in rotation machinery damage detection and condition monitoring [6][7][8], especially in cases with variable speed operations.For example, synchronous averaging can effectively separate the information synchronized with the shaft speed [9][10][11][12], such as the shaft speed response, the gear meshing response, and their high-order harmonics.For bearing fault analysis, the application of the virtual shaft concept [13] has been introduced with success.In the technique, a virtual shaft is introduced into signal processing, so that the bearing damage feature is synchronized with the virtual shaft.Thus, the effective synchronous averaging technique can then be performed in the virtual shaft domain, to improve the bearing damage signature signal-to-noise ratio (SNR) [14].Synchronous analysis is also frequently integrated with other bearing fault diagnosis methods, like the popular acceleration enveloping analysis (AEA) [15], and spectral kurtosis analysis [16,17].The AEA method uses a band-pass or high-pass filter to separate the dynamic responses induced by damage.Spectral kurtosis analysis can improve the band-pass filter design by figuring out an optimal frequency range through kurtogram analysis.These techniques will be utilized as a reference baseline in our research as well.More details regarding the utilizing of synchronous analysis technologies for bearing fault diagnosis can be found in the open literature, such as references [18,19].
The prerequisite for synchronous analysis is synchronous sampling.Modern analog-todigital conversion (ADC) is usually achieved by digitizing an analog signal to an equal time interval data series [20].The resulting digital series has an equal time interval, and the time interval is specified by the sampling rate.In the case of constant speed operations, equal delta time data acquisition is equivalent to the equal shaft rotation angle data acquisition.Nevertheless, the sampling is not necessarily the integer multiple of the shaft speed, and therefore a full-cycle sampling within a shaft cycle may not be guaranteed.As a result, in the spectrum, the shaft-related feature frequency may have energy leakage, although it has eliminated the smearing effect.In general, energy leakage can be minimized or even eliminated by increasing the sampling rate [21].For variable speed operations, however, if the equal delta time interval data series is directly converted into frequency domain, the shaft-related features will have the inevitable smearing phenomenon.To avoid it, a synchronization sampling is required.
The conventional synchronization sampling method relies on hardware devices to generate an equidistant sampling clock.This clock is subsequently employed to regulate the data acquisition process.Since such a sampling clock is determined by hardware, the sampling precision and resolution is limited.With the development of computer hardware and software technologies, synchronous resampling is currently achieved by reprocessing the equal delta time interval data series in the digital domain [22].At present, synchronous acquisition methods based on digital processing can be categorized into two groups.One is the hybrid synchronous analysis method, which uses a piece of hardware to record the shaft rotation pulse train.The recorded speed signal is then used to assist the synchronous resampling through a designed software [23][24][25].This process is known as computed order tracking (COT).The other group is based on the instantaneous frequency estimation from the vibration response, without a tachometer or a key-phasor installed on the shaft [26,27].
For a given function of time, mathematically, the synchronous resampling process is used to find the independent variable values for given dependent variable values.Usually, this is an iteration process, and therefore it consumes a considerable amount of CPU time.In advanced rotating machinery tests, such as the high-speed railway bench test [28,29], it takes a long time for the machine to experience many different operational modes during the testing procedure.Any uneconomical processing steps may prevent the application of the process to the entire time range.On the other hand, in quasi real-time analysis, it is always beneficial to minimize the processing time in each and every time step [30].Therefore, it is important to further reduce the processing time while maintaining high accuracy in data resampling for synchronous analysis.Towards improving the efficiency and accuracy of synchronous signal processing, this paper presents a fast yet accurate synchronization time-point calculation method based on the inverse function interpolation.This approach achieves the computation of timestamps corresponding to the equal rotation angles of the shaft through a process of evaluating dependent variables for given independent variables.Notably, this method avoids an iterative process, resulting in a significant reduction in computational time costs.Analysis of the numerical simulations and real examples demonstrated the feasibility, effectiveness, and practicality of the proposed approach.
The rest of this paper is structured as follows: Section 2 details the synchronous resampling method based on inverse function interpolation.Section 3 verifies the feasibility and advantages of this methodology through numerical simulation.Section 4 verifies the proposed methodology with two practical cases of railway bearing damage detection.Finally, Section 5 provides the conclusion.

Basic Principles of Synchronous Resampling
Modern synchronous resampling is mostly carried out in the digital domain, as shown in Figure 1.The vibration response signal and key-phasor signal are obtained from a rotating machine through sensors.They are then discretized by a high-frequency and highresolution ADC with an equal delta time interval simultaneously.For the convenience of explanation, the vibration response signal is represented by the shaft fundamental response only.A key-phasor signal is a signal correlated with the rotational phase of a mechanical shaft in a rotating system.The key-phasor signal is usually measured by a proximity probe with respect to the key on the shaft, and then converted to the format of transistortransistor logic (TTL) pulses.The pulses are generated at a rate of one pulse per shaft rotation cycle.Then, from the digitized key-phasor signal, the shaft speed 1/rev time mark and the synchronous timestamps for the equal delta shaft angle are determined numerically.In cases where the shaft speed remains constant, synchronous timestamps are uniformly distributed along the time axis.Conversely, for a variable-speed case, these synchronous timestamps will be uneven on the time axis and need to be evaluated accordingly.Following synchronous timestamps, the synchronous vibration response can be obtained through a numerical processing scheme such as interpolation.
always beneficial to minimize the processing time in each and every time step [30].Therefore, it is important to further reduce the processing time while maintaining high accuracy in data resampling for synchronous analysis.
Towards improving the efficiency and accuracy of synchronous signal processing, this paper presents a fast yet accurate synchronization time-point calculation method based on the inverse function interpolation.This approach achieves the computation of timestamps corresponding to the equal rotation angles of the shaft through a process of evaluating dependent variables for given independent variables.Notably, this method avoids an iterative process, resulting in a significant reduction in computational time costs.Analysis of the numerical simulations and real examples demonstrated the feasibility, effectiveness, and practicality of the proposed approach.
The rest of this paper is structured as follows: Section 2 details the synchronous resampling method based on inverse function interpolation.Section 3 verifies the feasibility and advantages of this methodology through numerical simulation.Section 4 verifies the proposed methodology with two practical cases of railway bearing damage detection.Finally, Section 5 provides the conclusion.

Basic Principles of Synchronous Resampling
Modern synchronous resampling is mostly carried out in the digital domain, as shown in Figure 1.The vibration response signal and key-phasor signal are obtained from a rotating machine through sensors.They are then discretized by a high-frequency and high-resolution ADC with an equal delta time interval simultaneously.For the convenience of explanation, the vibration response signal is represented by the shaft fundamental response only.A key-phasor signal is a signal correlated with the rotational phase of a mechanical shaft in a rotating system.The key-phasor signal is usually measured by a proximity probe with respect to the key on the shaft, and then converted to the format of transistor-transistor logic (TTL) pulses.The pulses are generated at a rate of one pulse per shaft rotation cycle.Then, from the digitized key-phasor signal, the shaft speed 1/rev time mark and the synchronous timestamps for the equal delta shaft angle are determined numerically.In cases where the shaft speed remains constant, synchronous timestamps are uniformly distributed along the time axis.Conversely, for a variable-speed case, these synchronous timestamps will be uneven on the time axis and need to be evaluated accordingly.Following synchronous timestamps, the synchronous vibration response can be obtained through a numerical processing scheme such as interpolation.Digital-domain synchronous resampling consists of three major steps.The first step is to determine the time mark of every full-cycle shaft rotation from the digitized keyphasor signal.The second step is to discern the timestamps corresponding to each equal delta shaft angle point.The last step is to utilize the calculated timestamps to resample the equal delta-time digitized signal to obtain the synchronous response data series.Digital-domain synchronous resampling consists of three major steps.The first step is to determine the time mark of every full-cycle shaft rotation from the digitized key-phasor signal.The second step is to discern the timestamps corresponding to each equal delta shaft angle point.The last step is to utilize the calculated timestamps to resample the equal delta-time digitized signal to obtain the synchronous response data series.
In general, the first step can be obtained by identifying the trigger position of the shaft key-phasor pulse train, or the shaft whole-revolution time location from the shaft instantaneous speed.While the third step can be achieved by a general one-dimensional interpolation algorithm (such as the interp1 function in MATLAB).The second step is more involved, especially if the shaft speed changes substantially.
As shown in Figure 2a, assuming that t 1 and t 2 are the time marks of two consecutive key-phasor pulses (or shaft whole-revolution time marks), the corresponding shaft speeds at t 1 and t 2 are f 1 and f 2 , respectively.The unit of shaft speed in this paper is revolution per second, r/s.For a variable speed operation, f 1 ̸ = f 2 .According to the definition of a complete shaft rotation, we have In the most general case, assuming the number of synchronous resampling points is N, we define the timestamps for equal delta angle points as t 1,i , i = 0, 1, . . ., N, and t 10 = t 1 , t 1,N = t 2 .The determination of the discrete time point t 1,i , i = 1, . . ., N − 1, can be solved one by one by the following equation For any f (t) > 0, it is difficult to obtain the analytical solution directly from Equation (2).To solve Equation ( 2), a time-consuming iterative process is usually needed.
If f (t) can be approximated by a linear function within the time interval [t 1 , t 2 ], as shown in Figure 2b, then , and Equation ( 2) can be rewritten as By simplifying Equation ( 4), we can get Solving Equation ( 5) and retaining the real solution within t ⊂ [t 1,i−1 , t 2 ], we can obtain the synchronous sampling discrete time point t 1,i , i = 1, 2, . . ., N − 1.
In general, the first step can be obtained by identifying the trigger position of the shaft key-phasor pulse train, or the shaft whole-revolution time location from the shaft instantaneous speed.While the third step can be achieved by a general one-dimensional interpolation algorithm (such as the interp1 function in MATLAB).The second step is more involved, especially if the shaft speed changes substantially.
As shown in Figure 2a, assuming that  and  are the time marks of two consecutive key-phasor pulses (or shaft whole-revolution time marks), the corresponding shaft speeds at  and  are  and  , respectively.The unit of shaft speed in this paper is revolution per second, /.For a variable speed operation,  ≠  .According to the definition of a complete shaft rotation, we have In the most general case, assuming the number of synchronous resampling points is , we define the timestamps for equal delta angle points as  , ,  = 0, 1, … , , and  =  ,  , =  .The determination of the discrete time point  , ,  = 1, … ,  − 1, can be solved one by one by the following equation For any () > 0, it is difficult to obtain the analytical solution directly from Equation (2).To solve Equation ( 2), a time-consuming iterative process is usually needed.

Synchronous Time-Point Calculation Based on Inverse Function Interpolation
Based on the above analysis, it is observed that mathematically Equation ( 2) has provided a solution for synchronous resampling.However, since numerically solving for a synchronous timestamp from Equation ( 2) involves an iterative process.It becomes time-consuming to apply Equation ( 2) directly in the case of long-time testing, or in the case of a quasi-real-time data processing requirement.With a constant shaft-speed assumption within a cycle, the solution for a synchronous timestamp changes from an iteration process, to solving for the roots of the first-order polynomial equation.This assumption has mitigated the computational time for synchronous resampling to a certain extent.However, a side effect is that the assumption of constant speed within a cycle may introduce calculated errors in synchronous timestamps, particularly in cases with significant variations in shaft speed.
In order to further speed up the resampling process, while retaining the shaft-speed variation fidelity in synchronous analysis, we have developed an approach to locating the synchronous timestamps based on the inverse function interpolation.In this technique, the instantaneous shaft speed is time integrated to obtain the instantaneous phase y = ϕ(t), as shown in Figure 3a.Assuming the shaft is rotating only in a single direction, y = ϕ(t) is monotonic in t ∈ [0, T].Here we assume ϕ(0) = 0 and ϕ(T) = Φ for simplicity of statement.The inverse function t = φ(y) exists, and is also monotonic in y ∈ [0, Φ] as shown in Figure 3b.

Synchronous Time-Point Calculation Based on Inverse Function Interpolation
Based on the above analysis, it is observed that mathematically Equation ( 2) has provided a solution for synchronous resampling.However, since numerically solving for a synchronous timestamp from Equation ( 2) involves an iterative process.It becomes timeconsuming to apply Equation (2) directly in the case of long-time testing, or in the case of a quasi-real-time data processing requirement.With a constant shaft-speed assumption within a cycle, the solution for a synchronous timestamp changes from an iteration process, to solving for the roots of the first-order polynomial equation.This assumption has mitigated the computational time for synchronous resampling to a certain extent.However, a side effect is that the assumption of constant speed within a cycle may introduce calculated errors in synchronous timestamps, particularly in cases with significant variations in shaft speed.
In order to further speed up the resampling process, while retaining the shaft-speed variation fidelity in synchronous analysis, we have developed an approach to locating the synchronous timestamps based on the inverse function interpolation.In this technique, the instantaneous shaft speed is time integrated to obtain the instantaneous phase  = (), as shown in Figure 3a In the instantaneous phase domain, the physical meaning of the synchronous timestamps is easily understood.As seen in Figure 3c, for a given full shaft rotation cycle,  and  are the time marks of two adjacent key-phasor pulses (relative whole-revolution time marks in the tacholess analysis).The corresponding instantaneous phases at time  and  are  and  , respectively.Synchronous resampling means dividing  segments equally between  and  to form  , ,  = 0,1, 2, … ,  , where  = y and In the instantaneous phase domain, the physical meaning of the synchronous timestamps is easily understood.As seen in Figure 3c, for a given full shaft rotation cycle, t 1 and t 2 are the time marks of two adjacent key-phasor pulses (relative whole-revolution time marks in the tacholess analysis).The corresponding instantaneous phases at time t 1 and t 2 are y 1 and y 2 , respectively.Synchronous resampling means dividing N segments equally between y 1 and y 2 to form y 1,i , i = 0, 1, 2, . . ., N, where y 10 = y 1 and y 1N = y 2 , and then determine t 1,i , i = 0, 1, . . ., N, where t 10 = t 1 and t 1,N = t 2 .In other words, the synchronous resampling process is to find the independent variable points for given Machines 2024, 12, 101 6 of 15 dependent variable points in the instantaneous phase function.The discrete timestamps, t 1,i , i = 1, . . ., N − 1, can be solved one by one according to the following equation Solving Equation (6) usually requires an iterative process, which is expected to be a time-consuming numerical operation.With the inverse function known, as shown in Figure 3d, the synchronous resampling timestamps are simply dependent variable data points evaluated at the inverse function by given independent variable data points y 1,i , i = 0, 1, 2, . . ., N. For given the synchronous resampling timestamps are evaluated at the inverse function t = φ(y); and therefore, In short, in numerical operations, evaluating the independent variable value for a given function value is often a time-consuming iterative process.Evaluating the function value for a given independent variable is a much simpler function valuation or interpolation process, and thus is much less time consuming.In rotating machinery signal processing, as long as the rotation direction of the shaft is unchanged, the instantaneous phase is a monotonic function with time.Its inverse function exists and is also monotonic.Therefore, we can always use the inverse function of the instantaneous phase to evaluate the synchronous resampling timestamps to improve data processing speed.In addition, since there is no approximation in the inverse function evaluation process, it is expected that the synchronous resampling based on the inverse function of the instantaneous phase will be more accurate.

Implementation of Synchronous Resampling Based on Inverse Function Interpolation
Synchronous resampling based on inverse function interpolation can be implemented as shown in Figure 4: 1.
Signal acquisition.It is assumed that both the vibration signal and the key-phasor signal are obtained simultaneously and are digitized at a high sampling frequency with high resolution, as shown in Figure 1.

2.
Instantaneous shaft speed recognition.If the shaft speed information is provided in the form of a pulse train, the corresponding instantaneous shaft speed can be obtained by conventional signal processing [25].The instantaneous shaft speed can also be extracted by processing the vibration signal if the shaft speed information is not available [26,27].

3.
Instantaneous phase recognition.Relative instantaneous phase function is obtained by numerical integration of the instantaneous shaft speed.Without losing the generality, the phase is set to 0 at initial time.The obtained relation between the phase of shaft rotation and time, y = ϕ(t), is as shown in Figure 3a.

4.
Inverse function conversion of the instantaneous phase.Find the inverse function of the instantaneous phase identified in step 3, as shown in Figure 3b.

5.
Shaft synchronous resampling clock evaluation.Using the interpolation algorithm on the inverse function t = φ(y), the synchronous resampling timestamps t i , are determined quickly by evaluation from the inverse function for given y i .6.
Vibration response resampling.According to the synchronous sampling clock, the vibration response signal is resampled synchronously to obtain the synchronous response signal in the equal shaft angle domain.

7.
Synchronous analysis and synchronous averaging.The synchronously resampled vibration response signal is analyzed and the results are used to diagnose the health condition of the rotation machinery.
Machines 2024, 12, x FOR PEER REVIEW 7 of 16 6. Vibration response resampling.According to the synchronous sampling clock, the vibration response signal is resampled synchronously to obtain the synchronous response signal in the equal shaft angle domain.7. Synchronous analysis and synchronous averaging.The synchronously resampled vibration response signal is analyzed and the results are used to diagnose the health condition of the rotation machinery.

Validation
This section verifies the feasibility and advantages of the proposed method through numerical simulations.Since data processing time is computer hardware-and softwaredependent, the following lists the computer environment used in the simulations: Windows 10 Professional 64-bit operating system; MATLAB R2021b; 16 11th Gen Intel(R) Core (TM) i7-11700K @ 3.60GHz; and 16 GB RAM.
In the numerical simulation, it was assumed that the shaft speed varies, as shown in Figure 5a.In actual situations, the shaft speed information is obtained from the keyphasor, or evaluated from the measured vibration signal through instantaneous frequency identification techniques [31].

Validation
This section verifies the feasibility and advantages of the proposed method through numerical simulations.Since data processing time is computer hardware-and softwaredependent, the following lists the computer environment used in the simulations: Windows 10 Professional 64-bit operating system; MATLAB R2021b; 16 11th Gen Intel(R) Core (TM) i7-11700K @ 3.60 GHz; and 16 GB RAM.
In the numerical simulation, it was assumed that the shaft speed varies, as shown in Figure 5a.In actual situations, the shaft speed information is obtained from the keyphasor, or evaluated from the measured vibration signal through instantaneous frequency identification techniques [31].
6. Vibration response resampling.According to the synchronous sampling clock, the vibration response signal is resampled synchronously to obtain the synchronous response signal in the equal shaft angle domain.7. Synchronous analysis and synchronous averaging.The synchronously resampled vibration response signal is analyzed and the results are used to diagnose the health condition of the rotation machinery.

Validation
This section verifies the feasibility and advantages of the proposed method through numerical simulations.Since data processing time is computer hardware-and softwaredependent, the following lists the computer environment used in the simulations: Windows 10 Professional 64-bit operating system; MATLAB R2021b; 16 11th Gen Intel(R) Core (TM) i7-11700K @ 3.60GHz; and 16 GB RAM.
In the numerical simulation, it was assumed that the shaft speed varies, as shown in Figure 5a.In actual situations, the shaft speed information is obtained from the keyphasor, or evaluated from the measured vibration signal through instantaneous frequency identification techniques [31].For simplicity, the vibration response simulated here only contained the first-order response of the shaft rotation.Both the speed signal and the vibration response signal were digitized at 10,000 Hz.The signals lasted for 1000 s.The first 10 s of the signal time history is shown in Figure 5b.
Before synchronous resampling, the instantaneous phase of the shaft was obtained by numerical integration of the identified shaft speed, as shown in Figure 5c.The corresponding inverse function is shown in Figure 5d.According to the highest shaft order required in the analysis, the synchronous resampling points within a shaft cycle were determined.In other words, the independent variable discretization points were specified for the inverse function of the instantaneous phase.An interpolation method can be applied to calculate the corresponding synchronous resampling timestamps.With the synchronous resampling timestamps, the equal delta time vibration response was converted to equal delta shaft angle synchronous vibration response by an interpolation algorithm.
The first four cycles of synchronous responses are shown in Figure 6.The performance of three synchronous resampling methods is listed in Table 1, where the accuracy was measured by the mean squared error (MSE): where  is the theoretical value of simulation, and  is the estimated value with different method.For simplicity, the vibration response simulated here only contained the first-order response of the shaft rotation.Both the speed signal and the vibration response signal were digitized at 10,000 Hz.The signals lasted for 1000 s.The first 10 s of the signal time history is shown in Figure 5b.
Before synchronous resampling, the instantaneous phase of the shaft was obtained by numerical integration of the identified shaft speed, as shown in Figure 5c.The corresponding inverse function is shown in Figure 5d.According to the highest shaft order required in the analysis, the synchronous resampling points within a shaft cycle were determined.In other words, the independent variable discretization points were specified for the inverse function of the instantaneous phase.An interpolation method can be applied to calculate the corresponding synchronous resampling timestamps.With the synchronous resampling timestamps, the equal delta time vibration response was converted to equal delta shaft angle synchronous vibration response by an interpolation algorithm.
The first four cycles of synchronous responses are shown in Figure 6.The performance of three synchronous resampling methods is listed in Table 1, where the accuracy was measured by the mean squared error (MSE): where Y i is the theoretical value of simulation, and Ŷi is the estimated value with different method.
Machines 2024, 12, x FOR PEER REVIEW 8 of 16 (c) (d) For simplicity, the vibration response simulated here only contained the first-order response of the shaft rotation.Both the speed signal and the vibration response signal were digitized at 10,000 Hz.The signals lasted for 1000 s.The first 10 s of the signal time history is shown in Figure 5b.
Before synchronous resampling, the instantaneous phase of the shaft was obtained by numerical integration of the identified shaft speed, as shown in Figure 5c.The corresponding inverse function is shown in Figure 5d.According to the highest shaft order required in the analysis, the synchronous resampling points within a shaft cycle were determined.In other words, the independent variable discretization points were specified for the inverse function of the instantaneous phase.An interpolation method can be applied to calculate the corresponding synchronous resampling timestamps.With the synchronous resampling timestamps, the equal delta time vibration response was converted to equal delta shaft angle synchronous vibration response by an interpolation algorithm.
The first four cycles of synchronous responses are shown in Figure 6.The performance of three synchronous resampling methods is listed in Table 1, where the accuracy was measured by the mean squared error (MSE): where  is the theoretical value of simulation, and  is the estimated value with different method.In this simulation, for all the methods tested, the accuracies of the synchronous resampling were very high as expected, due to the mild speed variation in this simulation.However, the data processing efficiency of the inverse function evaluation-based method has been significantly improved.As seen in Table 1, the processing time of the inversefunction-based method was about 37% of that used by the constant shaft-speed-based method, and was only about 17% of that used by the conventional iteration-based method.
In order to further verify the accuracy of the proposed method, a case with more drastic shaft speed variation was simulated.The frequency variation is assumed to be: where f min and f max are the frequency range and α is associated to the frequency change rate.The larger the α, the shorter the time required for the shaft rotation frequency to change from f min to f max .A case with f min = 20, f max = 100, and α = 100 is shown in the top portion of Figure 7, and the corresponding vibration is shown in the bottom portion of Figure 7.
Machines 2024, 12, x FOR PEER REVIEW 9 of 16 In this simulation, for all the methods tested, the accuracies of the synchronous resampling were very high as expected, due to the mild speed variation in this simulation.However, the data processing efficiency of the inverse function evaluation-based method has been significantly improved.As seen in Table 1, the processing time of the inversefunction-based method was about 37% of that used by the constant shaft-speed-based method, and was only about 17% of that used by the conventional iteration-based method.
In order to further verify the accuracy of the proposed method, a case with more drastic shaft speed variation was simulated.The frequency variation is assumed to be: where  Digital signal processing was carried out with the same procedures as the previous example.The accuracy comparison of different resampling methods is presented in Figure 8a.A zoomed version for the rapid speed change portion is shown in Figure 8b.According to the simulation results, we were able to observe that, in the range of rapid shaft speed change, the method proposed in this study indeed rendered higher accuracy than that of the approximation method with the assumption of constant speed within a shaft cycle.For the current simulation case, the maximum error rate (the maximum amplitude difference vs. the simulated amplitude) of the constant frequency assumption (red line) was as high as 55.5%, while the maximum error rate from the proposed inverse function-based processing (blue dashed line) was less than 9%.The maximum errors for the different changing rate  were also evaluated and listed in Table 2.As expected, the higher the speed change rate, the more maximum error was introduced in the constant speed-based Digital signal processing was carried out with the same procedures as the previous example.The accuracy comparison of different resampling methods is presented in Figure 8a.A zoomed version for the rapid speed change portion is shown in Figure 8b.According to the simulation results, we were able to observe that, in the range of rapid shaft speed change, the method proposed in this study indeed rendered higher accuracy than that of the approximation method with the assumption of constant speed within a shaft cycle.For the current simulation case, the maximum error rate (the maximum amplitude difference vs. the simulated amplitude) of the constant frequency assumption (red line) was as high as 55.5%, while the maximum error rate from the proposed inverse function-based processing (blue dashed line) was less than 9%.The maximum errors for the different changing rate α were also evaluated and listed in Table 2.As expected, the higher the speed change rate, the more maximum error was introduced in the constant speed-based method.The inverse function-based method accomplished consistently better results than those from the constant speed-based method.
method.The inverse function-based method accomplished consistently better results than those from the constant speed-based method.

Applications on Railway Bearing Damage Detection
In the above section, the computational efficiency and accuracy of the improved method was demonstrated with numerical simulations.In this section, the proposed method was used to process the data from a high-speed railway test stand to carry out embedded bearing damage detection.

Test Setup
A schematic diagram of a high-speed railway carriage and bearing locations is shown in Figure 9a, and the acceleration sensor locations are shown in Figure 9b.Amongst the eight bearings, some were preset with specific damages.During testing, the acceleration vibration response signals and the shaft key-phasor signals were recorded at the sampling rate of 25.6 kHz simultaneously for subsequent analysis.The data segment used for this demonstration was 20 s.
The test for the bearing geometry parameters was as follows: the number of rollers was 19; the pitch diameter was 183.929 mm; the roller diameter was 26 mm; and the contact angle was 10 deg.Based on these parameters, the bearing damage characteristic orders were calculated as follows: the inner race damage order was 10.822; the outer race damage order was 8.177; the roller damage order was 3.468; and the cage damage order was 0.429.The detailed results are shown in Table 3.

Applications on Railway Bearing Damage Detection
In the above section, the computational efficiency and accuracy of the improved method was demonstrated with numerical simulations.In this section, the proposed method was used to process the data from a high-speed railway test stand to carry out embedded bearing damage detection.

Test Setup
A schematic diagram of a high-speed railway carriage and bearing locations is shown in Figure 9a, and the acceleration sensor locations are shown in Figure 9b.Amongst the eight bearings, some were preset with specific damages.During testing, the acceleration vibration response signals and the shaft key-phasor signals were recorded at the sampling rate of 25.6 kHz simultaneously for subsequent analysis.The data segment used for this demonstration was 20 s.
The test for the bearing geometry parameters was as follows: the number of rollers was 19; the pitch diameter was 183.929 mm; the roller diameter was 26 mm; and the contact angle was 10 deg.Based on these parameters, the bearing damage characteristic orders were calculated as follows: the inner race damage order was 10.822; the outer race damage order was 8.177; the roller damage order was 3.468; and the cage damage order was 0.429.The detailed results are shown in Table 3.

Extraction of Bearing Outer Race Spalling Damage Features
In this case, the bearing outer race was artificially embedded with spalling damage as shown in Figure 9c. Figure 10a shows the shaft key-phasor signal and the acceleration response signal at the bearing seat when the test rig was running.Figure 10b shows the shaft speed time history identified from the key-phasor signal, where it was clearly shown that the shaft speed was a time-varying function.
Before carrying out the AEA, we used spectrum kurtosis analysis to determine the optimal frequency range [16,17].Under this operating condition, through spectrum kurtosis analysis, the frequency range of [4475, 5125] Hz had the highest kurtosis.This frequency range was used for AEA analysis.
The AEA results based on the spectrum kurtosis analysis are shown in Figure 10c.With identified shaft speed and theoretical damage characteristic orders, the outer race damage feature frequency was around [6.542, 10.072] Hz.In the AEA spectrum, as seen in Figure 10c, the damage signature was only barely visible due to the frequency smearing effect caused by the variable speed operation.After processing using the proposed method, the AEA order spectrum was obtained as shown in Figure 10d.The outer race damage feature was significantly enhanced.

Extraction of Bearing Outer Race Spalling Damage Features
In this case, the bearing outer race was artificially embedded with spalling damage as shown in Figure 9c. Figure 10a shows the shaft key-phasor signal and the acceleration response signal at the bearing seat when the test rig was running.Figure 10b shows the shaft speed time history identified from the key-phasor signal, where it was clearly shown that the shaft speed was a time-varying function.
Before carrying out the AEA, we used spectrum kurtosis analysis to determine the optimal frequency range [16,17].Under this operating condition, through spectrum kurtosis analysis, the frequency range of [4475, 5125] Hz had the highest kurtosis.This frequency range was used for AEA analysis.
The AEA results based on the spectrum kurtosis analysis are shown in Figure 10c.With identified shaft speed and theoretical damage characteristic orders, the outer race damage feature frequency was around [6.542, 10.072] Hz.In the AEA spectrum, as seen in Figure 10c, the damage signature was only barely visible due to the frequency smearing effect caused by the variable speed operation.After processing using the proposed method, the AEA order spectrum was obtained as shown in Figure 10d.The outer race damage feature was significantly enhanced.
To qualify the improvement, we defined an indicator based on the damage feature energy as: where ENV 2 (Ω i,damage ) denotes the energy of the ith order of the damage feature, and mean(ENV 2 (Ω)) denotes the averaged energy of the frequency (order) range of interest.The defined indicator is the ratio of the damage feature energy with respect to the averaged energy within the range of interest.The higher the energy ratio E amp is, the easier the damage can be identified from the background noise.To qualify the improvement, we defined an indicator based on the damage feature energy as: where  (Ω , ) denotes the energy of the  order of the damage feature, and ( (Ω) denotes the averaged energy of the frequency (order) range of interest.The defined indicator is the ratio of the damage feature energy with respect to the averaged energy within the range of interest.The higher the energy ratio  is, the easier the damage can be identified from the background noise.
In this example, we used the order range of [5,30], which approximately includes the first three orders of the outer race damage feature, excluding the influence of shaft order responses.The calculated results are listed in Table 4.As shown in Table 4, the results based on synchronous sampling had higher  value, which is consistent with the visual observations from the spectra as shown in Figure 10c,d.In this example, we used the order range of [5,30], which approximately includes the first three orders of the outer race damage feature, excluding the influence of shaft order responses.The calculated results are listed in Table 4.As shown in Table 4, the results based on synchronous sampling had higher E amp value, which is consistent with the visual observations from the spectra as shown in Figure 10c,d.

Extraction of Bearing Outer Race Indentation Damage Features
In this case, the outer race was imbedded with an indentation fault, as shown in Figure 9d. Figure 11a shows the shaft pulse signal and the acceleration response signal from the bearing seat when the test rig was in operation.Figure 11b shows the identified shaft speed time history.

Extraction of Bearing Outer Race Indentation Damage Features
In this case, the outer race was imbedded with an indentation fault, as shown in Fig- ure 9d. Figure 11a shows the shaft pulse signal and the acceleration response signal from the bearing seat when the test rig was in operation.Figure 11b shows the identified shaft speed time history.Similar to the previous case, the segment of vibration signal had the highest kurtosis at a frequency range of [7000, 9000] Hz.This frequency range was used for AEA analysis and the results are shown in Figure 11c.The proposed method was then employed to reduce the influence of speed variation, and the corresponding AEA order spectrum was obtained, as shown in Figure 11d.By comparing the results in Figure 11c,d, it can easily be visualized that the proposed method significantly improved the SNR of the damage features.Similar to the previous case, the segment of vibration signal had the highest kurtosis at a frequency range of [7000, 9000] Hz.This frequency range was used for AEA analysis the results are shown in Figure 11c.The proposed method was then employed to reduce the influence of speed variation, and the corresponding AEA order spectrum was obtained, as shown in Figure 11d.By comparing the results in Figure 11c,d, it can easily be visualized that the proposed method significantly improved the SNR of the damage features.
Similarly, the energy ratios were calculated and the results are listed in Table 5.From the results in Table 5, it can be seen that they are consistent with the visual observations from Figure 11c,d.Both the numerical simulation and experimental results indicate that the proposed synchronous processing method can significantly improve efficiency.When applied to a variable shaft speed situation, it can significantly improve the bearing damage feature extraction SNR.

Figure 1 .
Figure 1.A digital method for synchronous resampling.

Figure 1 .
Figure 1.A digital method for synchronous resampling.

Figure 2 .
Figure 2. Assumption of instantaneous shaft speed: (a) arbitrary variable shaft speed within a shaft cycle; and (b) locally linearized shaft speed within a shaft cycle.Figure 2. Assumption of instantaneous shaft speed: (a) arbitrary variable shaft speed within a shaft cycle; and (b) locally linearized shaft speed within a shaft cycle.

Figure 2 .
Figure 2. Assumption of instantaneous shaft speed: (a) arbitrary variable shaft speed within a shaft cycle; and (b) locally linearized shaft speed within a shaft cycle.Figure 2. Assumption of instantaneous shaft speed: (a) arbitrary variable shaft speed within a shaft cycle; and (b) locally linearized shaft speed within a shaft cycle.

Figure 6 .
Figure 6.Synchronous resampling result based on inverse function interpolation.

Figure 5 .
Figure 5. Simulation and analysis with mild speed variations: (a) simulated shaft speed; (b) shaft key-phasor and vibration signals; (c) instantaneous rotation cycle (phase); and (d) inverse function of instantaneous phase.

Figure 5 .
Figure 5. Simulation and analysis with mild speed variations: (a) simulated shaft speed; (b) shaft key-phasor and vibration signals; (c) instantaneous rotation cycle (phase); and (d) inverse function of instantaneous phase.

Figure 6 .
Figure 6.Synchronous resampling result based on inverse function interpolation.Figure 6. Synchronous resampling result based on inverse function interpolation.

Figure 6 .
Figure 6.Synchronous resampling result based on inverse function interpolation.Figure 6. Synchronous resampling result based on inverse function interpolation.
and  are the frequency range and  is associated to the frequency change rate.The larger the , the shorter the time required for the shaft rotation frequency to change from  to  .A case with  = 20,  = 100, and  = 100 is shown in the top portion of Figure 7, and the corresponding vibration is shown in the bottom portion of Figure 7.

Figure 7 .
Figure 7.A case with drastic change of speed and corresponding vibration response.

Figure 7 .
Figure 7.A case with drastic change of speed and corresponding vibration response.

Figure 8 .
Figure 8. Simulation error result with more drastic speed variations (a) resampling error; and (b) resampling error, zoomed version.

Figure 8 .
Figure 8. Simulation error result with more drastic speed variations (a) resampling error; and (b) resampling error, zoomed version.

Figure 9 .
Figure 9. Schematic of carriage test stand and bearing damage: (a) schematic of carriage and bearing locations; (b) sensor locations; (c) bearing damage-artificial spalling; and (d) bearing damageindentation.

Figure 9 .
Figure 9. Schematic of carriage test stand and bearing damage: (a) schematic of carriage and bearing locations; (b) sensor locations; (c) bearing damage-artificial spalling; and (d) bearing damageindentation.

Figure 10 .
Figure 10.AEA analyses based on proposed method: (a) shaft speed key-phasor and vibration signals; (b) identified shaft speed; (c) conventional AEA spectrum; and (d) synchronized AEA order spectrum.

Figure 10 .
Figure 10.AEA analyses based on proposed method: (a) shaft speed key-phasor and vibration signals; (b) identified shaft speed; (c) conventional AEA spectrum; and (d) synchronized AEA order spectrum.

Figure 11 .
Figure 11.AEA analyses based on proposed method: (a) shaft speed key-phasor and vibration signals; (b) identified shaft speed; (c) conventional AEA spectrum; and (d) synchronized AEA order spectrum.

Figure 11 .
Figure 11.AEA analyses based on proposed method: (a) shaft speed key-phasor and vibration signals; (b) identified shaft speed; (c) conventional AEA spectrum; and (d) synchronized AEA order spectrum.

Table 2 .
The maximum errors for different changing rate

Table 2 .
The maximum errors for different changing rate α.

Table 3 .
Test bearing parameters and calculated damage order.

Table 3 .
Test bearing parameters and calculated damage order.

Table 4 .
Improvement of damage feature amplitude on spectrum.

Table 5 .
Improvement in damage feature amplitude on spectrum.