Multiple instantaneous frequency ridge based integration strategy for bearing fault diagnosis under variable speed operations

Rolling element bearings are one of the key elements used in rotating machines. Their failure will result in system breakdowns and cause unexpected accidents. The health condition monitoring of bearings is, therefore, an emerging discipline to scientifically manage machine lifetime. Instantaneous frequency (IF) extraction and IF-based resampling consist of the major tasks used in conventional approaches for bearing fault diagnosis under variable speeds, where it is often desirable for the IF to be extracted based on the time-frequency analysis (TFA) of the vibration signal, instead of using a tachometer. However, an accurate IF extraction based on vibration signals from incipient bearing faults is often undermined by poor energy concentration and the weak readability of the TFA. Furthermore, the resampling procedure for bearing fault diagnosis is error-prone and vulnerable to noise. An approach based on multiple IF ridge integration is, therefore, proposed to address such problems. The proposed approach is dedicated to an accurate IF estimation and bearing fault diagnosis without tachometer utilization and resampling involvement. It is mainly comprised of four steps: (1) acquire multiple pre-IF ridges via a regional peak search algorithm (RPSA) from time frequency representations (TFRs); (2) integrate pre-IF ridges based on the probability density function (PDF) to gain an accurate IF estimation; (3) rectify the multiple pre-IF ridges; (4) diagnose bearing health condition according to the average ratios of any two rectified IF ridges, i.e. fault characteristic coefficient (FCC) or FCC-related numbers. Then, the IF can be accurately estimated based only on the collected vibration signals, and bearing fault diagnosis under variable speed conditions can be implemented. Numerical simulations and experimental signal analyses validate the effectiveness of the proposed method.


Introduction
Rolling element bearings are one of the most important elements in modern industrial systems such as wind turbines, thus their healthy condition is crucial to the entire mechanical system operation. However, bearing-related failures frequently appear in systems due to the complicated environment, resulting in large losses and even casualties [1,2]. Therefore, the prognostics and health management (PHM) of bearings is extremely important for maintaining the normal functions of wind turbines [3]. To achieve the PHM of bearings, their health condition and fault diagnosis are indispensable as they provide a health indicator for bearings [4,5]. When a local fault occurs in rolling bearings, a collision occurs between this fault and its mating surface; an impulse is, hence, excited every time [6]. Under constant shaft rotational speed, the impulse appears with a constant frequency referred to as fault characteristic frequency (FCF), which is a crucial feature related to the bearing fault type [7,8]. For different bearing components, i.e. outer race and inner race, there are different FCFs expressed as follows, respectively: where BPFO and BPFI stand for the frequency of the ball passing outer race and inner race, respectively, n is the number of rolling elements, f r is the shaft rotational frequency, ψ is the contact angle, D and d represent the pitch diameter and the diameter of the rolling element respectively, FCC 0 and FCC I represent the FCC for the outer race fault and inner race fault, respectively. The FCF could be detected for bearing fault diagnosis under a constant speed condition [9][10][11]. Nevertheless, in reality, bearings often operate at variable speed operation, where FCF accordingly changes over time and is referred to as instantaneous fault characteristics frequency (IFCF). Since IFCF varies with time, the capability of various approaches devised for bearing fault diagnostics under constant speeds is impaired [12,13]. Thus, technologies for bearing fault diagnosis under variable speed operation need to be developed [14][15][16][17].
Order analysis is one of the most common approaches for bearing fault diagnosis under variable speed by resampling [18,19]. To implement order analysis, the shaft rotational speed must be obtained in advance. Considering the design and cost, it is desirable for present order tracking methods to obtain the shaft rotation speed in the light of vibration signals instead of tachometers [20]. The time-frequency analysis of vibration signals is a popular way to estimate the shaft rotational speed, which is also referred to as instantaneous shaft rotational frequency (ISRF) for a time-varying speed case. However, due to the interference of noise, non-uniform energy distribution and other unfavorable factors, the instantaneous frequency (IF) components, including ISRF, IFCF or their harmonics, cannot always be accurately extracted. Therefore, much attention has been given to IF estimation. An algorithm using time frequency representation (TFR) with an adaptive window width has been proposed by Stankovic and Katkovnik for IF estimation [21]. Hussain and Boashash developed an adaptive IF estimation method for multi-comp onent FM signals based on tracking the component maxima in the time-frequency plane [22]. Orović et al [23] exploited the multi-window S-method for IF estimation by proposing a new distribution with high concentration in the time-frequency domain. Recently, Kwok and Jones developed an improved IF estimation technique using an adaptive short time Fourier transform (STFT) [24]. Rodopoulos et al applied a wavelet structure for IFCF estimation in rolling element bearings [25]. Among the various IF estimation methods, the peak search-based method receives an extensive application for IF extraction from TFR due to its simple yet computationally effective characteristics. It is, thus, widely used for rotating machine fault diagnosis under variable speed conditions. Gao et al employed the peak search algorithm to extract rotational speed-related IF ridges for gear fault detection [26]. Combet and Zimroz proposed an original method for IF estimation based on an instantaneous time-scaling factor estimation along the vibration signal [27]. Later, Zhao et al suggested searching the point which has a local maximum amplitude at each instant to estimate IF components after the use of a kurtogram and generalized demodulation [28]. By advancing STFT, the STFT-Viterbi algorithms were explored to obtain the IF estimation of distorted signals collected from train bearings [29,30]. These algorithms have enriched the literature of IF estimation based on peak value detection. However, the traditional peak amplitude search-based IF estimation is not sufficiently accurate in a high noise environment and complex working conditions, for instance, variable speed operations. To overcome the inaccuracy of IF estimation from the signal blurred by noise, one mainstream idea is to estimate IF ridges based on path functions or wavelet structures to reduce the impact of the signal itself or the noise. However, taxing computation, a complicated process, and low accuracy in the case of a low signal-to-noise ratio (SNR) limit the industrial applications of such methods. Therefore, an algorithm for IF estimation with computational efficiency and sufficient accuracy is demanding for bearing fault diagnosis under time-varying speed conditions. Furthermore, even if the IF is accurately estimated, the resampling would inevitably propagate errors in the result because resampling is achieved via polynomial interpolations, while vibration signals are generated by cyclic phenomena and thus are sinusoidal, not polynomial, in nature. Furthermore, order tracking might cause the carrier frequencies of the transient responses to extend to a wider scope as natural characteristics of the bearing system rarely vary, which is not conducive to bearing fault feature extraction [31,32]. Fortunately, the geometry parameters of a monitored bearing are all fixed; hence, the ratio, i.e. fault characteristic coefficient (FCC), between the IFCF and ISRF is a constant related to bearing structural parameters and can be applied for bearing fault identification [33].
According to the above statement, as long as both the IFCF and ISRF are accurately estimated, the FCC and its multiples can be obtained to perform bearing fault diagnosis under time-varying speeds. To extract the IFCF, ISRF as well as their harmonics, the lower and resonance frequency bands are chopped. Then, a novel integration strategy is proposed to integrate multiple ISRF and IFCF related ridges from the lower and resonance frequency bands to achieve an accurate IF estimation. In the previous standard deviation (STD)-based integration strategy, integration intervals are firstly determined by the difference between two IF ridges. Then the ridge with the smaller STD is directly selected as the final IF estimation result after the integration. This integration strategy can obtain an IF estimation with efficient computation; however, it has two main obstacles. Firstly, only two IF ridges are extracted for the integration, which may lead to insufficient IF-related information for an accurate IF approximation. Secondly, once two IF ridges in abnormal intervals have large errors simultaneously, the integration result would be inaccurate. To solve such a problem, an approach based on multiple IF ridge extractions is proposed in this paper to attempt to improve the IF estimation accuracy. Specifically, multiple IF ridges are firstly extracted from the vibration signal and then a new probability density function (PDF)-based integration strategy is developed to integrate IF ridges instead of simply choosing one IF segmentation via STD. Furthermore, to overcome the disadvantages of resampling [32,34], the bearing fault type is diagnosed based on the average ratios of any two rectified IF ridges, respectively, from the lower and resonance frequency bands. Ideally, the ratios should be the FCC or FCC-related numbers depending on which IF ridges are extracted featuring bearing fault types.
The rest of the paper is structured as follows. The proposed method is introduced in section 2. Section 2.1 presents the proposed integration strategy, which includes frequency band chopping, multiple IF ridge extractions and the PDF-based IF ridge integration strategy. Section 2.2 expounds the diagnosis strategy, including the rectification of multiple IF extractions and the FCC average-based fault diagnosis. The performances of the proposed approach via simulated and experimental signals are described in sections 3 and 4, respectively. Conclusions are drawn in section 5.

The proposed method
The proposed method is composed of two main modules, i.e. multiple ridge-based IF estimation and an FCC average-based fault diagnosis strategy, which are presented as follows.

IF estimation
To begin with, frequency bands of the signal are chopped to obtain a lower frequency band signal and resonance frequency band signal, paving the way for multiple IF ridge acquisition. Then a novel integration strategy based on PDF is developed to integrate multiple IF ridges for accurate IF estimation.

Frequency band chopping and IF extraction via
RPSA. To obtain multiple IF ridge candidates, the signal is segmented into the lower frequency band signal and resonance frequency band signal using band-pass filtering and the fast spectral kurtosis (FSK) algorithm, respectively [35,36]. With the separation of two frequency bands, their TFR can be obtained by STFT, which is the most popular analysis method developed to obtain the TFR with low computational complexity and a clear physical meaning [37]. Afterwards, the regional peak search algorithm (RPSA), which can effectively identify the frequency with the maximum spectrogram within a regional frequency band at each time instance, is utilized to extract IF estimations from the obtained TFRs [38]. Key steps of the RPSA are detailed as follows: Step 1: Pick the frequency with the maximum spectrogram at the starting time.
Step 2: Pick the neighboring frequency with the maximum spectrogram through the entire frequency band at time bin m, which can be formulated as f (m, n) = arg max SP(m, n), n = 1, 2, . . . , where N represents the N -point Fourier transform. m and n stand for time and frequency bins in the TFR, respectively, and SP(m, n) is the spectrogram obtained by the STFT.
Step 3: Evaluate the extracted f (m, n) via the criterion where α and β are user-specified parameters set to generate the frequency search range. The parameters α and β can be set according to the evaluation in the literature [38]. If neither equations (2) or (3) are satisfied, the corre sponding maximum spectrogram is set to zero and steps 2 and 3 are repeated until one of the equations is satisfied. The output frequency is the frequency estimator at the current time bin.
Step 4: Proceed with the time index m until all time bins with one IF ridge are extracted from the TFR.
Step 5: Set the corresponding maximum spectrogram of the extracted IF ridge to zero and repeat all steps above to extract another IF ridge.
By this means, l IF ridges can be received after l − 1 repetition extractions. In this paper, four IF ridges are taken into consideration-two extracted from the lower frequency band denoted by f 1 , f 2 and the other two in the resonance frequency band denoted by f 3 , f 4 . Then, at least two IF ridges contain information about the shaft rotational speed, and likewise, two IF ridges convey bearing fault features, which are sufficient for identifying bearing faults under time-varying speed operation. In this way, multiple IF ridges are acquired, paving the way for the integration. It is worth mentioning that the number of IF ridge extractions is not limited to four. The reason for selecting four IF ridges is that too few IF ridges may not be sufficient for carrying enough IF information, while too many IF ridges will not always improve the estimation accuracy as these IF ridge candidates themselves may have larger deviations from the true IF.

PDF-based IF ridges integration.
Although four IF ridges (i.e. pre-IF ridges) are extracted from the lower and resonance frequency bands, it cannot ensure that all IF ridges are accurate and effective enough for the PHM of bearings. Actually, it is most likely that IF ridges are inaccurate, at least at some time intervals, due to the noise interfering and nonuniform signal component energy distribution. Therefore, the PDF-based IF ridge integration strategy, which attempts to make four IF ridges supplement each other, is elaborated in this section to obtain an accurate IF estimation. Major steps of the integration strategy are presented as follows.
Step 1. Abnormal interval localization based on STD. As mentioned above, IF ridges are extracted from both the lower and resonance frequency bands. Obtained IF ridges are apparently different from each other. Specifically, they are differentiated by coefficients. In order to perform the integration, these obtained IF ridges need to be synchronized first. In this paper, three IF ridges ( f 1 , f 2 and f 4 ) are synchronized with the IF ridge ( f 3 ). The synchronization can be formulated as where f 1 stands for the synchronized shaft speed related IF f 1 , m is the synchronization coefficient and can be determined by finding the minimum absolute value between f 1 and f 3 , which can be expressed as where M is the maximum of m i and µ is the resolution of M. M is set to cover ten times of the fault characteristic order, which is wide enough to cover the main IF ridges with fault features. The resolution µ is set as 0.01 to ensure that the proper synchronization coefficient is not missed. The other synchronized IF f 2 and f 4 can be obtained in the same way. Then taking synchronized frequency bins f 1 , f 2 , f 4 and the based IF ridge f 3 at each time instance τ i into account, their STDs can be expressed as stands for the average of the four synchronized frequencies.
If the four pre-IF ridges are accurate enough, the STD of the four peak search results should be small, even zero. However, the interference of noise and non-uniform energy distribution in the TFR may cause inaccurate ridge peak search results to occur, as stated above. In this case, STD (τ i ) is likely to be larger. In other words, pre-IF values at this time instance appear to be abnormal. Abnormal IF values at consecutive time instances form abnormal intervals, which are also integration intervals and are determined by a rising edge τ l and a falling edge τ r . The rising and falling edge can be determined by a certain threshold: where ∆STD denotes the threshold. Apparently, too small a threshold produces too many abnormal intervals and leads to a high computational cost. In contrast, too large a threshold leads to an inadequate integration. Therefore, in this section, probability distribution statistics based on PDF are applied to determine the threshold, which can be expressed as Although some segments of the four pre-IF ridges extracted via the RPSA often inevitably deviate from the true IF ridge due to the adverse effects of noise and the non-uniform energy distribution, their frequency values at most time instances are correct or with errors at an acceptable level. The PDF of STD values is used to determine the threshold ∆STD. The peak of the PDF corresponds to the selected ∆STD, as illustrated in figure 1.
Once the threshold ∆STD is decided, abnormal intervals can be determined. The frequency segments, whose STDs are greater than the threshold, are taken as integration intervals. In contrast, the other frequency segments are considered accurate and do not need to be processed further. A schematic diagram of the abnormal intervals is illustrated in figure 2.
Step 2. Pretreatment of multiple IFs for abnormal intervals.
With the integration intervals determined, the average of four synchronized pre-IF candidates at each time instance is taken as the final frequency bins within normal intervals. For abnormal intervals, one suggests that two IF ridges are extracted for the integration. The STDs of these two pre-IF ridges are calculated to quantify the amount of variation. A lower STD indicates a small dispersion of a set of data. Therefore, the ridge with a STD STD lower STD is considered to be closer to the true IF and selected as the final IF estimation, while the other one with higher STD is abandoned. The above statement can be formulated as where f std indicates the pre-IF ridge with the smallest STD in abnormal intervals. However, when two estimations are inaccurate simultaneously, it is obvious that frequencies in abnormal intervals could be inaccurate as well. To address this problem, four IF ridges-instead of two-are extracted in this paper to obtain more information and improve the IF estimation accuracy. More importantly, a new integration approach based on PDF is proposed. Considering that the rotating machinery is an inertial object and its spinning speed is unlikely to experience sudden changes, the frequency accordingly does not experience an abrupt change during a short period. However, IF bins within abnormal intervals are usually subjected to relatively large fluctuations, and the STDs of the four IF ridges at each time instance could be larger. In order to reduce the frequency volatility and obtain more accurate results, four original frequencies are preprocessed at each time instance before the integration. The pretreatment of multiple IF ridges for abnormal intervals consists of two parts-namely the outlier elimination and frequency value redistribution. Firstly, a range to exclude outliers is specified. Frequency bins beyond the range are treated as outliers and frequency bins within the range are reserved. To deal with four frequency bins at time instance τ i in abnormal intervals, two ranges are considered and the smaller one is chosen as the final range, which can be formulated as follows: (12) where f i−2 and f i−1 represent frequencies at the first two time instances τ i−2 and τ i−1, respectively. f stdi represents the integration frequency value at time instance τ i , determined by the STD-based approach introduced above. With the elimination of outliers, frequency bins are redistributed in the following way: and k is the number of random parameters which is set to be four. Redistributed k frequency bins generate a random distribution, whose average is µ and STD is dyt. On one hand, k ensures that the number of frequency bins adopted for the redistribution is not reduced due to outliers. On the other hand, the average and STD ensure that redistributed frequency values remain characteristics of the original frequency bins. This pretreatment contributes to the improvement of the final integration results.
Step 3. Multiple processed IF integration for abnormal intervals using PDF. After the pretreatment of four synchronized IF ridge segments in abnormal intervals, PDF is applied to integrate the preprocessed IF ridges, which can be expressed as .., f k and f i−1 guarantee that frequency values do not change abruptly, mean f 1 , f 2 , f 3 , f 4 reserves the characteristics of the original extracted four frequency ridges, and f stdi enables the final frequency f i integrated by the PDF-based approach to outperform the original STD-based approach. The frequency value corresponding to the peak of the PDF is taken as the estimation of the true IF value. In this integration strategy, the distribution of PDF is based on k + 3 frequency parameters. To ensure the efficiency of the strategy, the number of these parameters contributing to the distribution should not be too small or large-5 to 10 is recommended [39].

Bearing fault diagnosis strategy
With the integration strategy, the final IF estimation (i.e. IF) after integration can be accurately obtained. However, the bearing fault type cannot be identified yet. The traditional bearing fault diagnosis is to resample the vibration signal using the IF estimation after integration. The disadvantage is that resampling is very sensitive to the IF estimation and the resampling accuracy is affected enormously by interpolation methods. Thus, this paper proposes a resampling-free diagnosis strategy based on the average ratios of every two rectified IF ridges from the lower and resonance frequency bands, respectively.

Pre-IF rectification.
As stated above, two IF ridges have been extracted from the lower frequency band, i.e. f 1 and f 2 , which should be the ISRF or/and its harmonic(s). The other two from the resonance frequency band, i.e. f 3 and f 4 , should be the IFCF or/and its harmonic(s). The ratios of any two IF ridges respectively from the lower and resonance frequency bands should be FCC-related, from which the fault type can be diagnosed. In order to rectify the four original pre-IF ridges which often inevitably deviate from true IFs at some time instances and obtain an accurate FCC, the IF estimation obtained by the integration is applied to rectify the remaining three synchronized pre-IF ridges, i.e. f 1 , f 2 and f 4 . The details are as follows.
Step 1. Abnormal interval localization based on PDF. As this integration is between two synchronized IF ridges, the absolute value between them is calculated, which can be obtained by where f L (τ i ) is one of the four synchronized frequency bins at time τ i . If the pre-IF ridge is ideally accurate, the difference between two IF estimations should be equal to zero. However, due to the interference of noise, non-uniform energy distribution or other adverse factors, ∆f (τ i ) is likely to present nonzero. Generally, the greater the difference, the less accurate the frequency estimation is. If the difference ∆f (τ i ) is greater than the threshold denoted by ∆f , the IF appears abnormal at the time instance τ i . Abnormal frequency bins form abnormal intervals, which are also integration intervals and are determined by a rising edge τ u and a falling edge τ d . The rising and falling edge can be decided by the threshold ∆f : The threshold ∆f is determined via PDF. The one corresponding to the peak of the PDF is selected as the threshold, which can be expressed as Once the threshold ∆f is decided, abnormal intervals can accordingly be determined. The ∆f greater than the threshold forms integrated intervals. In contrast, the ∆f which is not greater than the threshold constitutes non-integrated intervals, as illustrated in figure 3 Step 2. Dual IF integration for abnormal intervals using STD. The STD-based integration strategy is expounded in step 2 of section 2.1.2. The final integration result can be expressed as where fr (τ ) is the frequency rectified by the IF estimation after integration.

FCC average-based fault diagnosis.
For the bearing fault diagnosis, the fault type can be identified via ratios of any two rectified IF ridges from the lower and resonance frequency bands, respectively. Since four pre-IF ridges are extracted and rectified, i.e. f r 1 , f r 2 , from the lower frequency band and f r 3 , f r 4 from the resonance frequency band, four ratios can be obtained, namely These ratios θ 1 , θ 2 , θ 3 and θ 4 should be related to FCC. The FCC average can then be expressed as where η 1 , η 2 , η 3 and η 4 are synchronized coefficients which can be determined according to the ratios of each two rectified IF ridges from the lower frequency band and resonance frequency band, respectively. The fault diagnosis can then be performed based on the calculated average of the FCC. It is worth mentioning that the FCC average-based fault diagnosis approach is expected to improve the diagnosis accuracy under circumstances where IF estimation deviations remain after rectification. Even if individual ratios deviate, the average-based nature can eliminate the deviation to some extent. The flowchart of the proposed approach can be graphically represented in figure 4.

Simulated case study
In order to verify the effectiveness of the proposed integration and diagnostic approach, the simulated signal is constructed in this section. The model of the simulated signal of a faulty bearing vibration is expressed as where x 1 (t) represents fault feature-related components, which can be defined as where [1 + α cos(2πx 2 (t))] reflects the phenomena of the resonance being demodulated by the shaft rotating frequency x 2 (t), α is the amplitude of the modulation (α = 0.11), M is the signal length, A m is the amplitude of the mth impulse response, β is the coefficient related to damping (β = 800), ω r is the excited resonance frequency, u (t) is a unit step function, t m is the occurrence time of the mth impulse which is calculated as where µ m is the slippage ratio varying from 0.01-0.02 [40,6], x 2 (t) represents shaft rotational speed-related components, η (t) is additive white Gaussian noise (SNR = −12 dB). x 2 (t) can be defined as 3 are amplitudes and are set as 1.6, 1.2 and 1 respectively, the ISRF of the simulated signal is defined as φ = 2.49t 2 + 5t + 35 and the FCC is set to be 3.7. The sampling frequency is 20 kHz and the sampling time is 5.1 s.
The simulated signal is shown in figure 7(a) and the corresponding TFR is shown in figure 7(b). Three visible IF ridges in figure 7(b) show the IFCF and its harmonics; however, frequency components related to the rotational speed cannot be obtained. Therefore, the bearing fault type cannot be identified. Then, the proposed method is used to analyze the simulated signal.
Firstly, the lower frequency range is chosen as [0, f 0 ], where f 0 is 400 Hz, which is enough to cover the common ISRF and its harmonics. The resonance frequency band is selected via the FSK algorithm. The TFR of the lower frequency band is shown in figure 7(c). Subsequently, the RPSA is applied to extract two IF ridges from the lower frequency band. The two extracted ridges are presented in figure 7(d). Similarly, figures 7(e) and (f) show the TFR of the resonance frequency band and two extracted IF ridges from it via RPSA, respectively. It can be observed that the segment close to the end of the extracted IF ridges experiences abrupt and frequent changes, indicating that IF extractions are not accurate enough and can affect the accuracy of the IF estimation as well as the bearing fault diagnosis.
Then the proposed integration strategy is applied. The four original frequencies at each time bin should be synchronized prior to the integration. With the synchronization, the STDs of four synchronized frequency bins are calculated for the threshold determination. The threshold is 3.362 corre sponding to the peak of the PFD of STDs, as illustrated in figure 5.
With the determined threshold ∆STD, integration intervals are settled. Frequency segments with STDs larger than ∆STD form abnormal intervals, while the ones with STDs not greater than ∆STD are classified within normal intervals, as shown in figure 6. The plot between 0 and 2.5 s in figure 6 is zoomed in to clearly show the abnormal intervals.
Within normal intervals, the averages of four synchronized frequencies at each time bin are taken as the final frequencies. In terms of abnormal intervals, the traditional STD-based integration strategy is applied. The resulting IF estimation is shown in figure 7(g), showing that the IF segment close to the end still experiences abrupt changes. For comparison, the proposed integration results obtained by the PDF-based integration approach are shown in figure 7(h), where the IF   estimation after the integration is smooth and does not experience sudden variations. It can then be concluded that the IF estimation accuracy is improved by the PDF-based integration approach.
To assess the improvement, the mean relative error (MRE) is applied. The MRE is computed by 1 where f i and f reali are the integration frequency at time τ i and the theoretical frequency value at τ i , respectively, m stands for the number of the data. The MREs of the IF estimation integrated by the STD-based and PDF-based strategy are 4% and 2%, respectively, quantitatively showing that the IF estimation after the integration obtained by the PDF-based integration strategy is more consistent with real values.
With the integration result, four pre-IF ridges are rectified by the IF estimation after integration. The rectification results are shown in figure 8. The ratios of any two rectified IF ridges from the lower and resonance frequency bands are computed, as listed in table 1.
From table 1, the ratios of each of the two IF ridges can be revealed: 2f r 1 = f r 2 and 2f r 3 = f r 4 . Thus, η 1 = 1, η 2 = 2, η 3 = 1/2 and η 4 = 1 can be obtained for the calculation of the FCC average θ: According to equation (24), the FCC average θ is 3.848, approximately identical to the pre-set FCC. In reality, the fault type can then be diagnosed by matching the calculated FCC average to the FCCs of the inner race and outer race. The relative error between the calculated FCC and theoretical FCC is 4%, within the tolerable error range 0%-5% [41]. Additionally, with the ratios listed in table 1 and the fact that the IFCF and ISRF are more likely to be of a higher amplitude than that of their harmonics, it can be identified that the two lower-located ridge paths extracted from the lower and resonance frequency bands are the ISRF and IFCF, respectively. The other two ridges are the second harmonic of the ISRF and IFCF, respectively. The numerical analysis above demonstrates that the IF can be accurately estimated using the proposed integration strategy and the bearing fault diagnosis can also be implemented via the pre-IF rectification and average FCC calculation.

Bearing with an outer race fault
In this section, the vibration signal of bearings with an outer race fault is analyzed to validate the proposed approach. The experiment is performed on the Quest machine fault simulator (MFK-PK5M), as shown in figure 9.
Two bearings (type ER16K) placed inside the bearing holders are applied to support the shaft (1 inch in diameter) which is driven by a 3-hp AC motor through a coupling with a Hitachi controller (SJ200-022NFU). The left rolling bearing has a fault in the outer race and the load mass of 5.03 kg is mounted on this 1-inch steel shaft. The right end of the drive shaft is connected to a smaller sheave and the gearbox shaft is connected to a larger sheave. The belts resting on the two sheaves connect the drive shaft to drive the gearbox shaft. The diameter ratio of the smaller sheave to the larger is roughly 1: 2.6. The tachometer and accelerometer are mounted on the test rig to measure the shaft rotational speed and the vibration signal, respectively. To ensure an optimum vibration amplitude of the gearbox, the accelerometer should be mounted near the gearbox. Then, the measured signal is collected by a DAQ card and sampled using LABVIEW. The specifics of the test are listed in table 2.
The raw vibration signal and its corresponding TFR are presented in figures 12(a) and (b), respectively. From the TFR, frequency components related to the rotational speed cannot be easily observed and, accordingly, the frequency components of interest cannot be extracted accurately from the TFR of the original vibration signal. Therefore, the FCC cannot be obtained for the fault diagnosis. The proposed method is, therefore, applied.
Firstly, the lower frequency range is chosen as [0400] Hz. The normal ISRF and its harmonics can be covered within the range. The FSK algorithm is utilized to separate the resonance frequency band. Figure 12(c) shows the TFR of the lower frequency band. The RPSA is then applied to extract two IF ridges from the lower frequency band. Figure 12(d) reveals the two extracted IF ridges. The TFR of the resonance frequency band is shown in figure 12(e), and figure 12(f) reveals the two IF ridges extracted from the resonance frequency band via RPSA. Obviously, inaccurate IF extractions result in the mutation of the segment situated at the beginning of the estimated IF ridges, which may ultimately influence the bearing fault diagnosis results.
Then the proposed integration strategy is applied. Before the integration, abnormal intervals should be determined according to the threshold ∆STD. In this case, the threshold is 8.264 corresponding to the peak of the PDF of the STDs in figure 10.
Then abnormal and normal intervals can be distinguished according to the threshold ∆STD, as illustrated in figure 11. Similarly, abnormal intervals are defined as the ones with STDs all above the threshold. In addition to IF segments within the defined abnormal intervals, the other IF segments are considered to be accurate. More defined abnormal intervals can be identified in the zoomed-in plot in figure 11.    Then, the traditional STD-based integration strategy is used for integration within abnormal intervals. Figure 12(g) shows the resulting IF estimation while figure 12(h) presents the integration results obtained by the proposed PDF-based integration approach. The IF ridge integrated by the proposed method is smooth and does not reveal any abrupt changes, while the IF segment in figure 12(g) still experiences fluctuations near the beginning. Obviously, the accuracy of the IF estimation after integration is improved by the PDF-based integration approach.
Moreover, the MRE is also applied to evaluate the accuracy of the estimator. The MREs of the IF estimation obtained by the STD and PDF-based strategy are 9.1% and 3.4%, respectively. It can be concluded, therefore, that the IF estimation obtained by the PDF-based integration strategy is more accurate.
With the more accurate estimated IF, four pre-IF ridges are rectified. Figure 13 shows the rectification results. The ratios of any two rectified IF ridges in the lower and resonance frequency bands are shown in table 3.
Then the synchronization coefficients η 1 , η 2 , η 3 and η 4 can be obtained as 1, 2, 1/2 and 1, respectively, according to table 3. The FCC average θ can be calculated as θ = ( f r 3 /f r 1 + 2 × f r 3 /f r 2 +( f r 4 /f r 1 )/2 + f r 4 /f r 2 )/4 = 3.452. (25) With the calculated FCC average θ, the fault type can then be identified by matching θ with the theoretical FCCs, from which it can be confirmed that the fault is located on the outer race of the bearing. The relative error between the calculated FCC average and the theoretical one is 1.37%-less than the allowable maximum error 5%, indicating that the estimation result is reasonable. Similar to the simulated case, by examining table 3, it can be concluded that the lower-located two ridge paths extracted from the lower and resonance frequency bands are the ISRF and IFCF, respectively. The remaining two ridges are the second harmonic of the ISRF and IFCF, respectively.
In addition, it can be observed that both ratios θ 1 and θ 2 are larger than the theoretical FCC and both ratios θ 3 and θ 4 are less than the FCC. Fortunately, the proposed FCC averagebased strategy can effectively relieve the deviations between the calculated FCC and theoretical one, resulting in a more accurate average FCC for bearing fault diagnosis.

Bearing with an inner race fault
In this section, the vibration signal of a bearing with an inner race fault is investigated to validate the proposed approach. The experiment is carried out on a SpectraQuest Machinery fault simulator (MKF-PK5M) at the University of Ottawa lab. The experiment setup is shown in figure 14.
The shaft is driven by a motor controlled by an AC inverter. There are two bearings mounted to support the shaft and the load. The unhealthy bearing is positioned at the left side. The tachometer and accelerometer are mounted on the test rig to collect the shaft speed and vibration signal, respectively. Processed by a signal conditioner (PCB PIEZOTRONIC 482C15), the data is fed to an NI data acquisition module (NI USB-6212 BNC) and is recorded by a computer with Labview software. The specifics of the test are listed in table 4. Figure 17(a) presents the raw vibration signal and figure 17(b) shows the TFR of the signal. From the TFR of the original signal, no clear IF ridges are visible or can be extracted because of the noise interference and non-uniform energy distribution, which are adverse for the IF estimation and fault type diagnosis. Therefore, the proposed method is used to analyze the experimental signal.
Firstly, the lower and resonance frequency bands are separated via the bandpass filter and FK algorithm, respectively. The bandwidth remains as [0, f 0 ], where f 0 = 400 Hz, to cover the normal ISRF and its harmonics. The TFR of the lower frequency band is shown in figure 17(c). Subsequently, the RPSA is applied to extract two IF ridges from the lower frequency band. Two extracted IF ridges are presented in figure 17(d). Similarly, the TFR of the resonance frequency band and two IF ridge extractions in this frequency band are presented in figures 17(e) and (f), respectively. It is clear that the beginning  segment of the estimated IF ridges is extracted imprecisely and changes abruptly. The IF estimation and bearing fault diagnosis can be affected due to the frequency mutation.
With the synchronization of four IF extractions, the threshold ∆STD is then determined after the calculation of the STDs of four synchronized frequencies at each time bin. The threshold takes the STD corresponding to the peak PDF to distinguish abnormal intervals and normal ones, as shown in figure 15. Figure 16 reveals distinguished abnormal   and normal intervals. As marked in the figure, IF segments with STDs larger than the threshold ∆STD are recognized as abnormal. The zoomed-in plot in figure 16 exhibits details between 11.5 and 13.5, where two abnormal intervals can be found for the integration. Within determined abnormal intervals, the proposed integration strategy and the traditional STD-based integration approach are then applied to integrate pre-IF ridges to achieve a smoother IF estimation after integration. Integration results obtained by the two integration approaches are shown in figures 17(g) and (h), respectively. It can be seen that the IF estimation after integration is smoother and improved by the PDF-based integration approach.
To quantify the improvement, the MRE is applied to evaluate the dispersion between the estimator and real IF. The MRE of the IF estimation integrated by the STD-based strategy is 6.6%, while it is 4.6% when integrated by the PDF-based strategy. To conclude, the PDF-based integration strategy is more efficient than the traditional STD-based integration strategy and a more accurate IF estimation after integration can be obtained by the proposed integration method.
Then four pre-IF ridges are rectified by the IF estimation after integration. The rectification results are shown in figure 18. The ratios of any two rectified IF ridges in the lower and resonance frequency bands are computed, as exhibited in table 5. In the light of the ratios in table 5, the relationships between each two rectified IF ridges respectively from the lower frequency band and resonance frequency band are 3f r 1 = 2f r 2 , 2f r 3 = f r 4 . The coefficients η 1 , η 2 , η 3 , and η 4 can then be determined as 1, 3/2, 1/2, and3/4 . The FCC average θ is computed as θ = f r 3 /f r 1 + 3 2 (f r 3 /f r 2 ) + 1 2 ( f r 4 /f r 1 ) + 3 4 ( f r 4 /f r 2 ) /4 = 2.746 (26) which is identical to half of the FCC of the inner race, rather than FCC I . The reason is that the extracted shaft-related IF ridges are the second and third harmonic of ISRF, whereas the identified IFCF-related IF ridges are the fundamental IFCF and its second harmonic, as shown in figure 18. Given the calculated FCC average, it can be concluded that an inner race fault occurs on the bearing. Moreover, the relative error is calculated to be 1.14%, within the allowable error range 0%-5%, which means that the fault type can be accurately determined. It is likely that the error between the calculated average FCC and theoretical one can be reduced via the average strategy. Thus, the multiple ridge integration and FCC average-based bearing fault diagnosis strategies are demonstrated to be efficient, as per the experimental analysis above.

Conclusion
To perform the PHM of bearings, an effective method without a tachometer and resampling process is proposed in this paper for IF estimation and bearing fault diagnosis under time-varying speed conditions. The two main unique features of the proposed approach lie in the PDF-based ridge integration strategy for an accurate IF estimation and FCC   average-based fault diagnosis. The former contributes to the precise IF approximation without relying on a measurement device and the latter is devoted to bearing fault diagnosis without the involvement of order tracking. Major advantages of the proposed fault diagnosis method can be summarized as follows: (a) multiple IF ridge extractions instead of two can reveal more information regarding the IF and thereby improve the estimation accuracy; (b) the proposed PDF-based integration strategy is effective, even if all extracted pre-IF ridges are inaccurate in abnormal intervals; (c) the pre-IF estimation rectification-and FCC average-based diagnosis means the proposed method is no longer dependent on a single IF, which can tolerate estimation errors and would not affect the diagnosis results. With the proposed method, the health indicator of bearings can be obtained, which paves the way for the remaining useful lifetime prediction of bearings and health management in energy systems. It should be emphasised that the proposed integration strategy may become ineffective when all considered frequency bins at a time instance are inaccurate. However, it is unlikely that all frequency bin interests are invalid simultaneously. Moreover, more pre-IF ridge paths are taken into consideration, which will reduce the probability that all frequency bins at a time instance are inaccurate. Thus, the effectiveness of the proposed method would generally not be influenced.