A novel threshold segmentation instantaneous frequency calculation approach for fault diagnosis

1.   Introduction

Rolling bearings play a crucial role in the functioning of rotating machinery ^[1,2]. The maneuver of rotating machinery is entirely dependent upon the health state of the rolling bearings, which accounts for almost 45–55% of these equipment failures ^[3,4]. Therefore, fault monitoring of rolling bearings can increase the reliability of rotating machinery to a certain extent. Fault feature information extraction is one of the most important problems in fault monitoring of rolling bearing ^[5,6,7]. The traditional diagnosis technique is to extract the feature vector for fault recognition according to the time domain or frequency domain feature of vibration signal ^[8,9,10]. In fact, when the rolling bearing failure, the vibration signal of rolling bearings often shows non-stationary characteristics due to the nonlinear stiffness, friction, clearance, and external load ^[11,12,13]. Therefore, how to go from extract fault feature information of non-stationary signals in the fault diagnosis of rolling bearing is particularly important ^[14,15].

One of the common characteristics of the non-stationary signal is the varying frequency ^[16]. In general, frequency refers to the overall characteristics of the periodic signal in a certain period of time, which cannot describe the local characteristics of the signal ^[17,18]. Therefore, the concept of instantaneous frequency (IF) is introduced ^[19,20,21]. IF can well reflect the characteristics of the signal transient component, whereas the transient component is an important manifestation of the non-stationary signal local characteristics, which contains a large amount of effective information ^[22,23,24]. So IF can be used to judge the timing and type of mechanical equipment failure ^[25,26].

There are a vast number of commonly used IF calculation methods, among which the Hilbert-Huang transformation (HHT) is the most widely used adaptive time-frequency analysis in recent years ^[27,28,29]. This method was proposed by Norden E. Huang of NASA. It is a general method to solve the IF of the single component real signal ^[30]. The analytical signal is obtained by the Hilbert transform, and the IF is obtained by differentiating the phase of the analytical signal ^[31,32]. This is a widely accepted and used definition of IF in academic circles, but it is also controversial. First, the IF may not be one of the frequencies in the spectrum ^[33]. Second, if only a few distinct frequencies constitute a linear spectrum, then the IF can be continuous and can vary over an infinite range of values ^[34,35]. Third, although the spectrum of the analytic signal is zero for negative frequencies, the IF can be negative. However, this is meaningless in practice ^[36]. Fourth, for a band limited signal, its IF can be outside the band. Fifth, the Hilbert transform is a global operation, not an idea of local time-frequency analysis. In order to calculate the analytical signal at a certain moment, it is necessary to know the signal of all time ^[37]. Sixth, according to the physical nature of the signal, the signal can be divided into single component signals and multicomponent signals. The single component signal has only one frequency at any time, while the multicomponent signal can have multiple frequencies. The Hilbert transform only obtains one frequency value for any signal, so this definition only applies to single component signals, whereas multicomponent signals have no physical meaning for discussing a single frequency. In fact, not all single component signals can be the Hilbert transform, and the Bedrosian theorem and the Nuttal's theorem must be satisfied ^[38,39].

According to the problem of negative frequency, Norden E. Huang et al. proposed a method of estimating IF based on empirical automated mapping and facilities management (AMFM) decomposition, namely Normalized Hilbert transform ^[40]. The intrinsic mode function (IMF) is decomposed by empirical AMFM to obtain a pure FM signal $F\left(t \right) = \cos \varphi \left(t \right)$. At this time the instantaneous amplitude of $F\left(t \right)$ is 1, which is no longer restricted by the Bodrosian theorem. So Hilbert transform can be used to obtain the IF, which overcomes the defect that the Hilbert transform will appear negative frequency. There is a great improvement over direct the Hilbert transform. However, since the Hilbert transform is still adopted, energy leakage will be generated at the endpoint, and the endpoint effect is still inevitable.

Zheng et al. ^[41] proposed an empirical envelope method to obtain IF, for unexplained negative frequency and obvious the endpoint effects of the Hilbert transform. Firstly, the pure FM signal is obtained by empirical amplitude modulation and frequency decomposition, and then it is derivation. After that the empirically amplitude-modulated frequency decomposition is performed on the derivative result, and the envelope signal is extracted to obtain the IF of the original signal. The calculation is simple and convenient, only requires two empirical AMFM decompositions and one derivative. However, the use of cubic spline to fit the envelope in the calculation process and a certain number of iterations will also produce the endpoint effect.

Cicone A et al. compared the single component signal with the Second Ordinary Differential Equations (ODEs) in polar coordinates, and proposed a local definition of instantaneous phase and IF ^[42,43]. The $x$ of the single component function $f\left(x \right)$ is regarded as the $x$ coordinate in the second order differential equation, the coordinate $\left({x, y} \right)$ is mapped to $\left({r, \theta } \right)$ by $x = r\cos \theta, y = - r\sin \theta$, and the IF is the derivative of the phase angle $\theta '$ in the polar coordinate^[44,45,46]. The method is simple in calculation, strengthens the ability to quickly track the instantaneous frequency, eliminates the negative frequency, and makes the instantaneous frequency more physical. Moreover, it has no extreme value operation, and avoids the error caused by the endpoint effect ^[47,48,49]. But when the amplitude of the signal changes suddenly, the error of the calculated results of signal IF increases.

In summary, this paper proposes a method to deal with the IF of the non-stationary signal with amplitude abrupt change characteristics. This method not only has simple calculation principle, but also can reduce the unnecessary influence on the calculation of IF caused by the sudden change of amplitude in the non-stationary signal, and can effectively judge the fault of rolling bearing.

2.   Threshold segmentation method for instantaneous frequency

The signal whose statistical characteristics change with time is the non-stationary signal ^[50,51], and the amplitude mutation is also one of the manifestations of the non-stationary signal. When the amplitude of the signal changes greatly, the measured IF will generate some unnecessary fluctuations in the place where the amplitude changes greatly, and the IF of the original signal cannot be well reflected ^[52,53,54]. Based on the ordinary differential equation (ODE) method proposed by Cicone, this paper proposes a method of threshold segmentation which can deal with the IF of non-stationary signals with amplitude mutation characteristics ^[55,56]. When the original signal amplitude is abrupt, it needs to be decomposed into IMFs by empirical mode. The IMFs is the segmented and then the IF of each segment is calculated. The specific steps are as follows ^{[57,58,59,60]}:

(1) Let $f\left(x \right)$ is an IMF signal. Firstly, determine whether there is a sudden amplitude change of $f\left(x \right)$. Let a certain point ${x_i}$ have a corresponding amplitude value ${a_{\rm{i}}}$, then the amplitude value at point ${x_{i - 1}}$ is ${a_{i{\rm{ - }}1}}$, and the amplitude value at point ${x_{i + 1}}$ is ${a_{{\rm{i + 1}}}}$. Calculate the integral of the signal $f\left(x \right)$ in the interval $\left[{{x_{i- 1}}, {x_i}} \right]$ by using (1), that is, the area ${S_j}$ of the signal $f\left(x \right)$ in the corresponding interval with the $x$ axis, where $j = i - 1$:

The principle of integral formula used in this paper is the combined trapezoidal formula in numerical integration. The specific principle is as follows. Suppose that interval $\left[{a, b} \right]$ is divided into equidistant nodes ${x_k} = a + kh$, $k = 0, 1, \cdots, M$ and divided into $M$ subintervals $\left[{{x_k}, {x_{k + 1}}} \right]$ with width $h = \left({b - a} \right)/M$. The combined trapezoidal formula of $M$ subintervals can be expressed as any one of three equivalent ways:

Or:

Or:

This is an approximation of the integral of $f\left(x \right)$ in the interval $\left[{a, b} \right]$, written as:

(2) Differentiate ${S_j}$ in the interval $\left[{{x_{i- 1}}, {x_i}} \right]$, which is the change of area ${P_j}$. Moreover, the method of integration and differentiation is also equivalent to magnifying the change of signal, which is easy to observe. As shown in Figure 1, since the signal sampling interval is the same, the time interval of each group interval is the same as $\Delta t$, so the larger the area change in time $\Delta t$, the larger the amplitude change of signal $f\left(x \right)$ in this interval. Therefore, the change trend of the signal amplitude can be judged by the change of area, that is, the ratio $\rho$ of ${P_j}$ in two adjacent $\Delta t$ times can be determined. Moreover, the ratio is within the range of $\left[{0, 1} \right]$. As the value of $\rho$ approaches 1, the amplitude changes more gently, and the value of $\rho$ approaches 0, and the amplitude changes more rapidly.

The original signal is segmented on the basis of $\rho$. The smaller the value of $\rho$ is, the more segments the signal is segmented, and the better the instantaneous frequency calculation result will be, but the calculation amount will also increase accordingly. Therefore, this paper tested the relationship between the value of $\rho$ and the correlation coefficient of instantaneous frequency calculation results, as shown in Figure 2. Considering the real-time performance and accuracy of the calculation, $\rho$ in this paper is 0.3.

(3) Find all amplitude abrupt points. Assuming there are $n$ amplitude abrupt points, the original signal is divided into $n + 1$ segments based on these $n$ points, which are ${f_1}\left(x \right), {f_2}\left(x \right), \cdots, {f_{n + 1}}\left(x \right)$ respectively.

(4) Find envelope ${q_1}\left(x \right), {q_2}\left(x \right), \cdots {q_{n + 1}}\left(x \right)$ corresponding to ${f_1}\left(x \right), {f_2}\left(x \right), \cdots, {f_{n + 1}}\left(x \right)$, then:

Where $i = n + 1$.

(5) Calculate the derivative ${f'_1}\left(x \right), {f'_2}\left(x \right),$$\cdots, {f'_{n + 1}}\left(x \right)$of ${f_1}\left(x \right), {f_2}\left(x \right), \cdots, {f_{n + 1}}\left(x \right)$, and get the envelope signal ${r_1}\left(x \right), {r_2}\left(x \right), \cdots, {r_{n + 1}}\left(x \right)$ of ${f'_1}\left(x \right), {f'_2}\left(x \right), \cdots, {f'_{n + 1}}\left(x \right)$, then:

The envelope functions ${q_i}\left(x \right)$ and ${r_i}\left(x \right)$ are not unique. For example, the extreme values of ${f_i}\left(x \right)$ and ${f'_i}\left(x \right)$ can be obtained by cubic spline interpolation.

(6) Define:

Then ${F_i}\left(x \right)$ in the complex plane is a curve in $\left[{- 1, 1} \right] \times \left[{- 1, 1} \right]$.

(7) ${F_i}\left(x \right)$ is an orbiting unit circle, and its rotation angle is:

And ${\theta _1}\left(x \right), {\theta _2}\left(x \right), \cdots, {\theta _{n + 1}}\left(x \right)$ is the instantaneous phase of ${f_1}\left(x \right), {f_2}\left(x \right), \cdots, {f_{n + 1}}\left(x \right)$.

(8) Use equation (10) to find the IF ${\omega _1}\left(x \right), {\omega _2}\left(x \right), \cdots, {\omega _{n + 1}}\left(x \right)$ of each segment of signal ${f_1}\left(x \right), {f_2}\left(x \right), \cdots, {f_{n + 1}}\left(x \right)$:

Even though the envelope functions $q\left(x \right)$ and $r\left(x \right)$ are not unique, the instantaneous phase and instantaneous frequency obtained by the above definitions hardly depend on the method of solving the envelope functions $q\left(x \right)$ and $r\left(x \right)$.

(9) Connect the obtained IF ${\omega _1}\left(x \right), {\omega _2}\left(x \right), \cdots, {\omega _{n + 1}}\left(x \right)$ of each segment in time sequence, and the complete IF ${\omega _i}\left(x \right)$ of the original signal $f\left(x \right)$ can be obtained.

3.   Simulation calculation

In this chapter, three simulation signals are tested. The HHT method and the threshold segmentation method are respectively used to solve the IF.

Example 1. The test signal is given by (11) as a typical linear frequency modulation (LFM) signal and shown in figure 3. Figure 4 shows the IF obtained by the HHT method, and figure 5 shows the IF obtained by the threshold segmentation method.

From the mathematical derivation, the IF of (11) can be obtained by using the HHT and the threshold segmentation method as $\omega \left(x \right) = 2x + 5/2\pi$. Thus, the obtained IF should be a straight line in the time-frequency domain. It can be seen from figure 4 and figure 5 that the IF obtained by the threshold segmentation method is closer to the theoretical value than that obtained by the HHT method.

Example 2. A simulation signal with a low frequency and a low sampling frequency is given by (12), and the sampling frequency is 100 Hz. At the same time, add a Gaussian white noise to the simulation signal, as shown in figure 6.

The simulation signal given by (12) is a signal with relatively large amplitude change, and its ideal IF should be a straight line. Figure 7 shows the IF obtained by the HHT method. Although the overall trend is approximately a straight line, there are certain unnecessary fluctuations, especially in the two periods of 2s–4s and 6s–8s where the amplitude changes greatly. Figure 8 shows the process of calculating the IF of a signal using the threshold segmentation method. First, the signal is judged and divided into three sections according to calculation, as shown in figure 8(a), figure 8(c), and figure 8(e). According to the segmented signal, the IF of each segment is calculated separately, and is shown in figure 8(b), figure 8(d), and figure 8(f). The three IFs are connected in chronological order to obtain the complete IF of the original signal over the entire time period, as shown in figure 9. It can be seen that the obtained IF is relatively close to the ideal IF curve. By comparing these figures, it can be seen that the IF obtained by the threshold segmentation method is much smoother, which can reduce the influence of sudden amplitude change on IF calculation, and is closer to the ideal IF.

Example 3, a LFM simulation signal with a high sampling frequency is selected and given by (13) with a sampling frequency of 10kHz, as shown in figure 10.

Figure 10 shows a linear LFM signal. Similarly, two methods are used to obtain its IF, as shown in figure 11 and figure 13. Figure 12 shows the process of calculating the IF of a signal using the threshold segmentation method. First, the signal is judged and divided into three sections according to calculation, as shown in figure 12(a), figure 12(c), and figure 12(e). According to the segmented signal, the IF of each segment is calculated separately, and is shown in figure 12(b), figure 12(d), and figure 12(f). The three IFs are connected in chronological order to obtain the complete IF of the original signal over the entire time period, as shown in figure 13. It can be seen that the IF obtained by the segmentation method is relatively smooth, which can reduce the impact of the sudden change in amplitude on the IF calculation, remove unnecessary fluctuations, and get closer to the ideal IF value.

4.   Application in fault diagnosis of rolling bearings

In this experiment, the QPZZ-Ⅱ system experimental bench was used to measure the vibration signals of rolling bearings, as shown in figure 14. The QPZZ-Ⅱ series of rotating machinery vibration monitoring and fault diagnosis experiment platform is a widely used in bearing fault and gear box testing experiment platform. The platform is composed of variable speed drive motor, bearing, gear box, shaft and special rotating disk, etc. In recent years, it has been favored by the related mechanical vibration and signal analysis university laboratory, especially the related vibration research company in Japan also USES this platform as the new employee training platform. This system can quickly simulate multiple states and vibrations of rotating machinery, and can perform comparative analysis and diagnosis of various states.

The NI-USB6210 data acquisition card vibration model is used in the experiment, and the bearing model is N205. Bearing related parameters are shown in Table 1.

Equations (14) is the ball pass frequency outer race (BPFO) and equations (15) is the ball pass frequency inner race (BPFI), where ${f_r}$ represents the shaft rotation frequency, $n$ is the number of rolling elements, $\varphi$ is the contact angle, and $d$ is the diameter of the rolling element. $D$ is the diameter of the roller distribution circle. The above formula describes that the characteristic frequency of rolling bearings under different fault conditions varies with the speed.

Figure 15 is a schematic diagram of the rolling bearing fault diagnosis process, which is divided into two parts: signal acquisition, signal analysis and feature extraction.

Firstly, the continuous vibration signal of the rolling bearing is measured by the acceleration sensor, and then the signal is converted into digital vibration signal through A/D. The digital signal is pre-processed for noise reduction, and the processed signal is decomposed into several IMFs. The method proposed in this paper is used to calculate the instantaneous frequency and instantaneous amplitude (IA) of each IMF and generate the local Hilbert marginal spectrum. The peak value in the marginal spectrum is compared with the fault characteristic frequency of the rolling bearing to judge whether there is a problem in the bearing and the location of the problem.

4.1.   Outer ring fault

The vibration signal of the rolling bearing with outer ring fault is selected with a speed of 1200r/min and a sampling frequency of 25 kHz, as shown in figure 16. According to (14), the fault characteristic frequency of the outer ring is 131Hz. Using the method mentioned in the second chapter, the vibration signal is divided into several segments and the IF of each segment is calculated respectively. And then connect them in time order into a complete IF in the time domain. The results were compared with the traditional HHT method, as shown in figure 17(a) and figure 18(a). Moreover, its Hilbert marginal spectrum is obtained for comparison, as shown in figure 17(b) and figure 18(b).

It can be seen from figure 17(a) that the IF obtained by the HHT method has some negative frequencies, while the IF obtained by the threshold segmentation method in figure 17(a) eliminates this problem well. In addition, by comparing the two marginal spectra of figure 17(b) and figure 18(b), the peak value in figure 18(b) is more prominent and the outer ring fault of the rolling bearing can be judged more effectively. The result of fault detection can be achieved while simplifying the calculation principle.

4.2.   Inner ring fault

The vibration signal of the rolling bearing with inner ring fault is selected with a speed of 1200r/min and a sampling frequency of 25kHz, as shown in figure 19. According to equation (15), the fault characteristic frequency of the inner ring is 192.50Hz. In the same way, the IF of the signal is obtained by the HHT method and the threshold segmentation method, and then its Hilbert marginal spectrum is obtained for comparison, as shown in figure 20 and figure 21.

Similarly, comparing figure 20(a) and figure 21(a), the threshold segmentation method eliminates the phenomenon of negative frequency. The peak value in figure 21(b) is also more obvious than the peak value in figure 20(b), and the interference is reduced, which can effectively judge the inner ring failure of the rolling bearing. Moreover, the calculation principle is relatively simple, which can save time in fault detection.

5.   Conclusion

For fault diagnosis of rolling bearings, the IF calculation has many potential benefits, it can reflect the transient components of the non-stationary signal. In this paper, a single component signal is compared with ODEs in polar coordinates, and a relatively simple method for calculating the IF is developed. The formula is explained mathematically. Simulations, mathematical derivations and experimental tests are used to highlight the performance of the proposed method. It can eliminate meaningless negative frequency, and reduce the unnecessary impact of IF due to sudden changes in amplitude. Its calculation principle is simple, and it can be applied to the occasions with high real-time requirements in the actual production.

Acknowledgments

This work was supported in part by the Natural Science Foundation of Liaoning Province under Grant 2019ZD0112 and 2019ZD0099, and in part by the National Natural Science Foundation of China under grant 51475065, grant 51605068, grant 51879027, Traction Power State Key Laboratory of Southwest Jiaotong University under grant TPL2002, and Liaoning BaiQianWan Talents Program.

Conflict of interest

The authors declare no conflict of interest.