Signal processing - bearing fault diagnosis practice based on Hilbert Demodulation (envelope spectrum), which is explained in super detail through python code


Welcome to my official account "fault diagnosis and python learning"

1. Data introduction


Including inner ring, outer ring, rolling element and normal data, which are one respectively.
1730_12k_0.007-InnerRace: inner ring fault
1730_12k_0.007-OuterRace3: outer ring fault
1730_12k_0.014-Ball: rolling element failure
1730_48k_Normal: normal data

For students who do not understand CWRU bearing data set, see here:
CWRU dataset introduction phase 1
CWRU dataset introduction phase 2
CWRU dataset introduction phase 3
CWRU dataset introduction phase 3

2. Load CWRU inner ring fault data

Next, we will analyze the fault data of bearing inner race. The original data is a mat file, which is a matlab file. Define a function to read the data

def data_acquision(FilePath):
    """
    fun: from cwru mat File reading acceleration data
    param file_path: mat File absolute path
    return accl_data: Acceleration data, array type
    """
    data = scio.loadmat(file_path)  # Load mat data
    data_key_list = list(data.keys())  # The mat file is of dictionary type. Get all the keys of the dictionary and convert them to list type
    accl_key = data_key_list[3]  # Get'X108_DE_time'
    accl_data = data[accl_key].flatten()  # Get'x108_ DE_ The value corresponding to time 'is the vibration acceleration signal, and the two-dimensional array is expanded into a one-dimensional array
    return accl_data

data_ Acquisition (FilePath) input the parameter FilePath and output a 1-dimensional array data. Here is a demonstration

file_path = r'E:/graduate student/pytorch/CSDN code/fault_diagnosis_signal_processing/Part 4-Envelope spectrum/1730_12k_0.007-InnerRace.mat'
xt = data_acquision(file_path)
plt.figure(figsize=(12,3))
plt.plot(xt)
print(xt)

Output results:

[ 0.22269856  0.09323776 -0.14651649 ... -0.36125573  0.31138814
  0.17055689]

3. Hilbert Demodulation (envelope spectrum) analysis

Hilbert Demodulation, also known as envelope spectrum analysis.

3.1 Hilbert Huang Transform

set up x ( t ) x(t) x(t) is a real-time domain signal, and its Hilbert transform is defined as:
h ( t ) = 1 π ∫ − ∞ + ∞ x ( τ ) t − τ d τ = x ( t ) ∗ 1 π t h(t)=\frac{1}{\pi} \int_{-\infty}^{+\infty} \frac{x(\tau)}{t-\tau} \mathrm{d} \tau=x(t) * \frac{1}{\pi t} h(t)=π1​∫−∞+∞​t−τx(τ)​dτ=x(t)∗πt1​
Then the original signal x ( t ) x(t) x(t) and its Hilbert transform signal h ( t ) h(t) h(t) can construct a new analytic signal z ( t ) z(t) z(t):
z ( t ) = x ( t ) + j h ( t ) = a ( t ) e j φ t z(t)=x(t)+j h(t)=a(t) e^{j \varphi t} z(t)=x(t)+jh(t)=a(t)ejφt

# step1: do Hilbert transform
ht = fftpack.hilbert(xt)
print(ht)

Output results:

[-0.02520403 -0.28707983 -0.00610516 ...  0.1100125   0.22821944
 -0.11203138]

Output at this time h ( t ) h(t) h(t) is the analytic signal a ( t ) a(t) Imaginary part coefficient of a(t)

yes z ( t ) z(t) Take the mode of z(t) and get its amplitude a ( t ) : a(t): a(t):

a ( t ) = ∣ z ( t ) ∣ = x 2 ( t ) + h 2 ( t ) {a(t)=|z(t)|=\sqrt{x^{2}(t)+h^{2}(t)}} a(t)=∣z(t)∣=x2(t)+h2(t) ​

Note: a ( t ) a(t) a(t) is envelope signal, also known as analytic signal

3.2 obtaining envelope spectrum

sampling_rate = 12000
am = np.fft.fft(at)   # fft is used to obtain the amplitude of at after Hilbert transform
am = np.abs(am)       # Calculate the absolute value of the amplitude (at this time, the absolute value is very large)
am = am/len(am)*2
am = am[0: int(len(am)/2)]
freq = np.fft.fftfreq(len(at), d=1 / sampling_rate)  # Obtain fft frequency, which includes positive frequency and negative frequency
freq = freq[0:int(len(freq)/2)]  # Get positive frequency
plt.plot(freq, am)


Observations:
(1) Where the frequency is 0Hz, the amplitude is relatively large
(2) In the low-frequency part, it seems to see the first and second harmonic

3.3 remove DC component

The amplitude at 0Hz is relatively high, which makes the amplitude of other frequencies low and inconvenient to observe. This phenomenon is called DC component, the method of removing DC component, y = y-mean(y)

sampling_rate = 12000
at = at - np.mean(at)  # To DC component
am = np.fft.fft(at)    # fft is used to obtain the amplitude of at after Hilbert transform
am = np.abs(am)        # Calculate the absolute value of the amplitude (at this time, the absolute value is very large)
am = am/len(am)*2
am = am[0: int(len(am)/2)]
freq = np.fft.fftfreq(len(at), d=1 / sampling_rate)  # Obtain fft frequency, which includes positive frequency and negative frequency
freq = freq[0:int(len(freq)/2)]  # Get positive frequency
plt.plot(freq, am)

Output results:

sampling_rate = 12000
at = at - np.mean(at)  # To DC component
am = np.fft.fft(at)    # fft is used to obtain the amplitude of at after Hilbert transform
am = np.abs(am)        # Calculate the absolute value of the amplitude (at this time, the absolute value is very large)
am = am/len(am)*2
am = am[0: int(len(am)/2)]
freq = np.fft.fftfreq(len(at), d=1 / sampling_rate)  # Obtain fft frequency, which includes positive frequency and negative frequency
freq = freq[0:int(len(freq)/2)]  # Get positive frequency
plt.figure(figsize=(12,3))
plt.plot(freq, am)
plt.xlim(0,500)

Output results:

4. Calculate fault characteristic frequency

Inner ring fault characteristic frequency: F B P F I = n f r 2 ( 1 + d D cos ⁡ α ) F_{\mathrm{BPFI}}=\frac{n f_{r}}{2}\left(1+\frac{d}{D} \cos \alpha\right) FBPFI​=2nfr​​(1+Dd​cosα)

Outer ring fault characteristic frequency: F B P F O = n f r 2 ( 1 − d D cos ⁡ α ) F_{\mathrm{BPFO}}=\frac{n f_{r}}{2}\left(1-\frac{d}{D} \cos \alpha\right) FBPFO​=2nfr​​(1−Dd​cosα)

Rolling element fault characteristic frequency: F B S F = D f r 2 d [ 1 − ( d D cos ⁡ α ) 2 ] F_{\mathrm{BSF}}=\frac{D f_{r}}{2 d}\left[1-\left(\frac{d}{D} \cos \alpha\right)^{2}\right] FBSF​=2dDfr​​[1−(Dd​cosα)2]

n n n: Number of scrolls, f r f_{r} fr: shaft speed d d d: Diameter of ball (sub) D D D: Bearing pitch diameter

Bearing model: 6205-2RSL JME SKF deep groove ball bearing
d d d=7.94mm, D D D=39.04mm, α \alpha α=0, n n n=9

4.1 define a bearing fault characteristic frequency calculation function

For convenience, a calculation function of bearing fault characteristic frequency is defined, which only needs to input parameters d d d, D D D, α \alpha α, n n n and f r f_{r} fr is enough

def bearing_fault_freq_cal(n, d, D, alpha, fr=None):
    '''
    Basic description:
        Calculate the fault characteristic frequency of rolling bearings
    Detailed description:
        Enter 4 parameters n, fr, d, D, alpha
    return C_bpfi, C_bpfo, C_bsf, C_ftf,  fr
           Inner ring    Outer ring    Speed of needle cage

    Parameters
    ----------
    n: integer
        The number of roller element
    fr: float(r/min)
        Rotational speed
    d: float(mm)
        roller element diameter
    D: float(mm)
        pitch diameter of bearing
    alpha: float(°)
        contact angle
    fr:: float(r/min)
        rotational speed
    Returns
    -------
    BPFI: float(Hz)
        Inner race-way fault frequency
    BPFO: float(Hz)
        Outer race-way fault frequency
    BSF: float(Hz)
        Ball fault frequency
    FTF: float(Hz)
        Cage frequency
    '''
    C_bpfi = n*(1/2)*(1+d/D*np.math.cos(alpha))
    C_bpfo = n*(1/2)*(1-(d/D)*np.math.cos(alpha))
    C_bsf = D*(1/(2*d))*(1-np.square(d/D*np.math.cos(alpha)))
    C_ftf = (1/2)*(1-(d/D)*np.math.cos(alpha))
    if fr!=None:
        return C_bpfi*fr/60, C_bpfo*fr/60, C_bsf*fr/60, C_ftf*fr/60, fr/60
    else:
        return C_bpfi, C_bpfo, C_bsf, C_ftf, fr

Next, calculate the fault characteristic frequency of CWRU at 1730rpm

bpfi, bpfo, bsf, ftf, fr = bearing_fault_freq_cal(n=9, alpha=0, d=7.94, D=39.04, fr=1730)
print('Inner ring fault characteristic frequency',bpfi)
print('Outer ring fault characteristic frequency',bpfo)
print('Rolling element fault characteristic frequency',bsf)
print(ftf)
print(fr)

5. Verification of theoretical fault characteristic frequency and actual fault characteristic frequency

sampling_rate = 12000
at = at - np.mean(at)  # To DC component
am = np.fft.fft(at)    # fft is used to obtain the amplitude of at after Hilbert transform
am = np.abs(am)        # Calculate the absolute value of the amplitude (at this time, the absolute value is very large)
am = am/len(am)*2
am = am[0: int(len(am)/2)]
freq = np.fft.fftfreq(len(at), d=1 / sampling_rate)  # Obtain fft frequency, which includes positive frequency and negative frequency
freq = freq[0:int(len(freq)/2)]  # Get positive frequency
plt.figure(figsize=(12,3))
plt.plot(freq, am)
plt.xlim(0,500)
plt.vlines(x=156.13, ymin=0, ymax=0.2, colors='r')  # Frequency doubling
plt.vlines(x=156.13*2, ymin=0, ymax=0.2, colors='r')  # Double frequency

Output results:

Red is the theoretical inner circle fault characteristic frequency, and the blue line is the actual fault characteristic frequency. Although it is not completely coincident, this is within the allowable range. Therefore, in the actual situation, if the envelope spectrum appears, the inner ring fault of the bearing can occur.

6. Compared with fft

Can we see the fault characteristic frequency without Hilbert transform demodulation and direct fft? The following is a comparative analysis

sampling_rate = 12000
am = np.fft.fft(xt)   # fft is used to obtain the amplitude of at after Hilbert transform
am = np.abs(am)       # Calculate the absolute value of the amplitude (at this time, the absolute value is very large)
am = am/len(am)*2
am = am[0: int(len(am)/2)]
freq = np.fft.fftfreq(len(xt), d=1 / sampling_rate)  # Obtain fft frequency, which includes positive frequency and negative frequency
freq = freq[0:int(len(freq)/2)]  # Get positive frequency
plt.plot(freq, am)

Output results:

It can be seen that there are more fft high frequencies of the original inner ring fault signal

plt.plot(freq, am)
plt.xlim(0, 500)
plt.vlines(x=156.13, ymin=0, ymax=0.05, colors='r')  # Frequency doubling
plt.vlines(x=156.13*2, ymin=0, ymax=0.05, colors='r')  # Double frequency


It can be seen that the fault characteristic frequency amplitude is low and is covered by other frequency amplitudes in the case of direct fft. In turn, Hilbert Demodulation can make it easier to observe the low frequency of fault characteristic frequency.

7. Envelope spectral function

In order to use envelope spectrum more conveniently, an envelope spectrum function PLT is encapsulated here_ envelope_ spectrum

def plt_envelope_spectrum(data, fs, xlim=None, vline= None):
    '''
    fun: Draw envelope spectrum
    param data: Input data, 1D array
    param fs: sampling frequency 
    param xlim: Abscissa of picture xlim,default = None
    param vline: Picture vertical line, default = None
    '''
    #----To DC component----#
    data = data - np.mean(data)
    #----Do Hilbert transform----#
    xt = data
    ht = fftpack.hilbert(xt)
    at = np.sqrt(xt**2+ht**2)   # Obtain the analytic signal at = sqrt(xt^2 + ht^2)
    am = np.fft.fft(at)         # Obtain the amplitude by fft transformation of analytical signal at
    am = np.abs(am)             # Calculate the absolute value of the amplitude (at this time, the absolute value is very large)
    am = am/len(am)*2
    am = am[0: int(len(am)/2)]  # Take positive frequency amplitude
    freq = np.fft.fftfreq(len(at), d=1 / fs)  # Obtain fft frequency, which includes positive frequency and negative frequency
    freq = freq[0:int(len(freq)/2)]  # Get positive frequency
    plt.plot(freq, am)
    if vline:  # Whether to draw a vertical line
        plt.vlines(x=vline, ymax=0.2, ymin=0, colors='r')  # Height y 0-0.2, color red
    if xlim: # Whether xlim is set for the abscissa of the picture
        plt.xlim(0, xlim)  
    plt.xlabel('freq(Hz)')    # Abscissa label
    plt.ylabel('amp(m/s2)')   # Ordinate label

7.1 outer ring fault data test

file_path = r'E:/graduate student/pytorch/CSDN code/fault_diagnosis_signal_processing/Part 4-Envelope spectrum/1730_12k_0.007-OuterRace3.mat'
data = data_acquision(file_path)
plt_envelope_spectrum(data = data, fs=12000, xlim=300, vline=bpfo)

Output results:

It can be seen that the actual outer ring fault characteristic frequency is relatively obvious

7.2 test and analysis of rolling element fault data

file_path = r'E:/graduate student/pytorch/CSDN code/fault_diagnosis_signal_processing/Part 4-Envelope spectrum/1730_12k_0.014-Ball.mat'
data = data_acquision(file_path)
plt_envelope_spectrum(data = data, fs=12000, xlim=300, vline=bsf)

It can be seen that the actual rolling element fault characteristic frequency is not obvious

8. Summary

(1) Envelope spectrum can detect inner ring and outer ring faults, which is difficult for rolling elements
(2) Using fft directly is difficult to detect the fault characteristic frequency, and the fault characteristic frequency is easy to be covered up by high frequency

Tags: Python programming language

Posted by jmantra on Thu, 28 Jul 2022 02:17:42 +0930