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+Ddcosα)

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−Ddcosα)

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−(Ddcosα)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