# A Method for Determining Intrinsic Mode Function Number in Variational Mode Decomposition and Its Application to Bearing Vibration Signal Processing

Authors: Shoujun Wu; Fuzhou Feng; Junzhen Zhu; Chunzhi Wu; Guang Zhang
Journal: Shock and Vibration
 (2020)

---

## Abstract

Variational mode decomposition (VMD) method has been widely used in the field of signal processing with significant advantages over other decomposition methods in eliminating modal aliasing and noise robustness. The number (usually denoted byK) of intrinsic mode function (IMF) has a great influence on decomposition results. When dealing with signals including complex components, it is usually impossible for the existing methods to obtain correct results and also effective methods for determining K value are lacking. A method called center frequency statistical analysis (CFSA) is proposed in this paper to determine K value. CFSA method can obtain K value accurately based on center frequency histogram. To shed further light on its performance, we analyze the behavior of CFSA method with simulation signal in the presence of variable components amplitude, components frequency, and components number as well as noise amplitude. The normal and fault vibration signals obtained from a bearing experimental setup are used to verify the method. Compared with maximum center frequency observation (MCFO), correlation coefficient (CC), and normalized mutual information (NMI) methods, CFSA is more robust and accurate, and the center frequencies results are consistent with the main frequencies in FFT spectrum.


---

## Body

## 1. Introduction

Variational mode decomposition (VMD) has been used in various applications since it was proposed by Dragomiretskiy and Zosso [1]. VMD uses completely nonrecursive decomposition method to get intrinsic mode functions (IMF), which has a stronger antinoise ability than empirical mode decomposition (EMD) method. The number of IMF by EMD cannot be artificially set, while VMD can do so. The low-frequency component of VMD is easier to express the general fluctuation trend of the original signal, but it is difficult for EMD to observe this characteristic. In addition, VMD has great advantages in processing modal aliasing and high noise signal compared with EMD, ensemble empirical mode decomposition (EEMD) and empirical wavelet transform (EWT) [2, 3].When VMD is used to decompose a signal, there are six parameters (α,τ,K,DC,init,ε) that need to be determined. Here, α represents penalty factor, τ represents fidelity coefficient, K represents number of IMF, DC represents the updating parameter of first center frequency, init represents the initialization parameter of center frequency, and ε represents convergence threshold. Normally, the parameters τ,DC,init,ε can be assigned to the default value as τ=0,DC=0or1,init=1or2, ε=1e−7. However, the parameters K and α need to be optimized according to the actual signal. The greater the value of α is, the faster the attenuation is on both sides of the center frequency. α is mainly determined according to the principle of avoiding aliasing between mode functions and is generally 1/6∼2 times of sampling frequency Its specific value needs to be further determined according to the characteristics of the signal [4]. There is no unified criterion or method to determine the value of K so far [5]. Many researches have been done in recent years. A large number of methods have been proposed. At present, determination methods of IMF number in VMD can be divided into three categories: (1) methods based on center frequency observation, (2) methods based on threshold criteria, and (3) other methods.(1)
The first type of method forK value determination is to observe the center frequency or the IMF spectrum distribution. Liu [6] proposed maximum central frequency observation (MCFO) method for determining the value of K. When there are two center frequency sets similar to each other, the value of K is the best. Zheng et al. [7, 8] determined the K value according to the spectrum distribution of IMF components. The principle is that the spectrum of each IMF neither overlaps nor loses frequency information when the suitable K value is selected. This kind of method is practicable to some extent but lacks quantitative analysis.(2)
The second type of method ofK value determination is threshold criterion, such as correlation coefficient, mutual information, kurtosis, and distance measurement. In [9], the K value is determined by whether the correlation coefficient (CC) between the reconstructed signal and the original signal reaches the threshold or not. Zhang et al. [10] proposed the correlation and energy ratio to determine the number of IMF. Liu et al. [11] used the normalized mutual information (NMI) between the IMF and the original signal to judge whether the K value is appropriate or not. Wang et al. [12] determined K based on the number of IMF whose permutation entropy is greater than the threshold. Huang et al. [13] determined K based on proposed normalized distance (ND) indicator, which is used to describe the similarity degree between the reconstructed signal and original signal. The threshold of this kind of methods is difficult to determine and has strong subjectivity. Zhao and Li [14] used kurtosis to select sensitive IMFs. Kurtosis criterion is only suitable for processing signals with periodic big shock components. However, it is not suitable for the processing of signals without impact or signals with weak impact.(3)
There are other ways in the third type to determine the value ofK. Li et al. [15] proposed a method in which the value of K is equal to the number of IMF of EMD on the same signal. Liu et al. [16] selected the number of K based on detrended fluctuation analysis (DFA). Feng et al. [17] proposed a method to determine the IMF number, which is based on sampling frequency divided by two times meshing frequency in a gearbox. Zhang et al. [18] used the number of peaks in the envelope spectrum of product function (PF) components from local mean decomposition (LMD) to determine the value of K. Wang et al. [19] took the minimum average envelope entropy as the fitness function and used the particle swarm optimization algorithm to optimize the parameter K and penalty factor. The algorithm of this kind of method is more complex.In order to determine the IMF number simply and accurately while decomposing a complex signal by VMD, CFSA method is proposed in this work. The rest of the paper is organized as follows: Section2 introduces the algorithm of VMD, MCFO, CC, NMI, and CFSA method. In Section 3, the classic methods and CFSA method are studied by constructed simulation signals. In Section 4, the classic methods and CFSA method are applied to rolling bearing experimental signals processing. Conclusions are drawn in Section 5.


## 2. Theoretical Introduction

### 2.1. Variational Mode Decomposition

VMD can nonrecursively decompose a real-valued multicomponent signalf into a discrete number of quasi-orthogonal band-limited subsignals uk with specific sparsity properties in the spectral domain. Each mode is compact around a center pulsation ωk and its bandwidth is estimated using H1 Gaussian smoothness of the shifted signal. The VMD is written as a constrained variational problem:(1)minuk,ωk∑k=1K∂tδt+jπt∗ukte−jωkt22s.t.∑k=1Kuk=f,where uk and ωk are the kth intrinsic mode function and its center frequency, respectively.Equation (1) can be solved by introducing a quadratic penalty and Lagrangian multipliers. The augmented Lagrangian is given as follows:(2)Luk,ωk,λ=α∑k∂tδt+jπt∗ukte−jωkt22+ft−∑kukt22+λt,ft−∑kukt,where α denotes the balancing parameter of the data-fidelity constraint.Equation (2) is then solved with the alternate direction method of multipliers (ADMM). All the modes gained from solutions in spectral domain are written as(3)u^kn+1ω=f^ω−∑i≠ku^iω+λ^ω/21+2αω−ωk2,where ωk is a frequency corresponding to the center of gravity of power spectrum of the corresponding IMF. Thus, Wiener filtering is embedded in the VMD algorithm that makes it much more robust to sampling and noise. The update equation for the center frequency is expressed as(4)ωkn+1=∫0∞ωu^kω2dω∫0∞u^kω2dω,Complete algorithm of the VMD can be found in detail in [1]. The purpose of VMD for signals is to obtain several IMFs. How to determine the number of IMFs is crucial for the decomposition results. Next, several representative methods for determining the number of IMFs are introduced, and a novel method is proposed through research.


### 2.2. Traditional Methods to DetermineK

#### 2.2.1. Maximum Center Frequency Observation Method

The principle of maximum center frequency observation (MCFO) method is to observe the trend of maximum center frequencies. Maximum center frequency of each mode component is increasing gradually with increment of the mode number. When maximum center frequencies tend to be stable, the value ofK can be determined.


#### 2.2.2. Correlation Coefficient Method

Correlation coefficient (CC) between the mode components and the original signal can be obtained by (5). Decomposition will be stopped when the minimum correlation coefficient is less than the given threshold, and then the value of K can be determined:(5)ρuk,f=Euk⋅f−EukEfEuk2−E2ukEf2−E2f,where ρuk,f represents the CC between the IMF and original signal; f and uk represent the original signal and IMF obtained by VMD; E⋅ corresponds to the mathematical expectation.


#### 2.2.3. Normalized Mutual Information Method

Mutual information (MI) value of the IMF and the original signal is calculated by (6). Decomposition will be stopped when the minimum value of NMI is less than the given threshold, and thus the value of K can be determined:(6)MIX,Y=HY−HYX,where HY represents information entropy of Y; HYX represents the conditional entropy of Y when X is known. The stronger the correlation between X and Y, the smaller the conditional entropy value HYX, and the larger the mutual information MIX,Y. Normalized mutual information (NMI) is expressed as(7)δi=MIimaxMIi,where i represents the serial number of the IMF.


### 2.3. Center Frequency Statistical Analysis

Center frequency statistical analysis (CFSA) method is proposed in this work. Its main idea is to observe the number of frequencies higher than the mean value in the center frequency histogram, which is considered asK value. The detailed steps are as follows. (1) Initialize VMD parameters. (2) Decompose the signal by VMD to obtain IMFs. (3) Calculate the center frequency of IMFs using FFT. (4) Draw center frequency histogram. (5) Calculate average of center frequency and count number (N) of center frequencies that are higher than their mean value. K value plus 1, return to step (2). When the N value no longer increases, the decomposition is stopped, and the N value is the best number of IMF. The flowchart of CFSA method is shown in Figure 1.Figure 1
Flowchart of CFSA method.


## 2.1. Variational Mode Decomposition

VMD can nonrecursively decompose a real-valued multicomponent signalf into a discrete number of quasi-orthogonal band-limited subsignals uk with specific sparsity properties in the spectral domain. Each mode is compact around a center pulsation ωk and its bandwidth is estimated using H1 Gaussian smoothness of the shifted signal. The VMD is written as a constrained variational problem:(1)minuk,ωk∑k=1K∂tδt+jπt∗ukte−jωkt22s.t.∑k=1Kuk=f,where uk and ωk are the kth intrinsic mode function and its center frequency, respectively.Equation (1) can be solved by introducing a quadratic penalty and Lagrangian multipliers. The augmented Lagrangian is given as follows:(2)Luk,ωk,λ=α∑k∂tδt+jπt∗ukte−jωkt22+ft−∑kukt22+λt,ft−∑kukt,where α denotes the balancing parameter of the data-fidelity constraint.Equation (2) is then solved with the alternate direction method of multipliers (ADMM). All the modes gained from solutions in spectral domain are written as(3)u^kn+1ω=f^ω−∑i≠ku^iω+λ^ω/21+2αω−ωk2,where ωk is a frequency corresponding to the center of gravity of power spectrum of the corresponding IMF. Thus, Wiener filtering is embedded in the VMD algorithm that makes it much more robust to sampling and noise. The update equation for the center frequency is expressed as(4)ωkn+1=∫0∞ωu^kω2dω∫0∞u^kω2dω,Complete algorithm of the VMD can be found in detail in [1]. The purpose of VMD for signals is to obtain several IMFs. How to determine the number of IMFs is crucial for the decomposition results. Next, several representative methods for determining the number of IMFs are introduced, and a novel method is proposed through research.


## 2.2. Traditional Methods to DetermineK

### 2.2.1. Maximum Center Frequency Observation Method

The principle of maximum center frequency observation (MCFO) method is to observe the trend of maximum center frequencies. Maximum center frequency of each mode component is increasing gradually with increment of the mode number. When maximum center frequencies tend to be stable, the value ofK can be determined.


### 2.2.2. Correlation Coefficient Method

Correlation coefficient (CC) between the mode components and the original signal can be obtained by (5). Decomposition will be stopped when the minimum correlation coefficient is less than the given threshold, and then the value of K can be determined:(5)ρuk,f=Euk⋅f−EukEfEuk2−E2ukEf2−E2f,where ρuk,f represents the CC between the IMF and original signal; f and uk represent the original signal and IMF obtained by VMD; E⋅ corresponds to the mathematical expectation.


### 2.2.3. Normalized Mutual Information Method

Mutual information (MI) value of the IMF and the original signal is calculated by (6). Decomposition will be stopped when the minimum value of NMI is less than the given threshold, and thus the value of K can be determined:(6)MIX,Y=HY−HYX,where HY represents information entropy of Y; HYX represents the conditional entropy of Y when X is known. The stronger the correlation between X and Y, the smaller the conditional entropy value HYX, and the larger the mutual information MIX,Y. Normalized mutual information (NMI) is expressed as(7)δi=MIimaxMIi,where i represents the serial number of the IMF.


## 2.2.1. Maximum Center Frequency Observation Method

The principle of maximum center frequency observation (MCFO) method is to observe the trend of maximum center frequencies. Maximum center frequency of each mode component is increasing gradually with increment of the mode number. When maximum center frequencies tend to be stable, the value ofK can be determined.


## 2.2.2. Correlation Coefficient Method

Correlation coefficient (CC) between the mode components and the original signal can be obtained by (5). Decomposition will be stopped when the minimum correlation coefficient is less than the given threshold, and then the value of K can be determined:(5)ρuk,f=Euk⋅f−EukEfEuk2−E2ukEf2−E2f,where ρuk,f represents the CC between the IMF and original signal; f and uk represent the original signal and IMF obtained by VMD; E⋅ corresponds to the mathematical expectation.


## 2.2.3. Normalized Mutual Information Method

Mutual information (MI) value of the IMF and the original signal is calculated by (6). Decomposition will be stopped when the minimum value of NMI is less than the given threshold, and thus the value of K can be determined:(6)MIX,Y=HY−HYX,where HY represents information entropy of Y; HYX represents the conditional entropy of Y when X is known. The stronger the correlation between X and Y, the smaller the conditional entropy value HYX, and the larger the mutual information MIX,Y. Normalized mutual information (NMI) is expressed as(7)δi=MIimaxMIi,where i represents the serial number of the IMF.


## 2.3. Center Frequency Statistical Analysis

Center frequency statistical analysis (CFSA) method is proposed in this work. Its main idea is to observe the number of frequencies higher than the mean value in the center frequency histogram, which is considered asK value. The detailed steps are as follows. (1) Initialize VMD parameters. (2) Decompose the signal by VMD to obtain IMFs. (3) Calculate the center frequency of IMFs using FFT. (4) Draw center frequency histogram. (5) Calculate average of center frequency and count number (N) of center frequencies that are higher than their mean value. K value plus 1, return to step (2). When the N value no longer increases, the decomposition is stopped, and the N value is the best number of IMF. The flowchart of CFSA method is shown in Figure 1.Figure 1
Flowchart of CFSA method.


## 3. Simulation Studies

In this section, the four factors that influenceK were studied by simulation analysis. The four factors are components amplitude, components frequency, and components number and noise amplitude.


### 3.1. Simulation Signal

The simulation signalxt is synthesized by several sinusoidal signals with the following expression:(8)xt=x1t+x2t+x3t+x4t+x5t+x6t+xnt,x1t=A1⋅cos2π⋅f1⋅t,x2t=A2⋅cos2π⋅f2⋅t,x3t=A3⋅cos2π⋅f3⋅t,x4t=A4⋅cos2π⋅f4⋅t,x5t=A5⋅cos2π⋅f5⋅t,x6t=A6⋅cos2π⋅f6⋅t,xnt=An⋅η,where x1t,x2t,x3t,x4t,x5t,x6t are six sinusoidal signal components, A1,A2,A3,A4,A5,A6 are the amplitude of the corresponding components, respectively. f1,f2,f3,f4,f5,f6 are the frequencies of the components, respectively. xnt is the noise signal. An is the amplitude of the noise signal, and η∼0,1 is Gaussian white noise. Simulation signals under different influencing factors are shown in Table 1.Table 1
The value of the parameters related to the simulation signal.
Signal code
Factor
A
1
A
2
A
3
A
4
A
5
A
6
A
n
f
1
f
2
f
3
f
4
f
5
f
6
S1
Components amplitude
1
1/4
1/16
S2
1
5/12
25/144
0
0
0
0.1
2
24
288
0
0
0
S3
1
7/12
49/144
S4
1
1/3
1/9
2
12
72
S5
1
1/2
1/4
0
0
0
0.1
2
16
128
S6
Components frequency
0
0
0
S7
2
20
200
S8
2
24
288
S9
1
1/4
1/16
0
0
0
0.1
2
24
288
0
0
0
S10
Components number
1
1/4
1/16
1/8
0
0
2
24
288
100
0
0
S11
1
1/4
1/16
1/8
1/10
0
2
24
288
100
180
0
S12
1
1/4
1/16
1/8
1/10
1/12
2
24
288
100
180
240
S13
0
S14
Noise amplitude
1
1/4
1/16
0
0
0
0.05
2
24
288
0
0
0
S15
0.1
S16
0.15


### 3.2. Simulation Results Analysis

#### 3.2.1. Influence of Component Amplitude

Signals S1 to S4 were analyzed by methods mentioned above in this section. VMD parameters exceptK are initialized as follows: α=2000,τ=0,DC=0,init=1,ε=1e−7.According to MCFO method, the maximum center frequency trend under different component amplitudes is shown in Figure2. As can be seen, the maximum center frequency gradually stabilizes after the value of K becomes bigger than 6. Thus, six IMFs were obtained by the MCFO. However, the simulation signals only contain three components, and the decomposition result is inconsistent with the actual components number. So the method is not accurate in acquiring the IMF number of VMD.Figure 2
Maximum center frequency trend under different component amplitudes.According to CC method, correlation coefficients under different component amplitudes are shown in Table2. For a more intuitive analysis of the results, threshold ranges are plotted in Figure 3. It can be seen that when the component amplitudes are different, range of the correlation coefficient threshold body has different values, but there is an intersection (0.0649, 0.0823) area in each threshold range. When the correlation coefficient is greater than any value in the threshold intersection, three IMFs number can be obtained, which is consistent with the component number of the original signal.Table 2
CC under different component amplitudes.
Signal code
K
Correlation coefficient
Threshold range
S1
3
0.9599
0.2460
0.0837
(0.0627, 0.0823)
4
0.9599
0.2460
0.0823
0.0627
S2
3
0.9044
0.3786
0.1664
(0.0592, 0.1660)
4
0.9044
0.3786
0.1660
0.0592
S3
3
0.8224
0.4825
0.2856
(0.0571, 0.2849)
4
0.8224
0.4824
0.2849
0.0571
S4
3
0.9360
0.3158
0.1127
(0.0649, 0.1110)
4
0.9360
0.3158
0.1110
0.0649Figure 3
CC threshold range under different component amplitudes.According to NMI method, normalized mutual information under different component amplitudes is shown in Table3. It can be seen that there is no threshold range for the signal S4. Because when K = 4, the minimum value of NMI is 0.0658, which is larger than the minimum value of NMI 0.0579 when K = 3. In order to get a clearer understanding of the results, threshold ranges are plotted in Figure 4. It can be seen from analysis results of signal S4 that there is no threshold range. Therefore, there is no threshold intersection to determine the IMF number. So it is invalid to use NMI method for this signal.Table 3
NMI under different component amplitudes.
Signal code
K
Normalized mutual information
Threshold range
S1
3
1.0000
0.1692
0.0663
(0.0641, 0.0663)
4
1.0000
0.1685
0.0743
0.0641
S2
3
1.0000
0.2506
0.0930
(0.0605, 0.0930)
4
1.0000
0.2504
0.1002
0.0605
S3
3
1.0000
0.3442
0.1902
(0.1644, 0.1902)
4
1.0000
0.3439
0.1644
0.0756
S4
3
1.0000
0.1781
0.0579
No
4
1.0000
0.1930
0.0658
0.0833Figure 4
NMI threshold under different component amplitudes.According to CFSA method, center frequency histogram under different component amplitudes is shown in Figure5. It can be seen that the dominant center frequencies are about 12.6 Hz, 150.8 Hz, and 1809.6 Hz, respectively. There are three dominant frequencies. In this way, 3 IMFs can be determined, which is consistent with the component number of original signal. The results show that the proposed method is effective for this kind of signal.Figure 5
Center frequency histogram under different component amplitudes. (a) S1, (b) S2, (c) S3, and (d) S4.
(a)
(b)
(c)
(d)


#### 3.2.2. Influence of Component Frequency

Signals S5 to S8 were analyzed using four methods in this section. VMD parameters exceptK are initialized as follows.According to MCFO method, the maximum center frequency trend under different component frequencies is shown in Figure6. It can be seen that the maximum center frequency is different when the value of K is from 1 to 4. Center frequencies of the four signals tend to be stable after the value of K becomes bigger than 6 and no longer increases significantly. Therefore, the IMF number is 6; however, the results is inconsistent with the component number of original signal.Figure 6
Maximum center frequency trend under different component frequencies.According to CC method, correlation coefficients under different component frequencies are shown in Table4. The threshold range results are plotted in Figure 7. As can be seen, the intersection of the threshold ranges is 0.0591, 0.1007. When the correlation coefficient threshold is greater than any value in the intersection, it can be determined that the IMF number is 3. The results is consistent with the actual components number of simulated signal.Table 4
CC under different component frequencies.
Signal code
K
Correlation coefficient
Threshold range
S5
3
0.9250
0.2245
0.1007
(0.0578, 0.1007)
4
0.8680
0.4345
0.2189
0.0578
S6
3
0.8649
0.4436
0.2149
(0.0536, 0.2146)
4
0.8649
0.4436
0.2146
0.0536
S7
3
0.9502
0.2080
0.2160
(0.0591, 0.2080)
4
0.8658
0.4403
0.2159
0.0591
S8
3
0.8660
0.4370
0.2242
(0.0540, 0.2239)
4
0.8660
0.4370
0.2239
0.0540Figure 7
CC threshold range under different component frequencies.According to NMI method, normalized mutual information under different component frequencies is shown in Table5. It can be seen that there is no threshold range for signal S7. Because when K = 4, the minimum value of NMI is 0.0614, which is larger than 0.0599, the minimum value of NMI when K = 3. It can be seen from Figure 8 that there is no threshold range intersection, so the K value cannot be determined according to this method.Table 5
NMI under different component frequencies.
Signal code
K
Normalized mutual information
Threshold range
S5
3
1.0000
0.1677
0.1452
(0.1426, 0.1452)
4
1.0000
0.3034
0.2094
0.1426
S6
3
1.0000
0.2993
0.2295
(0.0956, 0.1900)
4
1.0000
0.3003
0.1900
0.0956
S7
3
1.0000
0.1097
0.0599
No
4
1.0000
0.2744
0.1185
0.0614
S8
3
1.0000
0.2864
0.1257
(0.0911, 0.1257)
4
1.0000
0.2675
0.1373
0.0911Figure 8
NMI threshold range under different component frequencies.According to the CFSA method, center frequency histogram under different component amplitudes is shown in Figure9. Three dominant center frequencies are obtained by signal S5 processing, which are approximately 12.6 Hz, 75.4 Hz, and 452.3 Hz, respectively. Another three frequencies (12.6 Hz, 100.5 Hz, and 804.3 Hz) can be obtained by processing signal S6. 12.6 Hz, 125.7 Hz, and 1255.3 Hz can be obtained by processing signal S7. 12.6 Hz, 151 Hz, and 1809.7 Hz can be obtained by processing signal S8. In summary, three IMF numbers can be obtained by the CFSA method, which is consistent with the actual component number of simulation signal.Figure 9
Center frequency histogram under different component frequencies, (a) S5, (b) S6, (c) S7, and (d) S8.
(a)
(b)
(c)
(d)


#### 3.2.3. Influence of Component Number

Signals S9 to S12 were analyzed by four methods in this section. VMD parameters exceptK are initialized as follows.According to MCFO method, the maximum center frequency trend under different components number is shown in Figure10. It can be seen that when the simulation signal contains 3, 4, and 5 components, the maximum center frequency trends are almost the same. However, when the simulation contains 6 components, the maximum center frequency trend is different from other situations. All maximum values of center frequency tend to be stable after the value of K becomes bigger than 7. Therefore, the IMF number is 7. But the results are inconsistent with the actual component number of the simulation signal.Figure 10
Maximum center frequency trend under different components number.According to the CC method, the correlation coefficients under different components number are shown in Table6. The threshold ranges are plotted in Figure 11. As can be seen, there is no threshold intersection for determining the IMF number. Therefore, when VMD is used to decompose the signals S9 to S12, it is difficult for CC method to determine the IMF number.Table 6
CC under different components number.
Signal code
K
Correlation coefficient
Threshold range
S9
3
0.9825
0.0843
0.0899
(0.0611, 0.0843)
4
0.9586
0.2496
0.0854
0.0611
S10
4
0.9750
0.1556
0.0846
0.0617
(0.0611, 0.0617)
5
0.9749
0.1547
0.0704
0.0831
0.0611
S11
5
0.9711
0.1479
0.1163
0.0821
0.0639
(0.0573, 0.0639)
6
0.9710
0.1477
0.1160
0.0798
0.0607
0.0573
S12
6
0.9673
0.1492
0.1183
0.0916
0.0821
0.0582
(0.0538, 0.0557)
7
0.9673
0.1491
0.1181
0.0910
0.0800
0.0538
0.0557Figure 11
CC threshold range under different components number.According to the NMI method, the normalized mutual information under different components number is shown in Table7. It can be seen from Figure 12 that there is no threshold range for the signal S9, because the minimum value of NMI is 0.0502, which is larger than 0.0436. Threshold ranges are plotted in Figure 13. As can be seen, there is no threshold intersection. So the NMI method cannot determine the IMF number when the signals S9 to S12 are processed by VMD.Table 7
NMI under different components number.
Signal code
K
Normalized mutual information
Threshold range
S9
3
1.0000
0.0508
0.0436
No
4
1.0000
0.0560
0.0543
0.0502
S10
4
1.0000
0.0772
0.0472
0.0577
(0.0429, 0.0472)
5
1.0000
0.0728
0.0688
0.0695
0.0429
S11
5
1.0000
0.0734
0.0610
0.0530
0.0670
(0.0458, 0.0469)
6
1.0000
0.0694
0.0607
0.0458
0.0469
0.0500
S12
6
1.0000
0.0740
0.0544
0.0530
0.0534
0.0594
(0.0500, 0.0511)
7
1.0000
0.0771
0.0654
0.0511
0.0678
0.0500
0.0544Figure 12
NMI threshold range under different components number.Figure 13
Center frequency histogram under different components number, (a) S9, (b) S10, (c) S11, and (d) S12.
(a)
(b)
(c)
(d)According to CFSA, center frequency histogram under different components number is shown in Figure13. From the center frequency histogram obtained by processing the signal S9, there are three dominant center frequencies (12.6 Hz, 151.1 Hz, and 1808.5 Hz). And there are four dominant center frequencies (12.6 Hz, 150.9 Hz, 627.9 Hz, and 1809.8 Hz) obtained by signal S10 processing. Five dominant center frequencies (12.6 Hz, 150.9 Hz, 627.5 Hz, 1129 Hz, and 1808.7 Hz) can be obtained by processing signal S11. Six dominant center frequencies (12.6 Hz, 150.8 Hz, 628.3 Hz, 1131 Hz, 1508 Hz, and 1809.6 Hz) can be obtained by processing signal S12. The simulation signals S9 to S12 contain exactly 3 to 6 components, respectively, so the CFSA is correctly verified.


#### 3.2.4. Influence of Noise Amplitude

Signals S13 to S16 were analyzed by four methods in this section. VMD parameters exceptK are initialized as follows.According to MCFO method, the maximum center frequency trend under different noise amplitudes is shown in Figure14. It can be seen that when there is no noise, the maximum center frequency tends to be stable after the value of K becomes bigger than 2. When the noise amplitude is 0.05 and 0.1, the maximum center frequency tends to be stable from K = 6. When the noise amplitude is 0.15, the maximum center frequency tends to be stable from K = 7. Therefore, when the noise amplitude is different, IMF number determined by MCFO method is different.Figure 14
Maximum center frequency trend under different noise amplitudes.According to the CC method, the correlation coefficients under different noise amplitudes are shown in Table8. Different threshold ranges are shown in Figure 15. There is no intersection in the four threshold ranges, so the value of K cannot be determined. Therefore, when VMD is applied to decompose signals S13 to S16, it is not suitable for CC method to determine the IMF number.Table 8
CC under different noise amplitudes.
Signal code
K
Correlation coefficient
Threshold range
S13
3
0.9684
0.2421
0.0604
(0.0237, 0.0596)
4
0.9684
0.2421
0.0596
0.0237
S14
3
0.9673
0.2392
0.0663
(0.0307, 0.0657)
4
0.9673
0.2392
0.0657
0.0307
S15
3
0.9825
0.0850
0.0871
(0.0612, 0.0835)
4
0.9588
0.2488
0.0835
0.0612
S16
3
0.9742
0.1121
0.1049
(0.0911, 0.1004)
4
0.9741
0.1104
0.1004
0.0911Figure 15
CC threshold range under different noise amplitudes.According to the NMI method, the normalized mutual information under different noise amplitudes is shown in Table9. Four threshold ranges are shown in Figure 16. There is no intersection of the four threshold ranges, so there is no mutual information threshold to determine the K. Therefore, the NMI method is ineffective in this case.Table 9
NMI under different noise amplitudes.
Signal code
K
Normalized mutual information
Threshold range
S13
3
1.0000
0.2450
0.1165
(0.0753, 0.1165)
4
1.0000
0.2450
0.0753
0.1513
S14
3
1.0000
0.1688
0.0452
(0.0431, 0.0452)
4
1.0000
0.1784
0.0431
0.0505
S15
3
1.0000
0.1716
0.0708
(0.0634, 0.0708)
4
1.0000
0.1672
0.0634
0.0552
S16
3
1.0000
0.1439
0.0586
(0.0532, 0.0586)
4
1.0000
0.0693
0.0729
0.0532Figure 16
NMI threshold range under different noise amplitudes.According to CFSA, center frequency histogram under different noise amplitudes is shown in Figure17. As can be seen, there are three dominant center frequencies in results of signals S13 to S16, which corresponds to the actual components number of signals.Figure 17
Center frequency histogram under different noise amplitudes, (a) S13, (b) S14, (c) S15, and (d) S16.
(a)
(b)
(c)
(d)In summary, the MCFO and NMI methods cannot accurately determine the IMF number in case of four influencing factors changes. CC method can only determine the IMF number accurately under first two influencing factors changes. The proposed CFSA method can obtain the IMF number accurately under the variation of four influencing factors, which verifies the advantages and effectiveness of the proposed method.


## 3.1. Simulation Signal

The simulation signalxt is synthesized by several sinusoidal signals with the following expression:(8)xt=x1t+x2t+x3t+x4t+x5t+x6t+xnt,x1t=A1⋅cos2π⋅f1⋅t,x2t=A2⋅cos2π⋅f2⋅t,x3t=A3⋅cos2π⋅f3⋅t,x4t=A4⋅cos2π⋅f4⋅t,x5t=A5⋅cos2π⋅f5⋅t,x6t=A6⋅cos2π⋅f6⋅t,xnt=An⋅η,where x1t,x2t,x3t,x4t,x5t,x6t are six sinusoidal signal components, A1,A2,A3,A4,A5,A6 are the amplitude of the corresponding components, respectively. f1,f2,f3,f4,f5,f6 are the frequencies of the components, respectively. xnt is the noise signal. An is the amplitude of the noise signal, and η∼0,1 is Gaussian white noise. Simulation signals under different influencing factors are shown in Table 1.Table 1
The value of the parameters related to the simulation signal.
Signal code
Factor
A
1
A
2
A
3
A
4
A
5
A
6
A
n
f
1
f
2
f
3
f
4
f
5
f
6
S1
Components amplitude
1
1/4
1/16
S2
1
5/12
25/144
0
0
0
0.1
2
24
288
0
0
0
S3
1
7/12
49/144
S4
1
1/3
1/9
2
12
72
S5
1
1/2
1/4
0
0
0
0.1
2
16
128
S6
Components frequency
0
0
0
S7
2
20
200
S8
2
24
288
S9
1
1/4
1/16
0
0
0
0.1
2
24
288
0
0
0
S10
Components number
1
1/4
1/16
1/8
0
0
2
24
288
100
0
0
S11
1
1/4
1/16
1/8
1/10
0
2
24
288
100
180
0
S12
1
1/4
1/16
1/8
1/10
1/12
2
24
288
100
180
240
S13
0
S14
Noise amplitude
1
1/4
1/16
0
0
0
0.05
2
24
288
0
0
0
S15
0.1
S16
0.15


## 3.2. Simulation Results Analysis

### 3.2.1. Influence of Component Amplitude

Signals S1 to S4 were analyzed by methods mentioned above in this section. VMD parameters exceptK are initialized as follows: α=2000,τ=0,DC=0,init=1,ε=1e−7.According to MCFO method, the maximum center frequency trend under different component amplitudes is shown in Figure2. As can be seen, the maximum center frequency gradually stabilizes after the value of K becomes bigger than 6. Thus, six IMFs were obtained by the MCFO. However, the simulation signals only contain three components, and the decomposition result is inconsistent with the actual components number. So the method is not accurate in acquiring the IMF number of VMD.Figure 2
Maximum center frequency trend under different component amplitudes.According to CC method, correlation coefficients under different component amplitudes are shown in Table2. For a more intuitive analysis of the results, threshold ranges are plotted in Figure 3. It can be seen that when the component amplitudes are different, range of the correlation coefficient threshold body has different values, but there is an intersection (0.0649, 0.0823) area in each threshold range. When the correlation coefficient is greater than any value in the threshold intersection, three IMFs number can be obtained, which is consistent with the component number of the original signal.Table 2
CC under different component amplitudes.
Signal code
K
Correlation coefficient
Threshold range
S1
3
0.9599
0.2460
0.0837
(0.0627, 0.0823)
4
0.9599
0.2460
0.0823
0.0627
S2
3
0.9044
0.3786
0.1664
(0.0592, 0.1660)
4
0.9044
0.3786
0.1660
0.0592
S3
3
0.8224
0.4825
0.2856
(0.0571, 0.2849)
4
0.8224
0.4824
0.2849
0.0571
S4
3
0.9360
0.3158
0.1127
(0.0649, 0.1110)
4
0.9360
0.3158
0.1110
0.0649Figure 3
CC threshold range under different component amplitudes.According to NMI method, normalized mutual information under different component amplitudes is shown in Table3. It can be seen that there is no threshold range for the signal S4. Because when K = 4, the minimum value of NMI is 0.0658, which is larger than the minimum value of NMI 0.0579 when K = 3. In order to get a clearer understanding of the results, threshold ranges are plotted in Figure 4. It can be seen from analysis results of signal S4 that there is no threshold range. Therefore, there is no threshold intersection to determine the IMF number. So it is invalid to use NMI method for this signal.Table 3
NMI under different component amplitudes.
Signal code
K
Normalized mutual information
Threshold range
S1
3
1.0000
0.1692
0.0663
(0.0641, 0.0663)
4
1.0000
0.1685
0.0743
0.0641
S2
3
1.0000
0.2506
0.0930
(0.0605, 0.0930)
4
1.0000
0.2504
0.1002
0.0605
S3
3
1.0000
0.3442
0.1902
(0.1644, 0.1902)
4
1.0000
0.3439
0.1644
0.0756
S4
3
1.0000
0.1781
0.0579
No
4
1.0000
0.1930
0.0658
0.0833Figure 4
NMI threshold under different component amplitudes.According to CFSA method, center frequency histogram under different component amplitudes is shown in Figure5. It can be seen that the dominant center frequencies are about 12.6 Hz, 150.8 Hz, and 1809.6 Hz, respectively. There are three dominant frequencies. In this way, 3 IMFs can be determined, which is consistent with the component number of original signal. The results show that the proposed method is effective for this kind of signal.Figure 5
Center frequency histogram under different component amplitudes. (a) S1, (b) S2, (c) S3, and (d) S4.
(a)
(b)
(c)
(d)


### 3.2.2. Influence of Component Frequency

Signals S5 to S8 were analyzed using four methods in this section. VMD parameters exceptK are initialized as follows.According to MCFO method, the maximum center frequency trend under different component frequencies is shown in Figure6. It can be seen that the maximum center frequency is different when the value of K is from 1 to 4. Center frequencies of the four signals tend to be stable after the value of K becomes bigger than 6 and no longer increases significantly. Therefore, the IMF number is 6; however, the results is inconsistent with the component number of original signal.Figure 6
Maximum center frequency trend under different component frequencies.According to CC method, correlation coefficients under different component frequencies are shown in Table4. The threshold range results are plotted in Figure 7. As can be seen, the intersection of the threshold ranges is 0.0591, 0.1007. When the correlation coefficient threshold is greater than any value in the intersection, it can be determined that the IMF number is 3. The results is consistent with the actual components number of simulated signal.Table 4
CC under different component frequencies.
Signal code
K
Correlation coefficient
Threshold range
S5
3
0.9250
0.2245
0.1007
(0.0578, 0.1007)
4
0.8680
0.4345
0.2189
0.0578
S6
3
0.8649
0.4436
0.2149
(0.0536, 0.2146)
4
0.8649
0.4436
0.2146
0.0536
S7
3
0.9502
0.2080
0.2160
(0.0591, 0.2080)
4
0.8658
0.4403
0.2159
0.0591
S8
3
0.8660
0.4370
0.2242
(0.0540, 0.2239)
4
0.8660
0.4370
0.2239
0.0540Figure 7
CC threshold range under different component frequencies.According to NMI method, normalized mutual information under different component frequencies is shown in Table5. It can be seen that there is no threshold range for signal S7. Because when K = 4, the minimum value of NMI is 0.0614, which is larger than 0.0599, the minimum value of NMI when K = 3. It can be seen from Figure 8 that there is no threshold range intersection, so the K value cannot be determined according to this method.Table 5
NMI under different component frequencies.
Signal code
K
Normalized mutual information
Threshold range
S5
3
1.0000
0.1677
0.1452
(0.1426, 0.1452)
4
1.0000
0.3034
0.2094
0.1426
S6
3
1.0000
0.2993
0.2295
(0.0956, 0.1900)
4
1.0000
0.3003
0.1900
0.0956
S7
3
1.0000
0.1097
0.0599
No
4
1.0000
0.2744
0.1185
0.0614
S8
3
1.0000
0.2864
0.1257
(0.0911, 0.1257)
4
1.0000
0.2675
0.1373
0.0911Figure 8
NMI threshold range under different component frequencies.According to the CFSA method, center frequency histogram under different component amplitudes is shown in Figure9. Three dominant center frequencies are obtained by signal S5 processing, which are approximately 12.6 Hz, 75.4 Hz, and 452.3 Hz, respectively. Another three frequencies (12.6 Hz, 100.5 Hz, and 804.3 Hz) can be obtained by processing signal S6. 12.6 Hz, 125.7 Hz, and 1255.3 Hz can be obtained by processing signal S7. 12.6 Hz, 151 Hz, and 1809.7 Hz can be obtained by processing signal S8. In summary, three IMF numbers can be obtained by the CFSA method, which is consistent with the actual component number of simulation signal.Figure 9
Center frequency histogram under different component frequencies, (a) S5, (b) S6, (c) S7, and (d) S8.
(a)
(b)
(c)
(d)


### 3.2.3. Influence of Component Number

Signals S9 to S12 were analyzed by four methods in this section. VMD parameters exceptK are initialized as follows.According to MCFO method, the maximum center frequency trend under different components number is shown in Figure10. It can be seen that when the simulation signal contains 3, 4, and 5 components, the maximum center frequency trends are almost the same. However, when the simulation contains 6 components, the maximum center frequency trend is different from other situations. All maximum values of center frequency tend to be stable after the value of K becomes bigger than 7. Therefore, the IMF number is 7. But the results are inconsistent with the actual component number of the simulation signal.Figure 10
Maximum center frequency trend under different components number.According to the CC method, the correlation coefficients under different components number are shown in Table6. The threshold ranges are plotted in Figure 11. As can be seen, there is no threshold intersection for determining the IMF number. Therefore, when VMD is used to decompose the signals S9 to S12, it is difficult for CC method to determine the IMF number.Table 6
CC under different components number.
Signal code
K
Correlation coefficient
Threshold range
S9
3
0.9825
0.0843
0.0899
(0.0611, 0.0843)
4
0.9586
0.2496
0.0854
0.0611
S10
4
0.9750
0.1556
0.0846
0.0617
(0.0611, 0.0617)
5
0.9749
0.1547
0.0704
0.0831
0.0611
S11
5
0.9711
0.1479
0.1163
0.0821
0.0639
(0.0573, 0.0639)
6
0.9710
0.1477
0.1160
0.0798
0.0607
0.0573
S12
6
0.9673
0.1492
0.1183
0.0916
0.0821
0.0582
(0.0538, 0.0557)
7
0.9673
0.1491
0.1181
0.0910
0.0800
0.0538
0.0557Figure 11
CC threshold range under different components number.According to the NMI method, the normalized mutual information under different components number is shown in Table7. It can be seen from Figure 12 that there is no threshold range for the signal S9, because the minimum value of NMI is 0.0502, which is larger than 0.0436. Threshold ranges are plotted in Figure 13. As can be seen, there is no threshold intersection. So the NMI method cannot determine the IMF number when the signals S9 to S12 are processed by VMD.Table 7
NMI under different components number.
Signal code
K
Normalized mutual information
Threshold range
S9
3
1.0000
0.0508
0.0436
No
4
1.0000
0.0560
0.0543
0.0502
S10
4
1.0000
0.0772
0.0472
0.0577
(0.0429, 0.0472)
5
1.0000
0.0728
0.0688
0.0695
0.0429
S11
5
1.0000
0.0734
0.0610
0.0530
0.0670
(0.0458, 0.0469)
6
1.0000
0.0694
0.0607
0.0458
0.0469
0.0500
S12
6
1.0000
0.0740
0.0544
0.0530
0.0534
0.0594
(0.0500, 0.0511)
7
1.0000
0.0771
0.0654
0.0511
0.0678
0.0500
0.0544Figure 12
NMI threshold range under different components number.Figure 13
Center frequency histogram under different components number, (a) S9, (b) S10, (c) S11, and (d) S12.
(a)
(b)
(c)
(d)According to CFSA, center frequency histogram under different components number is shown in Figure13. From the center frequency histogram obtained by processing the signal S9, there are three dominant center frequencies (12.6 Hz, 151.1 Hz, and 1808.5 Hz). And there are four dominant center frequencies (12.6 Hz, 150.9 Hz, 627.9 Hz, and 1809.8 Hz) obtained by signal S10 processing. Five dominant center frequencies (12.6 Hz, 150.9 Hz, 627.5 Hz, 1129 Hz, and 1808.7 Hz) can be obtained by processing signal S11. Six dominant center frequencies (12.6 Hz, 150.8 Hz, 628.3 Hz, 1131 Hz, 1508 Hz, and 1809.6 Hz) can be obtained by processing signal S12. The simulation signals S9 to S12 contain exactly 3 to 6 components, respectively, so the CFSA is correctly verified.


### 3.2.4. Influence of Noise Amplitude

Signals S13 to S16 were analyzed by four methods in this section. VMD parameters exceptK are initialized as follows.According to MCFO method, the maximum center frequency trend under different noise amplitudes is shown in Figure14. It can be seen that when there is no noise, the maximum center frequency tends to be stable after the value of K becomes bigger than 2. When the noise amplitude is 0.05 and 0.1, the maximum center frequency tends to be stable from K = 6. When the noise amplitude is 0.15, the maximum center frequency tends to be stable from K = 7. Therefore, when the noise amplitude is different, IMF number determined by MCFO method is different.Figure 14
Maximum center frequency trend under different noise amplitudes.According to the CC method, the correlation coefficients under different noise amplitudes are shown in Table8. Different threshold ranges are shown in Figure 15. There is no intersection in the four threshold ranges, so the value of K cannot be determined. Therefore, when VMD is applied to decompose signals S13 to S16, it is not suitable for CC method to determine the IMF number.Table 8
CC under different noise amplitudes.
Signal code
K
Correlation coefficient
Threshold range
S13
3
0.9684
0.2421
0.0604
(0.0237, 0.0596)
4
0.9684
0.2421
0.0596
0.0237
S14
3
0.9673
0.2392
0.0663
(0.0307, 0.0657)
4
0.9673
0.2392
0.0657
0.0307
S15
3
0.9825
0.0850
0.0871
(0.0612, 0.0835)
4
0.9588
0.2488
0.0835
0.0612
S16
3
0.9742
0.1121
0.1049
(0.0911, 0.1004)
4
0.9741
0.1104
0.1004
0.0911Figure 15
CC threshold range under different noise amplitudes.According to the NMI method, the normalized mutual information under different noise amplitudes is shown in Table9. Four threshold ranges are shown in Figure 16. There is no intersection of the four threshold ranges, so there is no mutual information threshold to determine the K. Therefore, the NMI method is ineffective in this case.Table 9
NMI under different noise amplitudes.
Signal code
K
Normalized mutual information
Threshold range
S13
3
1.0000
0.2450
0.1165
(0.0753, 0.1165)
4
1.0000
0.2450
0.0753
0.1513
S14
3
1.0000
0.1688
0.0452
(0.0431, 0.0452)
4
1.0000
0.1784
0.0431
0.0505
S15
3
1.0000
0.1716
0.0708
(0.0634, 0.0708)
4
1.0000
0.1672
0.0634
0.0552
S16
3
1.0000
0.1439
0.0586
(0.0532, 0.0586)
4
1.0000
0.0693
0.0729
0.0532Figure 16
NMI threshold range under different noise amplitudes.According to CFSA, center frequency histogram under different noise amplitudes is shown in Figure17. As can be seen, there are three dominant center frequencies in results of signals S13 to S16, which corresponds to the actual components number of signals.Figure 17
Center frequency histogram under different noise amplitudes, (a) S13, (b) S14, (c) S15, and (d) S16.
(a)
(b)
(c)
(d)In summary, the MCFO and NMI methods cannot accurately determine the IMF number in case of four influencing factors changes. CC method can only determine the IMF number accurately under first two influencing factors changes. The proposed CFSA method can obtain the IMF number accurately under the variation of four influencing factors, which verifies the advantages and effectiveness of the proposed method.


## 3.2.1. Influence of Component Amplitude

Signals S1 to S4 were analyzed by methods mentioned above in this section. VMD parameters exceptK are initialized as follows: α=2000,τ=0,DC=0,init=1,ε=1e−7.According to MCFO method, the maximum center frequency trend under different component amplitudes is shown in Figure2. As can be seen, the maximum center frequency gradually stabilizes after the value of K becomes bigger than 6. Thus, six IMFs were obtained by the MCFO. However, the simulation signals only contain three components, and the decomposition result is inconsistent with the actual components number. So the method is not accurate in acquiring the IMF number of VMD.Figure 2
Maximum center frequency trend under different component amplitudes.According to CC method, correlation coefficients under different component amplitudes are shown in Table2. For a more intuitive analysis of the results, threshold ranges are plotted in Figure 3. It can be seen that when the component amplitudes are different, range of the correlation coefficient threshold body has different values, but there is an intersection (0.0649, 0.0823) area in each threshold range. When the correlation coefficient is greater than any value in the threshold intersection, three IMFs number can be obtained, which is consistent with the component number of the original signal.Table 2
CC under different component amplitudes.
Signal code
K
Correlation coefficient
Threshold range
S1
3
0.9599
0.2460
0.0837
(0.0627, 0.0823)
4
0.9599
0.2460
0.0823
0.0627
S2
3
0.9044
0.3786
0.1664
(0.0592, 0.1660)
4
0.9044
0.3786
0.1660
0.0592
S3
3
0.8224
0.4825
0.2856
(0.0571, 0.2849)
4
0.8224
0.4824
0.2849
0.0571
S4
3
0.9360
0.3158
0.1127
(0.0649, 0.1110)
4
0.9360
0.3158
0.1110
0.0649Figure 3
CC threshold range under different component amplitudes.According to NMI method, normalized mutual information under different component amplitudes is shown in Table3. It can be seen that there is no threshold range for the signal S4. Because when K = 4, the minimum value of NMI is 0.0658, which is larger than the minimum value of NMI 0.0579 when K = 3. In order to get a clearer understanding of the results, threshold ranges are plotted in Figure 4. It can be seen from analysis results of signal S4 that there is no threshold range. Therefore, there is no threshold intersection to determine the IMF number. So it is invalid to use NMI method for this signal.Table 3
NMI under different component amplitudes.
Signal code
K
Normalized mutual information
Threshold range
S1
3
1.0000
0.1692
0.0663
(0.0641, 0.0663)
4
1.0000
0.1685
0.0743
0.0641
S2
3
1.0000
0.2506
0.0930
(0.0605, 0.0930)
4
1.0000
0.2504
0.1002
0.0605
S3
3
1.0000
0.3442
0.1902
(0.1644, 0.1902)
4
1.0000
0.3439
0.1644
0.0756
S4
3
1.0000
0.1781
0.0579
No
4
1.0000
0.1930
0.0658
0.0833Figure 4
NMI threshold under different component amplitudes.According to CFSA method, center frequency histogram under different component amplitudes is shown in Figure5. It can be seen that the dominant center frequencies are about 12.6 Hz, 150.8 Hz, and 1809.6 Hz, respectively. There are three dominant frequencies. In this way, 3 IMFs can be determined, which is consistent with the component number of original signal. The results show that the proposed method is effective for this kind of signal.Figure 5
Center frequency histogram under different component amplitudes. (a) S1, (b) S2, (c) S3, and (d) S4.
(a)
(b)
(c)
(d)


## 3.2.2. Influence of Component Frequency

Signals S5 to S8 were analyzed using four methods in this section. VMD parameters exceptK are initialized as follows.According to MCFO method, the maximum center frequency trend under different component frequencies is shown in Figure6. It can be seen that the maximum center frequency is different when the value of K is from 1 to 4. Center frequencies of the four signals tend to be stable after the value of K becomes bigger than 6 and no longer increases significantly. Therefore, the IMF number is 6; however, the results is inconsistent with the component number of original signal.Figure 6
Maximum center frequency trend under different component frequencies.According to CC method, correlation coefficients under different component frequencies are shown in Table4. The threshold range results are plotted in Figure 7. As can be seen, the intersection of the threshold ranges is 0.0591, 0.1007. When the correlation coefficient threshold is greater than any value in the intersection, it can be determined that the IMF number is 3. The results is consistent with the actual components number of simulated signal.Table 4
CC under different component frequencies.
Signal code
K
Correlation coefficient
Threshold range
S5
3
0.9250
0.2245
0.1007
(0.0578, 0.1007)
4
0.8680
0.4345
0.2189
0.0578
S6
3
0.8649
0.4436
0.2149
(0.0536, 0.2146)
4
0.8649
0.4436
0.2146
0.0536
S7
3
0.9502
0.2080
0.2160
(0.0591, 0.2080)
4
0.8658
0.4403
0.2159
0.0591
S8
3
0.8660
0.4370
0.2242
(0.0540, 0.2239)
4
0.8660
0.4370
0.2239
0.0540Figure 7
CC threshold range under different component frequencies.According to NMI method, normalized mutual information under different component frequencies is shown in Table5. It can be seen that there is no threshold range for signal S7. Because when K = 4, the minimum value of NMI is 0.0614, which is larger than 0.0599, the minimum value of NMI when K = 3. It can be seen from Figure 8 that there is no threshold range intersection, so the K value cannot be determined according to this method.Table 5
NMI under different component frequencies.
Signal code
K
Normalized mutual information
Threshold range
S5
3
1.0000
0.1677
0.1452
(0.1426, 0.1452)
4
1.0000
0.3034
0.2094
0.1426
S6
3
1.0000
0.2993
0.2295
(0.0956, 0.1900)
4
1.0000
0.3003
0.1900
0.0956
S7
3
1.0000
0.1097
0.0599
No
4
1.0000
0.2744
0.1185
0.0614
S8
3
1.0000
0.2864
0.1257
(0.0911, 0.1257)
4
1.0000
0.2675
0.1373
0.0911Figure 8
NMI threshold range under different component frequencies.According to the CFSA method, center frequency histogram under different component amplitudes is shown in Figure9. Three dominant center frequencies are obtained by signal S5 processing, which are approximately 12.6 Hz, 75.4 Hz, and 452.3 Hz, respectively. Another three frequencies (12.6 Hz, 100.5 Hz, and 804.3 Hz) can be obtained by processing signal S6. 12.6 Hz, 125.7 Hz, and 1255.3 Hz can be obtained by processing signal S7. 12.6 Hz, 151 Hz, and 1809.7 Hz can be obtained by processing signal S8. In summary, three IMF numbers can be obtained by the CFSA method, which is consistent with the actual component number of simulation signal.Figure 9
Center frequency histogram under different component frequencies, (a) S5, (b) S6, (c) S7, and (d) S8.
(a)
(b)
(c)
(d)


## 3.2.3. Influence of Component Number

Signals S9 to S12 were analyzed by four methods in this section. VMD parameters exceptK are initialized as follows.According to MCFO method, the maximum center frequency trend under different components number is shown in Figure10. It can be seen that when the simulation signal contains 3, 4, and 5 components, the maximum center frequency trends are almost the same. However, when the simulation contains 6 components, the maximum center frequency trend is different from other situations. All maximum values of center frequency tend to be stable after the value of K becomes bigger than 7. Therefore, the IMF number is 7. But the results are inconsistent with the actual component number of the simulation signal.Figure 10
Maximum center frequency trend under different components number.According to the CC method, the correlation coefficients under different components number are shown in Table6. The threshold ranges are plotted in Figure 11. As can be seen, there is no threshold intersection for determining the IMF number. Therefore, when VMD is used to decompose the signals S9 to S12, it is difficult for CC method to determine the IMF number.Table 6
CC under different components number.
Signal code
K
Correlation coefficient
Threshold range
S9
3
0.9825
0.0843
0.0899
(0.0611, 0.0843)
4
0.9586
0.2496
0.0854
0.0611
S10
4
0.9750
0.1556
0.0846
0.0617
(0.0611, 0.0617)
5
0.9749
0.1547
0.0704
0.0831
0.0611
S11
5
0.9711
0.1479
0.1163
0.0821
0.0639
(0.0573, 0.0639)
6
0.9710
0.1477
0.1160
0.0798
0.0607
0.0573
S12
6
0.9673
0.1492
0.1183
0.0916
0.0821
0.0582
(0.0538, 0.0557)
7
0.9673
0.1491
0.1181
0.0910
0.0800
0.0538
0.0557Figure 11
CC threshold range under different components number.According to the NMI method, the normalized mutual information under different components number is shown in Table7. It can be seen from Figure 12 that there is no threshold range for the signal S9, because the minimum value of NMI is 0.0502, which is larger than 0.0436. Threshold ranges are plotted in Figure 13. As can be seen, there is no threshold intersection. So the NMI method cannot determine the IMF number when the signals S9 to S12 are processed by VMD.Table 7
NMI under different components number.
Signal code
K
Normalized mutual information
Threshold range
S9
3
1.0000
0.0508
0.0436
No
4
1.0000
0.0560
0.0543
0.0502
S10
4
1.0000
0.0772
0.0472
0.0577
(0.0429, 0.0472)
5
1.0000
0.0728
0.0688
0.0695
0.0429
S11
5
1.0000
0.0734
0.0610
0.0530
0.0670
(0.0458, 0.0469)
6
1.0000
0.0694
0.0607
0.0458
0.0469
0.0500
S12
6
1.0000
0.0740
0.0544
0.0530
0.0534
0.0594
(0.0500, 0.0511)
7
1.0000
0.0771
0.0654
0.0511
0.0678
0.0500
0.0544Figure 12
NMI threshold range under different components number.Figure 13
Center frequency histogram under different components number, (a) S9, (b) S10, (c) S11, and (d) S12.
(a)
(b)
(c)
(d)According to CFSA, center frequency histogram under different components number is shown in Figure13. From the center frequency histogram obtained by processing the signal S9, there are three dominant center frequencies (12.6 Hz, 151.1 Hz, and 1808.5 Hz). And there are four dominant center frequencies (12.6 Hz, 150.9 Hz, 627.9 Hz, and 1809.8 Hz) obtained by signal S10 processing. Five dominant center frequencies (12.6 Hz, 150.9 Hz, 627.5 Hz, 1129 Hz, and 1808.7 Hz) can be obtained by processing signal S11. Six dominant center frequencies (12.6 Hz, 150.8 Hz, 628.3 Hz, 1131 Hz, 1508 Hz, and 1809.6 Hz) can be obtained by processing signal S12. The simulation signals S9 to S12 contain exactly 3 to 6 components, respectively, so the CFSA is correctly verified.


## 3.2.4. Influence of Noise Amplitude

Signals S13 to S16 were analyzed by four methods in this section. VMD parameters exceptK are initialized as follows.According to MCFO method, the maximum center frequency trend under different noise amplitudes is shown in Figure14. It can be seen that when there is no noise, the maximum center frequency tends to be stable after the value of K becomes bigger than 2. When the noise amplitude is 0.05 and 0.1, the maximum center frequency tends to be stable from K = 6. When the noise amplitude is 0.15, the maximum center frequency tends to be stable from K = 7. Therefore, when the noise amplitude is different, IMF number determined by MCFO method is different.Figure 14
Maximum center frequency trend under different noise amplitudes.According to the CC method, the correlation coefficients under different noise amplitudes are shown in Table8. Different threshold ranges are shown in Figure 15. There is no intersection in the four threshold ranges, so the value of K cannot be determined. Therefore, when VMD is applied to decompose signals S13 to S16, it is not suitable for CC method to determine the IMF number.Table 8
CC under different noise amplitudes.
Signal code
K
Correlation coefficient
Threshold range
S13
3
0.9684
0.2421
0.0604
(0.0237, 0.0596)
4
0.9684
0.2421
0.0596
0.0237
S14
3
0.9673
0.2392
0.0663
(0.0307, 0.0657)
4
0.9673
0.2392
0.0657
0.0307
S15
3
0.9825
0.0850
0.0871
(0.0612, 0.0835)
4
0.9588
0.2488
0.0835
0.0612
S16
3
0.9742
0.1121
0.1049
(0.0911, 0.1004)
4
0.9741
0.1104
0.1004
0.0911Figure 15
CC threshold range under different noise amplitudes.According to the NMI method, the normalized mutual information under different noise amplitudes is shown in Table9. Four threshold ranges are shown in Figure 16. There is no intersection of the four threshold ranges, so there is no mutual information threshold to determine the K. Therefore, the NMI method is ineffective in this case.Table 9
NMI under different noise amplitudes.
Signal code
K
Normalized mutual information
Threshold range
S13
3
1.0000
0.2450
0.1165
(0.0753, 0.1165)
4
1.0000
0.2450
0.0753
0.1513
S14
3
1.0000
0.1688
0.0452
(0.0431, 0.0452)
4
1.0000
0.1784
0.0431
0.0505
S15
3
1.0000
0.1716
0.0708
(0.0634, 0.0708)
4
1.0000
0.1672
0.0634
0.0552
S16
3
1.0000
0.1439
0.0586
(0.0532, 0.0586)
4
1.0000
0.0693
0.0729
0.0532Figure 16
NMI threshold range under different noise amplitudes.According to CFSA, center frequency histogram under different noise amplitudes is shown in Figure17. As can be seen, there are three dominant center frequencies in results of signals S13 to S16, which corresponds to the actual components number of signals.Figure 17
Center frequency histogram under different noise amplitudes, (a) S13, (b) S14, (c) S15, and (d) S16.
(a)
(b)
(c)
(d)In summary, the MCFO and NMI methods cannot accurately determine the IMF number in case of four influencing factors changes. CC method can only determine the IMF number accurately under first two influencing factors changes. The proposed CFSA method can obtain the IMF number accurately under the variation of four influencing factors, which verifies the advantages and effectiveness of the proposed method.


## 4. Experimental Validation

To further verify the superiority of the proposed method, four methods were tested by rolling bearing experimental data from Case Western Reserve University [20]. The experimental data parameters are set as follows: the load is 3 horsepower and the sampling frequency is 12 kHz. 12000 data points were used for analysis. Time domain signal and FFT spectrum of normal (Normal_3), inner race fault (IR028_3), outer race fault (OR021@6_3), and ball fault (B028_3) state are as shown in Figure 18.Figure 18
Bearing vibration signal and FFT spectrum, (a) Normal_3, (b) IR028_3, (c) OR021@6_3, and (d) B028_3.
(a)
(b)
(c)
(d)In the MCFO method, VMD parameters exceptK are initialized as α=2000,τ=0,DC=1,init=2,ε=1e−7. The maximum center frequency trend results are shown in Figure 19. As can be seen from Figure 19(a), the maximum center frequency fluctuates two times when K = 4 and 9, so the accurate K value cannot be obtained. As can be seen from Figure 19(b), the maximum center frequency tends to be stable from K = 5, so the number of IMF is 5. Similarly, it can be seen from Figures 19(c) and 19(d) that K = 3 in outer race fault and rolling ball fault situation.Figure 19
Maximum center frequency trend of bearing vibration signal, (a) Normal_3, (b) IR028_3, (c) OR021@6_3, and (d) B028_3.
(a)
(b)
(c)
(d)In the CC method, VMD parameters exceptK are initialized as α=2000,τ=0,DC=1,init=2,ε=1e−7. Correlation coefficients were shown in Table 10. As can be seen, there exists no obvious threshold, and the fixed numbers of correlation coefficients are always greater than this threshold from a certain K. So the accurate K value cannot be obtained according to correlation coefficients of four state signals.Table 10
CC of bearing vibration signal.
K
Normal_3
K
IR028_3
1
0.67
1
0.71
2
0.63
0.62
2
0.52
0.67
3
0.62
0.59
0.49
3
0.39
0.64
0.47
4
0.62
0.59
0.51
0.51
4
0.28
0.64
0.44
0.60
5
0.51
0.59
0.56
0.49
0.05
5
0.21
0.58
0.43
0.33
0.60
6
0.45
0.58
0.39
0.49
0.24
0.48
6
0.17
0.32
0.58
0.39
0.32
0.59
K
OR021@6_3
K
B028_3
1
0.36
1
0.46
2
0.26
0.89
2
0.16
0.92
3
0.25
0.23
0.89
3
0.09
0.84
0.60
4
0.17
0.22
0.44
0.87
4
0.07
0.56
0.48
0.78
5
0.13
0.21
0.44
0.17
0.87
5
0.05
0.48
0.49
0.71
0.45
6
0.12
0.19
0.17
0.87
0.44
0.15
6
0.05
0.71
0.49
0.45
0.48
0.18In the NMI method, VMD parameters exceptK are initialized as α=2000,τ=0,DC=1,init=2,ε=1e−7. The NMI results of four state signals were shown in Table 11. As can be seen, the K value cannot be determined according to the normalized mutual information values. Because starting from a certain layer, you cannot find a threshold that is less than a fixed amount of mutual information.Table 11
NMI of bearing vibration signal.
K
Normal_3
K
IR028_3
1
1.00
1
1.00
2
1.00
0.95
2
0.56
1.00
3
1.00
0.89
0.59
3
0.37
1.00
0.48
4
1.00
0.89
0.63
0.64
4
0.22
1.00
0.42
0.81
5
0.70
1.00
0.90
0.66
0.39
5
0.15
1.00
0.48
0.36
0.97
6
0.57
1.00
0.42
0.73
0.19
0.66
6
0.11
0.34
1.00
0.42
0.30
0.97
K
OR021@6_3
K
B028_3
1
1.00
1
1.00
2
0.29
1.00
2
0.08
1.00
3
0.29
0.32
1.00
3
0.09
1.00
0.43
4
0.25
0.34
0.47
1.00
4
0.12
0.48
0.34
1.00
5
0.26
0.33
0.47
0.33
1.00
5
0.14
0.43
0.47
1.00
0.39
6
0.29
0.24
0.33
1.00
0.47
0.20
6
0.14
1.00
0.46
0.37
0.43
0.11In the CFSA method, VMD parameters exceptK are initialized as α=2000,τ=0,DC=0,init=2,ε=1e−7. Center frequency histograms obtained by CFSA about four states of bearings were shown in Figure 20. As can be seen from Figure 20(a), the IMF number is 5 because there are 5 center frequencies whose counts are higher than the average. Similarly, the IMF number in the other three situations is 4.Figure 20
Center frequency histogram of bearing vibration signal, (a) normal_3, (b) IR028_3, (c) OR021@6_3, and (d) B028_3.
(a)
(b)
(c)
(d)In order to quantitatively evaluate the accuracy of the proposed method, the accuracy of the methods can be evaluated by errorce which can be expressed as(9)ce=∑i=1Nei⋅aiN,ei=fc,i−fo,ifo,i,where fo,i are the main frequencies of FFT spectrum. fc,i represents center frequencies obtained by four methods. ai represents the amplitude corresponding to fo,i. N is the number of center frequencies.According to theK values determined by four methods, the bearing vibration signals in four states are decomposed by VMD, and the center frequencies of IMFs in each state are obtained. The error of each method is calculated by formula (9),and the results are shown in Table 12.Table 12
Center frequency and comprehensive error.
Signal
Methods
Center frequency (Hz)
Error
Normal_3
FFT
90
234
602
1065
2102
—
CFSA (K = 5)
90
234
602
1065
2102
0
IR028_3
FFT
616
985
1481
2934
—
MCFO (K = 5)
616
763
1467
2782
0
2.2%
CFSA (K = 4)
616
739
1481
2934
0
OR021@6_3
FFT
618
1339
2724
3399
—
MCFO (K = 3)
599
3436
0
1.3%
CFSA (K = 4)
618
1339
2724
3399
0
B028_3
FFT
2714
2850
3143
3395
—
MCFO (K = 3)
2850
3311
0
0.7%
CFSA (K = 4)
2714
2850
3143
3395
0As you can see from the table, CC and NMI do not exist because these two methods do not get aK value. In Normal_3 state, the value of K is 5 obtained by CFSA method, while MCFO methods cannot obtain K value. In IR028_3 state, the value of K is 4 obtained by CFSA with no difference between center frequencies of FFT. However, the value of K is 5 obtained by MCFO with 2.2% error. In the state of outer race fault and ball fault, the K value obtained by CFSA method is 4, and the center frequencies results are completely consistent with the main components of FFT spectrum. However, the K values obtained by MCFO method are 3, and the error of outer race and ball fault vibration signals processing results are 1.3% and 0.7%, respectively. To sum up, the CFSA method can obtain the appropriate IMF number when VMD is used to process the bearing vibration signals.


## 5. Conclusion

This study proposes a novel method based on IMFs center frequency to determine IMF number of VMD. Compared with MCFO, CC, and NMI methods, the proposed method CFSA can accurately obtain theK value from complex simulation signal with variable components. The proposed method was demonstrated to be effective for processing rolling bearing vibration signal. Center frequencies obtained by processing the vibration signals of bearing in normal and fault state are consistent with the main components of FFT spectrum. Results show that the proposed method has high robustness and accuracy. In future work, the effectiveness of the proposed method needs to be further verified by processing other types of vibration signals, for example, vibration signal of planetary gearbox or wind gearbox.



---

*Source: 8304903-2020-05-07.xml*